Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

On nonlinear free surface potential flow by a Bubnov-Galerkin formulation in space and a semi-lagrangian… Allievi, Alejandro 1993

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1993 A1 A64.pdf [ 6.96MB ]
Metadata
JSON: 831-1.0302316.json
JSON-LD: 831-1.0302316-ld.json
RDF/XML (Pretty): 831-1.0302316-rdf.xml
RDF/JSON: 831-1.0302316-rdf.json
Turtle: 831-1.0302316-turtle.txt
N-Triples: 831-1.0302316-rdf-ntriples.txt
Original Record: 831-1.0302316-source.json
Full Text
831-1.0302316-fulltext.txt
Citation
831-1.0302316.ris

Full Text

ON NONLINEAR FREE SURFACE POTENTIAL FLOW BY A BUBNOV-GALERKIN FORMULATION IN SPACE AND A SEMI-LAGRANGIAN SEMI-IMPLICIT SCHEME IN TIME A L E J A N D R O A L L I E V I Dipl. Ing., Universidad de Buenos Aires, 1980 M.A.Sc., University of British Columbia, 1987 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R OF P H I L O S O P H Y in The Faculty of Graduate Studies Department of Mechanical Engineering We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A July 1993 © Alejandro Al l ievi , 1993 02 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the U N I V E R S I T Y O F B R I T I S H C O L U M B I A , 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. D E P A R T M E N T O F M E C H A N I C A L E N G I N E E R I N G U N I V E R S I T Y O F B R I T I S H C O L U M B I A 2075 Wesbrook Place Vancouver, B . C . Canada V 6 T 1Z4 Date: J U L Y , 1993 ABSTRACT The potential flow initial-boundary value problem describing fluid-structure interac-tion with fully nonlinear free surface boundary conditions has been studied using a mixed Lagrangian-Eulerian formulation. The boundary-value problem has been solved in the phys-ical domain by means of a Bubnov-Galerkin formulation of the Laplace equation. The init ial-value problem related to the behavior of some of the moving boundaries has been discretized using various numerical techniques. Among these is a series of predictor-corrector methods. These methodologies proved to require considerable numerical smoothing to maintain sta-bility of the numerical scheme. In turn, dissipation led to inaccuracies in the solution of the problem. In order to avoid this negative effect, a semi-implicit semi-Lagrangian two-time level iterative scheme that is almost free from smoothing has been developed. A Bubnov-Galerkin formulation of an elliptic system for the generation of boundary fitted curvilinear coordinates has been used. When solved iteratively, this method provides orthog-onal meshes of very good characteristics for both symmetric and non-symmetric domains. Previous publications concluded that the present system was inadequate for non-symmetric regions leading to lack of convergence in the iterative process. Solutions described in this work show that this limitation has been overcome. F lu id responses to periodic excitation of surface-piercing and submerged bodies have been calculated. Both linear and nonlinear cases show agreement with published results. Very low total energy/work error has been obtained which demonstrates accuracy, good stability and convergence characteristics of the numerical scheme. The impulsive response of tanks of various shapes has also been simulated. Resulting natural frequencies show good agreement with available data. A slender body representation of the flow around a hull advancing with forward speed in otherwise calm water has also been simulated. Numerical calculations of a number of quantities of engineering interest are presented for different length Froude numbers. Results compare favorably with experimental data. n T A B L E OF C O N T E N T S A B S T R A C T i i LIST O F T A B L E S v i LIST O F F I G U R E S v i i LIST O F S Y M B O L S x A C K N O W L E D G E M E N T x i i i 1. I N T R O D U C T I O N 1 1.1 S L E N D E R B O D Y T H E O R Y 2 1.2 A B R I E F R E V I E W O F F R E E S U R F A C E F L O W A N A L Y S I S 6 1.3 O B J E C T I V E S O F THIS I N V E S T I G A T I O N 11 2. G O V E R N I N G E Q U A T I O N S 14 2.1 P O T E N T I A L F L O W W I T H A N O N L I N E A R F R E E S U R F A C E . . . 14 2.2 O R D E R O F M A G N I T U D E A N A L Y S I S 17 2.3 L A G R A N G I A N D E S C R I P T I O N O F T H E P R O B L E M 24 3. W E I G H T E D R E S I D U A L F O R M U L A T I O N 26 3.1 B U B N O V - G A L E R K I N F O R M U L A T I O N 26 3.2 T H E I S O P A R A M E T R I C E L E M E N T 29 3.3 F I N I T E E L E M E N T F O R M U L A T I O N 32 3.4 N U M E R I C A L I N T E G R A T I O N 34 4. E L L I P T I C G R I D G E N E R A T I O N 35 4.1 E L L I P T I C G R I D G E N E R A T I N G S Y S T E M S 35 4.2 P O I S S O N S Y S T E M S 37 i i i 4.3 C O O R D I N A T E T R A N S F O R M A T I O N 37 4.4 B U B N O V - G A L E R K I N P R O C E D U R E 42 4.5 F I N I T E E L E M E N T F O R M U L A T I O N 44 4.6 C A L C U L A T I O N O F T H E S H A P E F A C T O R 46 4.7 N U M E R I C A L C O N S I D E R A T I O N S 47 4.8 N U M E R I C A L R E S U L T S 50 4.8.1 Case A : Complete Boundary Correspondence 51 4.8.2 Case B : Dirichlet and Neumann Boundary Conditions 55 4.8.3 Case C: Dirichlet and Neumann Boundary Conditions 57 N U M E R I C A L CONSIDERATIONS 62 5.1 B E R N O U L L I ' S C O N S T A N T 62 5.2 T R E A T M E N T O F T H E B O U N D A R Y C O N D I T I O N S 63 5.2.1 Free Surface Boundary Conditions 63 5.2.2 Body Boundary Condition 64 5.2.3 Treatment of the Right Side in Equation 3.28 65 5.3 S C H E M E S U S E D T O D I S C R E T I Z E T I M E D E R I V A T I V E S 66 5.3.1 Predictor-Corrector Methods 66 5.3.2 Semi-implicit Semi-Lagrangian Scheme 70 5.4 C A L C U L A T I O N O F F O R C E S , W O R K A N D E N E R G Y 72 5.5 A L T E R N A T I V E F O R M F O R T H E C A L C U L A T I O N O F V E L O C I T I E S 74 RESULTS A N D DISCUSSION 75 6.1 P E R I O D I C M O T I O N IN A T A N K : L I N E A R F R E E S U R F A C E . . . 75 6.1.1 Surface-Piercing Cylinder 76 6.1.1.1 Physical Aspects * 76 6.1.1.2 Computational Aspects 77 6.1.2 Submerged Cylinder 81 6.1.2.1 Physical Aspects 82 6.1.2.2 Computational Aspects 83 iv 6.2 P E R I O D I C M O T I O N IN A T A N K : N O N L I N E A R F R E E S U R F A C E . 84 6.2.1 Surface-Piercing Cylinder 84 6.2.1.1 Physical Aspects 86 6.2.1.2 Computational Aspects 89 6.2.2 Submerged Cylinder 90 6.3 I M P U L S I V E M O T I O N IN A T A N K . 91 6.3.1 Square tank 91 6.3.2 Circular tank 92 6.3.3 Parabolic tank 93 6.3.5 Triangular tank 93 6.3.5 El l ipt ic tank 94 6.4 S L E N D E R B O D Y W I T H F O R W A R D S P E E D IN A T A N K 95 C O N C L U D I N G R E M A R K S 135 7.1 F U T U R E R E C O M M E N D A T I O N S 137 B I B L I O G R A P H Y 140 C O N C E P T S O N D I F F E R E N T I A L G E O M E T R Y 150 LIST OF TABLES 4-1 Boundary conditions as given by Chikliwala and Yortsos 50 4-2 M D O , A D O and number of iterations for Case A 52 4-3 M D O , A D O and number of iterations for Case A 54 4-4 M D O , A D O and number of iterations for Case A 54 4-5 M D O , A D O and number of iterations for Case B 57 4-6 M D O , A D O and number of iterations for Case C 59 4- 7 M D O , A D O and number of iterations for Case C 61 5- 1 Explicit Formulae 68 5- 2 Implicit Formulae 68 6- 1 R M S error and number of iterations of conjugate gradient method 79 6-2 R M S error and number of iterations of conjugate gradient method 80 6-3 R M S error and number of iterations of conjugate gradient method 81 6-4 R M S error and number of iterations of conjugate gradient method 82 6-5 R M S error and number of iterations of conjugate gradient method 83 6-6 R M S error and number of iterations of conjugate gradient method 84 6-7 R M S error and number of iterations of conjugate gradient method 85 6-8 R M S error and number of iterations of conjugate gradient method 86 6-9 R M S error and number of iterations of conjugate gradient method 87 6-10 R M S error and number of iterations of conjugate gradient method 88 6-11 Natural frequencies for square tank 92 6-12 Natural frequencies for circular tank 93 6-13 Natural frequencies for parabolic tank 93 6-14 Natural frequencies for triangular tank 94 6-15 Natural frequencies for elliptic tank 95 6-16 Wigley H u l l Information 95 6-17 Computational Conditions for Wigley Hul l Calculations 97 6-18 Computational Demands for Wigley Hul l Calculations 97 vi LIST OF FIGURES 3- 1 Linear (4-noded) and quadratic (8-noded) isoparametric elements 30 4- 1 Typical Domain 51 4-2 16 X 16 grid, H — 0.15, ISH A — 4. Dirichlet boundary conditions with exponential distribution of points on left boundary 53 4-3 32 X 16 grid, H — 0.15, ISH A — 3. Dirichlet boundary conditions with equidistant distribution of points in all boundaries 53 4-4 Convergence characteristics for exponential and equidistant distribution of points. 55 4-5 16 x 16 grid, H — 0.25, ISH A — 2. Neumann boundary condition on left side. 56 4-6 16 x 16 grid, H — 0.15, ISH A — 1. Neumann boundary condition on left side. 57 4-7 Convergence characteristics for Figure 4.6 with H=0.1 58 4-8 16 x 16 grid, H — 0.15, ISH A = 2. Linear isoparametric elements 60 4-9 16 x 16 grid, H — 0.15, ISH A = 2. Quadratic isoparametric elements 60 4-10 Convergence characteristics for Figures 4.8 and 4.9 61 6-1 Typical mesh for surface-piercing cylinder 100 6-2 Linear free surface time record for surface-piercing cylinder (NI = 30, At — 0.1). 100 6-3 Linear free surface time record for surface-piercing cylinder (NI = 30, At = 0.1). 101 6-4 Linear free surface time record for surface-piercing cylinder (NI — 20, At = 0.1). 101 6-5 Linear free surface time record for surface-piercing cylinder (NI — 20, At = 0.1). 102 6-6 Hydrodynamic force due to sinusoidal excitation (NI = 20, At = 0.1) 102 6-7 Hydrodynamic force due to sinusoidal excitation (NI — 10, At = 0.1) 103 6-8 Energy (—) and work (- - -) (NI - 30, At = 0.1) 103 6-9 Energy (—) and work (- - -) (NI - 20, At = 0.1) 104 6-10 Energy error for linear elements 104 6-11 Energy error for quadratic elements 105 6-12 R M S of energy error versus At r 105 6-13 R M S of energy error versus number of nodes 106 6-14 Number of iterations of C G . method versus m (as in Section 5.5) 106 6-15 Number of iterations of C G . method versus m (as in Section 5.2.1) 107 6-16 Number of iterations of C G . method versus time-step size (as in Section 5.5). . 107 vi i 6-17 Number of iterations of C . G . method versus time-step size (as in Section 5.2.1). 108 6-18 Typical mesh for submerged cylinder 108 6-19 Linear free surface for submerged cylinder (NI = 20, A t — 0.1) 109 6-20 Linear free surface for submerged cylinder (NI — 20, A t = 0.1) 109 6-21 Hydrodynamic force due to sinusoidal excitation (NI — 20, A t = 0.1) 110 6-22 Energy (—) and work (- - -) (NI - 30, A t = 0.1) 110 6-23 Energy error for submerged cylinder I l l 6-24 Plan view of free surface particles for predictor-corrector methods I l l 6-25 Plan view of free surface particles for semi-implicit semi-Lagrangian method. . 112 6-26 Linear free surface evolution up to 20 cycles, NI=10, At=0.1 112 6-27 Nonlinear free surface evolution up to 20 cycles, NI=10, At=0.1 113 6-28 Mesh deformation at 8.25, 8.5, 8.75 and 9th cycles 114 6-29 Velocity vector plots at 8.25, 8.5, 8.75 and 9th cycles 115 6-30 Hydrodynamic forces for surface-piercing cylinder 116 6-31 Hydrodynamic forces for surface-piercing cylinder versus tank dimensions. . . . 116 6-32 Work (—) and total energy (- -) records, NI=10, At=0.1 , A=0.5, u =1.25. . . 117 6-33 Energy error versus time-step and number of grid divisions 117 6-34 R M S of energy error versus time-step (NI=10) and number of nodes (At=0.05). 118 6-35 Nonlinear free surface time record for submerged cylinder (NI = 10, A t = 0.05). 119 6-36 Nonlinear free surface time record for submerged cylinder (NI — 10, A t = 0.05). 119 6-37 Mesh deformation at 9.25, 9.5, 9.75 and 10th cycles 120 6-38 Velocity vector plots at 9.25, 9.5, 9.75 and 10th cycles 121 6-39 Hydrodynamic force due to sinusoidal excitation, A=1.5 122 6-40 Work (—) and total energy (- -) records, NI=10, At=0.1, A=0.5, u =1.25. . . 122 6-41 Energy error versus time-step and number of grid divisions 123 6-42 Free surface frequency spectra due to impulsive motion for square basin. . . . 124 6-43 Free surface frequency spectra due to impulsive motion for circular basin. . . . 125 6-44 Free surface frequency spectra due to impulsive motion for parabolic basin. . . 126 6-45 Free surface frequency spectra due to impulsive motion for triangular basin. . . 127 6-46 Free surface frequency spectra due to impulsive motion for elliptic basin. . . . 128 6-47 Linear free surface elevation for Wigley hull with forward speed, Fn — 0.267. . 129 6-48 Nonlinear free surface elevation for Wigley hull with forward speed, Fn = 0.267. 130 6-49 Wave profile at side of Wigley hull with forward speed, Fn — 0.267 131 6-50 Linear free surface elevation contours of Wigley hull with forward speed, vi i i Fn = 0.267 131 6-51 Nonlineax free surface elevation contours of Wigley hull with forward speed, Fn = 0.267 132 6-52 Wave resistance coefficient of Wigley hull versus Fn 132 6-53 Wave resistance coefficient of Wigley hull versus Fn 133 6-54 Sinkage of Wigley hull versus Fn 133 6-55 Tr im of Wigley hull versus Fn 134 1 Tangent to a coordinate line 151 2 Tangent to a coordinate line and normal to a coordinate surface 153 ix LIST OF SYMBOLS o* Contravariant base vectors, i— 1, 2, 3 o,i Covariant base vectors, i — 1, 2, 3 B(x,y) Body offset B (x, y) Nondimensional body offset Cb Block coefficient Cp Prismatic coefficient Cw Wave resistance coefficient Cwp Waterplane coefficient D Ship draft / Distortion function for mesh generation F(x, y, z) Function describing body surface Fn Length Froude number, Fn — —j= FP Draft Froude number, FP - -¥= Fv Force in y-direction g Acceleration of gravity g%j Contravariant components of metric tensor gij Covariant components of metric tensor II Longitudinal moment of inertia J Jacobian of the transformation from (x,y) to (£,77) coordinates J* Jacobian of the transformation from (e, q) to (£, 77) coordinates L Ship length Mz Moment about z-axis n Unit normal vector rii, 7i2, Ti3 Components of unit normal vector in x,y,z directions respectively Components of unit normal vector in x, y, z directions respectively qi Components of body velocity vector, i — 1,.., 6 aijl2}l3 Components of translational body velocity vector in x,y,z directions respectively 24>25)96 Components of rotational body velocity vector about x,y,z directions respectively x R Residual of partial differential equation Rw Wave resistance __F s Non-dimensional sinkage = Aw*g_ S Ship hull wetted surface at design waterline t Time t Non-dimensional tr im = T Wave encounter period Tw Wave period U Ship forward speed u, v Velocity components in z, y directions respectively v Position vector Vw Wave speed w Distance between body vertical side and tank wall Wi Weighting functions V(x, y, z) F lu id velocity vector x,y, z Cartesian coordinates x, y, z Nondimensional Cartesian coordinates x, y Approximation functions for mesh generation of x, y coordinates yi, Heave motion displacement function for surface piercing case Heave motion velocity function for surface piercing case Zb Heave motion displacement function for submerged case Zb Heave motion velocity function for submerged case a Parameter for modified Lax method 3mk Coefficients of predictor-corrector methods e Slender body parameter £, 77 Isoparametric coordinates e, g Boundary fitted curvilinear coordinates A Wave length A Function describing free surface elevation T Boundary of domain of integration Q w Body motion frequency Q Domain of integration (j> Velocity potential Nondimensional velocity potential xi Body velocity potential <f>b Nondimensional body velocity potential Wave velocity potential <f>VJ Nondimensional wave velocity potential $ Approximation function for velocity potential Linearly independent functions c(x,z,t) Free surface elevation Unsteady Bernoulli constant Nondimensional free surface elevation p Water density x i i A C K N O W L E D G E M E N T A number of individuals have been helpful during the course of this work. I extend my sincere gratitude to them: to Sander Calisal, my supervisor, for allowing me the freedom to search for answers to the problems addressed in this thesis; to Avrom Soudack from Electrical Engineering for his lectures in nonlinear systems as well as his encouragement, support and good nature throughout my P h D studies at U B C ; to Ricardo Foschi and Mervyn Olson from C i v i l Engineering for their teaching of the Finite Element Method and for their comments to improve the text of the thesis; to Bruce Bowen from Chemical Engineering for providing a key reference for the mesh generation scheme used in this thesis; and to Ian Gartshore whose affable nature made our discussions more enjoyable throughout my graduate studies at U B C . Assistance provided by the staff of the Department of Mechanical Engineering and by Tom Nicol and Freddy Echeverry from U C S U N I X G system is greatly appreciate it . Mention must be made of Rodolfo Bermejo, now at the Department of Applied Math-ematics of the Complutense University of Madrid, Spain, for his insightful suggestions and indispensable support in the implementation of some of the concepts presented in this disser-tation. I also thank Chung-Fa W u , from Brown & Root, U S A , for his assistance during our phone conversations. Deep appreciation is also given to Annette Johnston for her patience, support and proof-reading of the text. This thesis is dedicated to my parents Julia and Angel, and to all those who strive to make of this world a better place to live. From them I have learned my most important lessons. This study was partially supported by grants from the Natural Sciences and Engineering Research Council of Canada held by Sander Calisal. Much needed financial support in the form of teaching assistantships from the Faculty of Applied Science and the Department of Mechanical Engineering is also greatly appreciated. xm ONE INTRODUCTION In ship design, it is important to determine the dynamic behavior of a vessel at sea as a function of its form and dimensions. From this information, the designer is able • to select an optimal hull form that wil l satisfy predefined criteria such as type of cargo, speed, displacement and route of operation • to predict the behavior of the chosen configuration and to ascertain the consequences of this behavior on various effects such as safety, structural dimensioning, resistance and propulsion. A t the preliminary design stage, many ship parameters remain to be varied. Principal dimensions, body plan, and various hull coefficients are some of the variables involved in the primary selection of a ship hull . Testing of this initially large number of concepts is more easily executed by mathematical simulations than by physical model tests. Test facilities and costs may also be restrictive factors in the use of the latter models. In order to have reliability of the numerical simulation results, all possible phenomena that are related to the type of problem to be solved should be incorporated in the development of the computer code. Some of these phenomena are •Added mass and added inertia due to surrounding water •Hydrodynamic damping •Restoring forces •Inertia of the ship •Loading imparted by waves 1 •Steering mechanism The variables related to these phenomena can be considered as part of a three di-mensional nonlinear problem. This requires the solution of very large systems of algebraic equations at each of a large number of time-steps. In this case, excessive memory and C P U requirements may pose restrictions on the viability of the calculations. In order to make this problem more tractable in a numerical sense, certain simplifications in ship dimensions may be introduced to reduce the computational demands required for the solution of the problem. In this work, simplifying steps have been taken in l imiting the ratio of the transverse to lon-gitudinal physical dimensions of the ship. This is known as a slender body approximation[l]. 1.1 SLENDER BODY THEORY A marine vessel can be considered a slender body. A slender body is an elongated body wherein any tranverse direction is small compared to its length. The velocities at locations along the longitudinal body axis are extremely important. Certain valid approximations are then made near this axis so that the boundary conditions may be satisfied on the actual boundary of the body. The actual three-dimensional flow problem is then replaced by a series of two dimensional flow problems in transverse cross-flow planes. The original three dimensional equation governing potential flow is then reduced to the two-dimensional Laplace equation. For an elongated body, the underlying physical idea is that flow variations in the lon-gitudinal direction are much smaller than those occurring in other directions in the local vicinity of the body. One of the earliest attempts t© use this simplification was made by Munk[2] when calculating the transverse force on a rigid airship hull turning at angles of yaw and pitch in the absence of drag. Munk claimed that air gave way to the passing ship by flowing around the axis of the ship, not by flowing along it. He considered that the flow in ver-tical planes at right angles to the motion remained in that plane. The lateral force was then 2 calculated in two-dimensions as the rate of change of momentum M — M(U2, p, aB ,B„, jj^), where as and Bs are the pitch and yaw angles and S is the function defining the ship's sur-face. Integration along the length of the airship gave rise to an expression with the following three terms: a) a moment of a ship flying straight at a pitch angle b) a distributed transverse force c) transverse forces concentrated at the ends of the hull . Summation of b) and c) along the length of the airship gave a zero force. This left the airship with an unstable moment that was overcome by using stabilizing tai l fins. This work was extended to wings of low-aspect ratio by Jones[3]. More detailed mathematical treatments were presented by Ward[4] and Adams and Sears[5] for supersonic and subsonic flow respectively. The above-mentioned works were concerned with steady-flow problems. For unsteady flow, fundamental developments were performed by Miles[6,7] using a systematic perturba-tion expansion procedure. A consequence of using slender body theory is that the cross-flow at each cross-section is independent of the flow downstream of that section. The fluid flow at each cross-section is only determined by information conveyed from upstream conditions. The connection of these series of two dimensional flows with the actual three-dimensional flow is achieved through the body boundary condition that contains the longitudinal x-variable. The velocity potential is then found as (f> — <f>2D + /(a3)- 0 n e of the first attempts to use slender body theory for surface ships was made by Cummins[8] in the solution of the wave resistance problem. Later, Kaplan[9] studied the vertical force and pitching moment on a slender body moving normal to the crests of regular waves. 3 As a result of these investigations, a considerable number of researchers made use of slender body theory in works related to the marine environment. Lighthill[10] studied the thrust efficiency of a slender fish and concluded that this parameter would depend strongly on the wave like motion of the fish and its propagation speed along the spinal cord. Vossersfll] and Maruo[12] studied wave resistance of a ship with uniform forward speed. Ursell[13] approximated the flow around a stationary ship by an axial line distribution of known point singularities and the harmonic small motions by wave sources satisfying the free surface condition. Tuck[14] extended this work to the steady translation of a slender ship. Inner and outer solutions were used to give a systematic treatment of the linear boundary conditions by successive terms in the expansions. These expansions expressed the velocity potential ^ as a summation of <bn terms. The inner expansion was made to satisfy the boundary condition on the body, and the outer one was merely required to match the inner one at an intermediate region close to the body surface. Thus, the inner expansion had no outer boundary condition and the outer expansion had no inner one. Newman and Tuck[15] described a complete systematic linear theory of the motions of a slender ship in a seaway. Their work was divided into three parts that contained calculations of pitch and heave response at zero speed, a zero-speed theory for harmonic oscillations in oblique incident waves, and a general forward speed theory for harmonic oscillations in calm water. Newman[16] generalized the Haskind relations[17] for the six exciting forces and moments in waves to include the effect of forward speed as a function of the wave length, heading angle and forward speed. Tuck[18] solved the problem of the disturbance produced to a shallow water stream by an immersed slender body. This investigation was directed towards studying the steady motion of ships in restricted waters using the method of matched asymptotic expansions. This technique was then compiled in a textbook by Van Dyke[l]. Throughout these analyses, the fluid was considered inviscid and incompressible. Vis-cous effects cause separation of the flow and the existence of a boundary layer around the body. The former is present when analyzing the lateral motions of a body in a fluid. This 4 is part of the so-called 'end-effect'problem and is generally not considered in the literature. Newman[19] treated the lateral flow of an ideal fluid past a slender body when the flow was constrained by a pair of closely spaced walls parallel to the body axis. In the absence of walls, the flow field was nearly two-dimensional in the cross-flow plane normal to the body axis. The presence of the walls introduced effective blockage in the cross-flow plane that caused the flow field to become three-dimensional. Vorticity shed by the fins attached to a submerged body was introduced later by Newman and Wu[20]. Much of the progress in slender body theory as applied to submerged and floating marine vehicles was carried out in the 1960's. Newman[21] summarized the advances made during this period. Research until then had focused on analytical techniques for steady state problems. Numerical results related to slender-ship motions were then published. These results compared analytical and numerical calculations within the realm of linear boundary conditions compatible with small departures of the body from an equilibrium position and with infinitesimal deviations of the free surface from undisturbed conditions. Some of these works are by Chapman[22] and Jensen and Pedersen[23]. Troesch[24,25] obtained sway, roll and yaw motion coefficients for a slender ship with uniform forward speed. A simplified boundary value problem was later introduced by Maruo[26] in a new slender body approach to calculate resistance of a ship with uniform forward speed. Other slender-ship approxima-tions were presented by Noblesse[27] for the calculation of linear wave resistance problems using a new integro-differential equation for the velocity potential of the flow caused by the ship. Faltinsen[28] studied the flow around a ship bow at high Froude numbers and regular incident head sea waves with complete linear boundary conditions. A technique based on the distribution of two dimensional sources and dipoles on the body and the free surface was used. Chen and Noblesse[29,30] presented numerical results for wave resistance at a rela-tively wide range of Froude numbers for the theory developed by Noblesse in[27]. Diffraction of free surface waves using linear theory was presented by Sclavounos[31] using matched asymptotic expansions to match a two dimensional inner field with the three dimensional 5 so lut ion of the far f ie ld . Sclavounos[32] extended this work to address l inear di f f ract ion and r a d i a t i o n problems for heave and p i t c h motions of a ship w i t h forward speed i n waves. B r e i t and Sclavounos[33] used a l inear a p p r o x i m a t i o n for surface wave r a d i a t i o n by two adjacent slender bodies. M a r u o ' s approach[26] was numer ica l ly s tudied by Song et. al.[34] i n the ca lcu la t ion of wave patterns a n d wave resistance of a W i g l e y h u l l . Var ious c o m p u t a t i o n a l procedures m a k i n g use of slender b o d y approximat ions w i t h different degrees of nonl ineari t ies i n c l u d e d i n their schemes were also tackled by C h o i and Mei[35], Kashiwagi[36] , K a s h i w a g i and Ohkusu[37] , C a l i s a l and Chan[38] and Ando[39] . Further progress of the theory pre-sented b y Maruo[26] were later reported by this author[40]. Extens ive exper imenta l results for the W i g l e y h u l l were made available by Shearer[41] and N a m i m a t s u et. al.[42]. A p r e l i m i n a r y at tempt to use slender b o d y theory i n large ship motions at f o r w a r d speed was m a d e by Loeser, Y u e and Salvesen[43] w i t h a l inear free surface condi t ion and exact b o d y b o u n d a r y condi t ion . L inear boundary conditions are consistent w i t h a s m a l l departure of the free surface f r o m its e q u i l i b r i u m posi t ion . If a b o d y is made to undergo s m a l l osci l lat ions about a fixed e q u i l i b r i u m pos i t ion , then the smallness of the dis turbance at the free surface is ensured by the smallness of the osci l lat ion alone and the b o d y need not necessarily be s m a l l . Therefore, i f the oscil lations are of smal l a m p l i t u d e , the l inear ized free surface b o u n d a r y condit ions m a y be assumed as adequate and slenderness is not required[13]. It was c l a i m e d b y Tuck[14] that i f the f loating b o d y has a finite ve loc i ty of t rans la t ion slenderness is necessary for l inear izat ion . Slenderness is then used as a means of enabl ing the designer to s i m p l i f y this compl icated three dimensional prob lem and to obta in results that can be taken as acceptable f r o m an engineering perspective. 1.2 A B R I E F R E V I E W O F F R E E S U R F A C E F L O W A N A L Y S I S A s previously ment ioned, slender body theory assumes the va l id i ty of a two-dimensional cross-sectional flow that is carried along the remaining t h i r d dimension through appropriate boundary condit ions. In this case i t seems advisable to i n i t i a l l y develop a robust m o d e l that can describe two-dimensional flow with a nonlinear free surface in a reliable and accurate form. As a consequence of intensified research in relatively recent years, the importance of two-dimensional nonlinear free surface hydrodynamic problems has become apparent. A n early compilation of some of the methodology used to solve these types of problems during preliminary years has been outlined by Yeung[44]. The pioneer work by Longuet-Higgins and Cokelet[45] used an integral-equation formulation to simulate large amplitude steep waves and wave breaking. Later investigations by Vinje and Brevig[46] and L i n , Newman and Yue[47] added the presence of a moving body. A number of contributions in the area of integral-equation methods were also made by Isaacson[48], Baker, Meiron and Orzag[49], Dold and Peregrine[50] and Suzuki[51]. Sen et.al.[52] studied the propagation of steep two dimensional periodic waves and the large motions induced by these waves on free floating bodies. Cointe[53] studied the nonlinear forced motion of surface piercing bodies and the nonlinear motion of a rectangular barge in beam seas[54]. Finite difference techniques have also been employed in nonlinear wave-body problems by numerically mapping the physical domain onto a regular computational domain. Ghia, Shin and Ghia[55] and Haussling[56] used this approach to study the formation of breaking waves. Telste[57] and Yeung and Wu[58] solved nonlinear free surface problems of bodies moving with forced periodic motion in a tank. Other research in closed domains concentrates either on the response due to slosh-ing or the response due to impulsive motion. Sloshing phenomena have been studied by Faltinsen[59], Arai[60], Bridges[61], Yeung and Wu[62] and Kamiya and Miyazawa[63]. Washizu and Nakayama[64], Washizu[65] and Washizu et. al.[66] used a number of finite element and boundary element techniques to solve various nonlinear wave sloshing problems. Agreement between F E M and B E M results was observed. The response due to impulsive loading and natural frequencies for tanks of different shapes has been reported by Mclver and Smith[67], Chang and Wu[68] and Yeung and Wu[62]. 7 Finite element techniques have also been used in the context of free surface flow. A variational principle for the exact nonlinear free surface flow was introduced by Luke[69] but no numerical results were given. A modified form of Luke's formulation was used by Whitham[70,71] to study wave dispersion. Bai and Yeung[72] studied linear forced motion problems using the conventional variational formulation and a more efficient modified varia-tional method in order to truncate the original infinite domain of solution. Known analytic solutions in terms of eigenfunctions were used in the outer region of the domain, the co-efficients of which were determined as part of the solution. Finite elements were used to discretize the inner region giving rise to the so called localized finite element method. In ad-dition, the velocity potential and normal velocity of both regions were made to match at an imaginary juncture line connecting both subdomains. The first procedure presented by Bai and Yeung[72] was used by Bai[73] to study diffraction of oblique waves by an infinite cylinder with linear boundary conditions. From this time onwards, a number of hybrid formulations were used to circumvent the problem posed by finding a proper truncation to the infinite domain boundary. Yim[74] addressed the problem of nonlinear steady ship waves using an approach with linear theory in the outer region. Bai[75,76] calculated flow about two di-mensional bodies under a linear free surface using a weak form and variational formulations. Bai[77] also derived an approximate formula for blockage effects affecting the flow past ship models moving in a towing tank with a free surface. Wellford and Ganaba[78] implemented a pseudo-variational approach for unsteady water wave problems. Further contributions using finite element techniques in ship hydrodynamics free surface potential flow were presented at the 3rd International Conference on Numerical Ship Hydrodynamics. Oomen[79] used perturbation analysis of the velocity potential and free surface elevation. This technique produced a series of linearized two dimensional problems that were used to study the steady translation of a Series 60 model in calm water. Jami[80] solved the linearized flow past a submerged body with arbitrary transient motion of small amplitude. Numerical results were presented for two-dimensional transient forces acting on the body. Euvard et. al.[81] coupled finite element and singularity distribution procedures. They concluded that the lo-calized finite element method was the most efficient technique in two-dimensional problems when water depth was not too large. For infinite depth or for three-dimensional problems, they recommended coupling finite element and boundary element techniques. This recom-mendation was based on results obtained for wave-resistance calculations of a submerged cylinder. One of the first attempts to use three-dimensional analysis of free surface potential flow in ship hydrodynamics was made by Bai[82]. Linear three-dimensional steady flow problems were solved. Bai[83] later applied the localized finite element method to the solution of radiation and diffraction motion problems in three dimensions with linear boundary condi-tions. B a i , K i m and Kim[84] studied the steady translation of a ship i n a towing tank. The three-dimensional nonlinear free surface problem was solved using a variational formulation. Computations were conducted for finite element generation made in such a way that the projections on the horizontal plane remained constant, thus simplifying mesh generation and the computations significantly. Choi, Bai , K i m and Cho[85] solved the nonlinear three-dimensional free surface flow problem of a ship moving in shallow water. In this work, the original problem was replaced with an equivalent variational problem based on Hamilton's principle as derived by Miles[86]. The variational functional used in this case was slightly different to that defined by Luke[69]. The formal difference resided in the introduction of the velocity potential defined only on the free surface. This resulted in a more advantageous functional to treat the nonlinear free surface boundary conditions. Luke's formulation[69] had the additional requirement that the velocity potential be known at initial and final times. Considerable effort has been dedicated in producing results which include viscous effects in the calculation of the flow around a ship. As a first step, boundary layer theory in the immediate neighbourhood of the hull has been coupled with potential flow theory outside of the viscous layer. This approach has given satisfactory predictions of viscous and wave 9 resistance leading to a more self-reliant tool for ship design purposes. Pure potential flow theory can only predict wave-making resistance. Consideration of the free surface boundary layer and surface tension effects has lead to improvement in the calculation of wave-making resistance[87]. F l u i d flow past an advancing marine vessel is three-dimensional and viscous. Numerical simulation of this problem requires the solution of three coupled nonlinear partial differential equations and a continuity equation in time. This difficult problem is further complicated by having to impose boundary conditions on a free surface whose position is unknown. Present-day high speed computers make it possible to solve Navier-Stokes equations with a free surface. A number of contributions in this area were presented at the 5th International Conference on Numerical Ship Hydrodynamics. Masuko and Ogiwara[88] studied viscous flow around a cargo Series 60 model and tanker models using Reynolds-averaged Navier-Stokes equations for three-dimensional flow. For their finite-difference discretization of the Series 60 model, 50000 grid points were reported to be required to obtain converged solutions after about 300 iterations of the numerical procedure. Zhu, Miyata and Kajitani[89] also used a finite difference solution method to solve three-dimensional viscous flow past a Wigley hull . A fine grid system of 340000 points and a time increment of At = 10~ 5 were used. Hino[90] carried out computations using a finite-difference representation of the three-dimensional Navier-Stokes equations with nonlinear boundary conditions to simulate the flow past an advancing Wigley hull . In order to limit the appearance of strong wiggles of pressure on and around the hull , about 80000 nodes had to be used. Stabilization of the solution was then achieved after 11000 time-steps. A recent attempt to solve three-dimensional viscous flow with free-surface waves was made by Miyata, Zhu and Watanabe[91]. A finite volume method was developed producing very good agreement between measured and simulated wave profiles. It was stated in this publication that about 100000 grid points seemed to be the minimum required for resolution of the physical phenomena of interest. 10 It is clear that the enormous c o m p u t a t i o n a l demands required to solve these types of problems involve computer systems that m a y be unavailable i n most naval architecture design offices. F u r t h e r m o r e , efforts should be devoted to devising n u m e r i c a l techniques capable of o b t a i n i n g accurate, stable a n d convergent results w i t h larger g r i d sizes a n d t ime-steps. Cons iderat ion to this p r o b l e m , is given i n the next sub-section. 1.3 O B J E C T I V E S O F THIS I N V E S T I G A T I O N I n v iew of the above i n f o r m a t i o n , consideration has been given to produce a s impl i f ied m o d e l that can re l iab ly predict quantities required i n ship design w i t h affordable computer hardware. Viscous effects have been neglected based on the assumption that these are concentrated inside a t h i n b o u n d a r y layer surrounding the vessel. T h e b u l k of the fluid is then considered i n v i s c i d . A n objective of this work was to estimate f lu id flow w i t h f u l l y nonl inear b o u n d a r y condit ions around a slender ship. Therefore, i n a d d i t i o n to the previous assumption, a slender b o d y a p p r o x i m a t i o n of the three-dimensional potent ia l flow p r o b l e m has been made. In neglecting the three dimensional i ty of the flow, certain in format ion is lost. T h i s i n f o r m a t i o n is p r i m a r i l y associated w i t h the generation and propagation of the transverse wave system. For re lat ively h i g h Froude numbers i n a closed d o m a i n , the f o r m a t i o n of solitons ahead of the ship after a certain per iod of t i m e is another physica l phenomena that slender b o d y approximat ions cannot account for. It seems that these are the most significant l i m i t a t i o n s of the m e t h o d described i n this work w i t h respect to the three-dimensional fu l ly nonlinear p r o b l e m . Slender b o d y calculations require the solution of a series of two-dimensional problems along the length of the vessel. In order to obta in accurate, stable and convergent solutions, efficient a lgorithms have to be used to integrate the governing equations i n t i m e and space. 11 In this work, a Bubnov-Galerkin formulation of the boundary value problem describing two-dimensional potential flow has been developed for the spatial estimation of the velocity potential. It was foreseen that this formulation would require an efficient mesh generation system capable of adapting to the changes of the nonlinear moving boundaries. Again, a Bubnov-Galerkin formulation of the grid generation equations was performed. Significant computational advantages over previous solutions were realized[92,93]. Linear and quadratic isoparametric elements were used and their relative merits analyzed. The time derivatives describing the nonlinear boundary conditions were initially dis-cretized using various orders of explicit-implicit predictor-corrector methods. These methods exhibited excellent stability characteristics and improved accuracy for linear free surface ap-plications. However, for nonlinear free surface flow calculations these procedures displayed unstable behaviour. Numerical dissipation of the velocity potential and the free surface el-evation then had to be introduced. The resulting motion of the free surface fluid particles exhibited inaccuracies that negatively affected the calculation of quantities of engineering interest. In order to avoid erroneous estimations due to smoothing, efforts were directed towards devising a numerical scheme that would be stable with minimum numerical dissipa-tion. A n efficient two-time level semi-implicit semi-Lagrangian iterative integration scheme was then conceived. The method is almost free from so called smoothing. Iterative schemes of a similar nature were used by Robert [94,95] in problems associated with the integration of the primitive meteorological equations. Larger time-steps and higher stability character-istics with better accuracy were reported. Staniforth and Temperton [96,97] applied this technique to a finite-element barotropic model using the shallow water equations. Recently, Bermejo[98] reinterpreted the Galerkin-characteristicmethod to devise a computationally efficient scheme for the scalar advection equation within the framework of particle methods. Robert's iterative scheme was used to advance the solution in time, resulting in a convergent and unconditionally stable scheme. These works were concerned with meshes whose bound-aries remained unchanged as the solution progressed in time. In this work, the iterative 12 semi-implicit semi-Lagrangian scheme is used to locate the fluid particles in contact with the moving boundaries. In addition, values of the dependent variable are obtained on these material points. In order to validate the methodologies previously mentioned, a number of two-dimensio-nal problems have been solved. Results are given for fluid response due to forced oscillations of surface piercing and submerged bodies in a closed domain. F lu id response in tanks of different shapes due to impulsive loading is also studied. The two dimensional results are subsequently extended to study a mathematically defined slender body, the Wigley hull , advancing with forward speed in otherwise calm water. 13 C H A P T E R TWO GOVERNING EQUATIONS In this chapter we introduce the equations that govern three dimensional potential flow with a nonlinear free surface, section 2.1. Free surface and body boundary conditions are derived in terms of Eulerian derivatives. In section 2.2, an order of magnitude analysis is performed in order to reduce the general three dimensional problem to a form that is more amenable to numerical analysis. This is performed under the assumptions given by slender body theory. In section 2.3, we re-write the free surface boundary conditions in a Lagrangian form so as to trace the motion of the fluid particles. 2.1 P O T E N T I A L F L O W W I T H A N O N L I N E A R F R E E S U R F A C E We assume that the fluid is ideal (inviscid, incompressible and homogeneous), the flow motion is irrotational, and the only significant body force acting on the fluid is gravitational. A Cartesian coordinate system (x,y,z) is adopted. The plane of the undisturbed free surface is defined at y=0 with gravity acting in the negative y-direction. The fluid velocity V(x, y, z, t), where t denotes time, can be represented by the gradient of a scalar potential <b{x,y,z,i) that must satisfy Laplace's equation „ d2cf> d2cf> d2(f> V 2 < £ = — H T = 0 2.1 dx2 dy2 dz2 throughout the fluid domain. The above equation has to be solved subject to conditions applied at a number of boundaries. One of these boundaries is the free surface. The physical nature of a free surface requires that both a kinematic and a dynamic boundary condition be satisfied. The 14 kinematic boundary condition states that the normal velocity of the fluid at the free surface and that of the boundary surface must be equal. The dynamic boundary condition states that the pressure on the free surface must be constant and equal to the atmospheric pressure. These two conditions are satisfied on the free surface that may be described by a function A = y-q(z,z,t) = 0 2.2 The exact kinematic boundary condition may be derived by stating that the substantial derivative of y — c vanish on the free surface. In other words, fluid particles that belong to the free surface are assumed to remain on the free surface at all times. The differential operator that describes the time progression of a measurable quantity corresponding to a fluid particle can be expressed as + + + 23 dt dt dx dx dy dy dz dz Application of (2.3) to y — c(x,z,t) then gives d, x , 9 dtp d d(f> d d<j> d , n or d<j> dc dc/> 9? ^ 9 ? which is the exact kinematic boundary condition in terms of Eulerian derivatives. The second boundary condition to be applied on the free surface is Bernoulli's equation. This results in the so called dynamic boundary condition. For unsteady irrotational flow we can then write ^ + ^V<£V<A = - - ( p + pgy) + C{t) 2.6 ot 1 p 15 The term C(t) is independent of the space coordinates but it may depend on time. In general, it can be absorbed into the velocity potential since the difference between (p and (j> + C(t) is a function of time only. Thus the gradient is not affected and either function is an equally suitable velocity potential. It wil l be seen in later sections that this constant cannot be absorbed so freely into the velocity potential for certain cases. However, presently we write —(? - Patm) - ~kt + o V<£V<£ + gy - 0 ony -c(x,z,t) 2.7 p eft l which is the dynamic boundary condition. It states that the pressure is constant over the whole free surface, thus being independent of the spatial position of the surface and equal to the atmospheric pressure. Finally, a body boundary condition is applied at the surface of the moving body. This surface can be expressed as F(x,y,z) = z-B(x,y)-0 2.8 In the general case of an impermeable body, the boundary condition on its surface is dF — = riiqi 2.9 at where represents again the substantial differential operator. As before, this operator is applied to 2.8 resulting in deb dd> dB dcpdB — = mqi 2.10 dz dx dx dy dy where 16 ( « i i»2 , n 3 ) = 2.11 and WF ( n 4 , n B , n 6 ) = r x 2.12 (gi, g2, gs) = translational velocity along i direction, i = 1, 2, 3 (g4,gs,g6) = angular velocity about i direction 2.2 O R D E R O F M A G N I T U D E A N A L Y S I S The velocity potential function used to approximate the flow around the ship's surface can be written as the superposition of the following terms <b = Ux + (f>i + K 2.13 where the contributing factors are Ux: due to ship forward speed, cf>b'- due to ship motions, (bw: due to waves. The velocity potential given by 2.13 can be non-dimensionalized as *=4d 2 - 1 4 where the bar denotes non-dimensional quantities and D is the ship draft. The wave surface c and the body offset B(x,y) can be written as S = s{ztz) 2.15 17 9 z = B(x,y) 2.16 Similarly, in non-dimensional form ?(*,*) = 2.17 B = ?^- 2.18 D We then express our Cartesian coordinates as follows: _ x _ y _ z X = L V = D *=D 2.19 where L is the length of the ship. The dimensional derivatives can then be non-dimensionalized as d 1 d dx L dx a 1 d dy Ddy d _ 1 d Yz ~ d _ 1 d dt ~ Tdt The nondimensional form of the Laplace equation then becomes ,D,9d2A d2d> d2d> 18 A The ratio ^ has an order of magnitude given by the relative dimensions of the ship. It represents an overall estimate of the transverse to longitudinal dimensions of the vessel. For low values of jj, generally less than y, it is clear that the flow variations produced by a slender body placed in a uniform stream in the longitudinal x-direction, represented by wil l be of lower magnitude than that experienced in the remaining y and z directions. This condition wi l l inevitably lead to values of which are of considerable less importance than 2— 2— p ^ - and In addition to this, the higher order term ( ^ ) 2 in 2.20 is further indication that the first term in this equation can be neglected. Based on these considerations, we can write d2<f> dhf> dy2 + dz which is the Laplace equation in two dimensions. It should be noted that for a slender ship moving across incoming waves, both parameters y and should be considered. For wavelenghts A > L the above considerations are still valid. However, for y = 0(1) the first term in 2.20 is no longer a higher order term. It then seems appropiate to speculate that the application of slender body theory is restricted to wavelengths of similar length or longer than the ship length. We now introduce our potential (p given by 2.13 into our body boundary condition d<t>b | d<t>w f J J | dfo | dcpw dB ,dcf>b | d<f>w dB _ Q _ ^ dz dz dx dx dx dy dy dy which in non-dimensional form becomes _h + ?hL-e?l-e2f^!^ _ e2d--Ld£ _ ^dJL _ ______ - o 2.23 dz dz dx dx dx dx dx dy dy dy dy 19 ! I If we now keep terms of 0(1) and 0(e) we have dh 0 ^ _ dfadB _ d^dB_ _ dB = 0z 0z cty 0y 0y 01/ Sac which is the body boundary condition of highest magnitude in e The kinematic boundary condition was given in 2.5 as ^ _ 0 £ _ ^ £ _ ^ 1 = O 2 25 dy dt dx dx dz dz or in non-dimensional form 0 ? 3 _ _ j ^ 0 | _ 2 f ¥ # r _ < t y ^ = 0 2 2 6 dy TU dt € dxdx dz dz The period T corresponds to the period of wave encounter. Neglecting terms of 0(e 2 ) gives the kinematic boundary condition of highest order in e as dj> D 0? d<b 0? — = - — = 0 2.27 dy TU dt dzdz We can now make use of the nondimensional derivatives to determine again the nondi-mensional form of equation 2.25. In this case, we introduce 2.13 into 2.26 before dropping higher order terms. This gives ^ + A ^ _ e 3 f e - 1 + ^ + ^ ) — - f " ^ + ^ ) — = 0 2 28 dy dy TU dt K dx dx 1 dx K dz dz ' dz Neglecting terms of 0(e 2 ) gives 20 ^ i + ^ L _ ^ _ ^ _ e ^ i _ ^ l ^ l _ ^ ^ l - 0 2 29 dy dy TU dt dx dz dz dz dz which is the kinematic boundary condition of highest magnitude in e. In this case, unlike in 2.27 the wave slope in the x-direction is kept i n the analysis. The component of the ship velocity in the direction of the wave motion is U cos /j,. The relative speed of the vessel with respect to the waves is then Vw — U cos p. As the wavelength remains constant, the time of travel for the ship from one wave crest to the next is T = = P- 2.30 Vw — U cos fi 1 — ^ - cos p T = 2.31 A — Tw U cos fi The linear dispersion relationship for deep water waves can be stated as[99] T_ = ^ 2.32 9 which determines the period of wave encounter as g\f\ — U cos u.s/2ng For a given ship speed and heading angle of astern seas^  we have: 2.33 T — ^ 2.34 C 2 V A + C 3 which when introduced into 2.26 gives 21 d<b DC2V* + C3ds 2d<f>ds df_8f_ 2 3 5 e dy U C\\ dt ~ dx dx dz dz where C% are constants. From this last equation, the second term becomes of order e2 for wavelengths comparable to L 4 . This is equivalent to a calm water case and only then the term can be dropped. Therefore, for ship wave applications this term should not be neglected. We can now continue with the dynamic boundary condition written as 9 ? + ^ + iv<£V<£ = 0 2.36 or in non-dimensional form Keeping terms of O ( l ) and 0(e) we have For the same reasons stated above regarding 2.35, ^ should be kept in the analysis. The first term in 2.37 is the square of the inverse of the draft Froude number F%. In order for this term to be important in the analysis, it should be at least Ff = 0(i) This result suggests that the draft Froude number is higher than standard values encountered for normal merchant ships. For these vessels, F® may roughly range from 0.1 to 1.0. This makes the first term in 2.37 of greater importance. 22 We now substitute 2.13 into 2.37: gD D (dh d<pw U2 TUK dt dt ' Dropping terms of 0 (e 2 ) gives gD D d$b d$w U2 TU{ dt dt' + i [ 1 + 2e(§ + § ) + ( § + f ) 2 + ( § + f n = 0 2.40 As before, we notice from comparison of 2.38 and 2.40, that dropping terms in e2 before introducing $ causes the loss of some terms of 0(1) and O(e). Finally, we write the Laplace equation i n non-dimensional form as V 2 ( ^ + h) = 0 2.41 This is not quite sufficient, and it is generally necessary to add a so-called radiation condition that specifies that <j> shall approach <pw in a physically acceptable manner. At infinity one also requires the free surface to become the incident wave. In this work, we have considered a finite closed domain whose rigid boundaries S& and Sw correspond to its bottom and side walls. The slender body problem of highest order in e in an Eulerian form can then be written as d2~6 d2~d> —1 + -4 = 0 in 2 - D 2.42 dy2 dz2 23 ^ + ^ _ A ^ _ 3 _ ^ * _ * * w * £ = = o o » y = *(*,*) 2.43 dy dy TU dt dx dz dz dz dz v ; iR + A S + I[i + 2 c S + + = o on y = c(z, t) 2.44 U2 TU dt 2L dx ydy! ydz' 1 w / _ ^ + _ ? ^ _ ^ _ ^ _ _ e — = n i g i o » z = B(z,y,t) 2.45 dz dz dy dy dy dy dx ^ = 0 onSw , Sb 2.46 cm 2.3 L A G R A N G I A N D E S C R I P T I O N O F T H E P R O B L E M In this work, we have considered the Lagrangian counterpart of the equations described in the previous section. The initial-boundary value problem is then described in dimensional form as d2d> d2d> ^ 4 - ^ = 0 m2-D 2.47 | - f * = 0 on » = ,(*,*) 2.48 dz dip , . lH~dz ony = c(z,t) 2.49 ^7 + <^  - ^ V<£V<£ = 0 wi!/ = ?(z,t) 2.50 at 2 ^ — njjj = 0 on z = B(x,y,t) 2.51 dn — 0 on Sw , Sb 2.52 cm 24 Equation 2.47 describes a boundary-value problem that is solved subject to the moving-boundary conditions 2.48 through 2.52. 25 C H A P T E R T H R E E W E I G H T E D R E S I D U A L F O R M U L A T I O N 3.1 B U B N O V - G A L E R K I N F O R M U L A T I O N The Bubnov-Galerkin formulationflOO] is one of the so called methods of weighted resid-uals. The basic idea of these methods is to obtain an approximate solution to a differential equation -A4> - / = 0 in fi 3.1 where A is a differential operator, / is a given function and <fi is the unknown being sought. Equation 3.1 is solved subject to boundary conditions applied on the contour T surrounding i~l by introducing a set of functions in the form n $ = 3.2 where Cj are constants and $j are linearly independent functions chosen so that global boundary conditions may or may not be satisfied. Since 3.2 is an approximate solution for the unknown <fi, we note that when substituted into 3.1, it wil l not satisfy the governing differential equation. The result is -A$ — f — R 3.3 where R is an error or residual. Let us now introduce a set of weighting functions Wi with 26 i — 1, 2 , . . . , n and construct an inner product (i?, Wi). We now set the inner product equal to zero (R,Wi) = Q 3.4 This is equivalent to forcing the error of the approximate differential equation 3.3 to be zero in an average sense. The choice of the weighting functions Wi determines the weighted residuals' method. In a Bubnov-Galerkin procedure the weighting functions are made equal to the trial (basis, shape) functions <3?^ . Therefore, the Bubnov-Galerkin method is an orthogonal projection of the residual R onto a set of linearly independent functions as implied by the inner product 3.4. In order for <£ to be the exact solution of the given equation, it is necessary that R be identical to zero. If R is considered to be continuous, this is equivalent to requiring the orthogonality of the expression R to all the functions of the system = 1 ,2 , . . . , n.[101]. However, having at our disposal only n constants C i i ^ 2, • • •, C n , we can satisfy only n conditions of orthogonality. On the basis of these considerations, we arrive at the system of equations / R$idn = - / [iVCj$,- + / = 0 3.5 which serves for the determination of the coefficients Cj. Finding the Cj from this system of equations and substituting them in the expression for $ gives the approximate solution. However, if we set the Cj to be the unknowns, the problem is solved when solving the system 3.5. To illustrate this, let us look at the equation describing the flow field of our particular problem stated in Chapter 2, the Laplace equation. This is d24> d2<f> _ dz2 dy2 Let the approximate solution $ be given by 3.2. The residual is then 27 #2.1 a 2 * dz2 + dy2 3.7 and the Bubnov-Galerkin equation becomes (R, $ 0 = J J R$idzdy - j j d 2 4 a 2 $ + dz2 ' 6yJ $idzdy — 0 3.8 One of the most important concepts in mechanics of continua is the relationship con-necting the surface integral for the flux of a vector field with the volume integral of its divergence. The relationship is known as the Green-Gauss theorem and can be expressed as [99] / aVitidn= j aVimdr- f Vicc^dQ Jn Jv Ja 3.9 where the comma denotes a derivative. Applying 3.9 in two dimensions to 3.8 gives a 2 * a 2 * + dz2 dy2 $idzdy d$ d$i <9$ d<f>i dzdy+ J JIdz dz ' dy dy J / <9<£ / d$ + J §i—nzdT + j riydY - 0 3.10 Considering now an arbitrary boundary segment dT we can write dy dT dz dT nz — cos (ft, z) nv = cos (ft, y) then 3.10 becomes // d$d$i 5$ d$i] dz dz dy dy - L dzdy — / <3? L— cos (ft, z) + — cos (n, y) dT 3.11 The second and third integrals in the right side of 3.10 can be written as 28 &4 0<§ J a T = / *i[V#-n]dr= / *i[V-n]dr 3.12 then 3.11 becomes 0* a * a $ - dzdy — ( $j 3.13 Jr where the integral in the right side of 3.13 is only calculated on the boundary. So far, we have stated our problem in a global sense. The essential procedure of the finite element method is that the original continuous fluid domain is divided into small elements of convenient shapes. The differential equation to be solved is replaced by 3.13 which is expressed by a combination of appropriately selected interpolation functions. Since these functions can rarely be found for a whole domain, we find them for each element and satisfy some continuity conditions between elements. Furthermore, these shape functions are defined at suitably chosen nodes within these elements. Therefore, 3.13 can be considered as the equation for each individual element with the interpolation functions written for the corresponding element only. The equations written for each element are collected together to form a global system of algebraic equations. After introducing the boundary conditions, nodal values of the variables can be calculated. We now proceed to discuss the type of elements and the interpolation functions used in this work. 3.2 T H E I S O P A R A M E T R I C E L E M E N T The isoparametric element has been used extensively. Details on its applications can be found in the book by Zienkiewicz[102]. The name isoparametric derives from the fact that the same parametric functions that describe the element's geometry may be used for interpolating spatial variations of a variable within the element. 29 Figure 3-1 Linear (4-noded) and quadratic (8-noded) isoparametric elements. Consider now a quadrilateral 8-noded element as shown in Figure 3-1 . The isopara-metric coordinates £, rj whose values range from 0 to ± 1 are established at the centroid of the element. The reference Cartesian coordinates z, y for the two-dimensional quadratic element are related to £, r\ as follows z - ai + a2£ + a3rj + a 4£ 2 + a5£n + a6r)2 + a7(2rj + a^n 2 V = h + b2£ + hr] + b^2 + hiv + b6r)2 + b7£2n + bsCrj2 We can now write 3.14 in terms of the nodal values Zi = ai + a2£i + a3rji + a 4£ 2 + a5&nj + a&rjf + arifrji + ag&r]? yi - h + b2(i + b3r)i + b4£2 + b_£iTii + b6rif + b7(frii + bs£iTjf with i — 1 , . . . , 8. In matrix form equations 3.15 are {*} = [C]{a.} {yi} = [C]{bi} • with the coefficients of each row of the matrix C defined as Cij = [1 ii Tji iiT}i »j? (?r)i C.nf] for i zzz 1, 8. Then, we can write 30 3.14 3.15 3.16 {*} = [c]-l{*} 3.17 {bi} = [C\-l{Vi} which are the polynomial coefficients in terms of the nodal coordinates. We can see that the matrix [C], made up by the element's nodal coordinates in the £,n system, is the same for both {di} and {b{}. Inversion of [C] in 3.14 and substitution of the coefficients into 3.14 gives n x = y^^jXj i = 1 3.18 and[103] * i = - ( l -0 (1 - « ) ( 1 + ^  + TJ)/4 $ 2 - - ( l + C ) ( l - « ) ( l - ^ + ") /4 *3 = -(1 + 0(1 + - € - *7)/4 $ 4 = - ( 1 - 0 ( 1 + *?)(! + £ 3.19 * B = ( l - * a ) ( l - i j ) / 2 $ 6 = ( l - n 2 ) ( l + 0/2 $ 7 = ( l - a ( l + ^)/2 * a = (1 - «?*)(! - 0/2 for quadratic elements (n = 8) and * i = 4 ( i - o a - ' ? ) *2 = J (1+0 (1 - *? ) I 3.20 *3 = ^ 1 +0 (1 + »?) * 4 = ^ ( l - 0 ( l + «) 31 for linear elements (n = 4). As mentioned above, the isoparametric element uses the same parametric functions to define the relationship between global and local coordinates as to describe the behavior of a given variable within the element. We can then express our approximation function as n * = J ^ ^ $ < 3.21 i-l where $ i and are the same shape functions given by 3.19 (or 3.20) and the nodal values of the velocity potential respectively. These nodal parameters are then equivalent to the constants Cj in 3.2. 3.3 F I N I T E E L E M E N T F O R M U L A T I O N In section 3.1 we obtained an expression that was the result of a weighted residual procedure. This expression, 3.13, contained approximation functions that were to be defined on an elemental basis. In section 3.2, those functions were denned for linear and quadratic isoparametric elements. In this section, we combine equations 3.13 for the element with equations 3.18 and 3.21. The Bubnov-Galerkin procedure yielded Introducing 3.21 into this equation, our element equation 3.22 becomes or where the k and I subindexes refer, respectively, to the element and the side of the element 32 in contact with the contour of the domain. Equation 3.24 contains z and y variables, yet the interpolation functions were expressed in terms of the local £, r? coordinates. Therefore, we require a transformation between the two coordinate systems. Consider then the following chain rule expressions d$j _ d^j_dz_ d<$>j dy d$j _ d$j dz d$j dy dn dz dn dy dn 3.25 or in matrix form and after inversion we have - dij - ' dz - dij -»€ H dz dz_ dy - drj . drj drj - dy . - dij - - d$j - - dij -dz - T'1 Jll Jl2 dij — J dij h i J22 d$j . dy . . dr) • • drj -3.26 where J is the Jacobian matrix given by J = dz d]L dt dz d]t |_ drj drj The integration over the domain Qk refered to the cartesian coordinates must also be changed to the isoparametric coordinates. This transformation can be achieved as[104] J J dzdy = J + J + det[J]dfdn 3.27 Introducing 3.26 and 3.27 into 3.23 gives J11 d$ji + J '21- + J; d$j 22 dn J 2 i ^ r - r J 2 2 ^ j d( ' ' " dn J \det[J]d£dn = / *< — d-Yi > Jr dn 3.28 33 We can see that the left side in 3.28 gives an 8 X 8 symmetric algebraic system for each quadratic element and 4 x 4 system for linear elements. The resulting matrix wi l l be termed as the element coefficient matrix. 3.4 N U M E R I C A L I N T E G R A T I O N The integrals in equation 3.28 require numerical integration. This can be efficiently performed by means of Gaussian quadrature. For one and two-dimensional cases the formulae are + 1 m9 1=1 where Hi and Hj are the weight coefficients[105]. The functions / ( & ) and f(£i,rjj) are to be evaluated at the mg Gauss points. The value of mg depends on the order of the polynomial that can be integrated exactly using Gauss quadrature and is given by Op — 2mg — 1 However, inspection of equation 3.28 does not yield an order of an exact polynomial. This is because the determinant of the Jacobian is involved in the denominator of this expression. It was seen that for the elements used in this work, more than 3 x 3 integration points had little effect on the results. 34 C H A P T E R FOUR E L L I P T I C G R I D G E N E R A T I O N It has been pointed out by several researchers that regardless of the numerical method used for the solution of partial differential equations, the accuracy obtained is closely related to the 'quality' of the generated mesh. One component of this quality is the orthogonality characteristics of the grid. A thorough description of different procedures for the estimation of orthogonal grids has been given by Thompson[106]. Ell iptic grid generation systems were reported to provide orthogonal grids of good accuracy. One system that is particularly attractive is the one developed by Ryskin and Leal[107]. It requires the calculation of only one line control function for both grid generating equations in two dimensions. In this chapter, we use a Bubnov-Galerkin formulation to solve the weak constraint method[107] for the construction of boundary-fitted orthogonal curvilinear coordinates. N u -merical computations are carried out on a particular domain studied by Chikliwala and Yort-sos [108]. This domain has been used here in order to facilitate comparison with reported results and because the wave-like shape of one of its boundaries is commonly encountered in the field of marine/ship hydrodynamics. Although orthogonal coordinates are not required when applying the finite element method, it is known that the error in the calculation of a variable on the element depends on the angle formed by adjacent sides[109]. Therefore, the lowest departure from orthogonality is likely to give the best results. 4.1 E L L I P T I C G R I D G E N E R A T I N G S Y S T E M S Ell iptic systems have long been recognized for their advantages in generating a globally smooth grid even for complicated and/or discontinuous domains. This smoothness charac-35 teristic stems from the extremum principle possessed by Laplace-type elliptic systems. It can be measured in two dimensions by the following integral[110,lll,112]: V e |2 + | Vg ds where e and g are boundary fitted curvilinear coordinates. Using calculus of variations, it can be shown that minimization of Is is equivalent to the minimization over the entire domain of the weighted integral of a quantity P[113,114]. This integral is j J WiPids which is representative of the accumulation of some property of the grid. Here, the subindex i refers to summation over the repeated index. If the property P of the grid is chosen to be the gradient of the curvilinear coordinates and equal to the weighting function w, we obtain the same integral as for Is. We then seek a solution with Pi — Wi — Ve» that wi l l make Is a minimum over the whole physical domain. The actual solution of such integral for each curvilinear coordinate e% is[l 13,114] V 2 e = 0 4.1 V 2 p = 0 These are Laplace equations for each coordinate e, g. The above considerations clarify the notion of a globally smooth mesh as generated by an elliptic system. Equations 4.1 generate a system of coordinates that tends to be equally spaced provided that the boundaries are straight and that the physical domain has a uniform coordinate distribution. However, in the presence of curvature, lines wil l tend to be more closely spaced over convex boundaries. The opposite effect is obtained over concave boundaries[115]. 36 4.2 P O I S S O N S Y S T E M S The grid generated by system 4.1 only allows for control of the coordinate spacing on the boundaries of the physical domain where coordinates are specified. In order to exercise control of the spacing between curvilinear coordinate lines inside the physical region, equations 4.1 must be altered. A number of alternatives have been suggested in order to accomplish this effect [113,116,117,118,119]. A general form of a system where line control is possible can be written as [115] V 2 e = Fe 4.2 where the control functions Fe and Fs can be chosen so as to alter the orientation and spacing of the coordinate lines. In this case, there is a possibility of losing the extremum principle. Thus, the coordinate lines may cross each other and/or exceed the boundaries of the physical region. However, the extremum principle is a sufficient, but not necessary, condition for a one-to-one mapping[113]. 4.3 C O O R D I N A T E T R A N S F O R M A T I O N It has been shown [120] that if a curvilinear coordinate system that satisfies the Laplace equation V 2 X i = 0 (» = 1,2,3) is transformed to another coordinate system Ci [i — 1, 2, 3), the new curvilinear coordinates ei satisfy the Poisson or inhomogeneous elliptic system V2a-Fi ( i = l ,2 ,3) 37 4.3 where 3 3 ;'=1 k-l with gik being the components of the contravariant metric tensor and F j k _ d s m dxn d2ej 4 _ 1 ^ dey Set dx-ndxn m—1 n—1 J Warsi[120] has shown that a grid with lines concentrated by applying a subsequent transformation, also known as stretching transformation, to a grid generated as the solution of the Laplace system could have been generated directly as the solution of Poisson sys-tem 4.3 with the control functions F^k derived from the subsequent transformation using 4.4. Therefore, the Poisson system 4.3 seems appropriate as the generating system with control functions that wil l be arbitrarily specified rather than obtained from the subsequent transformation [115]. Thus, an appropriate generating system is given as [115] V 2 e ^ E E ^ ' f c 4.6 j-1k-l with F(k to be specified. The components g*k, defined by equation a-39 (all equations starting with a are described in Appendix A ) as the scalar product of the contravariant base vectors, are the components of the contravariant metric tensor whose determinant defines the Jacobian of the transformation. System 4.6 produces a coordinate system that is the result of applying a subsequent transformation to a system generated for maximum smoothness (solution of the Laplace equation). From equation 4.5, the control functions F£% are known as the control functions for one-dimensional stretching, i.e., the control functions that exercise concentration of the coordinate lines in only one direction. In practice, the other control functions are neglected thus Fik = 6)b\Fi 38 where 6j and 6%k are Kronecker delta functions. The generating system then becomes V2a = g"Fi 4.7 From equation a.39 in Appendix A , we have gV = gi* = cT • c? (a.39) that, when substituted into 4.6 , gives 3 3 j'=l k=l or from the definition of contravariant base vectors (a.3) j=ifc=i Equation 4.9 represents a system where ej are the dependent variables and x%, the Cartesian coordinates sought, are the independent variables. In order to solve for x%, equa-tion 4.9 is transformed to a domain known as transformed or computational, where ei are the independent variables and Xi the dependent ones. For computational convenience, the transformed domain is usually chosen as a unit square. The actual transformation of 4.9 is obtained using equation a.48, i.e., the non-conservative expression of the Laplacian operator 3 3 3 or va*= + £( v 2^)%- 4.10 i—Xj—X j—X 39 where v is the position vector of a point in Cartesian coordinates. It was previously mentioned that the appropriate generation system 4.6 was chosen because it was the result of a stretching transformation applied to a system governed by Laplace's equation designed for maximum smoothness. This last system is that obtained by 'back-transforming' our initial generating system 4.6 to obtain 4.10 . We then have W = 0 or 3 3 3 E E ^ ' ^ ; + E ( V 2 * = 0 4.11 t = l J = l k - l which is the transformation of the original Poisson system. At this point, a number of options[120] exist for the substitution of the term V2efc in 4.11. The option used in this thesis is given by equation 4.2. The resulting equation is[120] 3 3 3 E E + E ^ * = 0 4 - 1 2 i - l j - 1 k-l which is a generation system that has been widely used [121,122]. The present work is concerned with grid generation in two dimensions. The two dimensional grid generation equation 4.12 then becomes gnvS£ + 2g12viee + g22v,ee + F£€,e + Fev,e = 0 4.13 Equation 4.13 can be written in terms of the covariant metric tensor components. In two dimensions, the coordinate x3 is invariant and e3 coincides with x3. The tangent to the coordinate line e3, i.e., a3, and the normal to the plane where e3 is constant, i.e. a3, are the 40 same. In this case a3 = 03 = fc. We can then write g33 — 933 — a3 • a3 - of • a3 — k • k — 1 9i3 — ?3i — »i • 03 = Si • k — 0 923 — 932 — Q>2 • 0.3 — 0,2 • k — 0 From (a.42) we have 22 _ 9ii 9 9 1^1 _ 922 9 12 _ 21 _ 912 , -9 — -and finally, 922V,es - 2gi2v<es + guv,ee + g{FEv,6 + Fev,e) - 0 4.14 We now make a choice for Fe and Fs. Expressions for these control functions were derived by Ryskin and Leal[107] as 1 df /gii^/922 de Fe -9uy/9~22 de V/y f2 _ 922 J 9ii A necessary and sufficient condition for the coordinate lines to meet at right angles, that is, are orthogonal to one another, is that g^ — 0[123,124]. It should noted that Fe and Fs are specified control functions and not functions derived from the subsequent transformation 4.4 . Therefore, the lines could not always be perfectly orthogonal. Introducing 4.15 into 4.14 gives 41 1 -. df cf fi\_ or de + 1__ dg lfV'1 = 0 4.16 The vector v can be defined as v — zi + yj. Taking each directional component separately results in d_ dt dz f— J de. dy d de V de + 1 az-Cf dg ldy dgifdgl - 0 4.17 Equations 4.17 define an elliptic system first used by Ryskin and Leal[107] and applied to a specific geometry using a finite difference discretization by Chikliwala and Yortsos[108]. In the following developments a Bubnov-Galerkin formulation is applied to system 4.17 and compared with results reported in [108]. 4.4 B U B N O V - G A L E R K I N P R O C E D U R E In this section, we develop the method to solve equations 4.17 . This equation is valid in the whole transformed domain O surrounded by the contour curve T. Solutions are sought at discrete points selectively specified. The method used to attain this goal is, as in the previous chapter, a Bubnov-Galerkin procedure. For the z-coordinate we have Ye f dz~\ a r 1 dz de\ dg if dg = 0 Let the approximate function z for the variable z be given by the expression 42 z = j = l 4.18 then the residual for the z-coordinate equation in 4.17 becomes d_ de de d_ do ldz [fdg R 4.19 and the inner product (R, $ j ) gives {R,*i) = J JR$ided8= j J{4z\fdZ de I de d_\ldz dgifdg — $ i \dedg = 0 4. do\ J 20 Applying the Green-Gauss Theorem to 4.20 in a two-dimensional domain gives //{ dzd$i 1 05 3 $ n , , f—— h — — — - — fdeag de de f dg dg > 4.21 dz d$i f f 1 f Y e ^ d e d S - J J / dz 9$j dg dg dedg — 0 or //{ f°*?*i + ^^Xdedg J de de ^ fdg dg J * ;r <. \ _ de ) 4.22 Similarly, for the y-coordinate we write //{ dyd*j ldydij *de de fdg dg 4.23 43 with the approximate function y specified as n 4.5 F I N I T E E L E M E N T F O R M U L A T I O N In this section, we develop the integral equations to be solved for delivering nodal coordinates. The upper l imit in the summation signs wil l be left unspecified and understood as n = 4 and n — 8 for linear and quadratic elements respectively. The Bubnov-Galerkin procedure for the z-coordinate gives fff.dzd.i ld£d.i\__ f ( f r d ^ \ ^ (ldz\~\^ Introducing our approximate functions z and y gives an elemental equation as where the subindexes fc and I refer, respectively, to the element ilk and the boundary Ti in contact with the contour of the domain. Equation 4.25 contains the e and q variables, yet the interpolation functions are expressed in terms of the local (£, tj) coordinates. Therefore, we require a transformation between these two domains. This transformation corresponds to the mapping of an element on the computational domain (e, p) to an element defined by the isoparametric coordinates. For convenience, the transformed domain can be taken as a square with unit sides and our element equation 4.25 applied to a given element on this square. 44 Consider now the following chain rule expressions d ^ _ d $ y d e d$j dg d£ ~ l)idl, + dg dC d$j _ dij de d$j dg dr) de dr) dg dr) These equations can be written in matrix form as - d$j • ' ds - d$j • n n ds de dij . drj • . drl dr) . de . Inversion of this system gives - d$j • • d$j - - dij -ds Jll T* J12 H dij T* . J21 T* J22 . dij . de . - drj . - drj . where J* is the Jacobian matrix given by r = de 9£ de_ d_ drj drj 4.27 The integration over the domain i~lk refered to the (e, g) coordinates must also be changed to the isoparametric coordinates. From (a.15) we then have •+i r+i J J dedg = J J det[r)d£dr) 4.28 Introducing 4.27 and 4.28 into 4.25 gives i i fi ^ i i 9( '12 dr) •'11 „ j. T «' 12 d$i d£ dr) + J21 ~>  J2\ J J i ^ + ^ 2 2 ^ ] }«fc*[jr]d&*r, = 4.29 2 1 Similarly, for the ^-coordinate we have 45 }det[J]d(dri 4.30 We can see that the left sides of 4.29 and 4.30 give an n X n algebraic system for each element. Note that 4.29 and 4.30 wil l generally give a non-symmetric matrix. This inconvenience can be circumvented by obtaining an average value ff. for each element. The right side in 4.29 and 4.30 contains line integrals with derivatives of the Cartesian coordinates evaluated at any given boundary, the boundary condition is termed as Neumann boundary condition. When the right side of 4.29 or 4.30 is made equal to the z or y coordinate of the boundary and the row and column of the corresponding main diagonal element are zeroed, the boundary condition is termed as Dirichlet boundary condition. Equations 4.29 and 4.30 require the numerical estimation of surface and line integrals. These can be performed using Gauss quadrature as described in Chapter 3. In this chapter, we have used 2 x 2 points for linear elements and 3 x 3 for quadratic ones[103]. 4.6 C A L C U L A T I O N O F T H E S H A P E F A C T O R The shape factor / was defined in 4.15 as with respect to the (e, q) coordinates of the computational domain. When this integral is 2 _ 922 911 4.31 The components of the contravariant metric tensor were defined in (a.50) as 46 2 i 2 911 = * e + ytS 922 - z\ + y\ (a.50) Using the approximate functions for z and y gives de J ' de l ; ' n 4.32 In these equations, we have the derivatives of the parametric functions with respect to the computational coordinates (e, g). The same transformation defined by 4.27 is needed. We can then write 4.32 as t _ 1 8 - 1 4.33 »»= [ E <J« ^ + + l £ H j 2 ^ + j 2 2 ^ - J . 4.7 N U M E R I C A L C O N S I D E R A T I O N S Observation of 4.29 and 4.30 together with 4.33 shows that the matrix system to be solved for the generation of the coordinate system (z, y) is a nonlinear function of the Cartesian coordinates being sought. A n iterative procedure is then required to solve systems 4.29 and 4.30 . However, there is an initial difficulty. Equation 4.33 cannot be evaluated at the first iteration since only the (z, y) coordinates of the contour of the physical domain are known. The coordinates in the interior of the physical domain are the solution of 4.33 at this very first iteration. Therefore, the process followed in this chapter is the following: 47 a: b: • Define boundary points of the physical domain • Solve the following elliptic system using Dirichlet boundary conditions [107] d2f d2f — + — = 0 4.34 dz2 + dg2 In a Bubnov-Galerkin formulation this gives + 1 ^ 2 1 o/- + ^22 */21 n<? ' J22' d£ z z dr, J I * di ' " drj }det[J]dtdn - J $*|4*Ti 4.35 Dirichlet boundary conditions for 4.35 are defined using various functions for / = ri) a ^ ^ n e boundary. These functions are described in the following section, c: • W i t h calculated / values from step b solve system 4.29 and 4.30 using specified boundary conditions. At this stage, this was done using Cholesky's method[125]. d: • W i t h newly calculated z and y values, recalculate new / values using 4.33. e: • Recalculate z and y and go back to step d until convergence is obtained. The number of iteration cycles required for convergence is the number of times the program requires to loop through steps d to e. Although the convergence criteria utilized by Chikliwala and Yortsos[108] was considered somewhat arbitrary, it has been adopted here solely for comparison purposes. This convergence criteria was specified by an index termed MDO — \^ir — 6\ that measured the localized maximum deviation from orthogonality. Convergence was obtained when MDO at a fixed point stabilized within one decimal point to a constant value. The angle 9 is calculated from (a.52) as a 012 cos v — ——=. or 48 cos 6 — TXT2 + r 3 r 4 (T? + T*)1/2(T22 -rT*)1/* 4.36 with, the T% terms defined as i=i j-1 9t dr] A more general convergence index could be the average deviation from orthogonality (ADO). This value would give the grid's overall deviation from orthogonality since it takes into account all the corners of each element. The index ADO is defined for the elements used in this work as k - l i - l where the subindex k refers to the element number and the superindex i refers to the el-ement's corner - N being the number of elements. The index ADO was not used as a convergence criteria. However, values of ADO are presented together with those of MDO to give a global idea of the orthogonality characteristics of the whole grid. We now briefly mention the types of functions used to define the shape factor / on the contour of the domain used in step b to solve 4.35 . This definition wil l be identified by the parameter I SHA as follows 49 Table 4-1 Boundary conditions as given by Chikliwala and Yortsos Boundary Conditions Applied in z and y Directions N : Neumann, D: Dirichlet Boundary Case A Case B Case C Left y: D z: D y: N z: D y: N z: D Upper y: D z: D y: D z: D y: D z: D Right y: D z: D y: D z: D y: N z: D Lower y: D z: D y: D z: D y: D z: N I S H A =0. Finite element definition of / using 4.33. I S H A =1. Uniform shape factor distribution. ISHA =2. Exponential shape factor distribution. I S H A =3. Linear shape factor distribution. I S H A =4. Normal shape factor distribution given by the probability density function 1 rr e»-°- 6^ 2i 4.37 ao — a constant Note that ISHA = 1,2,3 are taken from [108]. 4.8 N U M E R I C A L R E S U L T S The following sections are dedicated to obtain results for a given geometry of a phys-ical domain presented by Chikliwala and Yortsos[108], Figure 4-1 , using the methodology developed in this chapter. Boundary conditions specified in[108] are also used. These are shown in Table 4-1 as cases A, B and C. 50 H H z » ed A Z z = d 6-1 Z = 0 E - 0 A = 0.75 d = 1.0 e = [0,1] Figure 4-1 Typical Domain. 4.8.1 Case A : Complete Boundary Correspondence It was reported by Chikliwala and Yortsos[108] that no satisfactory orthogonal grid could be obtained using complete boundary correspondence when solving system 4.17. The cause of this failure was attributed to the asymmetry of the physical domain with the en-suing suggestion that complete boundary correspondence could only be used for symmetric regions. The method presented in this work shows otherwise. Table 4-2 shows results for the asymmetric region used in [108] with equidistant boundary point distribution in all sides but the left one where an exponential boundary point distribution has been applied. The first quantity in each column refers to MDO, the second to ADO and the third to the number of iterations required to attain convergence. This format of presenting tabled information wil l be used throughout this chapter unless otherwise stated. Table 4-2 shows that good localized (MDO) and overall (ADO) accuracy can be obtained using complete boundary correspondence in a non-symmetric domain. Overall, the best values are obtained using a 51 Table 4-2 M D O , A D O and number of iterations for Case A Case A : Complete Boundary Correspondence Exponential Boundary Point Distribution On Left Side Linear Isoparametric Elements M D O - A D O - N u m b e r of Iterations Case H=0.05 H=0.10 H=0.15 H=0.20 H=0.25 ISHA=0 3.6-0.9-19 4.9-1.2-20 6.3-1.7-20 8.3-1.9-21 9.6-2.9-20 ISHA=1 1.7-0.6-6 3.3-1.1-6 4.9-1.5-6 6.9-1.7-6 6.9-2.1-9 ISHA=2 1.8-0.6-10 3.4-1.0-10 5.1-1.5-10 8.9-1.6-12 9.4-2.7-16 ISHA=3 1.6-0.6-8 3.1-1.1-8 4.8-1.4-7 8.0-1.6-12 10.7-2.6-12 ISHA=4 1.7-0.6-6 3.0-1.1-7 4.3-1.4-7 5.2-1.4-7 6.8-2.0-7 normal distribution for / as described by 4.37. However, the difference in MDO for each H does not show a value greater than about 3° among the five initial shape factor distributions. The corresponding difference in ADO is not greater than 0.9°. Assuming that the effect of these differences on subsequent engineering computation is not very significant, the method seems to yield acceptable orthogonal grids regardless of the initial shape factor distribution. Solutions in [108] were reported as strongly dependent on such distributions. The number of iterations required to attain convergence was lowest for ISH A — 1 and ISH A — 4, i.e., constant and normal initial distribution of the shape factor. Therefore, from a CPU time viewpoint these initial distributions seem to be more advisable. For ISH A = 0, 1, and 4, the number of iterations required for the convergence of MDO to be satisfied was almost independent of H. For ISH A — 2 and ISH A — 3 a slightly higher number of iterations is required for H = 0.2 and H - 0.25. Figure 4-2 shows an orthogonal grid using complete boundary correspondence for a 16 x 16 grid with H - 0.15, ISH A = 4. For this case, MDO was 4.3 and ADO was 1.4. A similar level of accuracy was obtained for a symmetric region with equidistant boundary correspondence in a 32 X 16 grid with ISH A — 3 and H — 0.15. For this case, Figure 4-3 , 52 Figure 4-2 16 x 16 grid, H — 0.15, I S H A — 4. Dirichlet boundary conditions with exponential distribution of points on left boundary. Figure 4-3 32 X 16 grid, H = 0.15, I S H A - 3. Dirichlet boundary conditions with equidistant distribution of points in all boundaries. M D O was 4.2 and A D O was 1.1. Greater values of H were used for I S H A — 4 in order to verify if the method would still converge for very large H values. Table 4-3 shows results for this case. Convergence as well as good overall orthogonality characteristics are obtained. The localized orthogonality is seen to deteriorate as H becomes extremely large. 53 Table 4-3 M D O , A D O and number of i terations for Case A Case A : C o m p l e t e B o u n d a r y Correspondence E x p o n e n t i a l B o u n d a r y Point D i s t r i b u t i o n O n Left Side L inear Isoparametric Elements M D O - A D O - N u m b e r of Iterations Case H = 0 . 3 H=0 .4 H=0 .5 I S H A = 4 12.8-6.1-5 14.6-5.8-7 22.7-11.8-4 Table 4-4 M D O , A D O and number of iterations for Case A Case A : Comple te B o u n d a r y Correspondence Equid is tant B o u n d a r y Point D i s t r i b u t i o n Linear Isoparametric Elements M D O - A D O - N u m b e r of Iterations Case H=0 .05 H=0.10 H=0.15 H=0.20 H=0.25 I S H A = 0 6.2-2.5-17 11.1-4.7-16 16.2-6.2-16 20.4-7.2-17 24.5-8.2-20 ISHA=1 4.2-2.8-6 8.6-5.6-6 13.1-8.1-7 17.5-10.5-11 21.6-12.6-12 I S H A = 2 4.2-2.8-5 8.5-5.6-11 12.9-8.1-15 17.2-10.5-13 21.2-12.6-13 I S H A = 3 4.6-2.8-8 9.2-5.6-7 14.1-8.1-11 18.5-10.8-14 23.1-12.6-15 I S H A = 4 4.4-2.8-6 8.3-5.6-7 12.6-8.1-10 16.9-10.5-13 20.9-12.6-14 Results for MDO and ADO are shown i n Table 4-4 using equidistant boundary point d i s t r i b u t i o n for the d o m a i n shown i n F igure 4-1. Values of MDO and ADO are roughly 2 and 4 t imes higher t h a n those for exponential boundary point d is t r ibut ion . However , very good convergence characteristics are kept. T h e number of iterations is again seen to be the lowest for ISHA — 1 and generally for ISHA - 4. It should be noted that results for MDO and ADO i n Table 4-4 correspond to only 3 i terations. T h e number of iterations stated i n this table remains the value necessary to a t ta in convergence. T h e reason for g iv ing results at i terat ion 3 can be seen f r o m curves shown i n 54 Figure 4-4 Convergence characteristics for exponential and equidistant distribution of points. Figure 4-4 . This figure shows typical convergence characteristics for MDO and ADO as a function of the number of iterations. Although the behavior of ADO has a consistently decreasing tendency, MDO does not. The index MDO tends to show a minimum at the third iteration but stabilization occurs at a later iteration value. A very high rate of convergence in the initial few iterations was observed throughout this study. This is particularly important for subsequent computations where orthogonality is not strictly required. A satisfactory global level of orthogonality can then be obtained to improve results after very few cycles of iteration without requiring convergence. 4.8.2 Case B : Dirichlet and Neumann Boundary Conditions In this section, equidistant boundary correspondence is applied to all boundaries except the left side where Neumann boundary conditions are applied. The same non-symmetric geometry as in case A is studied. 55 Figure 4-5 16 X 16 grid, H — 0.25, ISHA = 2. Neumann boundary condition on left side. Table 4-5 shows results of MDO and ADO for this case. Values of MDO for ISHA — 3 are slightly higher than those reported in [108], the difference roughly ranging from 0.2° to 0.8°. A similar difference occurs with ISHA — 2 for H < 0.15. However, for larger values of H, i.e. very steep wave forms, improved orthogonality characteristics are obtained especially for H — 0.25 where the scheme presented in [108] failed to converge. This last case is shown in Figure 4-5 . Also for ISHA — 1 (linear shape factor distribution), it was reported by Chikliwala and Yortsos[108] that orthogonal grids could not be obtained. This is not the case in this work. Figure 4-6 shows results for a 16 x 16 grid, H — 0.15, ISHA — 1. We are able to attain very good localized and overall orthogonality characteristics. Therefore, it can be concluded that the present method gives very good orthogonality characteristics regardless of the initial shape factor distribution. Unconditional convergence is also observed for all H values and shape factors used. In general, the number of iterations is only slightly affected as H increases. For the case shown in Figure 4-6 with H = 0.1, convergence rates for MDO and ADO as a function of the number of iterations are shown in Figure 4-7 . In the presence of Neumann boundary conditions both MDO and ADO exhibited a decreasing tendency until 56 Figure 4-6 16 X 16 grid, H — 0.15, I S H A = 1. Neumann boundary condition on left side. Table 4-5 M D O , A D O and number of iterations for Case B Case B : Dirichlet and Neumann Boundary Conditions Neumann Boundary Condition On Left Side Linear Isoparametric Elements M D O - A D O - N u m b e r of Iterations Case H=0.05 H=0.10 H=0.15 H=0.20 H=0.25 ISHA=0 12.2-3.7-9 10.8-3.3-10 10.6-3.1-11 11.1-3.3-11 12.1-2.9-8 ISHA=1 1.6-0.6-10 3.1-0.9-13 4.6-1.2-13 5.9-1.6-14 6.8-2.3-16 ISHA=2 1.5-0.5-15 2.9-0.8-15 4.3-1.0-15 5.6-1.5-15 8.5-2.5-15 ISHA=3 1.4-0.4-10 2.7-0.7-12 4.5-1.1-14 5.1-1.3-10 7.7-2.9-7 ISHA=4 1.3-0.3-10 3.0-0.8-13 4.3-1.1-14 5.2-1.4-10 6.7-2.1-13 convergence was obtained. However, for a similar level of accuracy the rate of convergence seems to be slower than for Case A using exponential boundary point distribution. 4.8.3 Case C : Dirichlet and Neumann Boundary Conditions In this case, a combination of Dirichlet and Neumann boundary conditions is used in all boundaries but the upper side where Dirichlet boundary conditions are specified for both 57 CASE B: NEUMANN BOUNDARY CONDITION ON LEFTSIDE 20 g= 12 4 0 ! ! ' : ! ' V j .A....; V \ : V : • - • i -MDO ADO \ : \ [ S • \ : \ \ ^ ^  ^  T.r..-,... \ i ; -i i 1 i ; i ; 0 2 4 6 8 10 12 14 NUMBER OF ITERATIONS Figure 4-7 Convergence characteristics for Figure 4.6 with H=0.1. z and y coordinates. Results for 4 and 8 noded isoparametric elements are presented. Table 4-6 shows results of MDO and ADO for the geometry of Figure 4-1 as a function of the number of elements and ISH A using linear elements. Values of MDO in this table are consistently higher than those reported by Chikliwala and Yortsos[108], the difference generally ranging from 1.5° to 5.0°. Values of ADO are very low securing overall orthogonality. This was observed to be close to 90° in most elements of the grid. The MDO parameter shows very little improvement, if any, over those given in Table 4-2 and Table 4-5, i.e., case A and case B. In other words, localized orthogonality is not improved by increasing the use of Neumann boundary conditions. However, significant improvement in the overall orthogonality characteristics (ADO) is observed compared to results shown in Table 4-2. In some cases the improvements were seen to be of the order of 50%. Although I SHA — 0 gives the poorest results in case B, some of the best results are obtained in Case C with this particular initial shape function distribution. Also, as N increases, improvements in MDO 58 Table 4-6 M D O , A D O and number of iterations for Case C Case C: Dirichlet and Neumann Boundary Conditions Distribution Of Boundary Conditions As Defined by Chikliwala et al. Linear Isoparametric Elements M D O - A D O - N u m b e r of Iterations Mesh Size ISHA=0 ISHA=1 ISHA=2 ISHA=3 ISHA=4 8 x 8 6.7-2.2-8 6.6-1.5-6 6.6-1.5-12 10.2-2.2-6 7.6-1.5-6 16 x 16 4.6-1.4-9 5.9-0.7-7 4.9-0.8-6 7.0-1.1-9 6.0-0.7-9 24 x 24 3.9-1.0-11 4.6-0.5-7 3.7-0.6-6 4.5-0.8-16 4.6-0.5-8 36 x 36 3.7-0.8-10 3.8-0.4-8 2.9-0.4-6 5.6-0.6-8 3.9-0.4-8 40 x 40 2.9-0.6-10 3.2-0.3-8 2.6-0.3-6 3.6-0.5-19 3.4-0.3-7 48 x 48 2.3-0.5-10 2.8-0.2-8 2.3-0.3-6 4.4-0.4-9 2.9-0.2-6 tend to be more noticeable. Table 4-7 shows results for the same conditions as Table 4-6 for quadratic elements. In this case, overall improved results are obtained for both MDO and ADO. For the latter, the improvement is very significant with reductions exceeding 50%. Figures 4-8 and 4-9 show typical results for a 16 X 16 grid using linear and quadratic elements respectively. The position of the mid-side nodes for quadratic elements is shown in Figure 4-9. A slight reduction in the number of iterations to attain convergence is evident from the comparison of Tables 4-6 and 4-7 suggesting a higher rate of convergence for quadratic elements. Convergence rates corresponding to the results shown in Figures 4-8 and 4-9 are shown in Figure 4-10 . The higher rate of convergence exhibited by quadratic elements is apparent in this figure. 59 Figure 4-8 16 x 16 grid, H = 0.15, ISHA — 2. Linear isoparametric elements. 60 Table 4-7 M D O , A D O and number of iterations for Case C Case C: Dirichlet and Neumann Boundary Conditions Distribution Of Boundary Conditions As Defined by Chikliwala et al. Quadratic Isoparametric Elements M D O - A D O - N u m b e r of Iterations Mesh Size ISHA=0 ISHA=1 ISHA=2 . ISHA=3 ISHA=4 8 x 8 6.3-1.8-19 4.3-1.0-13 3.7-0.8-9 3.3-0.5-4 3.3-0.5-5 16 x 16 4.9-1.0-7 2.3-0.4-6 1.4-0.2-5 1.9-0.4-19 3.9-0.5-5 24 x 24 3.3-0.6-9 2.0-0.3-6 1.9-0.3-4 1.8-0.3-17 4.8-0.4-3 36 x 36 2.9-0.4-11 1.7-0.3-6 1.9-0.3-4 2.8-0.4-11 2.9-0.4-11 40 x 40 2.7-0.3-11 1.4-0.3-6 2.0-0.3-4 4.3-0.3-5 4.6-0.3-5 48 x 48 3.4-0.4-8 1.6-0.2-6 2.2-0.3-4 1.7-0.2-8 1.0-0.1-8 CASE C: NEUMANN BOUNDARY CONDITIONS ON THREE SIDES 5? 40 ! I ! I 1 . I . \ \ :\ • \ • \ : \ : \ i MDOUNEAR ADOLINEAR MDOQUADRATIC ADO-QUADRATIC \\ ! IN . . . - - - - - - -i i • • 0 1 2 3 4 5 6 7 NUMBER OF ITERATIONS Figure 4-10 Convergence characteristics for Figures 4.8 and 4.9. 61 C H A P T E R F I V E N U M E R I C A L C O N S I D E R A T I O N S In this chapter, we present various numerical schemes that have been used in treating the boundary conditions involved in this work. Concepts and/or variables that are required in the solutions presented in Chapter 6 are also included in this chapter. 5.1 BERNOULLI'S CONSTANT It was mentioned in Chapter 2 that Bernoulli's constant should be taken into account when imposing the dynamic boundary condition in certain cases. One of these cases is a surface piercing body moving vertically in a closed domain. When considering the time derivative on (p, the term C(t) should be evaluated accurately since the pressure field is strongly dependent on this value. For a surface piercing body with vertical sides performing heaving yi in a tank of rectangular cross section, C[t) is[58] C(t) = c0(t) = ^ 5.1 w where w is the distance between the body's vertical side and the wall of the tank and y j is the vertical displacement of the cylinder. The value so represents the increase/decrease in water level due to a downward/upward quasi-static movement of the surface piercing body. Equation 5.1 vanishes for submerged bodies and for surface piercing bodies moving laterally in a tank. This is because the net flux from the body is zero[58j. The boundary value problem in Lagrangian form is then described by d2cf> d2d> 62 dy d(j> - - - = 0 ony = c(z,t) 5.3 dz d<f> dt dz 0 ony-c(z,t) 5.4 d(b 1 + - 2 V W = 0 on y = c(z, t) 5.5 00 — i%iqi = 0 on z — B(y} t) or z = B(x, y, t) (771/ 5.6 dn — 0 on Sw , Sb 5.7 5.2 T R E A T M E N T O F T H E B O U N D A R Y C O N D I T I O N S 5.2.1 Free Surface Boundary Conditions B y a p p l y i n g the definit ion given by 3.26 to the derivatives on cj> w i t h respect to z a n d y, we can w r i t e 5.3 a n d 5.4 as dz dt dy dt n 7=1 S ' 8 on y — c(z,t) any - c(z, t) 5.8 5.9 where the subscript n refers to a free surface node. T h e d y n a m i c boundary condit ion 5.5 can s imi lar ly be wr i t ten as 63 d(p] , ... 1 + + " A r d$j f 8$, ony- q(z, t) 5.10 5.2.2 Body Boundary Condition Equation 5.6 enters the problem as a Neumann-type boundary condition. It is the normal velocity on the surface of the moving body and is included in the integration of the right side of 3.28. It can be written for vertical motion as (V • n){ and for the submerged case as dip dn 1/ J i J \ nv qh\ft\ 5.11 (V • rt)? -with | rly |, | n z | and | n | defined as d$ dn J i J i 5.12 dB dy dB djj_ de v^8 2-j-i zi dij + Jh d$j " drj f* dij 9? + H2 d$f drj 5.13 nz 1=1 5.14 64 n = 5.15 The Jlm terms are the elements of the following matrix r = de 9£ de_ d± drj drj where J * is the Jacobian of the transformation from the computational domain to the domain with isoparametric coordinates. Details on this transformation were given in Chapter 4. 5.2.3 Treatment of the Right Side in Equation 3.28 Once the normal velocity has been evaluated as shown in Section 5.2.2, the right side of 3.28 must be obtained. This is where Y\ is the boundary of the element that coincides with the moving body contour in the physical domain. Since this contour is arbitrary, we again resort to certain concepts of differential geometry. A n increment of arc length described by a vector position f on a curvilinear coordinate line along which £j varies is given by [124,106] with gu being the components of the covariant metric tensor. In this work, the body surface 5.16 was made to coincide with the north side of the computational domain. On this side, the curvilinear coordinate e varies while Q is constant. We then write 5.17 65 where which when introduced into 5.17 gives or, using 4.33 5.18 5.19 r, d n dn JJ + 'j = l J l 1 ^ + J l 2 c?7? de 5.20 The differential on the straight side of the element in the computational domain can now be transformed into isoparametric coordinates as 5.21 where lB is the length of the element side in the computational domain. 5.3 S C H E M E S U S E D T O D I S C R E T I Z E T I M E D E R I V A T I V E S 5.3.1 Predictor-Corrector Methods Predictor-corrector methods are based on numerical quadrature[126]. They are con-sidered part of a family of Adams multipoint methods with explicit formulae for prediction of a variable at tn + 1 and implicit formulae for correction of the predicted variable. The 66 time derivatives given by equations 5.8, 5.9 and 5.10 are used to predict and correct values of the corresponding variable at t — tn + 1 based on values at previous time steps. We then write[127] • explicit /predicting formulae m y\+l = Vi + E^mfc/(si-k+i>!K-fc+i) 5 2 2 fc=l • implicit/correcting formulae m y\+l = yi + ^_^f{xi-h+i.yi-k+i) 5.23 k~o where m refers to the number of values involved in the prediction/correction process. Ta-bles 5-1 and 5-2 summarize the coefficients used in equations 5.22 and 5.23. Equations of equal m have been used. Use of equation 5.22 with the coefficients given in Table 5-1 (P) allows for the calculation of y\+l- This value is then used in the location corresponding to the coefficient /3^ 0 from Table 5-2 to find y)+[. Further correction is possible by introducing yl^l into 5.23 to obtain an updated corrected value of the variable being calculated. The equations to be solved succesively for the solution of our initial-boundary value problem are then U = I mk 0 ( y - « o ( * ) ) n + + dn J tn-(*- l ) A t ony - c(z,t) 5.24 67 Table 5-1 Explicit Formulae k = l k=2 k=3 k=4 Error 1 — — — 3 2 l 2 — — &h?fii{xi,yi) 23 12 16 12 5 12 — |fc4/Wf(».-,w) 55 24 59 24 37 24 9 24 Table 5-2 Implicit Formulae k=0 k = l k=2 k=3 k=4 Error Plk l 2 1 2 — — — -$h3f{*i,Vi) P2k 5 12 8 12 1 12 — — -T2h*r{xi,yi) 9 24 19 24 5 24 1 24 — 251 720 646 720 264 720 106 720 106 720 r m 8 a* a* 1 *n-(*-l)A* + " ( D « ? { 5 > + }„} »»= «('.*) 5.25 ( m 8 „ , „ , > < n - ( k - l ) A t + ^ X > - { £ * + } „ } - » = * ) 5.26 Using the predicted values obtained from 5.24, 5.25 and 5.26, we solve 3.28 as 68 — - + /•-• '12- + + 1^21-^/ + J r ; '22 *n + l J 2 i ~ + J22-^\} det[J]aXdr) = dr) J *=B{y,t) 5.27 Using values obtained from the solution of 5.27, we can correct the solution on the free surface using the following expressions (C7)Wl t n (C) mfe .fc=0 g(y - ?o(t))+ + + ~dJ + Ju-drT dr) tn-(k-l)At ony — q(z, t) 5.28 5.29 U=0 j=l <> 5.30 W i t h values from 5.28, 5.29 and 5.30, we now solve 5.27 again to obtain corrected values in the whole domain. The correction formulae can be re-applied until the free surface has 69 converged within a specified tolerance, i.e. JCimi - ^) mi - l <€ p s 5.31 where mi is the number of corrections performed. The value eps has been taken as , At eFS = 1 ( T 3 - 5.32 with Tmax a s the maximum simulation time. 5.3.2 Semi-implicit Semi-Lagrangian Scheme Equations 5.3 and 5.4 define the trajectories followed by fluid particles located on the fluid-air interface . In order to integrate these equations, an iterative scheme is implemented to obtain vectorial position and value of the velocity potential at each fluid particle on this boundary. We can express the integration of a quantity F as r t + A t At Ft+1 = Ft + J Fdt = F* + At Ft+^r 5.33 t . At In order to evaluate F + 2 we use a two-time level extrapolation scheme as[94,95, 96,97,98] Ft+ 2 — - i f * — -F 5.34 We then write this expression for the first increment of each Lagrangian derivative on the free surface in the following extrapolative form 70 ( 0 ) ' < + ^ = A t { - ^ (o),t+4A .2Lc9y 12 L<9z 3 \dcb^ * _ 1 i ~ 2 * _ I i _ 2 1 2 12 i dt J i L dt \ i / 5.35 5.36 5.37 where the subindex i refers to a fluid particle on the free surface. Now, with an initial value given by the equations above we iterate the following scheme k = 1, ...,ra (*),(+At _ t 1 (*_i),t +At *i — z i "r 2 * ( f c ) j t + A i ^ 1_ (k_1),t+At 5.38 + i ,+i 8 A* Q # <9$,-irf 0*< -+ •^ 21-^ 7- + ^ 2 2 " }det[mdn = j f * , ^ d r ; The iteration is continued until (*),i+At _ ( f t - l ) , i + A t / ? J * ) > t + A < _ / 3 ( * - l ) , t + A t < ^ 71 (k),t+At _ (*-l),t+At where e a , and eT are tolerances to be specified. In this work, we have used _ K At ea = eT = I O " " 6 — -When convergence has been achieved at the m iteration, the free surface is updated as vl+At = vl + ft )>t+" ^ = ^ + T H ^ A t In some instances numerical dissipation was applied to maintain stability of the nu-merical scheme. This was performed only on the velocity potential <pl+At using a modified Lax method[128] as ti+At = + d A < ) + ( i - < ^ + A t with a specified as the minimum value that would result in a stable run. Note that for a — 0 no dissipation is introduced and for a = 1 we have the traditional Lax method. 5.4 C A L C U L A T I O N O F F O R C E S , W O R K A N D E N E R G Y The hydrodynamic pressure is calculated using Bernoulli's equation f dd> 1 1 Phyd = p i - — - -V<pV<pj 5.39 72 1 The hydrodynamic force in heave for surface-piercing and submerged bodies as pre-sented in Chapter 6 is then Hi*) = Jphyd nydr 5.40 = J Phyd nzdT 5.41 and the work done by this force is Wi(t)- f F((t)ybdt 5.42 Jo WjJ(i) = f' F{(t)zhdt 5.43 Jo with prescribed motion functions given as yb — Acoswt 5.44 ib — Acosut 5.45 where A and u are the amplitude and frequency of motion respectively. The work done by the motion of the body up to time t is responsible for the change of energy in the fluid. It can be expressed as W(t) = L\KE[t) + L\PE(t) 5.46 where kinetic energy K E and potential energy P E can be calculated as [58,99] KE n 2 73 PE = \pg J(<;-<;o)2dz 5.48 5.5 A L T E R N A T I V E F O R M F O R T H E C A L C U L A T I O N O F V E L O C I T I E S A n alternative method for calculating velocities on the free surface and on the whole domain to that given by 5.8 and 5.9 can be stated as (u'*<) = (^*<) 5.49 = (|^i*0 5.50 dy where the left side yields the mass matrix and the right side is a vector composed of known quantities. Solutions of the systems of algebraic equations resulting from the discretizations pre-sented in this chapter have been performed using a preconditioned conjugate gradient method developed by Chandra[129]. Details on this procedure can also be found in the text by Golub and Van Loan[130]. Solutions of all matrices involved in the applications solved in Chapter 6 have been carried out in double precision using a convergence criteria of 10~15. 74 C H A P T E R S I X R E S U L T S A N D D I S C U S S I O N In this chapter, we present results for a number of applications using the concepts developed in previous chapters. Section 6.1 presents fluid response due to periodic motion of surface piercing and submerged bodies in a tank using linear free surface boundary conditions. Physical and computational aspects are studied using linear and quadratic isoparametric elements. In Section 6.2 these problems are analyzed with nonlinear free surface boundary conditions. Section 6.3 addresses the problem of fluid response due to impulsive motion of basins of various shapes. Nondimensional natural frequencies are obtained for linear and nonlinear free surface boundary conditions. Finally, in Section 6.4 we treat the problem of a mathematically defined slender ship advancing with forward speed in otherwise calm water. Linear and nonlinear free surface cases are studied and their relative merits are discussed. 6.1 P E R I O D I C M O T I O N IN A T A N K : L I N E A R F R E E S U R F A C E In this section, we present linear free surface results for forced periodic motion of a surface-piercing composite cylinder with a vertical side and a semi-circular bottom of unit radius. Results for a submerged cylinder of unit radius are also given. The equations presented in section 5.1 are linearized by removing all nonlinear terms from the boundary conditions on the free surface. Furthermore, all equations are solved at 75 the equilibrium position defined at t = 0. In other words, all calculations are carried out at all times on the mesh defined by the initial conditions. A l l results are non-dimensionalized using g (acceleration of gravity), p (density) and the radius of the circular sections. We have used two different amplitudes of body motion, A — 0.1 and A = 0.5, and a frequency u — 1.25, to facilitate comparison of results with those reported by Yeung and Wu[58]. A l l results are given for the entire domain unless otherwise stated. 6.1.1 Surface-Piercing Cylinder Figure 6-1 shows a typical configuration for the mesh corresponding to the surface-piercing composite cylinder. The variable NI appearing in this chapter corresponds to the number of divisions used to discretize the free surface at one side of the body and half'of the body surface. Apart from this, all graphs and tables are self-explanatory. A l l results in the following subsection are for m — 1 as defined in Section 5.3.1. 6.1.1.1 Physical Aspects Figures 6-2 and 6-3 show time histories of the free surface elevation produced by vertical forced excitation of the composite cylinder shown in Figure 6-1 using linear isopara-metric elements for A — 0.1 and A — 0.5. Figures 6-4 and 6-5 show corresponding results for quadratic elements. In both cases vertical elevations have been magnified 10 times to enhance possible differences. As noted by Yeung and Wu[58], the free surface in both cases remains regular. However, after the first few cycles of motion there are qualitative differ-ences in the free surface profiles. As shown in Figures 6-3 and 6-5, this effect is particularly noticeable towards the end of the simulation run. While quadratic elements tend to predict a smoother y — z contour of the free surface, linear elements show the appearance of a hump 76 that is clearly noticeable at the end of the run. Figures 6-6 and 6-7 show hydrodynamic force applied by the fluid on half the body for A — 0.5. These curves appear visually identical to those given in [58]. Both linear and quadratic elements essentially predict the same behavior and amplitude. Calculation of 5.46 provides a means of assessing time histories of total energy in the fluid and work done by the body. Figures 6-8 and 6-9 show such curves for linear and quadratic elements respectively. Results for both types of elements are practically the same indicating good conservation of energy. As pointed out by Yeung and Wu[58] there are two frequencies involved in the oscillation exhibited by energy-work curves. One of them takes place with a frequency equal to twice the frequency of motion. The other, typical of 'beating' phenomena, occurs with a frequency equal to the difference of the natural frequency of the tank and the frequency of motion. The amplitude of this particular oscillation depends on the amount of potential energy contributed by the free surface elevation. 6.1.1.2 Computational Aspects In order to assess the convergence and stability characteristics of the numerical scheme, we can define energy error as the difference between total energy in the fluid and the work done by the body. The result is normalized by the energy content in the fluid at t — 0+. Figures 6-10 and 6-11 show the effect on the energy error due to reductions in time-step and grid size. As expected, the difference between linear and quadratic elements is significant. Figure 6-10 shows that for linear elements the error, although low, tends to increase with time. This problem can be overcome by increasing NI. This is the case when NI=30 where the error remains low and bounded as time progresses. In the case of linear elements, decreasing the time-step has no effect in reducing the negative effect mentioned above. Figure 6-11 shows that for quadratic elements, the energy error is very low and bounded. This is the 77 case even for a small number of elements, NI=10, and regardless of the time-step used in this work. Stable behavior was observed even for very large Tmax. As a check, we used Tmax — 500 and At = 0.1 and the error time history still remained bounded. The R M S error for large Tmax tended to be smaller than for Tmax — 50. Figures 6-12 and 6-13 show root-mean square of the error time histories, such as those of Figures 6-10 and 6-11, as a function of time-step and number of nodes. Figure 6-12 shows the effect of decreasing At for a half-domain of 441 linear-element nodes and 341 quadratic-element nodes. For this case, Yeung and Wu used 441 nodes in their finite difference discretization. Figure 6-13 shows results for At = 0.1 as a function of the number of nodes in half of the domain. As expected, significant improvement is obtained for quadratic elements. This can be viewed as two-folded. The accuracy of the results is improved and obtained with time-steps and grid sizes significantly higher than for linear elements. In summary, we observe that quadratic elements lead to more accurate and stable solutions with reduced computational effort. The effect of the order of the predictor-corrector method on the R M S error was also studied. We have identified these cases by the variable m given in Tables 5-1 and 5-2. Ta-bles 6-1 and 6-2 show results for linear elements with velocities calculated as in Section 5.5. The variables LCG and LCGQ i n all tables refer, respectively, to the one dimensional array containing the global stiffness matrix with and without removal of all zeros in this matrix. It is observed that higher order predictor-corrector schemes do not improve the R M S error for single correction. For this reason, predictor-corrector methods of order higher than 5th (m > 4) were not used. In fact, for larger At ' s the best results are obtained for m = 1. As the time-step is reduced, R M S error for all values of m tends to the value given by m = 1. We also observe that, for velocities calculated as in Section 5.5, the average number of iter-ations per time-step of the conjugate gradient method, N I T C G , remains constant for all m, although it decreases slightly with the time-step. Tables 6-3 and 6-4 show a similar trend 78 R M S error and number of iterations of conjugate gradient method S U R F A C E - P I E R C I N G C Y L I N D E R NI=20 w=1.25 A = 0.5 Linear Isoparametric Elements Velocity Calculated as in Section 5.5 Case DT=0.025 DT=0.050 DT=0.075 DT=0.100 m = l N I T C G 2.08 109 2.17 112 2.28 115 2.36 116 m=2 N I T C G 2.09 109 2.22 112 2.45 114 2.76 116 m=3 N I T C G 2.08 109 2.19 112 2.35 114 2.52 116 m=4 N I T C G 2.08 109 2.19 112 2.35 114 2.51 116 LCG=4121 LCG 0 =19803 N O D E S = 861 E L E M E N T S = 8 0 0 for quadratic elements. R M S error values with this type of element are approximately 2.5 to 3 times smaller than those obtained using linear elements for a similar number of nodes. Tables 6-5 , 6-6 and 6-7 show results for quadratic elements for velocities calculated as in Section 5.2.1. This form gives higher R M S errors for low numbers of elements (NI=10) (see Table 6-3 and Table 6-5). For higher numbers of elements, R M S error values are slightly lower than those obtained with velocities calculated as in Section 5.5 (see Table 6-4 and Table 6-6). There is a computational advantage in the form of Section 5.2.1 with respect to that given in Section 5.5. The form given in Section 5.2.1 does not require the solution of the mass matrix at every time-step. Calculations are only performed at a very small percentage of the total number of nodes, i.e., those at the free surface. Further advantages can be visualized in Figures 6-14 , 6-15 , 6-16 and 6-17 . These are graphical representations of Tables 6-3 and 6-5. Figures 6-14 and 6-15 show that N I T C G decreases considerably as 79 R M S error and number of iterations of conjugate gradient method S U R F A C E - P I E R C I N G C Y L I N D E R NI=30 w=1.25 A = 0.5 Linear Isoparametric Elements Velocity Calculated as in Section 5.5 Case DT=0.025 DT=0.050 DT=0.075 DT=0.100 m = l 1.03 1.13 1.26 1.37 N I T C G 170 176 179 181 m=2 1.03 1.17 1.39 1.70 N I T C G 170 176 179 181 m=3 1.03 1.14 1.30 1.48 N I T C G 170 176 179 181 m=4 1.03 1.14 1.30 1.47 N I T C G 170 176 179 181 LCG=9181 LCG o =62403 N O D E S = 1891 ELEMENTS=1800 m is increased for results as in Section 5.2.1 while it remains constant when using the form given in 5.5. Figures 6-16 and 6-17 show N I T C G as a function of the time-step size. Greater reduction in the number of iterations can be seen when decreasing A t using equations in 5.2.1. A l l R M S error values in Tables 6-3 to 6-7 are within 1.5%. It can then be concluded that for linear free surface calculations acceptable results can be obtained with considerable computational savings when using equations as in Section 5.2.1 with the higher m values. The effect of multiple correction on the R M S error was also studied. Multiple correction was continued until condition 5.32 was satisfied. Table 6-8 shows results corresponding to some of the cases given in Table 6-5. It can be seen-that multiple correction practically removes the small differences in R M S arising from different values of m. The trend displayed by N I T C G in previous tables as a function of m and A t is essentially the same but more obvious than for single correction. This is due to the fact that as the number of corrections is incremented, the number of iterations for convergence of the conjugate gradient method 80 Table 6-3 R M S error and number of iterations of conjugate gradient method S U R F A C E - P I E R C I N G C Y L I N D E R NI=10 w=1.25 A = 0.5 Quadratic Isoparametric Elements Velocity Calculated as in Section 5.5 Case DT=0.025 DT=0.050 DT=0.075 DT=0.100 m = l 0.62 0.72 0.86 0.98 N I T C G 123 126 129 131 m=2 0.62 0.75 0.98 1.29 N I T C G 123 126 129 131 m=3 0.62 0.72 0.89 1.06 N I T C G 123 126 129 131 m=4 0.62 0.72 0.89 1.06 N I T C G 123 126 129 131 LCG=5151 LCG 0 =23135 N O D E S = 661 E L E M E N T S = 2 0 0 is reduced. Satisfaction of 5.32 is then approached more quickly. In general, not more than 3 corrections were necessary for 5.32 to be satisfied. 6.1.2 Submerged Cylinder Figure 6-18 shows a typical configuration for the mesh corresponding to the submerged cylinder. For this case, the cylinder radius is unity. The variable NI appearing in this section corresponds to the number of divisions used to discretize the free surface and a quarter of the body surface. Based on the experience gained in the case of the surface-piercing cylinder, in this section we have opted to use quadratic elements only. Velocities on the free surface are calculated as in Section 5.2.1. A l l results in the following subsection are for m = 1 81 Table 6-4 R M S error and number of iterations of conjugate gradient method S U R F A C E - P I E R C I N G C Y L I N D E R NI =20 w=1.25 A = 0.5 Quadratic Isoparametric Elements Velocity Calculated as in Section 5.5 Case DT=0.025 DT=0.050 DT=0.075 DT=0.100 m = l 0.18 0.28 0.43 0.56 N I T C G 251 260 265 269 m=2 0.18 0.32 0.57 0.90 N I T C G 251 260 265 269 m=3 0.18 0.29 0.47 0.66 N I T C G 251 260 265 269 m=4 0.18 0.29 0.47 0.65 N I T C G 251 260 265 269 LCG=20301 LCG 0=163865 N O D E S = 2521 E L E M E N T S = 8 0 0 6.1.2.1 Physical Aspects Figures 6-19 and 6-20 show free surface time histories for A — 0.1 and A — 0.5 respectively. In this case, vertical dimensions have also been magnified 10 times. The free surface remains smooth for both cases, although there seems to be an increase in elevation as time becomes relatively large for A — 0.5. As expected, the vertical elevation remains lower than that obtained for the case of the surface-piercing cylinder. The hydrodynamic forces acting on half the body are shown in Figure 6-21 . The same sinusoidal nature as that given by the surface-piercing body is observed. Figure 6-22 shows time records of the total energy in the fluid and the work done by the body. In this case, the 'beat' phenomena is greatly reduced with respect to that of Figures 6-8 and 6-9. This is partly due to a significantly lower content of potential energy arising from the free surface movement. 82 Table 6-5 R M S error and number of iterations of conjugate gradient method S U R F A C E - P I E R C I N G C Y L I N D E R NI=10 « = 1 . 2 5 A = 0.5 Quadratic Isoparametric Elements Velocity Calculated as in Equation 5.2.1 Case DT=0.025 DT=0.050 DT=0.075 DT=0.100 m = l 0.71 0.81 0.95 1.07 N I T C G 100 111 117 121 m=2 0.71 0.85 1.08 1.39 N I T C G 67 80 102 109 m=3 0.71 0.81 0.98 1.15 N I T C G 62 67 74 82 m=4 0.71 0.81 0.97 1.14 N I T C G 61 64 66 71 LCG=5151 LCG 0 =23135 N O D E S = 661 E L E M E N T S = 2 0 0 6.1.2.2 Computational Aspects Figure 6-23 shows error curves as a function of time-step and grid size. The error is very low throughout the run. However, there seems to be some increase in the error amplitude at relatively large time. We conjecture that this is related to the increase in wave amplitude mentioned above. This effect can be controlled by increasing the number of divisions N I on the free surface (NI=20, At=0.05). Tables 6-9 and 6-10 show the effect of m on the R M S error. In this case R M S error values can be kept below 1%. Trends regarding the effect of m on the R M S error and N I T C G observed for this case are similar to those obtained for the case of the surface-piercing cylinder. 83 R M S error and number of iterations of conjugate gradient method S U R F A C E - P I E R C I N G C Y L I N D E R NI =20 w=1.25 A = 0.5 Quadratic Isoparametric Elements Velocity Calculated as in Equation 5.2.1 Case DT=0.025 DT=0.050 DT=0.075 DT=0.100 m = l 0.17 0.27 0.42 0.55 N I T C G 187 224 238 246 m=2 0.17 0.32 0.57 0.89 N I T C G 136 158 186 218 m=3 0.17 0.28 0.46 0.65 N I T C G 126 135 150 162 m=4 0.17 0.28 0.46 0.64 N I T C G 126 130 135 144 LCG=20301 LCG 0=163865 N O D E S = 2521 E L E M E N T S = 8 0 0 6.2 P E R I O D I C M O T I O N IN A T A N K : N O N L I N E A R F R E E S U R F A C E In this section, we present nonlinear free surface results for forced periodic motion of the same cases treated in Section 6.1. A l l boundary conditions are imposed at the exact instantaneous position of all moving boundaries. 6.2.1 Surface-Piercing Cylinder In this section we present results due to prescribed periodic motion of the composite surface-piercing cylinder of Figure 6-1 with nonlinear free surface boundary conditions. Again all results have been non-dimensionalized using density, acceleration of gravity and body radius. For nonlinear boundary conditions, free surface fluid particle paths tended to drift 84 Table 6-7 R M S error and number of iterations of conjugate gradient method S U R F A C E - P I E R C I N G C Y L I N D E R NI =30 w=1.25 A = 0.5 Quadratic Isoparametric Elements Velocity Calculated as in Equation 5.2.1 Case DT=0.025 DT=0.050 DT=0.075 DT=0.100 m = l 0.09 0.20 0.36 0.49 N I T C G 265 332 352 366 m=2 0.09 0.25 0.51 0.83 N I T C G 203 234 265 320 m=3 0.09 0.21 0.40 0.59 N I T C G 189 202 223 240 m=4 0.09 0.21 0.39 0.58 N I T C G 189 196 202 215 LCG=45451 LCG 0=530195 N O D E S = 5581 ELEMENTS=1800 away from the body's contour. For the frequency and large amplitude motion parameters used here, particle paths tended to cross over each other at about the third cycle of oscil-lation leading to the break down of the numerical scheme. This negative effect could only be overcome by introducing numerical smoothing. Figure 6-24 shows a plan view of the free surface fluid particles using predictor-corrector schemes with numerical smoothing. In order to keep the scheme stable using these methods, both the free surface elevation and the velocity potential had to be smoothed. This resulted in a free surface elevation with un-realistically smooth characteristics. This is evidenced by the relatively straight trajectories described by the free surface fluid particles. Figure 6-25 shows a plan view obtained using the semi- implicit semi-Lagrangian scheme developed in Section 5.3.2. Free surface particles exhibit relatively periodic trajectories with periods of oscillation that are roughly equal to the prescribed motion period. For this method, numerical smoothing was only required on the velocity potential. Furthermore, this was applied in a very small proportion as described by the modified Lax method mentioned in Section 5.3.2 with a = 0.01. A l l results presented 85 Table 6-8 R M S error and number of iterations of conjugate gradient method S U R F A C E - P I E R C I N G C Y L I N D E R NI=10 w=1.25 A = 0.5 Quadratic Isoparametric Elements Velocity Calculated as in Equation 5.2.1 Multiple Correction Case DT=0.025 DT=0.050 m = l 0.71 0.82 N I T C G 74 80 m=2 0.70 0.81 N I T C G 45 56 m=3 0.70 0.81 N I T C G 41 45 m=4 0.70 0.81 N I T C G 41 42 LCG=5151 LCG 0 =23135 N O D E S = 661 E L E M E N T S = 2 0 0 from this point onwards have been calculated using this method. It was observed that for nonlinear boundary conditions, calculation of velocities on the free surface using the form given in Section 5.2.1 led to significant energy-work error and unstable behavior. A similar effect was observed using linear isoparametric elements. For this reason, all the results for nonlinear free surface boundary conditions were calculated using quadratic elements and velocities as in Section 5.5. 6.2.1.1 Physical Aspects Figure 6-26 and 6-27 show free surface elevations corresponding to w = 1.25 and A — 0.5 through 20 periods of oscillation. As reported in Yeung and Wu[58], the nonlinear free surface tends to have a more irregular nature than that predicted by linear theory. This effect can also be seen in the trajectory described by the fluid particles in contact with the 86 Table 6-9 R M S error and number of iterations of conjugate gradient method S U B M E R G E D C Y L I N D E R NI=10 « = 1 . 2 5 A = 0.5 Quadratic Isoparametric Elements Velocity Calculated as in Equation 5.2.1 Case DT=0.025 DT=0.050 DT=0.075 DT=0.100 m = l 0.35 0.46 0.64 0.81 N I T C G 142 165 190 205 m=2 0.35 0.47 0.65 0.84 N I T C G 119 136 147 156 m=3 0.35 0.46 0.64 0.82 N I T C G 114 123 131 139 m=4 0.35 0,46 0.64 0.82 N I T C G 114 119 122 128 LCG=10251 L C G 0 = 45535 N O D E S = 1301 E L E M E N T S = 4 0 0 side walls of the tank. There is a slight appearance of steeper crests and shallower troughs, typical of nonlinear waves. However, irregularities in the z — y direction for this work are not as pronounced as those obtained by Yeung and Wu[58]. It is conjectured that this may be due to smoothing. Figure 6-28 shows grids at approximately 8.25, 8.5, 8.75 and 9th cycles and can readily be compared with that given by Yeung and Wu[58]. It can be seen from Figure 6-28 that the mesh generation studied by All ievi and Cal-isal[92,93] adapts very well to the significant changes occurring in the physical domain due to the large displacements of the composite cylinder. It should be mentioned at this point that mesh regeneration throughout this work has been performed using the generating matrices for the y and z coordinates resulting from the initial conditions. In other words, calculation of the elements of the matrices defined by 4.29 and 4.30 for mesh regeneration purposes is performed once at t = 0. Solution of this matrix with moving boundary conditions delivers the moving meshes of Figure 6-28. The computational effort required for mesh regeneration 87 Table 6-10 R M S error and number of iterations of conjugate gradient method S U B M E R G E D C Y L I N D E R NI=20 w=1.25 A = 0.5 Quadratic Isoparametric Elements Velocity Calculated as in Equation 5.2.1 Case DT=0.025 DT=0.050 DT=0.075 DT=0.100 m = l 0.13 0.25 0.43 0.61 N I T C G 280 325 359 397 m=2 0.13 0.26 0.44 0.64 N I T C G 241 272 295 315 m=3 0.13 0.25 0.44 0.62 N I T C G 230 250 267 282 m=4 0.13 0.25 0.44 0.62 N I T C G 230 246 255 264 LCG=40501 LCG 0=325065 N O D E S = 5001 ELEMENTS=1600 is therefore greatly reduced to within 10% of the total computational effort per time step. As forseen in [93], good quality meshes for subsequent engineering computations can be ob-tained with very few cycles of iteration. In this work, no iterations were required to obtain the sequence of deforming meshes shown in Figure 6-28. These meshes are seen to lead to accurate calculations of physical phenomena with a significant reduction of the required computational effort. Figure 6-29 shows velocity vector plots for the configurations shown in Figure 6-28. Velocity fields display a smooth trend throughout the domain following the motion of the moving boundaries in a physically acceptable manner. As a check, we have monitored velocity values at the very bottom of the cylinder. The differences between prescribed velocities and those calculated from the solution were very small. Figure 6-30 shows hydrodynamic forces calculated using the concepts described in 88 Section 5.4. As pointed out by Yeung and Wu[58] for the dimensions and large oscillations described in this case, nonlinear theory predicts a larger negative peak and a smaller positive one with respect to the sinusoidal prediction of linear theory. However, for a given frequency of motion this effect is essentially a consequence of the relative dimensions of the tank and the moving body. Figure 6-31 shows the effect of increasing the tank dimensions by 2 and 3 times while keeping the body dimensions constant. The resulting effect is two-folded. First, a significant reduction of the amplitude of the hydrodynamic force is observed. Also, a decrease in the difference between positive and negative peaks of this force occurs. Figure 6-32 shows total energy-work records for amplitude of motion A — 0.5 and frequency w — 1.25. The difference between the two curves is barely noticeable displaying good energy conservation characteristics of the numerical procedure. Again, in this case we notice a more irregular trend with respect to that presented in Figure 6-9 for linear free surface boundary conditions. 6.2.1.2 Computational Aspects Figure 6-33 shows energy error curves for nonlinear boundary conditions. These curves show a very distinct trend compared to those obtained in Figure 6-11. It was observed that the energy error tended to increase towards the end of the run regardless of the time-step used. Again, this negative effect can be overcome by increasing the number of elements (NI=20, At=0.05). The convergence tendencies obtained for linear free surface boundary conditions are again observed for the nonlinear free surface case. Figure 6-34 shows R M S error of records such as those shown in Figure 6-33 as-a function of time-step (NI=10) and the number of nodes (At=0.05) for half the domain. Fine convergence is achieved in both cases. It should be noted here that the R M S error versus time-step curve is the only curve where a fair comparison with Yeung and Wu[58] can be made. In this case, R M S values are below those reported by Yeung and Wu[58]. For the figure corresponding to R M S error 89 versus the number of nodes, we have used At = 0.05. A time-step of At = 0.1 was used by Yeung and Wu[58]. Therefore, this curve should only be taken as a reference. For At = 0.1, N I — 30 and a = 0.01 the scheme broke down at about the seventh cycle of body motion. 6.2.2 Submerged Cylinder In this section, we present nonlinear free surface results corresponding to the config-uration shown in Figure 6-18 where the circular cylinder is made to move with prescribed heave motion. The variable NI corresponds, as before, to the number of divisions used to discretize the free surface and a quarter of the body surface. No smoothing was necessary for this case. Figures 6-35 and 6-36 show nonlinear free surface evolution for A — 0.1 and A — 0.5 respectively. Steeper peaks were observed to develop i n the regions next to the tank walls. Linear results predicted an almost sinusoidal behavior of the free surface across the width of the tank, Figure 6-20. Figure 6-37 shows mesh deformation corresponding to the 9.25, 9.5, 9.75 and 10th cycles of oscillation. Figure 6-38 shows corresponding velocity vector plots. Again, velocity fields behaved in a physically acceptable manner with very small differences between specified and calculated values. For A — 0.5, linear and nonlinear prediction of the hydrodynamic forces resulting from vertical motion of the circular cylinder displayed similar trends. Figure 6-39 shows linear and nonlinear hydrodynamic forces corresponding to an amplitude of motion A — 1.5. Even for this extremely large amplitude of motion, both formulations gave similar results. Some differences can again be seen in the force peaks, but these are significantly less noticeable than for the surface piercing case. Figure 6-40 provides time records of the total energy in the fluid and the work done 90 by the body. Again, the 'beat' phenomena is greatly reduced with respect to that of the surface-piercing cylinder. Very little difference is observed with respect to the prediction given by linear free surface. Figure 6-41 shows error curves as a function of time-step and grid size. Again, the trend of these curves is quite different to that shown in Figure 6-23. Although larger than that obtained with linear free surface boundary conditions, the energy error is low throughout the run. It can also be observed that the energy error tends to have a slightly increasing tendency, regardless of the time-step used. Again, this is removed by increasing the number of elements on the free surface with a considerable reduction in the energy error throughout the run. 6.3 I M P U L S I V E M O T I O N IN A T A N K In this section, we present results corresponding to free surface response due to im-pulsive excitation of tanks of different shapes. Excitation is prescribed as a step velocity at t — 0. When appropriate, we have opted to match the various types of excitations given by Yeung and Wu[58,62]. These are a. Square basin: vertical impulsive motion of cylinder. b. Circular, parabolic and triangular basins: symmetric 'squashing' of the tank walls. c. El l iptic basin: anti-symmetric roll motion of tank. Application of Fast Fourier Transform to the time record of the free surface elevation mea-sured at a fixed point yields the frequency spectrum. Nondimensional frequencies for linear and nonlinear free surface boundary conditions are presented. Smoothing in these cases was applied using a = 1 every 10th time-step. 6.3.1 Square tank Figure 6-42 shows frequency spectra corresponding to the tank of Figure 6-1 . Frequen-91 Table 6-11 Natural freque ncies for square tank N A T U R A L F R E Q U E N C I E S O F S Q U A R E T A N K Case L i n . Nonl. Y W . L in . Y W . Nonl W L 1st 1.0124 1.0124 1.0095 1.0095 1.0231 2nd 1.4112 1.4112 1.4072 1.4072 1.4472 3rd 1.6566 1.6873 1.6519 1.6825 1.7725 cies for linear and nonlinear cases do not differ significantly. Table 6-11 shows comparison of results with those given by Yeung and Wu[58] ( Y W ) and Wehausen and Laitone[131](WL). Note that the results presented by Yeung and Wu[58] are for a tank exactly half the size of that given by Figure 6-1. Results by Wehausen and Laitone[131] are for a rectangular basin of dimensions 3 in width by 4 in depth. Frequencies predicted due to symmetric excitation for the entire tank are exactly the same as those given for half the tank by Yeung and Wu[58]. This is because the vertical line extending from the bottom of the cylinder to the bottom of the tank is a line of symmetry. The flow is symmetric about this line. This then acts as a rigid wall, virtually splitting the whole tank into two independent tanks located at both sides of the contour defined by half the body and the line of symmetry. 6.3.2 Circular tank Figure 6-43 shows results for a circular tank. In this case, frequencies predicted by linear theory tend to be slightly higher than those giyen by nonlinear theory. Table 6-12 shows comparison with results given by Chang and Wu[68] (CW) using an integral equation formulation with linear boundary conditions. Their results also exhibit higher values than those predicted by nonlinear theory. 92 Table 6-12 Natural frequencies for circular tank N A T U R A L F R E Q U E N C I E S O F C I R C U L A R T A N K Case L i n . Nonl. C W 1st 1.7330 1.7180 1.7492 2nd 2.4697 2.4545 2.5632 3rd 3.0066 2.9759 Table 6-13 Natural frequencies for parabolic tank N A T U R A L F R E Q U E N C I E S O F P A R A B O L I C T A N K Case L i n . Nonl. Y W . L i n . Y W . Nonl 1st 1.6720 1.6720 1.6846 1.6846 2nd 2.4390 2.4390 2.4170 2.4170 3rd 3.0066 3.0066 3.0029 3.0029 6.3.3 Parabolic tank Corresponding results for a parabolic basin are given in Figure 6-44 . For this particular configuration, results tend to show no difference between linear and nonlinear predictions. Table 6-13 compares our results with those given by Yeung and Wu[62]. In general, very good agreement is achieved. 6.3.4 Triangular tank Figure 6-45 shows results for a triangular basin. For this case, no definite trend is shown between linear and nonlinear results. As noted in Table 6-14 , this effect is also shown i n the results of Yeung and Wu[62]. Results by Lamb[132] are also presented. The first natural frequency was given by Ghanimati and Naghdi[133] to be 1.548 using fluid-sheets 93 Table 6-14 Natural frequencies for triangular tank N A T U R A L F R E Q U E N C I E S O F T R I A N G U L A R T A N K Case L i n . Nonl. Y W . L in . Y W . Nonl Lamb 1st 1.5339 1.5646 1.5381 1.5381 1.5244 2nd 2.3930 2.3930 2.3439 2.3438 2.3447 3rd 3.0066 2.9912 2.9297 3.0029 2.9393 theory. 6.3.5 Elliptic tank Figure 6-46 shows results for an elliptic basin under anti-symmetric excitation. As noted by Yeung and Wu[58], only nonlinear free surface results display frequencies associated with symmetric modes. Linear results disclose only frequencies associated with the type of excitation used. Table 6-15 shows comparison of results with Yeung and Wu[58] and with Wehausen and Laitone[131] for a square tank of 2x0.75 width to depth ratio. Again, good agreement is obtained. Contrary to the results of Yeung and Wu[62], the energy associated with even frequencies is not small. It is at least of the same order as that predicted by linear theory for the energy associated with anty-symmetric modes. 94 Table 6-15 Natural frequencies for elliptic tank N A T U R A L F R E Q U E N C I E S O F E L L I P T I C T A N K Case L i n . Nonl. Y W . L in . Y W . Nonl W L 1st 1.1044 1.1044 1.0986 1.0986 1.140 2nd 1.7180 1.758 1.757 3rd 2.1322 2.1168 2.124 2.124 2.169 4th 2.4850 2.490 2.506 5th 2.7611 2.7611 2.783 2.783 2.802 Table 6-16 Wigley H u l l Information W I G L E Y P A R A B O L I C H U L L Length L[m] Beam B[m] Draft D[m] CB Cwp Cp 2 0.2 0.125 0.444 0.667 0.667 6.4 S L E N D E R B O D Y W I T H F O R W A R D S P E E D I N A T A N K In this section we present results corresponding to the forward motion of a mathe-matically defined ship hull advancing at zero angle of attack in otherwise calm water. The Wigley parabolic hull , whose main parameters are given in Table 6-16 , is defined as Coordinates z and y remained as defined in previous sections for surface-piercing cases. The coordinate x is directed along the length of the hull and considered positive towards the stern. Comparison wil l be made with numerical and experimental results obtained at Yoko-95 hama National University. These were reported by Maruo and Song[40]. The wave resistance coefficient in this work was defined as Data given by Aanesland[134] and experimental data from the 1st and 2nd Workshops on Wave Resistance Computations are also presented[135,136j. The wave resistance coefficient in these references was defined as where S is the wetted surface at the design waterline. The value R,j, for wave resistance was calculated by first integrating nodal values of the longitudinal component of pressure around each ship station. This value depends on the slope of the waterplane and the depth of the waterplane. This calculation gave sectional resistance forces. Integration of these forces along the length of the ship hull provided the value of Rw. Sinkage and tr im were not included in the calculation of wave resistance. Results presented for tr im, sinkage and wave-resistance coefficients have been obtained for the conditions shown in Table 6-17 . The tank width is in metres. The variable A x refers to the spatial increment used along the ^-coordinate in metres. For the length of the Wigley hull used in this study, this increment resulted in 2000 stations. C P U time estimates are presented in Table 6-18 . This is the time required for • Mesh generation after extrapolation given by 5.35, 5.36 and 5.37 • Iterative procedure given by 5.38. This involves solution of both stiffness and mass 96 Table 6-17 Computational Conditions for Wigley Hul l Calculations C O M P U T A T I O N A L C O N D I T I O N S F O R W I G L E Y H U L L Aa; Tank Width No of elements No of Nodes 0.001 3.0 1200 3741 Table 6-18 Computational Demands for Wigley Hul l Calculations C O M P U T A T I O N A L D E M A N D S F O R W I G L E Y H U L L Machine Type Mflops C P U [min/As] S U N MP/670 4.0 5.0 I B M RS/6000 31.0 0.75 matrices • Mesh generation at end of iterative procedure 5.38 • Solution of stiffness and mass matrices with converged mesh from the step above • Calculation of sectional forces Every solution of the stiffness matrix is carried out with an exact three-dimensional body boundary condition. Figures 6-47 and 6-48 show linear and nonlinear three dimensional views of the free surface elevation for the Wigley hull advancing with a forward speed of 1.182j^. This corresponds with a length Froude number of 0.267. It is clear that nonlinear theory predicts a more realistic shape of the free surface. This is evident in the aft-half of the hull where linear theory is not able to produce the recovery from the second crest of the wave next to the hull . This is more easily seen in Figure 6-49 where comparison has been made with results presented by Maruo and Song[40]. Nonlinear results show good agreement with experimental values for the wave past the midship section. However, there is a shift of the first crest and 97 trough of the fore-half of the wave. It is conjectured that one possible reason for this shift is due to the at rest (zero-valued) initial conditions used in this work. Figure 6-49 shows numerical results given by Maruo and Song[40] which clearly display non-homogeneous initial conditions. Free surface elevation contours are also shown in Figure 6-50 and Figure 6-51 for linear and nonlinear boundary conditions respectively. Nonlinear results show a 'jagged' nature due to small wiggles appearing on the free surface. These wiggles were typical of nonlinear free surface profiles and occasionally were the origin of unstable behavior. In order to maintain the size of these wiggles bounded, small time steps ( A x ) were used. In ship design it is important to estimate fluid resistance in order to determine the size of the power plant required to propel the vessel. Figure 6-52 shows the wave resistance coefficient as a function of the Froude number Fn for linear and nonlinear formulations. These are compared with mean values of the results presented i n the two workshops men-tioned above. A t low values of Fn nonlinear prediction shows significant improvement over results obtained using linear theory. In fact, linear theory gives very poor agreement with experimental results for this F n-range. For this set of experimental data, nonlinear results show superior performance throughout the Fn interval. Figure 6-53 shows the wave resis-tance coefficient compared with experiments from Yokohama National University. Again, at low values of the Froude number nonlinear results give a better prediction than linear ones. However, as the Froude number increases, linear results show improvement over their nonlinear counterpart. This is a clear example of the scatter encountered in experimental data for these type of physical problems. Figures 6-54 and 6-55 show non-dimensional sinkage J and trim t of the Wigley hull as a function of Fn. Improved predictions at the lower end of the Froude number interval used i n this work are again realized for nonlinear results. This improvement is more clearly seen i n the results for sinkage. These results suggest that nonlinear theory predicts a more accurate magnitude of the forces. This is consistent with results obtained for wave resistance 98 as shown in Figure 6-52. Tr im calculations show little discrepancy between linear and non-linear results. In turn, this suggests that the distribution of forces is predicted with similar accuracy in both formulations. 99 M E S H F O R S U R F A C E - P I E R C I N G C Y L I N D E R -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 Z - A X I S Figure 6-1 Typical mesh, for surface-piercing cylinder. Figure 6-2 Linear free surface time record for surface-piercing cylinder ( N I — 30, A t 0.1). 100 Figure 6-4 Linear free surface time record for surface-piercing cylinder (NI = 20, At — 0.1). 101 SURFACE-PIERCING CYLINDER A-0.6 QUADRATIC ELEMENTS Figure 6-5 Linear free surface time record for surface-piercing cylinder ( N I — 20, At 0.1). SURFACE-PIERCING CYLINDER, AMPLITUDE-0.5 LINEAR ELEMENTS Figure 6-6 Hydrodynamic force due to sinusoidal excitation ( N I — 20, At = 0.1). 102 SURFACE-PIERCING CYLINDER, AMPUTUDE-0.5 QUADRATIC ELEMENTS I " " l • 111 • • I i i i • 11 i I i 11 • 11 10 15 20 25 30 35 40 46 50 55 TIME Figure 6-7 Hydrodynamic force due to sinusoidal excitation ( N I = 10, At = 0. SURFACE-PIERCING CYLINDER, AMPUTUDE-0.5 LINEAR ELEMENTS ? 0.12 • l " " l ° 5 10 15 20 25 30 35 40 45 50 55 TIME Figure 6-8 Energy (—) and work (- - -) ( N I = 30, At = 0.1). 103 SURFACE-PIERCING CYLINDER, AMPLITUDE-0.5 QUADRATIC ELEMENTS | i i i i | i n i | 5 10 15 20 2B 30 33 40 46 50 55 TIME Figure 6-9 Energy (—) and work (- - -) ( N I = 20, A i = 0. SURFACE-PIERCING CYLINDER, AMPLITUDE-0.5 LINEAR ELEMENTS 10.00 8.00 •^e.oo IE o n cc LU 4.00 2.00 11 1111 i M ! 11 1 1 ! i • 1 1 ! 111 111111 11111 p i , , | 1 1 1 1 | 1 T l 1 | 1 1 1 1 N12O,0T.O5 • NI20.DT.10 NI30.DT.10 -.-•v\ . " 1 . . . . 1 . . . . I m l m i ' TIME Figure 6-10 Energy error for linear elements. 104 SURFACE-PIERCING CYLINDER, AMPLITUDE-0.5 TIME Figure 6-11 Energy error for quadratic elements. 3.5 • ' 1 • • 1 • | 1 1 I 1 441 N0D..L YHJNQ&WU ! 3.0 • 341 NOa.Q J 2.5 2.0 -1.5 1.0 0.5 0.0 " 0.000 0.025 0.050 0.075 0 1 00 DT Figure 6-12 RMS of energy error versus A i . 105 ' 1 i • > • • i • • • • ; DT=0.1.L YEUNG &WU i : DT=0.1,Q J : \ \ - \ \ \ \ ': . \ ^ . . i . . . . i . . . . j - • 1 .... ' 0 500 1000 1500 2000 2500 0000 NUMBER OF NODES Figure 6-13 R M S of energy error versus number of nodes. 2 d 6 so 25 -_ - - DT-0.025 - DT=0.0S0 DT-0.075 DT=0.100 .... i .... i .... i • 1 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 Figure 6-14 Number of iterations of C G . method versus m (as in Section 5.5). 106 o s s d d so DT-0.025 .- DT-0.050 _ DT=0.07S DT-0.100 i I i i i i I i i i i I i i . i l l °-° 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 Figure 6-15 Number of iterations of C G . method versus m (as in Section 5.2.1). Figure 6-16 Number of iterations of C G . method versus time-step size (as in Sec-tion 5.5). 107 Figure 6-17 Number of iterations of C . G . method versus time-step size (as in Sec-tion 5.2.1). -4.0 -3.0 -2.0 -1.0 0,0 1.0 2.0 3.0 4.0 Y-AXIS Figure 6-18 Typical mesh for submerged cylinder. 108 109 SUBMERGED CYLINDER, AMPLITUDE-0.5 QUADRATIC ELEMENTS 1 1 1 11 ' ' ' 1111 i i 111111 i i i i | | 11111 1 1 1 1 1 1 1 1 1 1.50 • ' • 1 • 1 • • • • 1 • • • I • . . . I I I • I I I I I . I 0 6 10 15 20 25 30 36 40 45 50 55 TIME Figure 6-21 Hydrodynamic force due to sinusoidal excitation ( N I = 20, At = 0. ul a. SUBMERGED CYLINDER, AMPLITUDE=0.5 QUADRATIC ELEMENTS • <l ' 1 1 • 0 5 10 15 20 25 30 35 40 45 50 55 TIME Figure 6-22 Energy (—) and work (- - -) ( N I = 30, At = 0.1). 110 Figure 6-23 Energy error for submerged cylinder. 4.0 3.0 2.0 1.0 3 0.0 -1.0 -2.0 -3.0 -4.0 P R E D I C T O R C O R R E C T O R M E T H O D P L A N VIEW O F NONLINEAR F R E E S U R F A C E FLUID P A R T I C L E S I j f WM •HI mm +--i- + i 1. 10 15 20 25 TIME i—r - i — i -30 35 40 45 50 Figure 6-24 Plan view of free surface particles for predictor-corrector methods 111 SEMI-IMPLICIT SEMI-LAGRANGIAN METHOD 4.0 3.0 2.0 1.0 CO 3 0.0 N -1.0 -2.0 -3.0 -4.0 PLAN VIEW OF NON . -LINEAR III i r f r \ FREE SURFACE FLUID PART l l l l l l CLES IIII -• i i i i i i i i i i i i i i j \ 1 I I I I I I I I I I I I -\ i i i 1HIHHI 1 — • • i i i mmm 10 15 20 25 TIME 30 35 40 45 50 Figure 6-25 Plan view of free surface particles for semi-implicit semi-Lagrangian method. Figure 6-26 Linear free surface evolution up to 20 cycles, NI=10, At=0.1. 112 113 Figure 6-28 Mesh deformation at 8.25, 8.5, 8.75 and 9th cycL 114 Figure 6-29 Velocity vector plots at 8.25, 8.5, 8.75 and 9th cycL 115 HYDRODYNAMIC FORCES. A-0.5. ' ' " 1 ' ' " 1 1 I • • " I i i • • I i i i i | i • • i | i i i i | i 10 15 20 25 30 35 40 45 TIME Figure 6-30 Hydrodynamic forces for surface-piercing cylinder. NONLINEAR HYDRODYNAMIC FORCES VS. TANK DIMENSIONS Figure 6-31 Hydrodynamic forces for surface-piercing cylinder versus tank dimensions. 116 NONLINEAR FREE SURFACE F i g u r e 6-32 W o r k (—) and tota l energy (- -) records, NI=10 , A t = 0 . 1 , A = 0 . 5 , u =1.25. NONLINEAR FREE SURFACE. FREQUENCY-1.25. AMPLITUDE-0 5 ' 1 1 1 I 1 1 1 1 I i i j i • • NI-20.DT-0.050 NI-10.DT-0.050 NI-10.DT-0.025 F i g u r e 6 - 3 3 Energy error versus time-step and number of gr id divisions. 117 NONLINEAR FREE SURFACE — i — i — i — i — i — i — i — p — £ 2.50 f O 341 QUAD. NODES ' YEUNG & WU 1.00 0.50 0.050 TIME STEP * 2.50 O £ 2.00 DT-0.05, QUAD. - YEUNG & WU 1000 1500 2000 NUMBER OF NODES Figure 6-34 RMS of energy error versus time-step (NI=10) and number of nodes (At=0.05). 118 SUBMERGED CYLINDER A-0.1 NONLINEAR FREE SURFACE Figure 6-35 Nonlinear free surface time record for submerged cylinder (NI = 10, At 0.05). SUBMERGED CYLINDER A-O.S NONLINEAR FREE SURFACE Figure 6-36 Nonlinear free surface time record for submerged cylinder (NI = 10, Af 0.05). 119 120 1  • ' 1 ' ' • " 1 1 1 ' ' 1 ' ' * 1 AT TIME- 46.75 • • • •' ' • • • • • . . . . i . . . . Figure 6-38 Velocity vector plots at 9.25, 9.5, 9.75 and 10th cycles. 121 AMPLITUDE 1.5 ' ' LINEAR FREE SURFACE NONLINEAR FREE SURFACE TIME Figure 6-39 Hydrodynamic force due to sinusoidal excitation, A=1.5. SUBMERGED CYLINDER, AMPLITUDE=0.5 NONLINEAR FREE SURFACE I ' • 1 1 I ' ' • • I ' • • ' I ' ' i • I i • i • | • • i i r Figure 6-40 Work (—) and total energy (- -) records, NI=10, At=0.1, A=0.5, u =1.25. 122 -10.0 SUBMERGED CYLINDER. NONLINEAR FREE SURFACE. NI-10, DT-0.050 - NI-10, DT-0.025 NI-20, DT-0.050 40 , M W m k igure 6-41 Energy error versus time-step and number of grid divisions. 123 1.0 0.0 •1.0 -2.0 -3.0 • I 1 1 1 1 1 1 1 ' 1 1 1 • I 1 1 I 1 1 1 1 1 1 I I 1 I I 1 1 SQUARE BASIN '-0_i IS— / l 1 1 / l I 1 /1 1 1 1 -2.0 -1.0 0.0 1.0 2.0 10.0 - - . . i . i i i i 1.0124 1 1 i i i i 1 • i i i 7.5 - LINEAR NONLINEAR " 5.0 -2.5 1.4112 U 1.6566 -0.0 • i • • » . i • ^ i l i L . i i 1:6873 0.5 1.0 1.5 2.0 FREQUENCY Figure 6-42 Free surface frequency spectra due to impulsive motion for square basin. 124 • 1 1 1 1 1 ' 1 1 ' 1 • 1 • 1 1 1 1 1 1 ; CIRCULAR BASIN — i — i i i ' • • — / • l A T - ~ i i 1 , , . 0.400 0.350 1 • • • ' 1 • i — | — • — I — • — I — ; 0.300 1.7333 1.7180 — LINEAR " NONLINEAR • 0.250 f ; 0.200 0.150 2.4697 I 2.4545 0.100 3.0066 0.050 1 2.6759 0.000 I 1 , 4fA' X J L. • : 2.0 3.0 FREQUENCY Figure 6-43 Free surface frequency spectra due to impulsive motion for circular basin. 125 Figure 6-44 Free surface frequency spectra due to impulsive motion for parabolic basin. 126 0.075 - ji 1.5339 • i • — — • — i — | — i — i — i — i — i 1.5646 • LINEAR NONLINEAR 0.050 -• 0.025 j -[ 2.3930 ,' 1 2.3930 J 3.0066 0.000 . . . ...i.._H 2.9912 . . . , . 2.0 3.0 FREQUENCY Figure 6-45 Free surface frequency spectra due to impulsive motion for triangular basin. 127 0.25 - ' ' ' ' 1 1 1 , 1 1 ' ' ' ) i i i j •0.00 E L L I P T I C A L B A S I N -0.26 -0.B0 -0.76 0.75 0.25 0.00 1 1 1 1 1 1 1 1 1 i 1 1 i { 1.1044 i LINEAR NONLINEAR -i i 2.1322 { i | 2.1168 / ; |\ 1.7180 r \ ft—/\ J,\ 2.4850 2 . 7 6 1 1 " 0.0 0.5 1.0 1.5 2.0 2.5 FREQUENCY Figure 6-46 Free surface frequency spectra due to impulsive motion for elliptic basin. 128 129 130 131 CONTOUR PLOT OF FREE SURFACE ELEVATION IN METRES 0.0 05 0.5 0.8 1.0 1.2 1.5 1.8 2.0 Figure 6-51 Nonlinear free surface elevation contours of Wigley hull with forward speed, Fn = 0.267. WAVE RESISTANCE COEFFICIENT OF WIGLEY HULL 1— MEAN OF EXPERIMENTS, MODEL FREE MEAN OF EXPERIMENTS, MODEL FIXED • FEM, LINEAR • FEM, NONLINEAR 1 i i I--/• /' w -• p ... ._ - 1 i' i i i 1 i 0.10 0.15 0.20 0.25 0.30 0.35 0.40 FROUDE NUMBER Figure 6-52 Wave resistance coefficient of Wigley hull versus Fn. 132 WAVE RESISTANCE COEFFICIENT OF WIGLEY HULL • MARUO & SONG, NUMERICAL MARUO & SONG, TOWING TEST • MARUO & SONQ, WAVE ANALYSIS > FEM, LINEAR <> FEM, NONLINEAR 0-15 0.20 0.25 0.30 0.35 FROUDE NUMBER Figure 6-53 Wave resistance coefficient of Wigley hull versus NONDIMENSIONAL SINKAGE OF WIGLEY HUU 2 -0.040 • AANESLAND, NUMERICAL AANESLAND, EXPERIMENTAL (UL) AANESLAND, EXPERIMENTAL (LL) > FEM, LINEAR • FEM, NONLINEAR 0.2 FROUDE NUMBER Figure 6-54 Sinkage of Wigley hull versus Fr 133 NON DIMENSIONAL TRIM OF WIGLEY HULL • AANESLAND, NUMERICAL AANESLAND, EXPERIMENTAL (UL) 0.00 0.05 0.10 0.16 0.20 0.25 0.30 0.35 0.40 FROUDE NUMBER Figure 6-55 Trim of Wigley hull versus Fn. 134 C H A P T E R SEVEN C O N C L U D I N G R E M A R K S en-The boundary value problem defined by the Laplace equation for the solution of pot tial flow in two dimensions has been formulated using a Bubnov-Galerkin weighted residual method to obtain the spatial variation of the velocity potential. The moving boundary prob-lem corresponding to the a-priori unknown dependent and independent spatial variables in the nonlinear boundary conditions has been solved using various numerical schemes to track the free surface in time. Predictor-corrector methods of different orders displayed good con-vergence and accuracy characteristics for linear free surface boundary conditions. However, in nonlinear motion these methods were found to require dissipation for the scheme to be sta-ble. Smoothing of the free surface elevation and the velocity potential led to unrealistically smooth free surface profiles. In order to avoid these inaccuracies, a semi-Lagrangian semi-implicit scheme was developed to track free surface fluid particles and obtain the value of the velocity potential. This method proved to be more accurate and stable than predictor-corrector methods. Furthermore, computations only required dissipation of the velocity potential. The typical number of iterations of the semi-implicit semi-Lagrangian scheme for the types of problems solved in this thesis ranged between 2 and 5. These values were seen to depend on a number of variables such as time-step, grid size, frequency and amplitude of motion, forward speed and tank dimensions. The amount of smoothing required to maintain stability is also dependent on these considerations. Because of the nonlinear nature of free surface flow, it is recommended that numerical experimentation is made for problems not considered in this investigation. The computational conditions given in the text should only 135 be used as general guidelines. A Bubnov-Galerkin procedure has also been used to solve the elliptic generation system developed by Ryskin and Leal[107] for constructing boundary-fitted curvilinear orthogonal coordinates. Numerical aspects have been studied using the geometries and some shape factor distributions proposed by Chikliwala and Yortsos[108] and some proposed here. A number of improvements have been obtained and essentially different conclusions have been reached. Orthogonal grids of very good overall/global accuracy for both Dirichlet and Neu-mann boundary conditions in symmetric and asymmetric regions were obtained. Localized accuracy is lowest for the case of equidistant boundary correspondence in an asymmetric domain. For all other boundary conditions, including complete boundary correspondence with exponential boundary point distribution on one side of the domain, localized accuracy is very good in both symmetric and asymmetric domains. The scheme showed to be stable for all degrees of asymmetry and shape factor distributions used here. Convergence was observed to be independent of these considerations. The number of iterations required for convergence was generally quite low indicating the fast convergence characteristics of the scheme. In most cases, independence of the number of iterations on the shape factor dis-tribution was observed. When applied to the types of geometries used in this work, the generated meshes adapted well to the significant changes of the physical domain. Therefore, the method presented by Ryskin and Leal[107] is considered highly reliable when solved using a Bubnov-Galerkin formulation. Hydrodynamic phenomena of two dimensional periodic excitation of surface-piercing and submerged cylinders with linear and nonlinear boundary conditions were studied. Non-linear free surface boundary conditions showed distinct behavior with respect to their linear counterpart in the case of a surface piercing cylinder. This effect was seen to be related to the relative size of the composite cylinder with respect to the tank dimensions. Impulsive excitation of tanks of different shapes was also studied. Linear and nonlinear results tend to 136 predict similar behavior for symmetric excitation. Vertical excitation of a composite cylinder i n a fu l l tank predicts frequencies of half tanks. The reason for this is that the line of flow-symmetry below the cylinder acts as a rigid wall, virtually splitting the tank in two. For non-symmetric excitation, results obtained by Yeung and Wu[58] were confirmed. Nonlinear theory is able to predict odd and even frequencies. Linear theory only predicts frequencies associated with the type of excitation used. Quantities of engineering interest pertaining to ship design were also calculated for a mathematically defined hull for a number of length Froude numbers. Nonlinear wave resistance coefficients showed good agreement throughout the Fn interval used in this work. Linear theory only predicted satisfactory results for the higher end of the Froude number interval. Linear and nonlinear t r im calculations showed good agreement with experimental results. This suggests that both theories give satisfactory prediction of the distribution of forces along the hull. Improved sinkage calculations given by nonlinear results suggest that the magnitude of forces is more accurately performed by this theory. This, in turn, is consistent with the improvement obtained for wave resistance coefficients. It can be concluded that coupling a weighted residual formulation with slender body theory within the realm of potential flow theory produces a reasonably accurate and inex-pensive tool for the calculation of some hydrodynamic phenomena arising in ship design. 7.1 F U T U R E R E C O M M E N D A T I O N S Results obtained with linear and quadratic elements showed that although the latter are superior in accuracy the former are preferable from CPU time considerations. For computational applications with limited computing power where quadratic elements are to be used, it is recommended that an interpolation scheme be devised from the mesh generated with linear/4-noded elements to locate the mid-side nodes of the 8-noded element. For applications where computing power is not a concern, it is speculated that Lagrangian 137 elements would give improved results with reduced memory requirements. In this work, homogeneous initial conditions were used for slender body calculations. The role of init ial conditions in nonlinear wave calculations has been reported by Calisal, Goren and Okan[137] to improve the results. Quantities of engineering interest such as wave resistance could be more accurately estimated when starting simulations with non-zero valued init ial conditions. This is an area where improvement of the results can be obtained. The role of viscosity in fluid-structure interaction problems with a free surface is a rel-atively new area. So far, most researchers have concentrated on potential flow calculations. Present trends in computer hardware developments predict that machines wil l be able to perform 10 6 mill ion operations per second by the year 2000. In the future, emphasis should be placed on developing efficient algorithms capable of handling viscous flow with a free surface. The initial stage could simply consider two dimensional cases of the type solved in this thesis to visualize the role of viscosity in the calculated quantities. Next, slender body calculations or three-dimensional calculations with linear free surface boundary condi-tions could be further researched. The latter could allow the use of a mesh defined at the equilibrium position of the free surface thus avoiding mesh regeneration. It has now become clear i n the marine hydrodynamics community that the effect of flow separation on conventional-ship motions is not only restricted to roll motions. Vortex methods are useful tools for cross-sectional shapes with sharp edges where the position of the separation point is easily known. For a more general class of ship-like cross-sections, in fact for virtually all smoothly varying ship-like stations, the position of the separation point is unknown. For this more realistic geometrical scenario, it would be better to solve Navier-Stokes equations It has been shown by many researchers that the finite element method is a powerful tool for a diverse range of engineering calculations. It is recommended that the formula-138 tions included in this thesis be extended to account for the important physical phenomena previously mentioned. 139 BIBLIOGRAPHY [1] Van Dyke M . D . , 'Perturbation Methods in F lu id Mechanics', Academic Press, 1964. [2] M u n k M . M . , 'The Aerodynamic Forces on Airship Hulls ' , National Advisory Commit-tee for Aeronautics, Report 184, 1924. [3] Jones R . T . , 'Properties of Low-Aspect Ratio Pointed Wings at Speeds Below and Above the Speed of Sound', N A C A , Report 835, 1946. [4] Ward G . M . , 'Supersonic Flow Past Slender Pointed Bodies', Quarterly Journal of Me-chanics and Applied Mathematics, vol. 2, 1949. [5] Adams M . C . and Sears W . R . , 'Slender Body Theory - Review and Extension', Journal of Aeronautical Science, vol. 20, 1953. [6] Miles J .W. , ' O n Non-Steady Motion of Slender Bodies', Aeronautical Quarterly, vol. 2, 1950. [7] Miles J .W. , 'Unsteady Supersonic Flow Past Slender Pointed Bodies', Naval Ordinance Test Station, Inyokern, Navord Report 2031, 1953. [8] Cummins W . E . , 'The Wave Resistance of a Floating Slender Body', Unpublished The-sis, American University, 1956. [9] Kaplan P., 'Application of Slender Body Theory to Forces Acting on Submerged Bodies and Surface Ships in Regular Waves', Journal of Ship Research, vol. 1, 1957. [10] Lighthil l M . J . , 'Note on the Swimming of Slender Fish' , Journal of F lu id Mechanics, vol. 9, 1960. [11] Vossers G . , 'Wave Resistance of a Slender Ship', Schiffstechnik, vol. 9, 1962. [12] Maruo H . , 'Calculation of the Wave Resistance of Ships, the draft of which is as small as the beam', Journal of the Society of Naval Architects of Japan, vol. 112, 1962. [13] Ursell F . , 'Slender Oscillating Ships At Zero Forward Speed', Journal of F lu id Mechan-ics, vol. 14, 1962. [14] Tuck E . O . , ' A Systematic Asymptotic Expansion Procedure for Slender Ships', Journal of Ship Research, vol. 8, 1964. [15] Newman J . N . and Tuck E .O. , 'Current Progress in the Slender Body Theory for Ship 140 Motions', Procceedings of the 5th Symposium on Naval Hydrodyanmics, 1964. [16] Newman J . N . , 'The exciting Forces on a Moving Body in Waves', Journal of Ship Research, vol. 9, 1965. [17] Haskind M . D . , 'The exciting Forces and Wetting of a Ship in Waves', Izvestia Akademii Nauk SSSR, 1957 (in russian). D T M B translation No 307, 1962. [18] Tuck E . O . , 'Shallow-Water Flows Past Slender Bodies', Journal of F lu id Mechanics, vol. 26, 1966. [19] Newman J . N . , 'Lateral Motion of a Slender Body Between Two Parallel Walls' , Journal of F lu id Mechanics, vol. 39, 1969. [20] Newman J . N . , and W u T . Y . , ' A Generalized Slender Body Theory for Fish-Like Forms', Journal of F l u i d Mechanics, vol. 57, 1973. [21] Newman J . N . , 'Applications of Slender Body Theory in Ship Hydrodynamics', Annual Review of F l u i d Mechanics, vol. 2, 1970. [22] Chapman R . B . , 'Free Surface Effects for Yawed Surface Piercing Plates', Journal of Ship Research, vol. 20, 1976. [23] Jensen J.J . and Pedersen P.T. , 'Wave-Induced Bending Moments in Ships', Transac-tions of the Royal Institute of Naval Architects, vol. 121, 1979. [24] Troesch A . W . , 'Sway, Rol l and Yaw Motion Coefficients Based on a Forward Speed Slender Body Theory: Part 1', Journal of Ship Research, vol. 25, 1981. [25] Troesch A . W . , 'Sway, Rol l and Yaw Motion Coefficients Based on a Forward Speed Slender Body Theory: Part 2', Journal of Ship Research, vol. 25, 1981. [26] Maruo H . , 'New Approach to the Theory of Slender Ships with Forward Speed', Bulletin of the Faculty of Engineering, Yokohama University, vol. 31, 1982. [27] Noblesse F. , ' A Slender Ship Theory of Wave Resistance', Journal of Ship Research, vol. 27, 1983. [28] Faltinsen O . M . , 'Bow Flow and Added Resistance of Slender Ships at High Froude Numbers and Low Wave Lengths', Journal of Ship Research, vol. 27, 1983. [29] Chen C Y . and Noblesse F. , 'Preliminary Numerical Study of a New Slender Ship Theory of Wave Resistance', Journal of Ship Research, vol. 27, 1983. [30] Chen C Y . and Noblesse F. , 'Comparison Between Theoretical Predictions of Wave 141 Resistance and Experimental Data for the Wigley H u l l ' , Journal of Ship Research, vol. 27, 1983. [31] Sclavounos P . D . , 'The Diffraction of Free Surface Waves by a Slender Ship', Journal of Ship Research, vol. 28, 1984. [32] Sclavounos P .D. , 'The Unified Slender Body Theory: Ship Motions in Waves', Pro-ceedings of the 15th Symposium on Naval Hydrodynamics, 1984. [33] Breit S.R. and Sclavounos P .D. , 'Wave Interactions Between Adjacent Slender Bodies', Journal of F lu id Mechanics, vol. 165, 1986. [34] Song W . , Ikehata M . and Suzuki K . , 'Computation of Wave Resistance and Ship Wave Pattern by the Slender Body Approximation', Journal of the Kamsai Society of Naval Architects of Japan, vol. 209, 1988, in Japanese. [35] Choi H.S. and Mei C . C . , 'Wave Resistance and Squat of a Slender Ship Moving Near the Crit ical Speed in Restricted Waters', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989. [36] Kashiwagi M . , 'Theoretical prediction of Tank Wall Effects on Hydrodynamic Forces Acting on an Oscillating and Translating Slender Ship', Extended Abstracts of the 4th International Workshop on Water Waves and Floating Bodies, 1989. [37] Kashiwagi M . , 'Side Wall Effects on Hydrodynamic Forces Acting on a Ship with Forward and Oscillating Motions', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989. [38] Calisal S . M . and Chan J . L . K . , ' A Numerical Model of Ship Bow Waves', Journal of Ship Research, vol. 33, 1989. [39] Ando S., 'Forward Speed Effect in Head Wave Diffraction over the Forebody of A Slender H u l l ' , Extended Abstracts of the 4th International Workshop on Water Waves and Floating Bodies, 1989. [40] Maruo H . and Song W . , ' A Numerical Appraisal of the New Slender Ship Formulation in Steady Motion' , Proceedings of the 18th Symposium on Naval Hydrodynamics, 1990. [41] Shearer J .R., 'The Experimental Determinatiorr of the Components of Ship Resistance for a Mathematical Model ' , Royal Institution of Naval Architects, 1965. [42] Namimatsu M . , Ogiwara S., Tanaka H . , Hinatsu M . and Kajitani H . , ' A n Evaluation of Resistance Components on Wigley Geosim Models', Journal of the Kamsai Society of Naval Architects of Japan, vol. 197, 1985, in Japanese. 142 [43] Loeser D . , Yue D . and Salvesen M , 'Slender-Body Calculations of Large-Amplitude Motions', Proceedings of the 15th Symposium on Naval Hydrodynamics, 1984. [44] Yeung R. W . , 'Numerical Methods in Free Surface Flows', Annual Review of F lu id Mechanics, vol. 14, 1982. [45] Longuet-Higgins M . S . and Cokelet E . D . , 'The Deformation of Steep Surface Waves on Water I: A Numerical Method of Computation'. Proceedings of the Royal Society of London, Series A 350, 1976. [46] Vinje T. and Brevig P., 'Nonlinear Ship Motion' , Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics, 1981. [47] L i n W . M , Newman J . N . and Yue D . K . , 'Nonlinear Forced Motions of Floating Bodies', Proceedings of the 15th International Symposium on Naval Hydrodynamics, 1985. [48] Isaacson, M . de St. Q., 'Nonlinear Wave Effects on Fixed and Floating Bodies', Journal of F lu id Mechanics, vol. 120, 1981. [49] Baker G .R . , Meiron D . J . and Orzag S.A., 'Generalized Vortex Methods for Free-Surface Flow Problems', Journal of F luid Mechanics, vol. 123, 1982. [50] Dold J . W . and Peregrine D . H . , 'Steep Unsteady Water Waves: A n Efficient Compu-tational Scheme', Proceedings of the 19th International Conference on Coastal Engi-neering, A S C E , 1984. [51] Suzuki K . , 'Calculation of Nonlinear Water Waves Around a Two-Dimensional Body in Uniform Flow by Means of a Boundary Element Method', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989. [52] Sen, D . , Pawlowsky J.S., Lever J . and Hinchey M . J . , 'Two-Dimensional Numerical Modelling of Large Motions of Floating Bodies in Waves', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989. [53] Cointe R., 'Nonlinear Simulation of Transient Free Surface Flows', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989. [54] Cointe R., Geyer P., King B . , Mol in B . and Tramoni M . , 'Nonlinear and Linear Mo-tions of a Rectangular Barge in a Perfect F lu id ? , Proceedings of the 18th International Symposium on Naval Hydrodynamics, 1990. [55] Ghia U . , Shin C T . and Ghia K . N . , 'Analysis of a Breaking Free Surface Wave Using Boundary Fitted Coordinates for Regions Including Reentrant Boundaries', Proceed-ings of the 3rd International Conference on Numerical Ship Hydrodynamics, 1981. 143 [56] Haussling A . J . , 'Solution of Nonlinear Water Wave Problems Using Boundary Fitted Coordinate Systems', Numerical Grid Generation: J . Thompson, editor, 1982. [57] Telste J . G . , 'Calculation of F luid Motion from Large-Amplitude Forced Heave Motion of a Two-Dimensional Cylinder in a Free Surface', Proceedings of the 4th International Conference on Numerical Ship Hydrodynamics, 1985. [58] Yeung R . W . and W u C , 'Nonlinear Wave-Body Motion in a Closed Domain' , Com-puters & Fluids, vol. 17, 1989. [59] Faltinsen O . M . , ' A Numerical Nonlinear Method of Sloshing in Tanks with Two-Dimensional Flow', Journal of Ship Research, vol. 22, 1978. [60] A r a i M . , 'Experimental and Numerical Studies of Sloshing Pressure in Liquid Cargo Tanks', Society of Naval Architects of Japan, Collected Papers, 1984. [61] Bridges T . J . , ' A Numerical Simulation of Large Amplitude Sloshing', Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics, 1981. [62] Yeung R . W . and Wu C. , 'On Nonlinear Wave Motion in a Closed Domain', Jahrbuch der Schiffbautechnischen Gesellschaft, vol. 83, 1989. [63] Kamiya N . and Miyazawa H . , 'Sloshing in Two Dimensional Container with Inclined Side Walls: Boundary Element Approach', Communications in Applied Numerical Methods, vol. 6, 1990. [64] K . Washizu and T. Nakayama, 'Application of the Finite Element Method to Some Free Surface F lu id Problems', Finite Element Methods in Water Resources, 1976. [65] K . Washizu, 'Some Applications of Finite Element Techniques to Nonlinear Free Sur-face F lu id Flow Problems', Proceedings of the 4th International Symposium on Finite Elements in Flow Analysis, 1982. [66] K . Washizu, T. Nakayama, M . Ikegawa, Y . Tanaka and T. Adachi, 'Some Finite Ele-ment Techniques for the Analysis of Nonlinear Sloshing Problems', Finite Elements in Fluids, vol. 5, 1984. [67] Mclver P. and Smith S.R., 'Free Surface Oscillations of Fluids i n Closed Domains', Journal of Engineering Mathematics, vol. 21, 1987. [68] Chang S.C. and Wu S.T., 'On the Natural Frequencies of Standing Waves in a Canal of Arbitrary Shape', Journal of Applied Mathematics and Physics, vol. 23, 1972. [69] Luke J . C , ' A Variational Principle for a F luid with a Free Surface', Journal of F lu id Mechanics, vol. 27, 1967. 144 [70] Whi tham G . B . , 'Nonlinear Dispersion of Water Waves', Journal of F luid Mechanics, vol. 27, 1967. [71] Whi tham G.B . , 'Two-timing, Variational Principles and Waves', Journal of F lu id Me-chanics, vol. 44, 1970. [72] K . J . Bai and R . W . Yeung, 'Numerical Solutions to Free Surface Flow Problems', Pro-ceedings of the 10th Symposium on Naval Hydrodynamics, 1974. [73] Ba i K . J . , 'Diffraction of Oblique Waves by an Infinite Cylinder', Journal of F lu id Mechanics, vol. 68, 1975. [74] Y i m B . , ' A Variational Principle Associated with a Localized Finite Element Technique for Steady Ship-Wave and Cavity Problems', Proceedings of the 1st International Con-ference on Numerical Ship Hydrodynamics, 1975. [75] Ba i K . J . , ' A Localized Finite Element Method for Steady, Two Dimensional Free Sur-face Flow Problems', Proceedings of the 1st International Conference on Numerical Ship Hydrodynamics, 1975. [76] Ba i K . J . , ' A Localized Finite Element Method for Two-Dimensional Steady Potential Flows with a Free Surface', Journal of Ship Research, vol. 22, 1978. [77] B a i K . J . , 'Blockage Correction with a Free Surface', Journal of F lu id Mechanics, vol. 94, 1979. [78] Wellford O L . and Ganaba T., 'Finite Element Procedures for F luid Mechanics Prob-lems Involving Large Free Surface Motion' , Proceedings of the 3rd International Con-ference in Finite Elements in Flow Problems, 1980. [79] Oomen A . , 'Free Surface Potential Flow Computation Using a Finite Element Method', Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics, 1981. [80] Jami A . , 'Numerical Solving of Transient Linear Hydrodynamics Problems by Coupling Finite Elements and Integral Representation', Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics, 1981. [81] Euvard D . , Jami A . , Lenoir M . and Martin D.,-'Recent Progress Towards an Optimal Coupling between Finite Elements and Singularity Distribution Procedures', Proceed-ings of the 3rd International Conference on Numerical Ship Hydrodynamics, 1981. [82] Bai K . J . , ' A Localized Finite Element Method for Steady Three-Dimensional Free Surface Flow Problems', Proceedings of the 2nd International Conference on Numerical Ship Hydrodynamics, 1977. 145 [83] B a i K . J . , ' A Localized Finite Element Method for Three-Dimensional Ship Motion Problems', Proceedings of the 3rd International Conference on Numerical Ship Hydro-dynamics, 1981. [84] Ba i K . J . , K i m J . W . and K i m Y . H . , 'Numerical Computations for a Nonlinear Free Surface Flow Problem', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989. [85] Choi H.S . , B a i K . J . , K i m J . W . and Cho I .H. , 'Nonlinear Free Surface Waves due to a Ship Moving Near the Critical Speed in a Shallow Water', Proceedings of the 18th International Symposium on Naval Hydrodinamics, 1990. [86] Miles J .W. , 'On Hamilton's Principle for Surface Waves', Journal of F lu id Mechanics, vol. 83, 1977. [87] Campana E . , Lal l i F . and Bulgarelli U . , 'Some Numerical Computations about Free Surface Boundary Layer and Surface Tension Effects on Nonlinear Waves', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989. [88] Masuko A . and Ogiwara S., 'Numerical Simulation of Viscous Flow Around Practi-cal H u l l Forms', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989. [89] Zhu M . , Miyata H . and Kajitani H . , ' A Finite-Difference Simulation of a Viscous Flow About a Ship of Arbitrary Configuration', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989. [90] Hino T. , 'Computation of a Free Surface Flow Around an Advancing Ship by the Navier-Stokes Equations', Proceedings of the 5th International Conference on Numer-ical Ship Hydrodynamics, 1989. [91] Miyata H , Zhu M . and Watanabe 0. , 'Numerical Study on a Viscous Flow with Free Surface Waves About a Ship i n Steady Straight Course by a Finite-Volume Method', Journal of Ship Research, vol. 36, 1992. [92] A . Al l iev i and S . M . Calisal, 'Finite Element Application to Grid Generation for Sub-merged and Floating Bodies', International Symposium on Hydro &; Aero Dynamics In Marine Engineering, 1991. [93] Al l ievi A . and Calisal S. M . , ' A Bubnov-Galerkin Formulation For Orthogonal Grid Generation', Journal of Computational Physics, vol. 98, 1992. [94] Robert A . J . , ' A Stable Numerical Integration Scheme for the Primitive Meteorological Equations', Atmosphere-Ocean, vol. 19, 1981. 146 [95] Robert A . J . , ' A Semi-Lagrangian and Semi-Implicit Numerical Integration Scheme for the Primitive Meteorological Equations', Journal of the Meteorological Society of Japan, vol. 60, 1982. [96] Staniforth A . and Temperton C , 'Semi-Implicit Semi-Lagrangian Integration Schemes for a Barotropic Finite-Element Regional Model ' , Monthly Weather Review, vol. 114, 1986. [97] Staniforth A . and Temperton C , ' A n Efficient Two-Time-Level Semi-Lagrangian Semi-Implicit Integration Scheme', Quarterly Journal of the Royal Meteorological Society, vol. 113, 1987. [98] Bermejo R., ' A n Analysis of an Algorithm for the Galerkin-Characteristic Method' , Numerische Mathematik, vol. 60, 1991. [99] Newman J . N . , 'Marine Hydrodynamics', M I T Press, 1977. [100] Finlayson B . A . , 'The Method of Weighted Residuals', Academic Press, 1972. [101] Galerkin B . G . , 'Rods and Plates, Series Occuring in Various Questions Concerning the Elastic Equil ibrium of Rods and Plates', Vestn. Inghenevov, 1915. [102] Zienkiewicz O . C . , 'The Finite Element Method In Engineering Science', McGraw H i l l , 1971. [103] Olson M . D . , 'Advanced Finite Element Methods - Class Notes', 1989. [104] Kreyzig E . , 'Advanced Engineering Mathematics', John Willey &; Sons, Inc., 1962. [105] Press W . H . , Flannery B.P. , Teukolsky S.A. and Vetterling W . T . , 'Numerical Recipes: The Art of Scientific Computing', Cambridge University Press, 1989. [106] Thompson J .F . , ' G r i d Generation Techniques in Computational F luid Dynamics', A . L A . A . Journal, vol. 22, 1985. [107] Ryskin G . and Leal L . G . , 'Orthogonal Mapping', Journal of Computational Physics, vol. 50, 1983. [108] Chikliwala E . D . and Yortsos Y . C . , 'Application of Orthogonal Mapping to Some Two-i Dimensional Domains', Journal of Computational Physics, vol. 57, 1985. i» [109] Babuska I., 'Finite Element Method for Domains with Corners', Computing, vol. 6, 1970. \ [110] Brackbill J . U . , 'Adaptive Zoning for Singular Problems in Two Dimensions', Journal 147 of Computational Physics, vol. 46, 1982. [ I l l ] Brackbil l J . U . , 'Coordinate System Control: Adaptive Zones', Numerical Grid Gener-ation: J . Thompson,editor. Elsevier Science Publishing Co, 1982. [112] Saltzman J.S. and Brackbill J . U . , 'Applications and Generalizations of Variational Methods for Generating Adaptive Meshes', Numerical Grid Generation: J . Thompson, editor. Elsevier Science Publishing Co, 1982. [113] Mastin C. W . and Thompson J .F. , 'Ell iptic Systems and Numerical Transformations', Journal of Mathematical Analysis and Applications, vol. 62, 1978. [114] Mast in C . W . and Thompson J .F. , 'Transformation of Three Dimensional Regions onto Rectangular Regions by Ell iptic Systems', Numerische Mathematik, vol. 29, 1978. [115] Thompson J .F. , Warsi Z . U . A . and Mastin W . C , 'Numerical Grid Generation, Foun-dations and Applications', Elsevier Science Publishing Co., 1985. [116] Thames F . C , Thompson J .F. , Mastin C. W . and Walker R . L . , 'Numerical Solutions of Viscous and Potential Flow About Arbitrary Two-Dimensional Bodies Using Body Fitted Coordinate Systems', Journal of Computational Physics, vol. 24, 1977. [117] Barfield W . D . , ' A n Optimal Mesh Generator For Lagrangian Hydrodynamic Calcula-tions in Two Space Dimensions', Journal of Computational Physics, vol. 6,1970. [118] Godunov S.K. and Prokopov G.P. , 'The Use of Moving Meshes in Gas-Dynamical Computations', U.S.S.R. Comp. Math. Phys., vol. 12, 1972. [119] Amsden A . A . and Hirt C . W . , ' A Simple Scheme for Generating General Curvilinear Grids' , Journal of Computational Physics, vol. 11, 1973. [120] Warsi Z . U . A . , 'Basic Differential Models For Coordinate Generation', Numerical Gr id Generation: J.Thompson, editor. Elsevier Science Publishing Co., 1982. [121] Thompson J .F. , Warsi Z . U . A . and Mastin C . W . , 'Boundary-Fitted Coordinate Systems For Numerical Grid Solution of Partial Differential Equations - A Review', Journal of Computational Physics, vol. 47, 1982. [122] Thompson J .F. , Thames F . C . and Mastin C . W . , 'Automatic Numerical Generation of Body-Fitted Coordinate Systems For Fields Containing Any Number of Arbitrary Two-Dimensional Bodies', Journal of Computational Physics, vol. 15, 1974. [123] Kreyszig E . , 'Introduction to Differential Geometry and Riemannian Geometry', Uni-versity of Toronto Press, 1968. 148 [124] Eisenhart L . P . , ' A n Introduction to Differential Geometry W i t h Use of Tensor Calcu-lus', Princeton University Press, 1940. [125] James M . L . , Smith G . M . and Wolford J .C . , 'Applied Numerical Methods for Digital Computation', Harper & Row Publishers, 1977. [126] Young D . M . and Gregory R .T . , ' A Survey of Numerical Mathematics', Addison- Wesley Publishing Co., 1973. [127] Bowen B. , 'Mathematical Operations in Chemical Engineering', Unpublished class notes. [128] Vreugdenhil C . B . , 'Computational Hydraulics', Springer-Verlag, 1989. [129] Chandra R., ' P h . D . Dissertation', Yale University, 1977. [130] Golub G . H . and Van Loan C . F . , 'Matrix Computations', 1989. [131] Wehausen J . V . and Laitone E . V . , 'Surface Waves', Handbuch der Physik 9, Springer, 1960. [132] Lamb H . , 'Hydrodynamics', Cambridge University Press, 1932. [133] Ghanimati G .R . and Naghdi P . M . , 'Oscillations over Basins of Variable Depth' , Journal of F lu id Mechanics 164, 359-384, 1984. [134] Aanesland V . , ' A Hybrid Model for Calculating Wave Making Resistance', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989. [135] 'Proceedings of the First Workshop on Ship Wave Resistance Computations', David Taylor Naval Ship Research and Development Center, 1979. [136] 'Proceedings of the Second Workshop on Ship Wave Resistance Computations',David Taylor Naval Ship Research and Development Center, 1983. [137] Calisal S . M . , Goren 0 . and and Okan B. , 'On an Iterative Solution for Nonlinear Wave Calculations', Journal of Ship Research, vol. 35, 1991. [138] Springer C . E . , 'Tensor and Vector Analysis Wi th Applications To Differential Geome-i try' , The Ronald Press Company, 1962. j [139] Sokolnikoff I.S., 'Tensor Analysis - Theory and Applications to Geometry and Mechan-ics of Continua', John Wiley, 1951. 149 A P P E N D I X A C O N C E P T S O N D I F F E R E N T I A L G E O M E T R Y The following concepts are only meant to be an aid to the development of the equations described in Chapter 4 on elliptic grid generation and to make this thesis as self-contained as possible. Excellent references on the subject such as[122,123,138,139] are cited. Parts of the following notes are based on these books. B A S E V E C T O R S The curvilinear coordinate system resulting from the solution of an equation such as 4.12, consists of space curves formed by the intersection of planes where one of the curvilinear coordinates is constant, as shown in Figure 1. The tangent to a coordinate line and the normal to a coordinate surface are known as the covariant and contravariant base vectors respectively. Covariant Base Vectors Consider the coordinate line along which only ei varies passing through point P located on the surface 5 , Figure 2. The tangent to this line at P is given by v(ei + det) - v(ei) _ i ^ o dTt = v** W where v = v(x}y, z) is a vector defining the position of point P. The covariant base vectors are then defined as [119] 150 Figure 1 Tangent to a coordinate line. o_=iT l6,. i = 1,2,3 (a.2) where the subscript i corresponds to the coordinate line along which only a varies. Contravariant Base Vectors Consider now a vector at a point P normal to a plane where ei is constant, Figure 2. The three vectors normal to the three surfaces where the coordinates e» are constant are the contravariant base vectors. They are defined as [119] a? = Va (o.3) In general, the normal to a coordinate surface does not coincide with the tangent to a coordinate line. Therefore, covariant and contravariant base vectors are not parallel in the most general case. Parallelism occurs when the system of coordinates is orthogonal. ] 151 D I F F E R E N T I A L E L E M E N T S We now derive the expressions for differential elements corresponding to the length on a curve, area surface and volume. These concepts are needed for the evaluation of integrals related to the Divergence Theorem and the definition of the covariant metric tensor. Element of Length of a Curve Consider a vector v defining the position of a point P located on any line in space. A n increment of this vector position v is given by 3 3 Q£ 3 dv = det = J7:dei - X/ ®id€l (°-4) i = l i = l * »=1 The magnitude ds of the arc length expressed by dv, also known as linear length or element length of the curve, is obtained as 3 (ds)2 = | dv |2= dv • dv (a.5) 3 3 (ds)2 — y ] di • djdeidej (a.6) The term 5j • a;- represents the nine scalar products of the covariant base vectors and it defines a second order symmetric tensor known as the covariant metric tensor. The covariant base vectors are then a tensor of first order, in either case the order is equal to the number of subscripts. The covariant metric tensor components are defined as g%j — di • oij - gji i = l , 2 , 3 j - 1,2,3 (a.1) and then the increment of the arc length is 152 a1 Figure 2 Tangent to a coordinate line and normal to a coordinate surface. 3 3 (ds)2 - ] T ^ 9ijdeidej (a.8) i = l ;' = 1 From (a.6) and (a.8), we can see that the increment of arc length on a coordinate line along which only £i varies is dsi =\cii \dei — sfcjiidzi (a.9) Element of Surface Area The area of the quadrilateral defined by two vectors vi and if2 is defined as S —\ vi X v2 | A differential element of area is then defined as 153 dSi —\ dvj x dilk | i— 1,2,3 i,j,k : cyclic (a.10) or using (a.4) dSi =| tTEj. x v ) 6 f c | dejdek —\cij x Cjt | dejdck i — 1,2,3 i,j,k: cyclic (a. 11) Geometrically, (a. 11) represents an increment in area on a coordinate surface where Ci is constant due to an increment in length of the position vector in the directions of the other two coordinate lines that vary on the coordinate surface. Using the vector identity [yi x v2) • ($3 x - (vi • v3){v2 • Vi) - (vi • Vi)(v2 • v3) (a.12) we h. ave dj x dk | = (dj • dj)(dk • d k ) - (o;- • dfc)(o7- • sk) (0.13) and from (a.7) dj x ak \2- 9jj9kk - 9jk (o.l4) or dSi - \Jgjj9kk - gfkdcjdek i= 1,2,3 i,j,k: cyclic (o.l5) Volume Element A volume generated by three vectors is defined by V — vi • (v2 x v3) 154 then using (o.4), a differential volume can be expressed as dV — viBi • (viBj x vjek)deidejdek - a* • (OJ x ak)deidejdek (0.I6) i—1,2,3 i,j,k: cyclic The scalar triple product is[112] [3i • (a2 x a3)]2 - (01 • ax)[(a2 x a3) • (a2 x a 3)]- | a x x (a2 x o3) |2 (a.17) and «i X (v2 x v3) - (vi • v3)i)2 - (vi • v2)v3 (a.18) Using (o.l2) for the first term of (a.17) and (0.I8) for the second gives [ai • (a2 x a3)]2 = (ax • ai)[(a2 • a2)(a3 • a3) - (a2 • o3)2] - [(ax • a 3)a 2 - (01 • a 2)a 3] 2 Using the definition of the components of the covariant metric tensor (a.7) gives [a r(a 2 X 0 3 ) ] 2 = 0n[ff22033 - ff23] _ 1013022 - 2ffl 3 0 2 3ffl2 + 012033] — (a.19) - ffll[022033 - 0 2 3] - 012[012033 - 013023] + 013L012023 ~ 013022] Expression (a.19) is the determinant of the covariant metric tensor expanded by cofactors. Therefore, [ai • (o2 x a 3 )] 2 = det | j= g (a.20) From (0.I6) we then have dV — ^Jgdeide2de3 (o.21) 155 where ^/g is the Jacobian of the subsequent transformation. D E R I V A T I V E O P E R A T O R S In order to develop expressions for the estimation of quantities involving derivatives, the Divergence/Gauss Theorem is applied to a differential element of volume dV bounded by a differential surface dS formed by coordinate surfaces J J y V • fdV = J J T-ndS (o.22) with T representing any tensor and n the unit normal to the surface S surrounding V. We have seen in (a . l l ) that the cross-product of the covariant base vectors determined the area of the surface differential element. We can then write rtdSi - ±(ctj x dk)dcjdek (a.23) where the sign depends on the location of the surface relative to the volume. Using (o.21) and (o.23) into (a.22) then gives / / / (V • f )Vstfenfe3efe3 = 2 / / T • (a-j x dk)dejdek (a.24) J J JAV i = 1 J JAS or / / / (V • T)y/gdeidemde3 = J J JAV = >~Z I I T-{dj x dk)dejd€k -__ f f T-(Sj X ak)dejdek (°'25) ~[J J*sf i_XJ JAS7 i — 1,2,3 i , j, k : cyclic where the surface integrals in (a.24) are applied on faces opposite to each other of the 156 elemental volume A V . These are defined as A S + and AS{ corresponding to the larger and smaller values of the curvilinear coordinate e» respectively. Divergence Expanding (a.25) gives III ( V • f ^ c f c i e f e j d c s = J J J AV / / T • (a 2 x S3)de2de3 + / / T • (S2 x S3)de2de3 J JAS+ J JAS~ f f f f _ (°- 26) / / T • (a 3 X ai)de3<fei + I T • (a 3 x oi)de 3de:i y 7AS+ ./ -/AS-/ / f • (ai x a2)deide2 + / / T • (a x x a 2 )deide 2 Taking limits at both sides of (a.26) as Av —> 0 or, in other words, taking derivatives with respect to ex,e2 and e 3 on both sides of a.26 gives V • f = [ f • (oa x a 3 )] , £ l + [ f • (a 3 x Si) ] | 6 3 + [ f • (a x x a 2 ) ] , £ 3 (a.27) since the quantity T • (Si x ay) only varies in the k direction (i, j, k : cyclic) we have V • f = — ] T [ T • (a,- X Sk)]tH i, j, k : cyclic (a.28) Equation (a.28) is fundamental in the estimation of derivative operators. It is also respon-sible, as it wil l be shown below, for the two major forms of elliptic grid generation systems; i.e., conservative and non-conservative forms. Expanding (a.28) gives 1 V • T = — £ [ T i e < • & x a fc) + f • (Sj x Sk)tH) (a.29) 157 and with (a.2) in the second term of (a.29) we have 3 3 «=1 8=1 3 = E f • K«,-.i x W + E ^ • K.,- x j=l t=l i , j , A : cyclic (o.30) Since the indexes are cyclic, (o.30) can be written without loss of generality as E ^ - ( o V xa f c), e i i=l - E f • ft.,-., x vEk] + E f • E-* x if,,.,.] i=l i=l i , j , k : cyclic (a.31) but from vector calculus we know that vi X v2 — —v2 x u i . Therefore, (o.31) becomes 3 E ^ -(ay x o f c ) , E i = 0 (a.32) i=l Equation (a.29) then becomes 1 3 -V • f - — J2[TEi • (ay x afc)] i , : q / d i c (a.33) We have shown that equations (a.28) and (a.33) are entirely equivalent in an analytical sense; however, in a numerical sense they are not. In equation (a.28) the quantity T • (ay X a*) represents the flux of the quantity T through the incremental area (ay X a*). In (a.26) that led to (a.28), the two incremental surfaces at opposing sides of the elemental volume are considered. Thus, the quantity T • (ay X a*) is estimated on both of these faces. This form of estimating V • T is known as the conservative form. On the other hand, (a.33) expresses 158 the flux of the quantity T through an average incremental surface element (Sj X oj.) inside the elemental volume. This form of evaluating V • T is known as the non-conservative form. Gradient The gradient of a scalar quantity T defined as V T can be obtained from the equation derived for the divergence by replacing the dot product by simple operation on the left hand side and multiplication on the right hand side. The conservative form of the gradient is 1 V T = — Y_[T(Sj x ak)\H i, j, k : cyclic (a.34) V ? i = i and the non-conservative form 3 -» 1 . ^ V T = — Y^[TCi{dj x ak)] i, j, k : cyclic (a.35) Non-Conservative Gradient of a Curvilinear Coordinate Base Vectors and Metric Tensors Having derived the gradient of a scalar function T, (a.35), we can apply the operator V to each of the curvilinear coordinates a. Thus from (a.35) V c ; = — *Y_\?l>Hx 1 - 1> 2> 3 *'3>  k : cVclic (a.36) V ? t=l but since each of the z% coordinates are mutually independent, we have that eitH - 6\. Then 1 Ve; = x ak) i— 1, 2,3 i, j, k : cyclic (a.37) s/9 but (a.37) is nothing but the definition, (a.3), of the contravariant base vectors, thus 159 if — Vej = —p(aj x dk) i— 1,2,3 i, j, k : cyclic (a.38) Equation (a.38) is fundamental since it expresses the relationship between the contravari-ant and the covariant base vectors. In other words, it expresses the relationship between the derivatives of the curvilinear coordinates with respect to the Cartesian coordinates ei,a>ii {contravariant) and the derivatives of the Cartesian coordinates with respect to the curvilinear coordinates x t ) e i , (covariant). We are now able to define the components of the contravariant metric tensor in a fashion similar to that of the covariant counterpart. Thus, as in (a.7), we have g(i -b? bl -gji i = 1,2,3 j = 1,2,3 (a.39) Introducing (a.38) into (a.39) gives 1 gtl - b? - a 1 = -(a,- X a-kYiim x 9 i - 1,2,3 i,j,k: cyclic (a-40) £ = 1,2,3 l,m,n: cyclic Using (a. 12) we have 9 U - -[(<*7 " <*m)(»* • Om) ~ (dj • an)(dk • am)] (tt.41) 9 and introducing (a.7) which defines the components of the covariant metric tensor into (a.41) gives 9 %l — -[9jm9km - 9jn9km] 9 i— 1,2,3 i,j,k: cyclic (a.42) 1 — 1,2,3 l,m,n: cyclic 160 The right hand side of (a.42) is the signed cofactor of the li component of the covariant metric tensor divided by the determinant g of this tensor. The determinant of the contravariant metric tensor | g%l | is then equal to the inverse of the covariant metric tensor | gu | det | gil |= , ] . = - (o.43) det \ gu \ g Non-Conservative Laplacian Operator Using (a.38) into the non-conservative form of the gradient (a.35) gives V T = £ [ T £ ; . c F ] (a.44) thus the non-conservative form of the gradient operator V becomes V = £ o > A (a.45) Applying (a.45) to (a.44) gives V • ( V T ) - V 2 T = E E ^ - * T „ . . , + E E ° " • ^ U T ' i («-46) i=l j = l i=l}=1 Using (a.45), the Laplacian of a curvilinear coordinate is defined as 3 V 2 e,- = V • (Ve,-) = V • cP = ' * (°-47) i=l Introducing (a.47) into the last term of (a.46) produces 3 3 3 V 2 r = E E ^ T > ^ + E ( V % ) T > V (°-48) ! = 1 7=1 J=l which is the non-conservative form of the Laplacian operator of the quantity T. 161 A N G L E B E T W E E N C O O R D I N A T E LINES In order to quantify the orthogonality characteristics of the generated mesh, the angle between coordinate lines can be calculated as[123] cos 6=—^= (o.49) Using definitions given by (a.2) and (a.7) on a vector v — zi + yj gives 2 1 2 9u = 2 , + y% 922 = z\ + y% (a.50) 012 = z,ez,e + y,sy,e and with (a.50) into (a.49) we have cos 9 = 9 s 9 s 9 E 9 e (a.51) /^[(ff )2 + (|^ )2][(|f)2 + (If)2] with 9 the angle formed by adjacent sides of an element. 162 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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-0302316/manifest

Comment

Related Items