NUMERICAL PROCEDURE FOR POTENTIAL FLOW PROBLEMS WITH A FREE SURFACE By JOHNSON LAP-KAY CHAN M.A.Sc. The Un i v e r s i t y of B r i t i s h Columbia, 1984 B.A.Sc. The U n i v e r s i t y of Toronto, 1982 A THESIS SUBMITTED IN PARTIAL FUFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n THE FACULTY OF GRADUATE STUDIES ( DEPARTMENT OF MECHANICAL ENGINEERING ) We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1987 © Johnson Lap-Kay CHAN, 1987 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia 1956 Main Mall Vancouver, Canada Department V6T 1Y3 DE-6(3/81) ABSTRACT A numerical procedure based upon a boundary i n t e g r a l method f or gra v i t y wave making problems i s studied i n the time domain. The free-surface boundary conditions are combined and expressed i n a Lagrangian notation to follow the free-surface p a r t i c l e ' s motion i n time. The corresponding material d e r i v a t i v e term i s approximated by a f i n i t e d i f f e r e n c e expression, and the v e l o c i t y terms are extrapolated i n time f o r the completion of the formulations. The fluid-body i n t e r s e c t i o n p o s i t i o n at the free surface i s predicted by an i n t e r p o l a t i o n function that requires information from both the free surface and the submerged surface conditions. Solutions corresponding to a l i n e a r free-surface condition and to a non-linear free-surface c o n d i t i o n are obtained at small time increment values. Numerical modelling of surface wave problems i s studied i n two dimensions and i n three dimensions. Comparisons are made to l i n e a r a n a l y t i c a l solutions as well as to published experimental r e s u l t s . Good agreement between the numerical solutions and measured values i s found. For the modelling of a three dimensional wave d i f f r a c t i o n problem, r e s u l t s at high wave amplitude are r e s t r i c t e d because of the use of q u a d r i l a t e r a l elements. The near cy l i n d e r region of the free surface i s not considered to be well represented because of the coarse element s i z e . Wave forces c a l c u l a t e d on the II v e r t i c a l c y l i n d e r are found to be af f e c t e d by the modelled tank length. When the simulated wave length i s comparable to the wave tank's dimension, numerical r e s u l t s are found to be less than the experimental measurements. However, when the wave length i s shorter than the tank's length, solutions are obtained with very good p r e c i s i o n . I l l TABLE OF CONTENTS ABSTRACT II TABLE OF CONTENTS IV NOMENCLATURE V LIST OF FIGURES VIII ACKNOWLEDGEMENT XI CHAPTER I INTRODUCTION 1 CHAPTER II -- BASIC THEORY 5 II.1 Boundary Conditions 7 I I . 2 Time Stepping Formulation 12 CHAPTER III - NUMERICAL METHOD 17 I I I . l Boundary Integral Method 18 I I I . 2 Extrapolation i n Time 22 CHAPTER IV -- APPLICATIONS 24 IV. 1 Wave Tank Simulation 25 IV.2 Deformation of High Amplitude Waves 40 IV.3 Predictions of Ship's Bow Wave .... 48 IV. 4 Wave D i f f r a c t i o n of a C i r c u l a r Cylinder ( 3-D ) 66 CHAPTER V --- DISCUSSIONS AND CONCLUSIONS 85 V. 1 General Discussions 85 V.2 Conclusions 91 REFERENCES 92 FIGURES 95 APPENDIX 127 IV N O M E N C L A T U R E Symbol d> V e l o c i t y P o t e n t i a l Function r) Wave Elevation p F l u i d Density V D i f f e r e n t i a l Operator ^ M a t e r i a l Derivative St Time Increment value or Time Step w Angular Frequency A Wave Length e Small Parameter f o r Dimensional Analysis a Half Angle of Wedge Model A Wave Amplitude a Radius of V e r t i c a l Cylinder Am Maximum Wave Height Measured Along Ship's Bow A<£=B System of Linear Equations Aw/Am Wave Amplitude to Wave Paddle Motion Amplitude Ratio B Geometric Function for Ship H u l l c Phase V e l o c i t y of Free-Surface Wave c Geometric Center of Qu a d r i l a t e r a l Elements d Mean Water Depth F Draft-Froude Number 5 Extrapolation Function i n Time V Longitudinal Force c a l c u l a t e d on V e r t i c a l Cylinder Non-dimensional, maximum F^ value i n a Wave Cycle Green's Function Gravity Constant Wave Height to Wave Length Ratio Wave Height to Water Depth Ratio Wave Number Ship Length [ n , n , n ] Normal Unit Vector x y z Point of Interest Pressure Control Point on the Control Surface,S Distance Measured Between Two Points Total Control Surface Submerged Surface of V e r t i c a l Cylinder Control Boundary at the Free Surface Control Surface of Ship's Bow Control Surface at the Wave Paddle Impermeable Control Surface V e r t i c a l Imaginary Boundary at one Wave Length Separation Wave Period Ship's Draft Time for High Amplitude Wave to Break Time Variable Model Speed VI v [ v , v , v ] V e l o c i t y Vector x y z J X Motion Amplitude of Wave Paddle m x Displacement Function of Wave Paddle Ymax Non-dimensional Value of Am a,b,c,d,..Unknown Constants i n Quadratic Functions F , F ,F , 2 2b 3 and K Known Functions Defined i n the Boundary Integral Equation of associated Problems I,J,K,L Element Numbers P,Q,R,S Four Corners of Qu a d r i l a t e r a l Elements VII LIST OF FIGURES Figure Page 1. A T y p i c a l Wave Making Problem 95 2. Numerical Model with Mirror Image at y=-d .. 95 3. Linear and Non-linear Wave Tank Models 96 4. Int e r p o l a t i o n of V e l o c i t y Between Elements . 96 5. . Wave Generations by Numerical Wave Tank .... 97 6. P r e d i c t i o n of Wave Length by Linear Numerical Model 98 7. Comparison of Phase V e l o c i t y Between Solution and Wave Theory 98 8. Wave Amplitude Predictions by Numerical Solution 99 9. Wave P r o f i l e ' s Smoothness vs Element Number 99 10. Comparison Between Linear and Non-linear Wave Tank r e s u l t s 100 11. Non-linear Wave Modelling by Numerical Wave Tank 101 12. Motions of Free - Surface's P a r t i c l e s i n Time 101 13. Single Wave Length Model f or S p a t i a l l y P eriodic Wave 102 14. Simulation of Free-Surface Wave 103 15. Deformation of High Amplitude Wave 103 VIII Figure Page 16. Deformation of High Amplitude Wave 104 17. Deformation of High Amplitude Wave 104 18. Waves Generated by a Wedge-shaped Model .... 105 19. Sectional View of the Considered Domain .... 105 20. Bow Wave Pre d i c t i o n Calculated Along Wedge-shaped Model 105 21. Bow Wave Predictions at D i f f e r e n t Draft .... 106 22. Non-dimensional Bow Wave Amplitude Comparisons at D i f f e r e n t Drafts ( a = 7.5° ) 107 23. Non-dimensional Bow Wave Amplitude Comparisons at D i f f e r e n t Drafts ( a = 15° ) 108 24. Longitudinal P o s i t i o n of Wave Peak 109 25. Bow Wave Modelling at Low Model Speed 110 26. Bow Wave Modelling at Intermediate Model Speed I l l 27. Secondary Bow Wave Obtained at Low Model Speed 112 28. Formation of Water Spray at High Model Speed 112 29. Numerical Modelling of Wave Tank with a Surface P i e r c i n g Cylinder ( 3-D ) 113 30. Q u a d r i l a t e r a l Element 113 IX Figure Page 31. H a l f - c y l i n d e r Modelling by Qu a d r i l a t e r a l Panels 113 32. Free-Surface of Numerical Tank Represented by Q u a d r i l a t e r a l Elements 114 33. Interpolation of V e l o c i t y Components Between Elements ( 3-D ) 114 34. Numerical simulations of D i f f r a c t e d Wave i n a Tank 115 35. Deformation of Free-Surface Elements 118 36. D i f f r a c t e d Waves Studied i n a Short Tank ... 119 37. Deformation of Free-Surface Elements Affected by Shorter Tank Length 121 38. Simulations of D i f f r a c t e d Waves at Higher Frequency 122 39. F a i l u r e of Non-linear Model at High Wave Amplitude 125 40. Non-linear Forces Calculated at the Cylinder 126 X ACKNOWLEDGEMENT This thesis i s written under many unforgettable contributions from many of my best friends, family members, and f a c u l t y members. I would l i k e to single out a few persons who have a s s i s t e d and supported me throughout these years. Thanks to the Lam's family, f o r feeding me with t h e i r d e l i c i o u s food when I was starving at the end of each month. Their continuous encouragement never run out when I was desperate. I wish them always good luck i n t h e i r future. Thanks to Mr. Ken Law, who i s one of my best friends i n Vancouver. Without h i s help, the three dimensional graphical p l o t s i n t h i s thesis could never have been accomplished. I wish him and h i s company a very b r i g h t future. Thanks to Annie, my loveable g i r l f r i e n d , who has been so patient and understanding f o r these years. I would l i k e to t e l l her how much I appreciated her encouragement when I was depressed, and I enjoy being her boyfriend f o r a l l these years. Thanks f o r my eldest brother Lap-Ngai, who p a r t i a l l y supported my l i v i n g f o r these years, e s p e c i a l l y on h i s taking care of our parents independently. I wish him continuous success i n h i s business. Thanks f or Professor Isaacson who gave me the in s t r u c t i o n s at the f i n a l stage of my thesis work when my supervisor was absent on h i s sabbatical leave. I would l i k e to apologize f o r my numerous interruptions to him at h i s work. F i n a l l y , I would l i k e to devote t h i s thesis to my supervisor, Dr. C a l i s a l , f o r h i s enormous patience and encouragement to guide me throughout a l l these years. I am not a good student, and I never was. Thanks f or h i s t o l e r a t i o n within the l a s t f i v e years. I wish he has a very b r i g h t and successful future i n h i s research career, and good luck to h i s family. with love and s i n c e r i t y Johnson L.K. Chan XI I. INTRODUCTION In recent years, ocean engineers have become more and more interes t e d i n the development of numerical procedures that p r e d i c t free-surface wave phenomena. This i s understandable since today's computer technology i s advancing at such a tremendous rate, the need for engineering software to be developed so as to u t i l i z e the computer power has become very strong. Moreover, a well developed numerical procedure can reduce the cost of conducting numerous exp ens ive exp e r iment s. Although wave making problems have been studied for almost a century, achievement i n t h i s area i s s t i l l quite l i m i t e d . The involved numerical s o l u t i o n procedure i s considered d i f f i c u l t because of the moving boundaries. Numerical models that simulate t h i s transient, non-linear problem i n the time domain seem to be unstable. This numerical i n s t a b i l i t y remained as an obstruction to the progress of ocean engineering technology for years. I t was not u n t i l 1976 that Longuet-Higgins and Cokelet s u c c e s s f u l l y simulated the deformation of high amplitude waves, and proposed a s o l u t i o n to the i n s t a b i l i t y problem. Numerical i n s t a b i l i t y associated with free-surface p o t e n t i a l flow problems i s believed to be i n i t i a t e d by the 1 formation of l o c a l short gravity waves when the long wave crest i s compressed i n the h o r i z o n t a l d i r e c t i o n . In r e a l l i f e , these l o c a l short gravity waves are damped by the f l u i d ' s v i s c o s i t y as well as the surface tension e f f e c t . However, i n the numerical model, the presence of these short gr a v i t y waves exaggerates the computation of the free surface, and leads to an unstable s o l u t i o n . In order to suppress the numerical i n s t a b l i t y i n the sol u t i o n , Longuet-Higgins and Cokelet have proposed the use of a smoothing function to correct the r e s u l t every few time step i n t e r v a l s . Solutions obtained with such a smoothing procedure then become very stable. Another e x i s t i n g problem that causes d i f f i c u l t i e s i n the numerical model i s the determination of the p o s i t i o n where the free-surface i n t e r s e c t s an impermeable body's surface, such as the wave paddle. This point was in d e n t i f e d by L i n (1984) as a s i n g u l a r i t y , and an a n a l y t i c s o l u t i o n could not be obtained. However, i t can be accommodated by applying both the v e l o c i t y p o t e n t i a l and the streamline function as prescribed conditions to reinf o r c e a s o l u t i o n numerically. Results obtained by following t h i s approach are very promising. Unfortunately, formulations are accomplished by using complex va r i a b l e s which are applicable for two dimensional problems only. Since i t i s more desirable to e s t a b l i s h a three dimensional numerical procedure for r e a l 2 l i f e a p p l ications, a free-surface formulation that contains only r e a l v a r i a b l e s i s , therefore, introduced and-investigated. In the following chapters, a s o l u t i o n method based upon the Boundary Integral Equation i s studied. The presented method d i f f e r s from Longuet-Higgins' approach by coupling a f i n i t e difference formulation onto the free-surface boundary conditions. A d d i t i o n a l assumptions are made to obtain a more appropriate form to march the s o l u t i o n i n time. The s i n g u l a r i t y problem at the fluid-body i n t e r s e c t i o n point i s handled by a simple i n t e r p o l a t i o n function. Four d i f f e r e n t kinds of wave making problems are examined, and comparisons are made between numerical solutions and published experimental r e s u l t s . The s o l u t i o n procedure i s f i r s t applied to simulate a two dimensional wave tank. A p i s t o n type wave maker i s modelled for the wave generating process, and numerical solutions are compared with l i n e a r wave theory. The next a p p l i c a t i o n i s on the modelling of a two dimensional, s p a t i a l l y p e r i o d i c wave problem. Deformation of high amplitude waves at a f i n i t e water depth i s considered. Results s i m i l a r to those by Longuet-Higgins and Cokelet (1976) are obtained. Then, ship's bow waves which are generated by a wedge-shaped model at d i f f e r e n t d r a f t s are 3 investigated. Predictions using a slender ship assumption are compared with measurements by O g i l v i e (1972). F i n a l l y , the three dimensional a p p l i c a b i l i t y of the free-surface formulations i s i l l u s t r a t e d by modelling a simple wave d i f f r a c t i o n problem. Forces experienced by a surface-piercing c y l i n d e r located i n a wave tank are ca l c u l a t e d and compared to published experimental works. Numerical r e s u l t s obtained i n the following chapters do not require any smoothing procedure. The free-surface formations with a small time step assumption are considered to be adequate to y i e l d a numerically stable s o l u t i o n . I t i s believed that the presented procedure provides a subsequent way to model wave making problems other than the Longuet-Higgins' approach. Moreover, the introduced free-surface formulations enable p r a c t i c a l problems to be solved i n three dimensions. 4 II. BASIC THEORY When an object i s moving at or near the water surface, the air-water i n t e r f a c e i s disturbed from i t s equilibrium, and the free surface i s set into motion. However, since a g r a v i t a t i o n a l force i s acting on the f l u i d , there i s a tendency f o r the air-water interface to return to i t s equilibrium p o s i t i o n . As a r e s u l t , gravity waves are formed at the free surface as a balance of the k i n e t i c and the p o t e n t i a l energy. These waves which are commonly known as grav i t y waves, propagate away from the source of disturbances i n the r a d i a l d i r e c t i o n . In general, c e r t a i n parameters such as water depth, and v e l o c i t y of disturbance are considered important for wave formation, while some other f l u i d properties l i k e surface tension, v i s c o s i t y , and compre s s i b i l i t y of the f l u i d medium are considered less important for the study of gravity waves. In the following chapters, wave making problems are formulated by assuming the f l u i d as incompressible and i n v i s c i d . Surface tension e f f e c t i s considered as n e g l i g i b l e compared to gravity e f f e c t and f l u i d i n e r t i a l forces. These assumptions are generally accepted when there i s no flow separation. For gravity wave problems, flow separation w i l l occur when the c h a r a c t e r i s t i c dimension of the structure i s not large r e l a t i v e to the o r b i t a l path length of f l u i d p a r t i c l e s . That i s , the i n v i s c i d assumption i s acceptable when the dimension of the structure 5 to the wave length r a t i o i s at about the order of unity. With these assumptions, the problem can be treated as a p o t e n t i a l flow problem. A v e l o c i t y p o t e n t i a l function, <j>, i s defined within the flow f i e l d . This p o t e n t i a l function, whose gradients are the f l u i d v e l o c i t y components, s a t i s f i e s the Laplace Equation : V 2 <(> = 0 or <j> + <j> + <f> = 0 (1) . ,xx ^,yy r , z z In the following chapters, a right-hand Cartesian coordinate system with y pointing up as the p o s i t i v e d i r e c t i o n i s used ( Figure 1 ). The o r i g i n of the coordinate axis i s located at the undisturbed free-surface, and the wave ele v a t i o n i s designated by the symbol r\. Equation (1) , which i s the governing equation for the flow f i e l d , can be solved by many well established numerical or a n a l y t i c a l methods. Although the Boundary Integral Method i s used i n the following studies, there i s no reason that other numerical methods cannot be used to solve the problem. In the next chapter, the Boundary Integral Method f o r so l v i n g 6 p o t e n t i a l flow problems w i l l be discussed. However, boundary conditions must f i r s t be arranged i n an appropiate form for the use of the selected numerical method. II.1 Boundary Conditions Although the Laplace Equation i s a l i n e a r , second order d i f f e r e n t i a l equation, the associated problem can be a non-linear one, depending on the kind of boundary conditions involved. In f a c t , there e x i s t d i f f e r e n t kinds of boundary conditions f o r free-surface wave problems. In order to derive the boundary conditions i n forms suitable f o r the d i f f e r e n t applications investigated i n l a t e r chapters, they are discussed under two main categories : i . The Wetted Surface Condition i i . The Free-surface Condition The following formulations are arranged s p e c i f i c a l l y f o r the use with the Boundary Integral Method. For some other numerical methods, c e r t a i n rearrangements should be considered. I l . l . i . The Wetted Surface Conditions The wetted surface i s defined as the submerged impermeable surface of s o l i d bodies that i s e i t h e r moving or 7 stationary. In wave making problems, the wave maker, the tank wall, and the submerged surfaces of any f l o a t i n g or standing structures belong to t h i s category. This type of boundary condition can be formulated as follows : 4> n = v - Von (2) . , n n 4> ^ i s the f l u i d v e l o c i t y component i n the normal d i r e c t i o n to the considered surface, and i s equal to the dot product of V and n. V i s the v e l o c i t y vector at any point on the moving surface, and n i s the normal u n i t vector pointing outward from the considered domain ( Figure 1 ). In equation (2) , f l u i d p a r t i c l e s at the wetted surface are forced to have a v e l o c i t y component, normal to the surface, equal to the normal v e l o c i t y component of the surface because of the impermeable boundary requirement. As for the tangential v e l o c i t y component of f l u i d p a r t i c l e s on th i s surface, there i s no r e s t r i c t i o n because of the i n v i s c i d f l u i d assumption. The no-slip boundary condition f o r viscous f l u i d s does not apply here. When equation (2) i s applied to the submerged surface of stationary objects, such as the wall of a wave tank, v i s i d e n t i c a l l y equal to zero. Therefore, <f> i s also equal to zero, which means f l u i d p a r t i c l e s on t h i s surface have a zero 8 normal v e l o c i t y component, and penetration of f l u i d p a r t i c l e s into the surface i s not allowed. I f the surface under consideration i s on the wave piston, then v i s a known function of time. F l u i d p a r t i c l e s on any point of t h i s surface are moving at a v e l o c i t y equal to that of the surface i n the normal d i r e c t i o n . However, i f the submerged surface i s on a f l o a t i n g object and the object's motion ( v ) i s the unknown to be solved as part of the solut i o n , the boundary condition, <j> , cannot be defined by following the form as i n i n equation (2). A d d i t i o n a l arrangements or assumptions have to be made to avoid s o l v i n g f o r V d i r e c t l y . This can be done by a pr e d i c t o r - c o r r e c t o r kind of procedure f o r the value of V i n the time domain. Since f l o a t i n g objects are not included i n the following studies, the formulation of t h i s kind of boundary condition i s , therefore, excluded from the discussions. In Figure 1, a phys i c a l wave tank problem i s i l l u s t r a t e d by a schematic diagram. A pi s t o n type wave maker i s located on the left-hand side of the tank. A displacement function, x , which i s p e r i o d i c i n time, i s assigned to t h i s p i s t o n f o r m wave generation. Further down the other end, the considered domain i s completed by a v e r t i c a l wall. In frequency domain problems, an a r t i f i c a l boundary that simulates a r a d i a t i o n condition i s necessary. However, the r a d i a t i o n boundary condition i s not required f o r a time stepping s o l u t i o n 9 scheme. Although a non- r e f l e c t i n g boundary condition can be used to reduce the computational e f f o r t , no such attempt was made i n t h i s study. Instead of using a no n - r e f l e c t i n g boundary, a standard wave tank was considered as a ph y s i c a l model. In order to ensure that the s o l u t i o n i s free from the interference of r e f l e c t e d waves, a reasonable tank siz e i s necessary, and the simulation time i s l i m i t e d . A discussion of the r a d i a t i o n boundary condition can be found i n Sarpkaya and Isaacson (1981) as well as i n Newman (1977). I l . l . i i . The Free-Surface Condition The free surface i s i d e n t i f i e d as the a i r - f l u i d i n t e r f a c e which i s free to move under any disturbance. One of the d i f f i c u l t i e s involved i n solving the problem i s to pre d i c t the p o s i t i o n of the free surface before the s o l u t i o n i s obtained. At the free surface, the flow has to s a t i s f y two conditions, namely the kinematic and the dynamic boundary conditions. Kinematically, f l u i d p a r t i c l e s on the free surface must remain there. That i s , at the free surface, f l u i d p a r t i c l e s are moving with a v e l o c i t y component equal to the normal v e l o c i t y component of the free surface i n a d i r e c t i o n normal to the free surface. Dynamically, the flow condition must obey the B e r n o u l l i Equation. These two conditions are defined as i n equations (3) and (6) res p e c t i v e l y i n Lagrangian notation. 10 Dx Dt Dy Dt = 4> on y=r? (3) Dz Dt and <t>r + i(<l>l + <f>l + <t> I ) + ^ + g y = 0 at y=r? (4), where ^ i s the material de r i v a t i v e given as : D T " a t + *,* a*. + +,y ay" + +,z a l ( 5 ) -, <f> , and ^ are the f l u i d v e l o c i t y components i n the • x > y i z x , y , and z d i r e c t i o n s , ry i s the free-surface e l e v a t i o n measured i n the y d i r e c t i o n , P i s the pressure taken as zero at rj, p i s the f l u i d density, and g i s the g r a v i t a t i o n a l constant. In order to trace the motion of f l u i d p a r t i c l e s at the free surface, equation (4) i s rearranged as : T>6 1 , , 2 , 2 , 2 , Dt " 2 , x + , y + , z + g y = 0 at y=„ When t h i s above equation i s combined with equation (3), the 11 following equation can be obtained : 2 Equation (6) which i s a combined form of the kinematic and dynamic free-surface conditions, has to be s a t i s f i e d at y=r?. A d e t a i l e d discussion of the free-surface boundary conditions i s given i n Newman (1977). Although, the kinematic and dynamic conditions are combined as i n equation (6), most of the terms are not known before the problem i s solved. In order to proceed and rearrange equation (6) into a useable form, two a d d i t i o n a l assumptions must be made. A f i n i t e difference expression i s coupled into the equation to approximate the material d e r i v a t i v e term, and a second order extrapolation function i s used i n time to p r e d i c t unknown v e l o c i t y terms. They are discussed i n sections II.2 and III.2 r e s p e c t i v e l y . II.2 Time Stepping Formulation Upon applying the Boundary Integral Method, e i t h e r <j> , d> , or t h e i r combinations must be defined at the c o n t r o l surface as boundary conditions. Therefore, the following rearrangements of equation (6) are made. 12 On y=r/ , <j> i s defined by the dot product of V(j> and n. I f n^, n^, n^ are the components of the normal u n i t vector, n, located at the free- surface, then, 4> = < f > T l + < f > T \ + ( j ) T l (7), ,n x r , y y \ z z where <j> can be obtained from equation (6). Upon the > y s u b s t i t u t i o n of 4> from equation (6) into (7) , the following > y form i s obtained. <f> = d> n + <j> n ,n , x x , z z 2 1 , 1 D / J 2 . 2 , 2 . D i i , , n K + I { 2 bT ( ^ ) X + ^ ,y + ^ , z > " ^ 1 ) n y a t ^ ( 8 >-Two a d d i t i o n a l assumptions are made here to approximate the unknown values f o r computing the time stepping s o l u t i o n . These assumptions are considered as follows : i . In time stepping s o l u t i o n problems, a time increment value ( time step ), St, must be chosen. This time increment value has to be small f o r the s o l u t i o n to be accurate. I f St i s small enough, unknown v e l o c i t y terms on the right-hand-side of equation (8) can be extrapolated i n time ( from t h e i r previous known values) without introducing large errors into the so l u t i o n . In general, the smaller the time step, the more precise the extrapolated values w i l l be. The order of t h i s extrapolation function can be chosen a r b i t r a r i l y . In the 13 following studies, a second order extrapolation function i s used. The dicussion of the extrapolation procedure can be found i n section I I I . 2. A prime notation i s used to i d e n t i f y the extrapolated terms i n equation (10). i l . A f i n i t e difference form i s used to approximate the second material d e r i v a t i v e term i n equation (8). A four points function with constant time i n t e r v a l , fit , i s written as : where the p o t e n t i a l functions with a negative superscipt represent solutions from previous time steps, and <f>° i s the unknown value to be solved at the current time step. There i s no r e s t r i c t i o n on the number of terms used i n the expression. With these two assumptions, equation (8) can be D2<f> 2<t>° - 5(f)'1 + kfZ - <{> - 3 (9), Dt 2 ( fit ) 2 rewritten as 1 J l -1 - 2 - 3 (St)2 + 4>' n + </>' n + ,X X , z z 2g at y=r? (10), 14 where <f> i s now the only unknown v a r i a b l e on the right-hand side of equation (10). The l i n e a r i z e d form of equation (6) applied at the undisturbed free surface and used to obtain the l i n e a r numerical s o l u t i o n (see Newman, 1977), i s given as : <t> r r + g<f> v = 0 at y=0 (11) where a l l the non-linear terms are dropped. Moreover, since the l i n e a r i t y assumption i s made, the following expressions can be accepted when formulating the l i n e a r free-surface boundary condition : n = 1 ( n = 0 , ii « 0 ) y x z Dt 2 Therefore, equation (11) can be rearranged i n a form s i m i l a r to equation (10) as : , 1— ? + ± J £ L I J ± L L £ L a t y = 0 ( 1 2 ) . n 2 2 g (St)2 g (St)2 F i n a l l y , i t must be emphasized that the non-linear free-surface boundary condition i s applied at r), while the l i n e a r i z e d condition i s s a t i s f i e d at y=0. Equation (10) and 15 (12), which are now arranged i n an su i t a b l e form, can be included i n the Boundary Integral Equation f or d i f f e r e n t a p p l i c a t i o n s . 16 I I I . NUMERICAL METHOD In t h i s chapter, the Boundary Integral Method i s b r i e f l y introduced. This p a r t i c u l a r numerical method i s chosen because i t i s considered to be an e f f e c t i v e s o l u t i o n method fo r modelling p o t e n t i a l flow problems with a moving boundary. The g r i d generation procedure involved with t h i s method i s not as complicated as i t i s i n the other numerical methods. When applying the Boundary Integral Method, the numerical model i s constructed by representing the c o n t r o l surface using small facets ( elements ). Only values at the control surface are stored and calculated, whereas i n other methods such as the F i n i t e Element Method, d i s c r e t i z a t i o n of the e n t i r e control volume i s necessary, and a l l the unknowns inside the c o n t r o l volume must be calculated. This surface g r i d generation requirement of the Boundary Integral Method o f f e r s a r e l a t i v e l y easy management of the r e s u l t s for the moving boundary i n comparison with other numerical methods e s p e c i a l l y i n three dimensional a p p l i c a t i o n s . With a boundary i n t e g r a l equation derived from Green's Theorem, a system of l i n e a r equations can be written with the p o t e n t i a l function and t h e i r d e rivatives at d i s c r e t e points on the control surface as unknowns. Upon so l v i n g t h i s set of equations, the p o t e n t i a l d i s t r i b u t i o n on the c o n t r o l surface i s obtained, and values at any point in s i d e the 17 control volume can be calculated. Although the Boundary Integral Method may not be the most e f f i c i e n t numerical method, i t i s considered to be a very e f f e c t i v e numerical approach for moving boundary problems. One of the disadvantages associated with the Boundary Integral Method, c r i t i c i z e d by most F i n i t e Element Method users, i s on i t s generated matrix. This matrix, which i s neither 'band' nor 'sparse', must be solved by the very time consuming Gaussian Elimination Method.Therefore, not only more computer storage w i l l be required, but the computational e f f i c i e n c y w i l l be poor. However, since problems solved i n t h i s thesis are used as an i l l u s t r a t i o n of the formulation's v a l i d i t y , the comparison between d i f f e r e n t numerical methods i s of no concern here and i s l e f t to someone with the i n t e r e s t and adequate knowledge i n numerical an a l y s i s . I l l . 1 Boundary Integral Method The Boundary Integral Method was f i r s t introduced by Fredholm (1903) who proved the uniqueness of a s o l u t i o n by using a f i n i t e number of elements on the boundary of the considered domain. I t was not u n t i l 1963 when Jawson and Ponter confirmed the use of Green's function i n the Boundary Integral Equation that the method became popular. With the 18 continuous developments of the computer's capacity and performance i n the l a s t two decades, the method has become an e f f e c t i v e numerical t o o l , and i s now broadly applied on a great v a r i e t y of engineering problems ( see Banerjee and B u t t e r f i e l d , 1981). The Boundary Integral Equation which can be derived through the Weighted Residual Method for so l v i n g the Laplace Equation has a general form given as : 0(P) + I <f>(Q) G dS = <f>(Q) G dS (13), ,n where G i s the Green's function defined between two points, P and Q. Point P i s the point of i n t e r e s t inside the flow f i e l d , and Q i s a control point on S ( Figure 2 ). The above Green's function has a d i f f e r e n t form for two dimensional and three dimensional problems. I t i s given as : G = - \ - ln( \ ) ( 2D ) ( 3D ) (14), where r i s the distance measured between P and Q. 19 • / ( x ~ - x . ) 2 + ( y « - y . ) 2 ( 2 D ) Q P Q P / 2 2 2 ( x - x ) + ( y - y ) + ( z - z ) (3D) Q P Q P Q P The term G ^ can be cal c u l a t e d by taking the dot product of VG and n , where n i s the normal u n i t vector defined according to the l o c a t i o n of Q. In some conditions, when the problem involves a f l a t impermeable surface ( such as the bottom of the wave tank ), the Green's function can have an a d d i t i o n a l image term to represent the image of the source point. This permits a reduction on the number of elements or unknowns involved i n the problem ( Figure 2 ). The Green's function then has the following form : G = { ln( i ) + l n ( p ) } (2D) (3D) (15), 47rr 47rr' where r remains the same as i n equation (14), and r ' = /( X - X ) 2 + ( y + y + 2d )2 (2D) / Q p o p 2 2 2 ( x - x ) + ( y + y + 2d ) + ( z - z ) (3D) Q P Q P Q P d i s the mean water depth measured from the undisturbed free-surface ( y=0 ). 20 A physical i n t e r p r e t a t i o n of equation (13) i s that, i f the v e l o c i t y p o t e n t i a l and i t s normal d e r i v a t i v e ( i e . normal v e l o c i t y component ) are known on the control boundary, S, the p o t e n t i a l value at any point, P, inside the control domain can be evaluated. Therefore, i f S i s represented by N elements, where on each of these elements, e i t h e r the p o t e n t i a l value or i t s normal d e r i v a t i v e i s known, a single equation of N unknowns can be written with point P on one of the elements. Since there are N elements, N equations of N unknowns can be written, and the unknown p o t e n t i a l values on each element can then be solved. I t must be emphasized that i n such a case, P i s a point on S, a s p e c i a l consideration should be made i n the c a l c u l a t i o n of the i n t e g r a l when P i s at Q. A d e t a i l e d discussion on the d e r i v a t i o n of equation (13) and i t s applications can be found i n Brebbia (1978). There i s also a choice on the order of the boundary element being used. When zeroth order elements are adopted, the p o t e n t i a l value as well as i t s normal d e r i v a t i v e are assumed as constant along each element. Solutions are c a l c u l a t e d at elements' mid points. For f i r s t order elements, p o t e n t i a l values are assumed to vary l i n e a r l y within each element. Results are computed at the node points between adjacent elements. In general, i t i s true that the higher the 21 order of the elements, the better the r e s o l u t i o n and accuracy of the r e s u l t s w i l l be, and more complicated numerical i n t e g r a t i o n procedures w i l l be involved. Therefore, zeroth order elements are used here to s i m p l i f y the programming procedure of the studied problems. The c a l c u l a t i o n s of matrix entries based on zeroth order elements by equation (13) i s b r i e f l y discussed i n Appendix I. A d e t a i l e d discussion on the choice of element orders as well as the involved numerical integrations can be found i n Banerjee and B u t t e r f i e l d (1981) . In the following chapters, most of the two dimensional applications are based on a Gaussian 5 t h order numerical i n t e g r a t i o n procedure. Whereas the three dimensional problem i s solved by considering the Green's function at the geometric node point of the q u a d r i l a t e r a l elements rather than a complicated numerical i n t e g r a t i o n procedure. I l l . 2 E x t r a p o l a t i o n i n Time As mentioned i n section II.2, extrapolation of v e l o c i t y terms i n equation (10) i s required to e s t a b l i s h the time stepping s o l u t i o n scheme. These v e l o c i t y terms are approximated by a second order function of time. I t must be emphasized that t h i s order i s a r b i t r a r i l y chosen, and solutions should not be affected by using higher order functions provided the time increment value i s reasonably small. The extrapolation function being used here i s 22 9 = a t 2 + bt + c (17) where ? i s the function of i n t e r e s t , t i s the time v a r i a b l e , and a, b, and c are unknowns d i f f e r e n t at each element. I f the value of 5 i s known at time t-5t, t-25t, and t-35t, a, b, and c can be obtained from the s o l u t i o n of a 3x3 matrix. 9 can then be approximated at time t, and — i s given as 2at+b. Equation (17) i s used i n the following chapters to approximate the v e l o c i t y of f l u i d p a r t i c l e s on free-surface elements before the s o l u t i o n i s achieved at each time step. These extrapolated values are found to be good approximations to the computed solutions except at the zero s t a r t up stage. However, once the s o l u t i o n proceeds f o r some time steps, the extrapolated values follow the trend of the solutions and pr e d i c t the numerical r e s u l t s quite accurately. Extrapolation using a higer order function than equation (17) has never been used, and whether i t can r e s u l t i n a better s o l u t i o n cannot be t o l d . 23 IV. APPLICATIONS In t h i s chapter, the non-linear free-surface formulations and the s o l u t i o n procedure defined i n previous chapters are tested by modelling four d i f f e r e n t problems. They are : i ) Wave Tank Modelling i i ) Deformation of High Amplitude Wave i i i ) P r e d i c t i o n of Ship's Bow Wave iv) Wave D i f f r a c t i o n by C i r c u l a r Cylinder Comparisons are made to a v a i l a b l e experimental and a n a l y t i c a l r e s u l t s i n each case to evaluate the s o l u t i o n procedure. The developed s o l u t i o n procedure i s f i r s t applied to simulate a two dimensional wave tank of f i n i t e water depth. A p i s t o n type wave paddle i s assigned a sinusoidal motion for wave generation. The numerical r e s u l t i s compared to l i n e a r wave theory to c a l i b r a t e the procedure. Then, the deformation of high amplitude waves i s studied by a two dimensional, s p a t i a l l y p e r i o d i c wave model. D i f f e r e n t breaking waves are obtained as solutions to confirm the introduced free-surface formulations as well as the assumptions. Although the same study was reported by Longuet-Higgins i n 1976, modelling of such high amplitude waves i s believed to be a good te s t to e s t a b l i s h the effectiveness of the free-surface formulations. A f t e r the s o l u t i o n scheme i s confirmed, i t i s applied to p r e d i c t a ship's bow wave. A wedge-shaped model of constant 24 d r a f t i s used f o r studies. Numerical r e s u l t s are obtained at d i f f e r e n t draft-Froude numbers, and are proved to be within the acceptable l i m i t s when compared to experimental measurements. F i n a l l y , the three dimensional a p p l i c a b i l i t y of the s o l u t i o n procedure i s i l l u s t r a t e d by a wave d i f f r a c t i o n model. A surface p i e r c i n g c i r c u l a r c y l i n d e r i s considered i n a wave tank. Forces experienced by the cy l i n d e r are obtained at d i f f e r e n t wave frequencies. Comparisons are made to some published experimental r e s u l t s . Each of these problems i s reported i n d i v i d u a l l y i n a separate section. Solutions obtained i n each case are numerically stable. Smoothing of the r e s u l t i s never applied and i s considered as unnecessary. IV.1 Wave Tank Simulation The wave tank ( or towing tank ) i s probably the most common f a c i l i t y f o r ocean engineers and naval a r c h i t e c t s to test the performance of t h e i r designs. Although the main purpose of the tank i s to tow ship models i n calm water and measure the resistance, i t also has f a c i l i t i e s to generate waves to simulate a regular sea state or a random sea state during the experiment. This wave making process i s u s u a l l y operated by a wave maker located at one end of the tank. E l e c t r o n i c signals are usually input to the wave maker to 25 govern i t s motion at a desired frequency and amplitude. Waves are then generated and t r a v e l down the tank. IV.1.1 Formulation In t h i s section, the motion of the wave maker i s governed by the following equation : x = X ( 1 - e" C t ) sin wt (18), m m where X i s the piston's motion amplitude, and u> i s the m - C t angular frequency. The term ( 1 - e ) i s a ramp function and i s used to ensure that the wave maker s t a r t s i t s motion smoothly at time ( t ) equal to zero, c i s chosen as 2.303 to recover 99% of the piston's motion within 2.00 seconds. The v e l o c i t y of the wave pi s t o n i s , therefore, simply given by taking the deri v a t i v e of equation (18) as : -ct -ct v = X { ( 1 - e ) w cos u>t + c e s i n wt } (19) . X ro At time t=0, the wave maker s t a r t s at zero v e l o c i t y . In Figure 3, the considered problem i s i l l u s t r a t e d by a diagram. The wetted surface of the wave pi s t o n i s designated by a symbol . The free-surface i s denoted by , and the 26 v e r t i c a l wall on the r i g h t as S q . Since a mirror image assumption i s used as i l l u s t r a t e d , the tank bottom i s not included i n the numerical model. With the control surface, S, given as : S = S + S^ + S m f o equation (13) can be rewritten as 4>(?)+j<f>G^dS + S o I* G dS + ' 4> G dS ,n ,n S f m I* G dS + [• G dS n n s f m G dS + S o where the Green's function, G, i s given i n equation (15). Then, by su b s t i t u t i n g the corresponding boundary conditions into the above equation, the following form i s obtained : 2 n <4(P) + I 4> G d S + [ < £ ( G + — G ) dS + j <f> G dS J r ,n J ,n . 2 J Y ,n S S^ g ( 5 t ) S o f m F G dS + 2 f v n G dS (20), J x x S_ S f m where 27 n 5<f> 1 - k4> 2 + 4> 3 F ^ 2 g ( S t ) 2 n + 0 ' n + -5Z- -£-( </> 2 + <f> 2 ) >x x 2g Dt v r,x r , y (21) F i s obtained by dropping the z d e r i v a t i v e terms from 2 2 n equation (10) , and moving [ - — G <f> ] from the g ( 6 t ) 2 right-hand-side to the left-hand-side of equation (20) . <j> ,n on the wave maker, S , i s given by v n m " J x x For l i n e a r applications, equation (12) i s used instead of equation (10) i n equation (13). The corresponding i n t e g r a l equation becomes : 4>(.v) + dS + \ <f> ( G + > n g ( f i t ) 2 G ) dS + f ^ G d S + f v n G dS J J x x m dt G dS .n m (22), where 54> -1 / .-2 ,-3 (sty (23) An e s s e n t i a l difference between equation (20) and equation (22) i s on the d e f i n i t i o n of the c o n t r o l surface. In equation (20) , S and S,, are considered to be the m t instantaneous locations of the wave maker and the free 28 surface. However, i n equation (22), and only designate the mean p o s i t i o n of the wave pi s t o n and the undisturbed free-surface. Boundary conditions are not s a t i s f i e d at the exact boundary p o s i t i o n . The control surface i s considered moving i n time for the non-linear cases but remains stationary i n the l i n e a r numerical model. When these two equations are applied for wave tank simulations, the ph y s i c a l problem i s represented by a numerical model of N elements as i n Figure 3 . Then, depending on whether the non-linear or the l i n e a r free-surface boundary condition i s used, a set of equations obtained by these two equations i s considered with point P on each element. The procedure i s very s i m i l a r to the discussion i n Appendix I. Upon sol v i n g the r e s u l t i n g set of equations, the p o t e n t i a l value on each i n d i v i d u a l element i s obtained, and the v e l o c i t y components can be cal c u l a t e d . In non-linear problems, the free-surface i s redefined by considering the motion of f l u i d p a r t i c l e s with the resolved v e l o c i t y components. As for l i n e a r problems, the free-surface el e v a t i o n i s computed by the l i n e a r i z e d B e r n o u l l i Equation which has the following form : <f> + gr, = 0 at y=0 A s i m i l a r procedure i s used to proceed the s o l u t i o n at the 29 s p e c i f i e d time increment except that the corresponding i n t e g r a l equation using a l i n e a r free-surface condition has fewer terms. IV.1.2 Interpolation In t h i s example, the considered problem i s i n two dimensions, and the control boundary i s defined by l i n e segments as shown i n Figure 3. Although the element siz e can be chosen very small to increase the r e s o l u t i o n of the p o t e n t i a l value along the control surface, p o t e n t i a l values can only be obtained at d i s c r e t e points along the boundary. In order to define the free-surface p o s i t i o n at the next time step, i n t e r p o l a t i o n s and extrapolations of the p o t e n t i a l values between elements are usually required i n the c a l c u l a t i o n s . In chapter I I I , i t was mentioned that zeroth order elements are applied here to model the problems for s i m p l i c i t y . Therefore, an i n t e r p o l a t i o n function i s used to determine the v e l o c i t y of f l u i d p a r t i c l e s at the free surface. I t must be emphasized that when the free surface i s considered moving i n time, i t s motion i s simulated by a time-step displacement of the free-surface elements. This time-step displacement of a surface element i s the resultant motion by considering the motions of i t s two ends and not i t s 30 mid point movement. That i s , the element's length i s not a constant but changes through each time increment. Before v e l o c i t y values are inte r p o l a t e d between elements, the normal and tangential v e l o c i t y components at the mid point of each element are determined. The normal v e l o c i t y component at the free-surface element's mid point i s cal c u l a t e d by equation (10) without the z associated terms. As f o r the tangential component, i t i s evaluated by the change of <f> value along the free surface. Since these two v e l o c i t y components are obtained i n d i r e c t i o n s normal and tangential to the free- surface, i t s x-y components can be resolved by a vector transformation procedure. A f t e r the v e l o c i t y components on each element are known, t h e i r values between adjacent elements are inte r p o l a t e d by the following function : v = ax + by + c x (24), v = dx + ey + f y where and v^ are v e l o c i t y components between elements, x and y are the coordinates of elements, and a, b, c,...etc are unknown constants ( d i f f e r e n t from those i n equation (17) ). With equation (24), values on three successive elements are required f o r the inte r p o l a t i o n s ( see Appendix II for 31 d e t a i l s ). The unknown constants, a , b, c,...etc, are found i n each case so as to e s t a b l i s h the i n t e r p o l a t i o n functions, and the v e l o c i t y components between two elements can then be evaluated ( Figure 4 ). At the fluid-body i n t e r s e c t i o n point, the same c r i t e r i o n i s applied. Three elements are chosen with two of them at the free-surface and one located on the wave piston. In t h i s case, ( Figure 4 ) the i n t e r s e c t i o n point has a value equal to the v e l o c i t y of the piston, and v^ i s the only quantity to be determined. When v^ of the i n t e r s e c t i o n point i s known, i t s l o c a t i o n at the next time step can be predicted. Although i t i s very c r u c i a l to evaluate t h i s i n t e r s e c t i o n point for a proper representation of the free surface near the moving body, i t i s observed from the numerical studies that the approximation of the i n t e r s e c t i o n point by equation (24) i s quite r e a l i s t i c . A numerical problem r e l a t e d to t h i s s i n g u l a r i t y as reported i n L i n (1984) was never encountered. I t i s also believed that improvements can be achieved by a higher order function than equation (24). However, the e f f e c t of using a higher order function as well as the l i m i t a t i o n s of such an i n t e r p o l a t i o n scheme are not studied here. As long as the element s i z e i s small enough, predictions of the fluid-body i n t e r s e c t i o n point by 32 equation (24) seem to be a good approximation f o r the use of the zeroth order element. I V . 1 . 3 Results and Discussions F i r s t order solutions are ca l c u l a t e d by using the l i n e a r i z e d free-surface condition : where > + gf? - 0 at y=0 (25) , t Dt 2 fit The c a l c u l a t i o n of <j> i s expressed i n a f i n i t e d i f f e r e n c e form since d> i s ca l c u l a t e d at a constant time increment, fit. Solutions obtained by equation (22) are l a b e l l e d as l i n e a r solutions while computations through equation (20) are defined as non-linear solutions. In t h i s example, comparisons are mainly between numerical solutions and the f i r s t order a n a l y t i c a l theory. The f i r s t set of r e s u l t s which are presented i n Figure 5 i s on wave generations by a numerical wave tank. Linear free-surface conditions are used, and motion of the wave maker i s at 0.50 Hz. A motion amplitude of 0.3 m i s considered. The tank i s chosen a r b i t r a r i l y to be 40.00 m i n length, and 4.00 m i n depth. Time domain solutions are obtained at a 5t increment of 0.05 second. Results are i l l u s t r a t e d i n Figure 5. Quantities l i k e wave amplitude, phase v e l o c i t y , and wave length can then be obtained d i r e c t l y or i n d i r e c t l y from the p l o t s . Computations are repeated at d i f f e r e n t frequencies and the corresponding wave lengths are measured. In Figure 6, the ca l c u l a t e d wave lengths are p l o t t e d against the dispersion r e l a t i o n s h i p at d i f f e r e n t frequencies. Agreements between the l i n e a r numerical solutions and the dispersion r e l a t i o n s h i p are very good over a wide range of frequency. This d i s p e r s i o n r e l a t i o n s h i p obtained from l i n e a r theory i s defined by the following expressions : 2 k tanh kd % g c 2 = tanh kd (26) , 2 k k where k = (27) w i s the angular frequency, k being the wave number, d i s the water depth, and A i s the wave length. When w and d are chosen, the value of k can be solved by i t e r a t i o n of equation (26), and the wave length, A, can be obtained through 34 equation (27). Comparisons of the phase v e l o c i t y , c, are also made with equation (26). In Figure 7, the computed value i s denoted by symbols at d i f f e r e n t frequencies. Similar to the case i n Figure 6, the numerical predictions ( which are l i n e a r ) agree well with equation (26) . In addition to the wavelength and phase v e l o c i t y computations, wave amplitudes are also c a l c u l a t e d and compared to the wave maker theory by Dean and Dalrymple (1984). A deep water s i t u a t i o n ( kd = 4.0 ) i s considered i n Figure 8, and r e s u l t s are expressed i n wave amplitude to piston's motion amplitude r a t i o s ( Aw/Am ). Small discrepancies ( les s than 5% ) are found between the numerical solutions and the theory, which i s probably explained by a numerical truncation error. ( S i t u a t i o n i s considered as deep water when kd > n ) . With the above comparisons, r e s u l t s obtained by the introduced s o l u t i o n scheme with a l i n e a r free-surface condition are considered to be s a t i s f a c t o r y . However, i t was observed that numerical r e s u l t s are p a r t i a l l y a f f e c t e d by the number of elements per wave length. In Figure 9, a si n u s o i d a l curve i s represented by d i f f e r e n t numbers of l i n e segments. In order to maintain a smooth p r o f i l e , i t i s found that at 35 l e a s t 24 elements per wave length should be used. This observation, unfortunately, can only give some ind i c a t i o n s to the r e l a t i o n s h i p between the r e s o l u t i o n of a curve and the number of elements per wave length. When the simulated wave amplitude to length r a t i o becomes very high, more elements should be used. Although i t i s true that the smaller the element s i z e , the better the free-surface d e f i n i t i o n w i l l be, there i s always a l i m i t a t i o n on how small an element can be used. I f the size of elements i s too small, not only does the computation become very i n e f f i c i e n t , but numerical problems can possibly a r i s e i n the solutions when adjacent elements are too close to one another. This can be r e a l i z e d i f the v e l o c i t y i s changing r a p i d l y between adjacent free-surface elements. Any improper choice of the time step ( too large ) can r e s u l t i n elements crossing each other and causing the numerical model to f a i l . Solutions of wave generation with l i n e a r and non-linear free-surface condition are found i n Figure 10. The tank's dimensions used i n t h i s case are 80.00 m i n length, and 4.00 m i n depth. A wave p i s t o n with motion amplitude of 0.15 m i s moving at 0.25 Hz. Solutions are computed at a 0.05 second increment f o r 20.00 seconds. Free-surface elevations are calc u l a t e d f o r both solutions, and presented i n graphical form. 36 The non-linear s o l u t i o n which i s p l o t t e d on the r i g h t h a l f of Figure 10 i s very s i m i l a r to the l i n e a r s o l u t i o n . Non-linearity e f f e c t i s i n s i g n i f i c a n t i n the s o l u t i o n since the simulated wave amplitude i s small. However the c a l c u l a t e d wave p r o f i l e s i n the l i n e a r s o l u t i o n are not very smooth a f t e r computations have been c a r r i e d on for some period of time. The same problem i s not observed i n the non-linear s o l u t i o n . A possible explanation f o r t h e i r behaving d i f f e r e n t l y i s believed to be i n the l i n e a r i z a t i o n assumption that the l i n e a r boundary condition i s s a t i s f i e d at the undisturbed free surface. This error that i s t o l e r a t e d i n the l i n e a r model accumulates as the s o l u t i o n proceeds i n time. As for the non-linear s o l u t i o n , though some unce r t a i n t i e s e x i s t i n the free-surface predictions, r e s u l t s are generally smooth and remain numerically stable i n time. Since i t i s known that the difference between l i n e a r and non-linear s o l u t i o n i s s i g n i f i c a n t when wave amplitude i s high, and that the non-linear wave w i l l break when the wave slope i s steep, the comparison can only be c a r r i e d out i n a very narrow range of wave height to length r a t i o . However, the main purpose of Figure 10 i s to show the reader that the non-linear s o l u t i o n i s c a l i b r a t e d and better than the l i n e a r s o l u t i o n . A vigorous comparison of the two solutions i s beyond the scope of t h i s study There e x i s t s a r e l a t i o n s h i p between the wave's breaking 3 7 mode and the i n i t i a l wave height to wave length r a t i o . When the wave amplitude i s very high, i t i s observed i n nature that the wave w i l l deform and eventually f a i l at d i f f e r e n t breaking modes. This i s the next topic to be discussed i n the following section. F i n a l l y , the d r i f t motion of f l u i d p a r t i c l e s at the free surface i s studied. Figure 11 i s the p l o t of the non-linear free surface at time t equal to 5.00 and 6.00 seconds. The simulated wave tank's dimensions are 10.00 m i n length and 6.00 m i n depth. The motion amplitude of the wave maker i s at 0.10 m, and waves are generated at a frequency of 0.5 Hz. 50 elements are used at the free surface, and the time step i s chosen as 0.05 second. Solutions e x h i b i t a wave length of approximately 6.00 m, and a wave height at about 0.36 m. The equivalent wave height to length r a t i o (H/A) i s 0.06, which i s s t i l l within the 1/7 limit ( t h e o r e t i c a l l i m i t before breaking ). A sharper crest and shallower trough are observed i n the so l u t i o n . F l u i d p a r t i c l e s at the free surface are also traced i n time, and t h e i r l o c i are recorded i n Figure 12. P a r t i c l e s at three locations inside the tank are chosen at 0.20 m, 2.00 m, and 6.00 m from the wave maker. I t can be observed i n Figure 12 that these p a r t i c l e s a l l experience a non-linear d r i f t motion. At 0.2 m from the wave maker, f l u i d p a r t i c l e s , though 38 mainly following the piston's motion, are slowly d r i f t i n g to the r i g h t . Whereas at 2.00 m from the piston, f l u i d p a r t i c l e s are d r i f t i n g down the tank at a mean d r i f t v e l o c i t y of 0.1 m/s. According to the p r e d i c t i o n by Newman (1977), a f l u i d p a r t i c l e ' s second order motion i s given as : dx , , - 3 — = oA e ^ cos ( kx-wt ) + cokA e 2 ^ + 0( A 3) dt where A i s the wave amplitude, k i s the wave number, and co i s the frequency. X q i s given i n the Lagrangian notation, and y i s taken at the undisturbed free- surface. In t h i s equation, the only term that i s responsible f o r the non-linear d r i f t 2 2ky motion i s the second order term ( wkA e J ) . This second order mean d r i f t v e l o c i t y i s c a l c u l a t e d as 0.1066 m/s, which indicates that the numerical p r e d i c t i o n has an differ e n c e of about 7% . At 6.00 m down the tank, f l u i d p a r t i c l e s do not experience any d r i f t motion u n t i l a well developed wave cr e s t a r r i v e s . The motion due to the f i r s t wave cycle d r i f t s the p a r t i c l e s at a greater distance than t h e i r migration i n l a t e r wave cycles. Since these p a r t i c l e s are close to the tank wall, t h e i r motions i n l a t e r cycles are believed to be affec t e d by the r e f l e c t e d waves. With the above observations i n the numerical solutions, 39 i t i s believed that the natural d r i f t phenomenon i s obtained by the proposed free-surface formulations. However, although the numerical solutions are stable within the considered simulation time, f a i l u r e of the numerical model i s expected for a longer simulation time. This f a i l u r e i s r e l a t e d to the improper s i z i n g of free-surface elements near the wave maker, and to the use of the i n t e r p o l a t i o n function ( equation (24)) at the fluid-body i n t e r s e c t i o n point. Moreover, as the s i z e of the elements near the wave maker i s growing i n time ( a consequence of elements d r i f t i n g downstream ), r e s o l u t i o n of the computed wave p r o f i l e near the p i s t o n w i l l become so poor that the s o l u t i o n down stream w i l l be a f f e c t e d and become imprecise. Improvement to the s o l u t i o n i s possible by breaking an element which becomes too long into two elements during the computation. However, t h i s involves not only a reconstruction of the matrix, but the i n i t i a l condition and the p o s i t i o n h i s t o r y of the newly created elements must be defined. This i s considered very d i f f i c u l t and i s not being c a r r i e d out i n t h i s study. IV.2 Deformation of High Amplitude Waves The s o l u t i o n procedure outlined above i s applied to the same problem investigated by Longuet-Higgins and Cokelet i n 1976. This study on the deformation of high amplitude waves i s considered a straight-forward way to t e s t how well the 40 introduced free-surface formulations can model a highly non-linear s i t u a t i o n . A two dimensional, s p a t i a l l y p e r i o d i c wave model i s used to simulate a high amplitude wave u n t i l i t becomes unstable. I n i t i a l conditions are obtained from l i n e a r wave theory to s t a r t the computation with a well developed wave p r o f i l e . Although t h i s numerical model by Longuet-Higgins seems a r t i f i c i a l , i t i s s t i l l considered as an e f f e c t i v e way to study breaking wave phenomenon. In r e a l l i f e , wave breaking i s expected to occur when the v e l o c i t y of f l u i d p a r t i c l e s exceeds the phase v e l o c i t y at the free surface. In other cases, i t can also happen when the p a r t i c l e s ' a c c e l e r a t i o n i s larger than the g r a v i t a t i o n a l a c c e l e r a t i o n . I t can be observed l a t e r i n Chapter IV.2.2 that p a r t i c l e s at the wave crest are moving at a r e l a t i v e l y higher v e l o c i t y than those at the trough. As a r e s u l t , the cr e s t t r a v e l s f a s t e r and eventually forms a j e t . This wave form cannot be described with l i n e a r theories as the function that defined the wave p r o f i l e i s multi-valued. The numerical computation stops when p a r t i c l e paths d e f i n i n g the free surface i n t e r s e c t each other. This condition i s taken to represent numerical wave breaking. IV.2.1 Formulation In Figure 13, numerical models f or l i n e a r and non-linear 41 applications are i l l u s t r a t e d . These numerical models are at exactly one wave length so as to take advantage of the s p a t i a l l y p e r i o d i c assumption. F i n i t e water depth i s considered, and a mirror image i s assumed about the bottom surface i n order to reduce the number of unknowns involved i n the problem. The modified Green's function ( equation (15) ) i s used i n the Boundary Integral Equation. The control surface i s divided into S,. , S , and S' . f r r S_ represents the free- surface, while S and S' denote two f r r r a r t i f i c i a l , v e r t i c a l boundaries f o r the completion of the cont r o l domain. According to the s p a t i a l l y p e r i o d i c assumption, and S^ , are one wave length apart, and the p o t e n t i a l d i s t r i b u t i o n as well as v e l o c i t y values should also be i d e n t i c a l at corresponding points ( Figure 13 ) . Therefore, i f N elements are used on the t o t a l boundary and element i i s on , the following conditions should be true: 1 H + l - l v x ' i = V x L + i - i (28), where the subscripts are element numbers. These a d d i t i o n a l conditions provide an extra set of equations to solve the problem. With 42 s = s + r S f + S' r the corresponding i n t e g r a l equation can be written as : (f>l?) + \ 4> G dS + ^ ( G + 2 - G ) dS J ,n J ,n G ( 5 T ) 2 S + S' r r f G dS + 2 ( v x n x ) G dS (29), S_ S + S' f r r where F i s as given i n equation (21). However, since v on 2 X and i s s t i l l not known, i t s associated i n t e g r a l terms could be moved to the left-hand-side of the equation to give the following form : 2n 4>{V) + \ 4> G d S + 4> ( G + 2L G ) dS + I v n G dS J ,n J Y ,n , . , 2 J x x S + S* Sc &K° } S + S' r r f r r = | F 2 G dS (30) . S f 2 Now, with the v associated terms as extra unknowns, i t x i s evident that there are more than N unknowns to be solved though only N equations can be written by equation (30) . Fortunately, equation (28) provides an a d d i t i o n a l set of conditions to complete the. matrix with the r i g h t number of equations. Thus, the v e l o c i t y on both and are solved simultaneously with the p o t e n t i a l values. A complete form of the established matrix can be found i n C a l i s a l and Chan 43 (1987) . For l i n e a r r e s u l t s , the l i n e a r i z e d free-surface conditions are applied, and equation (30) i s rewritten as : <(><.?) + [ <f> G dS + \ <f> ( G + —-— G ) dS + f v n G dS J ,n J ,n ( c r \ 2 J x x S + S' S., g^t.; S + S' r r f r r = J K G dS (31) , S f where K i s as given i n equation (23). This equation, because of l i n e a r i z a t i o n , has a d i f f e r e n t d e f i n i t i o n f o r the con t r o l surface. Boundary conditions are applied at the undisturbed free surface, while S^and are considered as not moving i n time. In the non-linear model, , , and are moving f r e e l y i n time. Since the two v e r t i c a l boundaries are at the same v e l o c i t y , the single wave length r e s t r i c t i o n can always be maintained. In t h i s problem, since there e x i s t s no external d r i v i n g mechanism ( such as the wave maker ) to govern the change of the free surface i n time, function and K i n equation (30) and (31) are considered as the only d r i v i n g functions f o r the solutions to proceed i n time. In fa c t , during the moment of s t a r t up, these two functions only contain information of a 44 well developed free-surface wave, and wave theory of any order can be used to ca l c u l a t e F and K at time t=0 . For 2 s i m p l i c i t y , the l i n e a r wave theory i s chosen here. The a n a l y t i c p o t e n t i a l function from l i n e a r theory i s given i n Newman (1977) as : Ag cosh k(y+d) s i n ( fcx-wt ) (32) , co cosh kd where A i s the wave amplitude, w i s the angular wave frequency, d i s the mean water depth, and k i s the wave number defined as i n equation (27). At time t = 0, wave elevation rj i s ca l c u l a t e d by equation (25) as : = A cos kx at y=0, t=0 (33). The free-surface boundary f or non-linear simulation i s then defined by t h i s expression. As for the functions and K, they can be evaluated by s u b s t i t u t i n g equation (32) into equations (21) and (23). IV.2.2 Results and Discussions The numerical s o l u t i o n with a l i n e a r free-surface 45 boundary condition i s obtained and i l l u s t r a t e d i n Figure 14. The considered wave length i s 10.00 m with a depth at 4.00 m. I n i t i a l wave amplitude i s assigned as 0.50 m, and 50 elements are used at the free surface. The corresponding wave height over length r a t i o ( H/A ) i s 0.10, and the height to depth r a t i o ( H/d ) i s 0.25 . The s o l u t i o n i s allowed to proceed for 20.00 seconds at a time increment of 0.05 second. Wave p r o f i l e s c a l c u l a t e d by equation (25) are i l l u s t r a t e d at 19.00, 19.50, and 20.00 second. Results are observed to be smooth and propagate at constant speed. At 20.00 seconds, the wave amplitude has l o s t 2% of i t s o r i g i n a l magnitude, which i s believed to be a consequence of numerical truncation errors. These numerical r e s u l t s obtained with a l i n e a r free-surface condition are, therefore, considered as numerically stable and compatible to l i n e a r theory within 2% error. As f o r non-linear solutions, the change of wave shape i s so dramatic that i t eventually f a i l s at d i f f e r e n t breaking modes. In the following r e s u l t s ( Figures 15, 16, 17 ), wave length i s chosen at 10.00 m, and 1.00 m of water depth i s considered. 50 elements are used to represent the free surface, and solutions proceed at 0.01 second increments u n t i l breaking waves are obtained. In Figure 15, the H/A r a t i o i s selected at 0.10 . The 46 corresponding H/d r a t i o i s 1.00 . Wave p r o f i l e s are calc u l a t e d u n t i l a s p i l l i n g breaker i s obtained which stops the computation at 1.25 second. The breaking time to wave period r a t i o ( T /T ) i s therefore c a l c u l a t e d as 0.3686 . The b numerical model i s found unstable at the wave cr e s t where the f l u i d flow condition changed r a p i d l y . Elements i n t h i s region become so close to each other that eventually, they cross each other and f a i l to represent the free-surface anymore. At higher H/A r a t i o , the wave possesses more energy per wave length, and formation of a plunging breaker at a much shorter time i s expected. This i s i l l u s t r a t e d i n Figures 16 and 17 . In these two cases, the i n i t i a l H/A r a t i o s are at 0.125 and 0.150 , and the respective H/d r a t i o s are 1.25 and 1.50. Plunging breakers are formed at 0.87 and 0.75 second, equivalent to T /T r a t i o s of 0.2565 and 0.2212 re s p e c t i v e l y . Although a c a r e f u l study of the r e s u l t s i n Figures 15, 16, and 17 shows that the wave crests are not as smooth as those reported by Longuet-Higgins and Cokelet (1976) , they behave very s i m i l a r l y , and f a i l at s i m i l a r breaking modes. This imperfection i s not sur p r i s i n g , since smoothing i s never applied to the cr e s t region where the flow condition i s changing so ra p i d l y . With the provided r e s u l t s , there i s enough evidence to 47 believe that the introduced free-surface formulations are adequate for modelling the free-surface behavior, and the s o l u t i o n procedure provides an a l t e r n a t i v e to the method by Longuet-Higgins (1976). IV.3 Predictions of Ship's Bow Wave In t h i s section, a more p r a c t i c a l problem w i l l be used to demonstrate the s o l u t i o n procedure, as well as to t e s t the body free-surface i n t e r p o l a t i o n function. The p r e d i c t i o n of ship's bow waves has been considered as one of the important topics f o r Naval A r c h i t e c t s . The understanding of such a phenomenon i s considered important because i t greatly a f f e c t s a ship's performance. A major portion of the ship's wave resistance i s believed to be contributed by the flow conditions near the bow. This problem, which i s a high Froude number problem ( as indicated by O g i l v i e (1972) ), i s formulated within a slender ship assumption. Such an assumption not only means a small parameter, e, i s found ch a r a c t e r i z i n g the beam/length or draft/length r a t i o , i t also implies that the s i z e and shape of the c r o s s - s e c t i o n a l h u l l form changes gradually i n the l o n g i t u d i n a l d i r e c t i o n . In t h i s section, a wedge-shaped model i s used to simulate the ship's bow. A wide range of the draft-Froude number i s considered i n the c a l c u l a t i o n s . Solutions are 48 obtained with l i n e a r and non-linear free-surface conditions, and comparisons are made to the experimental and a n a l y t i c a l work by O g i l v i e (1972). IV.3.1 Formulation When a ship i s moving at a constant speed, U , a steady wave pattern moving with the ship i s formed at the bow. This pattern, which i s three dimensional, i s affe c t e d by the ship's speed as well as the form of the ship's bow. In order to s i m p l i f y the problem, a wedge-shaped bow of constant d r a f t i s considered. A r i g h t hand car t e s i a n coordinate system i s used with i t s o r i g i n located at the bow. The x axis i s pointing i n the l o n g i t u d i n a l d i r e c t i o n of the ship, and y i s chosen p o s i t i v e i n the upward d i r e c t i o n . In Figure 18, the model i s moving at a constant speed, U , i n the - x d i r e c t i o n . Water depth i s considered to be large compared to the model's d r a f t . With the slender ship assumption, the three dimensional wave pattern can be treated as a two dimensional problem and solved along the ship length at constant x increments. The same argument was made by O g i l v i e (1972), that i f the ship's h u l l form i s defined by a function z = ± B ( x , y ) , the slope of the waterlines i n the x d i r e c t i o n can be assumed small since the change i n the lo n g i t u d i n a l d i r e c t i o n i s gradual. 49 Therefore, B - 0 ( e ) 0 < x < L s where B i s a geometric function that described the h u l l form, L g i s the ship's length, and e i s any small parameter such as the beam/length r a t i o . Whereas i n the bow near f i e l d , since the rate of change of flow condition i n the l o n g i t u d i n a l d i r e c t i o n i s higher than that i n the far f i e l d , the order of x should be smaller than the usual t h i n ship assumption ( i e . 0(x) = 1 ), but remain less than e. Therefore, The order of d i f f e r e n t i a l operators i s , therefore, given as : The associated order of terms for equation (1), then becomes: x = 0( e 1/2 ) and R = ( y 2 + z2 ) 1 / 2 - 0 ( « ) 50 + <k + <t> = o ,xx r,yy r , z z This indicates that the x associated terms are much smaller than the other terms, and solutions can be obtained by t r e a t i n g the problem as i n the y-z plane. The above order of magnitude analysis follows the arguments given by O g i l v i e (1972). As a consequence, the governing equation becomes : Now, since the model i s moving at constant speed, U, the t o t a l v e l o c i t y p o t e n t i a l , $ , can be assumed to have the following form : where <f> i s a function of y and z only. When equation (35) i s substituted into equation (4), and the pressure P , i s taken as zero at atmospheric pressure, the B e r n o u l l i Equation can be rewritten as : .yy + 4> zz - 0 (34). $ = Ux + 4>(y,z) (35), at y=r7 (36) . 51 I t must be c l e a r that, the term $ i s equal to zero since * t the wave pattern i s steady at the ship's reference. However, i t i s also true that at a f i x e d reference, the ship i s moving i n the -x d i r e c t i o n at v e l o c i t y -U. Tracing the ship's motion i n time at a reference f i x e d i n space i s given by : and the change along the ship at the f i x e d reference can be written as : Equation (36) i s therefore re-expressed at a f i x e d reference as : where the change of 4> i n the x d i r e c t i o n at the ship's reference i s expressed as a change i n time. Following the approach i n chapter I I , equation (39) i s combined with equation (3) to y i e l d the combined free-surface boundary condition at y=r). This i s given as : x = Ut (37) 5x = U St (38). 52 2n * T, = y J o n ,-1 y 5^ 4^~2 + <ft~3 (*t)' n \ z z 2g Dt v \ y *\z ' at y=r; (40), where n i s defined on the y - z plane. As for the impermeable boundary that represents the ship h u l l , </> = U ( f r ^ n + f r ^ n ) on z = B(x,y) ,n 3x y dx z J However, since a wedge-shaped model of constant d r a f t i s 3y considered, -J- i s taken as zero on the ship h u l l . Moreover, ox i f the wedge h a l f angle i s a , the above condition can be written as : <f> - U B n >n ,x z = U tan a n^ on z = B(x,y) (41). A schematic drawing for the y - z plane i s shown i n Figure 19. In Figure 19, the i l l u s t r a t i o n i s very s i m i l a r to a wave tank problem. The only difference between Figure 19 and Figure 3 i s that the wave maker i s now being replaced by the cross section of the ship h u l l . This cross section of the ship h u l l i s moving i n the z d i r e c t i o n at constant v e l o c i t y instead of o s c i l l a t i n g back and f o r t h . A more precise way to 53 explain Figure 19 i s to consider the i l l u s t r a t e d y-z plane moving along the ship length, so that the ship's section appears on the figure as i f a v e r t i c a l p a r t i t i o n i s moving to the r i g h t at constant speed. Therefore, a s i m i l a r approach to the modelling of a wave tank can be applied i n t h i s problem, except that the time v a r i a b l e , t , i s used instead of x to proceed the s o l u t i o n i n the y-z plane along the ship's length d i r e c t i o n . F i n a l l y , the boundary condition on the v e r t i c a l wall at the r i g h t end of the considered domain i s simply defined as zero. As i l l u s t r a t e d i n Figure 19, only the r i g h t h a l f of the considered problem i s modelled because of symmetry about the ship's mid plane. The v e r t i c a l wall i s located very f a r away from the ship to reduce any possible interference on the r e s u l t by r e f l e c t e d waves. Water depth i s also taken to be deep enough to avoid any shallow water e f f e c t i n the s o l u t i o n . I t must be c l e a r that, because of the slender ship assumption, wave ele v a t i o n ahead of the bow i s not formulated i n the s o l u t i o n . Since i t i s assumed not a f f e c t i n g the s o l u t i o n by a s i g n i f i c a n t amount, the flow condition at the bow i s considered as undisturbed to s t a r t the computation. With the control surface S, divided into S , S , . , and 54 , as i l l u s t r a t e d i n Figure 19, the associated Boundary Integral Equation can be written as : c A ( P ) + J c 6 t 7 n d S + 2n * ( G , n + G ) dS + g(5t) "o " f I F G dS + [ U tan a n G J 2b J Z dS <l> G dS .n (42) where F 2b g y 5<p - 4c6 + 0 n + ^' n + -=2_ JL( ^ 2 + ^ 2 )• *\z z 2g Dt v r , y r , z (43) As for l i n e a r a p plications, the i n t e g r a l equation i s , of course, very s i m i l a r except simpler. The l i n e a r i z e d boundary conditions are applied at the undisturbed free surface, and the ship's thickness i s ignored. This i n t e g r a l equation i s written as : cA(P) + f 6 G d S + f < 4 ( G + — " — G ) dS + f d> G J ,n J ,n , 2 J ,n s g ( 5 t ) su o f h K G dS + | U tan a n^ G dS S f S h dS where K i s defined as i n equation (23). 55 IV.3.2 R e s u l t s and Discussions The following r e s u l t s are obtained by considering wedge models of 15° and 30°. The model i s at d i f f e r e n t d r a f t s and moving at a wide range of v e l o c i t i e s . The chosen tank width i s taken as 10.00 feet to reduce the possible interference on the r e s u l t s by r e f l e c t e d waves. Water depth i s also taken at 5.00 feet to avoid any shallow water e f f e c t f o r the considered v e l o c i t y range. Linear and non-linear r e s u l t s are obtained and compared to the measurements by O g i l v i e (1972). An example of the wave elev a t i o n c a l c u l a t e d along the model i s shown i n Figure 20. The wedge's h a l f angle ( a ) i s 7.5°, and the wedge i s moving at 8.00 f t / s e c . The d r a f t i s considered to be 8.00 inches. Solutions are ca l c u l a t e d at increments of 0.05 feet i n the l o n g i t u d i n a l d i r e c t i o n of a 5.00 feet model. The obtained non-linear wave c r e s t i s observed to be higher than the l i n e a r r e s u l t , and i t i s located c l o s e r to the bow. A secondary crest, which i s absent i n the l i n e a r s o l u t i o n , i s observed to show a d i f f e r e n t c h a r a c t e r i s t i c i n the non-linear s o l u t i o n . The peak of the bow wave ( Am ) and the l o n g i t u d i n a l p o s i t i o n of i t s occurrence are cal c u l a t e d at d i f f e r e n t d r a f t s and v e l o c i t i e s . Both l i n e a r solutions and non-linear 56 solutions are p l o t t e d against the experimental r e s u l t s by Og i l v i e (1972), as shown i n Figure 21. In order to avoid confusion, Am i s p l o t t e d against v e l o c i t y U, and presented separately according to t h e i r d r a f t s ( T ). In Figure 21, comparisons are made i n dimensional units to give some i n d i c a t ion of the scale range being used. The l i n e a r predictions of the wave peak obtained on the wedge model behave l i n e a r l y with the model speed. The slope of t h i s s t r a i g h t l i n e increases as model d r a f t increases. In Figure 21, at T = 4 inches, the non-linear s o l u t i o n though over p r e d i c t i n g O g i l v i e ' s measurements, agrees with the experimental r e s u l t s much better than the l i n e a r solutions. Moreover, the trend of the non-linear predictions are s i m i l a r to that of the experimental values. However, the hump-hollow behavior of the experimental r e s u l t s i s absent i n the non-linear p r e d i c t i o n s . This i s probably due to an inadequacy of the slender body assumption that i s made i n the formulations. The hump-hollow behavior of the bow wave amplitude at increasing speed i s known to be a cont r i b u t i o n of the transverse wave. As the problem i s modelled and s i m p l i f i e d to be solved i n two dimensions, i t i s expected that t h i s hump-hollow behavior w i l l be absent i n the sol u t i o n . As the model d r a f t increases, the over-prediction of the 57 non-linear s o l u t i o n becomes more and more serious. At T = 16 inches, the non-linear s o l u t i o n can no longer p r e d i c t the measured values p r e c i s e l y . In order to study the differ e n c e between the experimental r e s u l t s and the numerical solutions, the same figure i s reproduced based on a non-dimensional scale. In Figure 22, the same r e s u l t s from Figure 21 are reproduced i n non-dimensional units. Y i s used instead of max Am, and a draft-Froude number, F, i s used instead of U. Their r e l a t i o n s h i p s are given by the following equations : F — (44), and Y = ^ y. fJ/T Am (45) , max 2aU / °' From Figure 22, the l i n e a r p r e d i c t i o n of Y i s found to be ° r max at a constant value, independent of the draft-Froude number, F. This constant Y value from l i n e a r solutions appears to max be quite i n s e n s i t i v e to the model d r a f t , T. I t va r i e s s l i g h t l y from 1.5 at T = 4.00 inches to 1.7 at T = 16.00 inches. The a n a l y t i c Y value obtained by O e i l v i e has a J max constant value at 1.65, and i s shown as a h o r i z o n t a l s t r a i g h t l i n e . 58 The non-linear solutions of Y i n Figure 22 agree max & ° better i n the shallow d r a f t case than i n the deep model d r a f t s i t u a t i o n . Although the case where T = 16 inches i s compared at a lower draft-Froude number, i t i s obvious that the l i n e a r solutions p r e d i c t the experimental r e s u l t s more accurately than the non-linear solutions. As one can see i n Figure 22, because of the over-prediction i n the non-linear solutions, the non-linear s o l u t i o n can no longer be claimed as a better p r e d i c t i o n than the l i n e a r s o l u t i o n . However, i r r e s p e c t i v e of i t s inaccuracy, the non-linear s o l u t i o n s t i l l e x h i bits a trend that i s followed by the experimental data. Such a trend on Y i s never shown i n the l i n e a r s o l u t i o n , max This discrepancy between the non-linear solutions and the experimental r e s u l t s i s believed to a r i s e by d i f f e r e n t f a c t o r s . F i r s t l y , although the slender body assumption s i m p l i f i e s the problem to be solved i n two dimensions, i t has induced some error into the s o l u t i o n . In Figure 22, i t i s observed that the hump-hollow behavior does not appear i n the solutions, which means that the transverse wave i s not included i n the s o l u t i o n . This transverse wave which contributes to ship resistance i s expected to counteract and a f f e c t the o v e r a l l wave pattern of a ship moving i n a calm sea. As a r e s u l t , the f a i l u r e of the non-linear s o l u t i o n at deep model d r a f t i s believed to be p a r t i a l l y a f f e c t e d by an 59 inadequacy of the slender body assumption. Anothor reason f o r the inaccuracy appearing i n the non-linear s o l u t i o n i s believed to be i n the i n v i s c i d assumption. For wave generation by moving objects, the p o t e n t i a l flow s o l u t i o n i s known as an adequate s o l u t i o n f o r the outer flow region. The inner s o l u t i o n ( i e . very near the fluid-body i n t e r s e c t i o n ) requires a more rigorous approach. Although an i n t e r p o l a t i o n scheme used to take care of the s i n g u l a r i t y at the fluid-body i n t e r s e c t i o n point i s derived i n chapter IV.1.2, a rigorous v e r i f i c a t i o n of t h i s scheme has never been c a r r i e d out. However, i t i s believed e s p e c i a l l y i n the bow wave problem, that the bow wave amplitude i s damped by the v i s c o s i t y of f l u i d . Since an i n v i s c i d assumption i s made i n the formulation, over-prediction i s expected i n the so l u t i o n . F i n a l l y , the numerical model has a d i f f e r e n t i n i t i a l c ondition than i n the r e a l l i f e phenomenon. In the r e a l l i f e s i t u a t i o n , there i s a free surface elevation ahead of the ship's bow, which i s not a small magnitude. However, i n the numerical model, because of s i m p l i c i t y , a zero i n i t i a l c ondition i s used at the wedge's leading edge. Therefore, with a d i f f e r e n t i n i t i a l condition, a difference i n the p r e d i c t i o n of the bow wave amplitude i s expected. Perhaps i t i s more sui t a b l e to declare that the assumptions, the 60 formulations, and the slender body approach have r e s t r i c t e d the a p p l i c a b i l i t y of the s o l u t i o n to a shallow d r a f t problem. When model d r a f t i s deep, some corrections or remedies must be made. At low draft-Froude numbers, both solutions f a i l to pr e d i c t the measured Y values. In f a c t , t h i s i s not max s u r p r i s i n g since upon the acceptance of the draft-Froude number as an independent parameter, i t i s r e a l i z e d that the problem under i n v e s t i g a t i o n i s e s s e n t i a l l y a high Froude number problem, and the low Froude number solutions may not be v a l i d . Therefore, as the draft-Froude number approaches zero, predictions w i l l become more and more u n r e l i a b l e . Moreover, the slender ship assumption i s also being contradicted at the bow, which was explained by O g i l v i e (1972) i n d e t a i l . From a general point of view, although the non-linear solutions always over p r e d i c t the measurements, the corrections ( to the l i n e a r r e s u l t s ) with a non-linear free-surface condition are i n the r i g h t d i r e c t i o n . Since surface tension and viscous e f f e c t s are not included i n the formulations, over predictions i n the s o l u t i o n can be expected. In Figure 23, the same comparisons are made with a 30° 61 wedge model. The numerical predictions are observed to be less accurate than the r e s u l t s f o r the 15° wedge. A possible explanation f o r t h i s i s the considered wedge angle being too large, making the slender ship assumption no longer v a l i d . The l o n g i t u d i n a l positions ( X ) where wave peaks ° max r occurred are also c a l c u l a t e d i n Figure 24. From most of the ca l c u l a t i o n s , the f i r s t peak measured along the wedge model i s higher than i t s successive one. I t i s again observed that the l i n e a r predictions of X ex h i b i t a l i n e a r r e l a t i o n s h i p r max r with the model's v e l o c i t y . The slope of the s t r a i g h t l i n e s drawn through the computed X values increases at deeper d r a f t . Although the non-linear solutions over p r e d i c t the measured values, they have corrected the l i n e a r r e s u l t s i n the r i g h t d i r e c t i o n . At T = 16.00 inches, the non-linear predictions tend to agree with O g i l v i e ' s measurements quite well i n the high Froude number range. In order to v i s u a l i z e the computed solutions, some of these r e s u l t s are presented i n graphical form. In Figure 25, a model at 12.00 inches d r a f t and with a h a l f wedge angle of 7.5° i s moving at 5.00 f t / s e c . Linear and non-linear r e s u l t are computed f or graphical i l l u s t r a t i o n s . The model thickness i s ignored i n the l i n e a r numerical s o l u t i o n , and no s i g n i f i c a n t differences are found between the two solutions since the v e l o c i t y i s considered low. 62 At a shallower d r a f t ( 8.00 inches ) and higher model speed ( 8.00 f t / s e c ), the differences observed i n the two solutions become quite obvious. As i n Figure 26, the non-linear wave generated by the wedge exhibits a higher crest and tra v e l s at a higher speed i n the transverse d i r e c t i o n . As a consequence, the bow wave angle observed i n the non-linear s o l u t i o n i s larger than that i n the l i n e a r s o l u t i o n . In order to show the formation of a secondary bow wave, the computation i s repeated f o r a model with a larger wedge angle. In Figure 27, a 30° wedge at 12.00 inches d r a f t i s moving at 4.00 f t / s e c . A secondary bow wave i s e a s i l y observed i n the so l u t i o n . I t can be seen that i n t h i s case, the f i r s t wave decays very f a s t i n the l o n g i t u d i n a l d i r e c t i o n . At approximately 2.50 feet down the wedge's leading edge, i t becomes almost i n s i g n i f i c a n t . Whereas the secondary wave picks up i t s amplitude at about 3.00 feet down the wedge's leading edge. F i n a l l y , a 15° model at 12.00 inches d r a f t i s considered at 14.00 f t / s e c . As the model speed i s very high, spray formation along the model may be expected. I t i s remarkable that t h i s natural phenomenon i s obtained i n the s o l u t i o n without any numerical problem. In Figure 28, the spray 63 formation can be seen c l e a r l y along the model. Although t h i s i s not reported by O g i l v i e , such a natural phenomenon i s not unexpected when the model speed i s high enough. Before entering the next section, i t i s worth while to mention the weakness involved i n the use of the slender ship assumption. When a ship i s assumed to be slender, not only are i t s beam and d r a f t considered as small compared to i t s length, i t also c a r r i e s the meaning that the ship's geometry i s changing gradually i n the length d i r e c t i o n . Moreover, the rate of change of f l u i d flow conditions i s assumed to be very large i n the transverse d i r e c t i o n i n comparison with the l o n g i t u d i n a l d i r e c t i o n . However, near the bow, because of the very sudden change of flow around the bow, the assumption of the i n e r t i a force i s comparable to the gravity force can be l o c a l l y i n v a l i d . Upon the acceptance of the draft-Froude number, i t i s already r e a l i z e d that the bow wave problem i s a high Froude number problem, and the low Froude number so l u t i o n cannot be expected to be very accurate. It i s also true that, at some distance along the wedge, the wave peak occurs at a region where model thickness can no longer be assumed small. There i s always a question on how t h i n the model should be so as to v a l i d a t e the slender ship assumption . From the solutions, the wave height predictions i n Figure 23 are obviously worse than the p r e d i c t i o n s i n 64 Figure 22, simply because the wedge angle i s larger. Therefore, when one considers the factors that cause the discrepancies between the c a l c u l a t i o n s and the measurements, i t i s more appropriate to include the p o s s i b i l i t y of the slender ship assumption being v i o l a t e d by using a model too thick. Another p o s s i b i l i t y f o r the discrepancies can be due to an inadequacy i n the introduced i n t e r p o l a t i o n function ( Equation (24) ) which over predicts the fluid-body i n t e r s e c t i o n p o s i t i o n . This function, which i s l i n e a r i n order, probably requires higher order terms to improve the sol u t i o n . However, from a general point of view, the proposed procedure as well as the free-surface formulations are considered as acceptable f o r the applications i n bow wave modelling. The corrections to l i n e a r numerical solutions by including the non-linear terms i n the free-surface formulations, p r e d i c t the trend of the bow wave amplitude at increasing draft-Froude number more accurately than the an a l y t i c s o l u t i o n provided by O g i l v i e . Moreover, the formation of water spray along the model at high speed found i n the s o l u t i o n i s very encouraging. I t i s believed that with some minor adjustments, t h i s s o l u t i o n approach o r i g i n a t e d by Og i l v i e can become a very e f f i c i e n t modelling technique f o r studying the performance of a ship. 65 IV.4 Wave D i f f r a c t i o n of a C i r c u l a r Cylinder ( 3D ) In t h i s f i n a l example, the three dimensional a p p l i c a b i l i t y of the introduced free-surface formulations i s tested by modelling a wave d i f f r a c t i o n problem. A surface-p i e r c i n g c y l i n d e r i s located i n a wave tank of constant water depth. Waves are generated by assigning s i n u s o i d a l motions to the wave maker. The d i f f r a c t i o n of incident waves by the cyl i n d e r i s cal c u l a t e d i n the time domain and presented i n graphical form. Both l i n e a r and non-linear r e s u l t s are obtained and p l o t t e d at constant time i n t e r v a l s . Forces experienced by the cyl i n d e r are compared to published experimental r e s u l t s by Hogben and Standing (1974) as well as Mogridge and Jamieson (1975). IV.4.1 FormulatIon A basic difference between sol v i n g a two dimensional problem and a three dimensional problem i s purely i n the d e f i n i t i o n of the control surface. In equation (13), which i s the basic form of the Boundary Integral Equation, the con t r o l boundary, S, can no longer be represented by l i n e elements i n the following a p p l i c a t i o n . I t i s probably one of the most p a i n f u l procedures to define S by using surface patches when applying equation (13) i n three dimensions. As one can see i n the following d e r i v a t i o n of the s p e c i f i c form of equation 66 (13) , the procedure i s very s i m i l a r to those i n previous applications except an a d d i t i o n a l mesh generation procedure i s required. I t i s therefore considered to be more appropriate to discuss the d i s c r e t i z a t i o n of three dimensional surfaces under a separate topic heading. In t h i s section, discussions are concentrated on the problem's formulation. An i l l u s t r a t i o n of the problem i s given i n Figure 29 . A surface p i e r c i n g c i r c u l a r c y l i n d e r i s considered numerically i n a small wave tank. The coordinate axis i s located as shown i n the diagram with y pointing up. In order to cut down the number of unknowns involved i n the problem, Green's function ( G ) defined as i n equation (15) i s applied to simulate a mirror image e f f e c t about the bottom of the numerical tank. Moreover, because of symmetry about z = 0 , modelling of ei t h e r h a l f of the tank i s considered adequate. As i n previous studies, the c o n t r o l surface that encloses the considered domain i s divided into S , S , . , and o f . S q i s used to represent impermeable surfaces such as the tank wall, and the cylinder's submerged surface. The free-surface i s designated by , and the wave p i s t o n by S . Therefore, equation (13) can be rewritten as : 67 4>(. P > + <l> G dS + f 4> G ,n J ^ ,n dS + ci G dS m d> G dS + f d> G dS + ,n J r , n ^ G dS ,n (46), m where G i s defined i n equation (15). Since the impermeable surface i s not moving, f l u i d p a r t i c l e s on t h i s surface have a zero normal v e l o c i t y , that i s , * „ = 0 on S (47). o As for the free surface, equation (10) i s considered as appropriate. 2n n c , - i ,,-2 , ,-3 ^ z_ / + _i J l -Ja± + * ' n g ( 5 t ) 2 g ( 6 t ) 2 n + 0 ' n + ^ ' n + - ^ £ - ( c 6 2 + ^ 2 + ^ 2 ) ' r,x x *\z z 2g Dt v r,x ^,y r , z .at y=r? (10) For convenience, a function F 3 i s used to represent the known terms i n the above expression, and i t i s rewritten as : 2n g(fit) <f> + F 2 r 3 at y=r? (48) , where 68 n r - i " i / i " 2 , - 3 F 4 * + * + r n + * ' n 3 g ( g t ) 2 ,X X ,Z Z n + _ z D g Dt^ \ x v , y v , z (49) I t can be r e c a l l e d from Chapter II.2 that the prime notation i s used to i d e n t i f y the terms that are approximated by an extrapolation function i n time. As for the l i n e a r i z e d free-surface condition, 4 — + * at y=0 (50) , g ( S t ) 2 where K i s given i n equation (23) F i n a l l y , f o r completion of the problem, a displacement function should be assigned to the wave maker. Although t h i s can be a r b i t r a r i l y chosen, i t i s very desirable to have a function that s t a r t s the pi s t o n smoothly at time equal to zero. For t h i s reason, the following function i s used. x = X ( 1 - cos tot ) on S (51), m m m where X i s the motion amplitude, and w i s the angular m frequency. The corresponding v e l o c i t y function i s , therefore, given as : 69 v = X co s i n cot X m on S (52) m With equations (47), (48), and (52) substituted into equation (46), the corresponding form of the i n t e g r a l equation i s obtained as follows : 2n <!>(?) + \ <t> G dS+ <f> ( G + G ) dS + | <f> G dS s g ( 6 t ) s o f m = [ F G dS + [ ( v n ) G dS J 3 J x x (53) S f m This equation i s then used to e s t a b l i s h a system of equations and y i e l d solutions at designated time i n t e r v a l s . I f l i n e a r solutions are desired, equation (50) should be used instead of equation (48) upon the s u b s t i t u t i o n , and r e s u l t s are ca l c u l a t e d at the undisturbed free surface. Since there i s nothing s p e c i a l about the l i n e a r i z e d form of the i n t e g r a l equation, i t s d e r i v a t i o n i s not included here. IV.4.2 Mesh Generation Before the system of equations can be constructed, the control surface must be defined. For three dimensional ap p l i c a t i o n s , the boundary surface i s d i s c r e t i z e d and divided into small elements. These small patches with t h e i r s i z e , 70 p o s i t i o n , and o r i e n t a t i o n known, enclose the considered domain. Certain numeric codes are then developed to keep track of the p o t e n t i a l and v e l o c i t y values on p a r t i c u l a r elements during the computations. The choice on the shape of elements i s a r b i t r a r y . Under normal conditions, t r i a n g u l a r patches are the most popular to use. However, i n the case of i n s u f f i c i e n t computer memory, q u a d r i l a t e r a l elements can be used provided that the leakage problem found between elements i s treated c a r e f u l l y . Unfortunately, there i s no general rule or guideline for the choice of what kind of elements should be used. Therefore, i t i s up to the programmer to fi g u r e out the most e f f i c i e n t way to represent the problem numerically. At the U n i v e r s i t y of B r i t i s h Columbia ( UBC ), the best current computer f a c i l i t y a v a i l a b l e i s the Michigan Terminal System ( MTS ). I t i s linked to a FPS-164/MAX array processor located i n the computer center. With the main processor rated at 11 MFLOPs ( M i l l i o n FLOating point Operations per second), and the MAX boards rated at 55 MFLOPs, solutions of large matrices can be c a l c u l a t e d at high speed. However, the si z e of the computed matrix i s r e s t r i c t e d by i t s accessible memory. Its 1 Mwords of main memory with two 135 Mbyte drives, can only handle a system of equations up to approximately 900 unknowns. For t h i s reason, the following problem can only be solved i n the present system with poor re s o l u t i o n . In order to maintain solutions at reasonable 71 r e s o l u t i o n with small enough element numbers, q u a d r i l a t e r a l elements are r e l u c t a n t l y used. A disc u s s i o n on the advantages of using t r i a n g u l a r patches over q u a d r i l a t e r a l elements can be found i n Webster (1975). In h i s paper, he pointed out that there i s a danger of source leakage and d i s c o n t i n u i t y found between elements when a three dimensional surface i s represented by q u a d r i l a t e r a l elements. Since larger memory computers are cur r e n t l y inaccessible, q u a d r i l a t e r a l elements are used under the following assumptions : i . Wave amplitude i s assumed to be small compared to wave length, such that the slope of the free surface i s small i n a l l d i r e c t i o n s , and a l l changes are gradual i n time. i i . In order to ensure that there i s no source leakage between elements, each q u a d r i l a t e r a l element i s considered to be composed of 4 t r i a n g u l a r planes. In Figure 30, the 4 corners of an element are l a b e l l e d as P, Q, R, and S. Its geometric center, c, i s ca l c u l a t e d by : 72 X + X + X + X P Q R S X = -. c 4 y + y + y + y z + z + z + z P Q R S Z c 4 There i s no necessity f o r these 4 corners to l i e on the same plane. In other words, each q u a d r i l a t e r a l element can twist i n space. Its area i s given by the sum of the 4 t r i a n g u l a r facets' area, and i t s normal unit vector i s considered by taking the average of the normal u n i t vectors on these t r i a n g l e s . When computing the sol u t i o n , the respective p o t e n t i a l value and v e l o c i t y on an element are cal c u l a t e d at c . Therefore, source leakage and di s c o n t i n u i t y problems found between elements are considered to be reduced. Provided the small amplitude wave assumption i s not v i o l a t e d , the above assumption for the q u a d r i l a t e r a l element can be accepted. Although these two assumptions have greatly s i m p l i f i e d the programming work, there i s no in t e n t i o n to i n s i s t on the use of q u a d r i l a t e r a l elements. In fa c t , the reason why they are used here i s c l e a r l y explained. When bigger memory computers are accessible, t r i a n g u l a r patches should be used. Therefore, t h i s example serves mainly as a numerical experiment f o r the three dimensional a p p l i c a b i l i t y of the free-surface formulations. A f t e r the q u a d r i l a t e r a l element i s chosen, the dimensions of the considered domain have to be determined. In most of the following computations, the wave tank i s chosen to be 5 times the cylind e r ' s diameter i n length and 2.5 times i n h a l f width. As a r a d i a t i o n condition i s not involved i n the formulations, t h i s chosen tank length i s only good enough for studying the f i r s t two waves without s u f f e r i n g from the disturbance of r e f l e c t e d waves. The simulation time i s ca l c u l a t e d by l i n e a r wave theory f o r the f i r s t wave to t r a v e l back and f o r t h along the tank. Computations should be terminated before the flow near the cy l i n d e r i s disturbed by r e f l e c t e d waves. It i s also important to maintain a f i n e enough siz e of elements on the free-surface i n order to provide good r e s u l t s . In most of the following c a l c u l a t i o n s , the free surface i s represented by 450 elements. A minimum of 12 elements are used to define the submerged surface of the h a l f c y l i n d e r . They are arranged i n layers of s i x , depending on water depth ( see Figure 31 ). The generated mesh of the free surface i s i l l u s t r a t e d i n Figure 32. The wave maker i s on the left-hand side, while on the other end, a r i g i d wall i s placed. 74 IV.4.3. I n t e r p o l a t i o n (3D) The material discussed i n t h i s section serves mainly as an extension of section IV.1.2 . In the three dimensional numerical model, a f t e r the p o t e n t i a l d i s t r i b u t i o n i s known on the free surface, v e l o c i t y components can be resolved since the free - surface's gradient i s known. These v e l o c i t y components are then used to pr e d i c t the free-surface motion within the next time increment, so as to advance the s o l u t i o n i n time. The motion of an element i s ca l c u l a t e d by considering the resultant motions at i t s 4 corners. Its geometric center, c, i s always given by equation (54). In order to define the v e l o c i t y at an element's corner, a simple i n t e r p o l a t i o n function ( equation (55) ) i s used. In Figure 33a, point P i s considered as the common point surrounded by element I , J , K, and L . Since the v e l o c i t y i s known at the center of these element, the following function i s used f o r the in t e r p o l a t i o n s of v e l o c i t y at P. ax + by + cz + d (55), where a, b, c, and d are constants d i f f e r e n t from element to element. 75 When t h i s equation i s applied at the centers of 4 adjacent elements I , J , K, and L , 4 equations with 4 unknowns can be written. The v e l o c i t y component,for example v , on the four elements contributes to the column vector on y the right-hand-side of the matrix equation, and the constants ( a, b, c, and d ) are the unknowns to be solved. When these unknowns are found, v at point P can be calculated. The y procedure i s very s i m i l a r to the discussion i n Appendix I I . It i s found that equation (55) can be applied to a point which i s located at the fluid-body i n t e r s e c t i o n boundary. In Figure 33b, element I and J are on the free- surface, while K and L are on the cyl i n d e r w a l l . The normal v e l o c i t y component on the cylinder wall i s taken as zero to prevent f l u i d p a r t i c l e s from penetrating into the surface. Other components are computed as discussed i n Appendix I I . Although there was a h e s i t a t i o n i n applying equation (55) to the fluid-body i n t e r s e c t i o n boundary as i t s l i m i t a t i o n i s not quite c l e a r , reasonable r e s u l t s are obtained i n most of the studies. I t only happens to f a i l i n modelling a high amplitude wave case, where the wave run-up on the cy l i n d e r i s improperly represented by coarse elements. This computer deficiency s i t u a t i o n can be improved by using more free-surface elements near the cy l i n d e r . I t i s , of course, possible to use a second order function f o r better 76 i n t e r p o l a t i o n s , except such an idea i s not very p r a c t i c a l f o r the following studies. Computer cost i s a prime f a c t o r to be considered upon the acceptance of a second order i n t e r p o l a t i o n function. IV.4.4 Results and Discussion In the following r e s u l t s , the diameter of the v e r t i c a l c y l i n d e r i s chosen to be 50.00 m. The wave tank i s determined to be 250.00 m i n length, and i t s h a l f width i s taken as 125.00 m. 450 elements are used at the free surface. Water depth i s at 10.00 m, and 12 elements, arranged i n two layers, are used to represent the cy l i n d e r . Motion amplitude ( X ) of the wave maker i s assigned at 1.00 m. Results are obtained at two frequencies. The f i r s t set of solutions i s c a l c u l a t e d by moving the wave maker at 0.061 Hz. The corresponding ka value i s equal to 1.0 , where k i s the wave number and a i s the radius of the v e r t i c a l c y l i n d e r . Computations proceeded i n time at 0.2 second r e a l time increments f o r 36.00 seconds. Results are pl o t t e d at every 4.00 second i n t e r v a l . In Figure 34, both l i n e a r and non-linear solutions are i l l u s t r a t e d i n t h e i r graphical form. Results obtained by a l i n e a r free-surface condition are on the left-hand-side of the diagram. The v e r t i c a l p l o t t i n g scale i s enlarged 50 times f o r better 77 i l l u s t r a t i o n . At t=0, the pi s t o n located at the l e f t end of the tank s t a r t s i t s motion. Waves are generated and propagate down the tank at constant speed. Since the simulated wave amplitude i s small, the difference between the two solutions i s not very obvious. The second wave crest i s ca l c u l a t e d at 0.20 m by the l i n e a r formulation, while i t i s at 0.21 m from the non-linear s o l u t i o n . Wave length i s observed at about 157 m which corresponds to a height over length r a t i o ( H/A ) of 0.00255 (l i n e a r ) and 0.00268 (non-linear) r e s p e c t i v e l y . The non-linear waves are observed to t r a v e l s l i g h t l y f a s t e r than the l i n e a r waves. At 36.00 seconds, the non-linear s o l u t i o n exhibits a crest high enough to d i f f e r from the l i n e a r r e s u l t . I t i s believed that the solutions at t h i s moment are disturbed by r e f l e c t e d waves from the end wall. Although at t=36 seconds, these r e f l e c t e d waves are almost i n s i g n i f i c a n t , the force c a l c u l a t e d on the c y l i n d e r i s thought to be influenced and become u n r e l i a b l e . Figure 35 i s a h o r i z o n t a l p r o j e c t i o n of the free-surface elements d i s t r i b u t i o n at t=0 and t=36 seconds. B a s i c a l l y , there i s no major difference between the two, except at the near cylinder region, elements seem to s u f f e r from very great d i s t o r t i o n s . I t i s believed that, at the forward and a f t e r 78 stagnation points of the c y l i n d e r , more elements should be used to define the wave run-up s i t u a t i o n . Moreover, i t i s desirable to model the tank length as long as possible, so as to delay the disturbance from r e f l e c t i o n u n t i l l a t e r cycles. In order to i l l u s t r a t e how the r e s u l t s are a f f e c t e d by the chosen tank length, a shorter length i s used to repeat the computations. In Figure 36, a tank length of three times the c y l i n d e r diameter i s considered ( i e . 150.00 m ). Other parameters remain the same as i n Figure 34 . The free surface now consists of 270 elements ( see Figure 37 ), and the simulation continues for 32.00 seconds. Results are b a s i c a l l y s i m i l a r to Figure 34 , except that the generated waves are s l i g h t l y higher than i n the previous case. However, when the free-surface elements i n Figure 37 are examined at time t=32, i t shows that some of these elements have p a r t i a l l y penetrated into the c y l i n d e r . This observation indicates that the tank length has a very strong influence on the solutions, and probably there e x i s t s a numerical weakness i n the constructed model. I f the chosen tank dimension i s too short, numerical solutions w i l l be exaggerated and consequently a f a i l u r e within a very short period of time. This i s observed by comparing the c a l c u l a t e d wave amplitude i n Figure 34 and Figure 36. As the same p l o t t i n g scale i s used, the wave amplitude obtained i n the shorter tank i s obviously higher i n 79 magnitude.Therefore, a tank with a length 5 times the cyli n d e r diameter i s taken to be the minimum acceptable l i m i t i n t h i s study. The motion frequency of the wave maker i s doubled to 0.122 Hz i n Figure 38 i n order to model a shorter wave length. In t h i s case, the v e r t i c a l dimension i s scaled up only 25 times f o r p l o t t i n g . Wave length i s 1.5 times the cyl i n d e r diameter. Computation i s designed to terminate at 36.00 seconds. Linear solutions are shown on the left-hand side and non-linear solutions are shown on the right-hand side of Figure 38. No major difference i s observed between the two. The non-linear wave cr e s t i s s l i g h t l y higher and tra v e l s s l i g h t l y f a s t e r . In order to te s t the l i m i t a t i o n s on the constructed model, the piston's motion amplitude i s increased to 2.5 m to simulate a high amplitude wave s i t u a t i o n . The frequency remains at 0.061 Hz. With a v e r t i c a l p l o t t i n g scale 25 times the c a l c u l a t e d dimension, r e s u l t s are shown i n Figure 39. The f i r s t wave, which i s highly transient i n nature, has an amplitude of about 1.00 m before i t reaches the c y l i n d e r . Although the H/A i s equal to 0.0127 ( which i s considered as acceptably small ), computations of the non-linear solutions f a i l e d at 16.00 seconds. This i s explained by an improper modelling of the wave run-up s i t u a t i o n on the c y l i n d e r with 80 the coarse si z e elements. In Figure 39, the free surface exhibits a depression at the forward stagnation point of the cyl i n d e r . Similar f a i l u r e i s not found i n the l i n e a r s o l u t i o n . I t i s believed that the f a i l u r e i s caused by an improper modelling of the wave run-up at the c y l i n d e r by q u a d r i l a t e r a l elements. A f i n e r g r i d i s probably required to represent t h i s region where flow changes r a p i d l y . Moreover, as the H/A r a t i o i s higher, the q u a d r i l a t e r a l element w i l l be twisted badly at the near c y l i n d e r region. However, i r r e s p e c t i v e of t h i s imperfection, the constructed numerical model behaves very well at low amplitude wave simulations. There i s enough evidence to support that the introduced free-surface formulation and so l u t i o n procedure are s u f f i c i e n t l y accurate f o r modelling transient, non-linear, three dimensional wave making problems i n the time domain. F i n a l l y , the wave force exerted on the c y l i n d e r i s computed and p l o t t e d i n Figure 40. The continuous l i n e appears i n the diagram i s the frequency domain s o l u t i o n based on Chan (1984). The non-linear time domain c a l c u l a t i o n s are compared to measurements made by Hogben and Standing (1974) as well as Mogridge and Jamieson (1975). Since the non-linear solutions are ca l c u l a t e d i n a time domain wave tank, forces due to the f i r s t wave cycle are ignored. The presented 81 r e s u l t s are considered to be the maximum force c a l c u l a t e d a f t e r the f i r s t wave has passed by the c y l i n d e r . The non-dimensional form was suggested by Sarpkaya and Isaacson (1981) as : max pgHda [ tanh kd / kd ] and J P dS (57) i s the submerged surface of the c y l i n d e r . P i s the pressure c a l c u l a t e d i n equation (4). p i s the f l u i d density, H i s the incident wave height, d i s the water depth, and a i s the radius of the c y l i n d e r . F c a l c u l a t e d at d i f f e r e n t ka values i s shown i n max Figure 40. At a ka value greater than 1.00, the agreement between the time domain solutions ( non-linear ) and the experimental r e s u l t s i s excellent. Not only the forces are accurately predicted, and the trend of the non-linear s o l u t i o n i s observed very s i m i l a r to the measurements by Hogben (1974). As for the frequency domain s o l u t i o n , which i s l i n e a r , the p r e d i c t i o n i s not as good as the non-linear s o l u t i o n . In f a c t , the trend of the experimental r e s u l t s cannot be observed i n the l i n e a r s o l u t i o n when ka i s greater 82 than 1.00. However, at a ka value less than 1.00, the experimental r e s u l t s follow the frequency domain s o l u t i o n while the time domain s o l u t i o n diverge from the measured values. The smaller the ka value i s , the bigger the discrepancy i s found. Since the error i s behaving systematically, there must be an explanation f o r i t . These discrepancies are believed to be associated with the chosen tank dimension. At ka equal to 1.00, the wave length i s ca l c u l a t e d to be 157 m, which i s about 0.628 times the chosen tank length. Below t h i s value, the wave length becomes very long and eventually exceeds the tank's dimension. When t h i s happens, since F i s the force max cal c u l a t e d a f t e r the f i r s t wave cycle, i t i s very possible that the f i r s t wave being r e f l e c t e d from the end wall i n t e r f e r e s with the second wave before i t reaches the cyl i n d e r . Although, the time domain p r e d i c t i o n i s very poor at low ka values, the f a c t that t h i s r e s u l t ( which i s non-linear ) at ka values greater than 1.0 i s much better than the l i n e a r s o l u t i o n that t h i s argument i s strongly supported. As for the frequency domain s o l u t i o n , the same problem does not e x i s t since a r a d i a t i o n boundary i s used ( see Chan 1984 ). I t only happens at ka value greater than 1.00 that the l i n e a r s o l u t i o n becomes not as r e l i a b l e as the 83 non-linear so l u t i o n . I f the modelled tank length can be increased to about 10 times the cylinder's diameter, the time domain s o l u t i o n of ^max ° a n P o s s i b l y be c a l c u l a t e d more p r e c i s e l y without being disturbed by r e f l e c t e d waves. From the experience gained i n t h i s study, the free-surface representation i s very c r i t i c a l to the numerical model's s t a b i l i t y . More elements should be used at the near cyl i n d e r region. A longer tank i s thought to be h e l p f u l to eliminate the interference from r e f l e c t e d waves i n the time doamin so l u t i o n . Moreover, i f i t i s possible, t r i a n g u l a r patches should be used to obtain a much better d e f i n i t i o n of the wave run-up s i t u a t i o n on the cy l i n d e r . Results obtained i n t h i s s ection are regarded as encouraging, and the three dimensional a p p l i c a b i l i t y of the proposed free-surface formulations i s shown with some r e s t r i c t i o n s . These r e s t r i c t i o n s seem to depend on the free-surface element d e f i n i t i o n rather than the methodology. 84 V. DISCUSSIONS AND CONCLUSIONS In t h i s chapter, discussions are made from a general point of view. Although, some of the discussions seem to overlap the ideas from previous chapters, conclusions cannot be established without an overview of a l l the accomplishments. The following discussions are mainly on three kinds of topics. They are on the formulation of problems, numerical r e s u l t s , and recommendations. Since there i s c e r t a i n unavoidable overlapping between the ideas, they are discussed simultaneously i n the following section. V.1 General Discussions The proposed s o l u t i o n procedure i s , i n f a c t , p a r a l l e l to Longuet-Higgins' approach, except that the free-surface boundary conditions are expressed i n r e a l v a r i a b l e s only. The time or material d e r i v a t i v e term i s replaced by a f i n i t e difference form based on known values at previous time steps. A predictor method i s used to extrapolate future values f o r v e l o c i t y terms and used to ca l c u l a t e the free-surface boundary condition to advance the computations i n time. As a consequence, a very systematic procedure i s established to solve gravity wave problems i n the time domain. 85 The f i r s t two problems presented i n chapter IV a c t u a l l y serve as t e s t cases for the study of the non-linear free-surface boundary conditions. The free-surface condition i s emphasized here because i t has remained as the source of problems for many years. As mentioned i n the introduction, t h i s moving boundary causes d i f f i c u l t i e s , not because of the uncertainty involved i n the free-surface d e f i n i t i o n , but because of the numerical i n s t a b i l i t y found i n the time domain so l u t i o n . One of the contributions achieved by Longuet-Higgins was h i s i ntroduction of a smoothing procedure during computations to i n h i b i t the growth of i n s t a b i l i t y i n the time stepping s o l u t i o n . From then on, many transient free-surface problems were modelled numerically. I t must be admitted by the author that the rearrangement of the free-surface boundary conditions presented i n t h i s t hesis i s f i n a l i z e d through a series of numerical t r i a l and error studies. The idea of extrapolating the unknown v e l o c i t y terms and maintaining them i n an unexpanded form are believed to be the c r u c i a l assumption made i n these studies. The wave tank study and the modelling of breaking waves enabled the author to study the d i f f e r e n t algebraic forms of the free-surface boundary conditions to r e s u l t i n a numerically stable s o l u t i o n without any smoothing requirement. 86 The choice of t h i s unexpanded form given i n equation (10) i s based purely on numerical t e s t experience gained by the author. Other possible arrangements, a l l based on various expansions of the non-linear terms i n equation (10), have been t r i e d without any success. One possible explanation f o r t h i s chosen form to y i e l d a numerically stable s o l u t i o n i s probably r e l a t e d to the nature of the B e r n o u l l i Equation. In equation (4) , the v e l o c i t y square terms can be considered associated with the k i n e t i c energy of the f l u i d p a r t i c l e . Its material d e r i v a t i v e ( as i n equation (6) ) i s thought to be re l a t e d to the change of k i n e t i c energy of the same p a r t i c l e being traced i n time. I t i s possible that f l u i d p a r t i c l e s ' k i n e t i c energy, which i s a scalar quantity, changes gradually i n time and can be approximated momentarily by a second order function of time. As f o r an expanded form of the material d e r i v a t i v e of the v e l o c i t y square terms, numerical error t o l e r a t e d i n the approximation of the change of k i n e t i c energy could be exagerated by increasing the number of unknown terms i n the representation. Instead of t r e a t i n g the associated k i n e t i c energy terms as a single quantity, each v e l o c i t y component as well as i t s material derivative have to be approximated by a time extrapolation scheme. As a r e s u l t , the error involved i n the free-surface formulation could become i n t o l e r a b l e , and 87 the s o l u t i o n could become unstable. I t i s also remarkable how the i n t e r p o l a t i o n function works at the fluid-body i n t e r s e c t i o n point. As studied by L i n (1984), t h i s point requires s p e c i a l treatment. I t can be proved that, when t h i s i n t e r s e c t i o n point i s extrapolated by values at the free-surface only, i t behaves reasonably for low amplitude wave cases. The idea of including information on the wave p i s t o n to interpolate the v e l o c i t y at t h i s i n t e r s e c t i o n point i s believed to be an acceptable approach. In f a c t , the wave splashing phenomenon obtained i n the bow wave simulation at high model speed has proven to a c e r t a i n extent that the idea i s correct. However i t i s also a disadvantage that t h i s i n t e r p o l a t i o n s k i l l only works on the zeroth order elements. For higher order boundary elements, since end points are node points and defined the i n t e r s e c t i o n point, another i n t e r p o l a t i o n c r i t e r i o n must be developed. The d r i f t of free-surface elements, though bringing some excitement since the natural phenomenon i s modelled, i s a bad i n d i c a t i o n on the s t a b i l i t y of the numerical model. I f a l l the elements at the free surface are d r i f t i n g down the tank, i t i s only a matter of time for the numerical model to f a i l because of the improper si z e d elements representing the free surface. Although, there are many remedies to t h i s problem, such as breaking an element which i s too long into two, i t i s 88 considered to be more appropriate to leave t h i s study to the future. Fortunately, i t takes many wave cycles to worsen the problem, and i t can be ignored for the f i r s t few cycles. From the studied problems, c e r t a i n l i m i t a t i o n s are found. T h e o r e t i c a l l y , better solutions can always be obtained by using a f i n e r g r i d , increasing the tank length, adopting t r i a n g u l a r patch elements, or even reducing the si z e of the time step. Most of these recommendations require more computer memory. Since i n r e a l l i f e , there i s always a hardware l i m i t a t i o n on computers, i t i s very unwise to r e l y on expanding computer memory for better solutions. Some actions should be taken to reduce the chances of having a bad s o l u t i o n . For time stepping s o l u t i o n problems of t h i s kind, i t would be very h e l p f u l to formulate some sort of beach condition. As i t was shown i n the three dimensional modelling problem, a great obstruction found i n programming the problem i s the l i m i t e d accessible computer memory. I f such a boundary condition can be formulated for time domain problems, at l e a s t the s o l u t i o n can be free from the disturbance of r e f l e c t e d waves. However, at t h i s present moment, such a boundary condition i s not ready, so upon constructing a numerical model, s p e c i a l care must be taken on the d e c i s i o n of the dimensions of the control domain, the s i z e of elements, as well as the problem's computer simulation time. 89 P r i o r to considering the conclusions, the e f f i c i e n c y of the time stepping s o l u t i o n procedure i s worth some discussion. Although solutions are obtained f o r various problems, the computer cost i s very high. In two dimensional problems, solutions can be obtained at quite a low cost since the number of unknowns i s usu a l l y small. However, f or three dimensional problems, i t i s extremely easy to come up with a numerical model with a few thousand unknowns. Access to very huge and super computers seems to become a necessity. It is very unfortunate that the established matrix i s unsuitable for i t e r a t i o n s . Since t h i s huge matrix must be solved by a d i r e c t method ( such as the Gaussian Elimination Method ), a major p o r t i o n of the computer time i s consumed i n solv i n g the matrix. Moreover, s t a r t i n g the simulations at the zero i n i t i a l stage consumes a large amount of computer time. I t would be id e a l i f the i n i t i a l conditions were given at a time j u s t p r i o r to the moment of in t e r e s t , and the transient behavior were studied from that point on. I f t h i s can be done, free-surface problems can be studied very e f f i c i e n t l y . However, at t h i s present moment, such an idea i s not very p r a c t i c a l . For example, many of the steady state solutions are c a l c u l a t e d f o r an i n f i n i t e domain or with a r a d i a t i o n condition. In order to apply these a n a l y t i c solutions as i n i t i a l conditions, the boundary conditions between the time 90 domain model and the frequency domain model must be compatible within some degree. V.2 Conclusions The presented numerical procedure i s considered as a r e l i a b l e approach f o r solving p o t e n t i a l flow problems with a free surface. The time domain solutions obtained f o r two and three dimensional problems are numerically stable and do not require any smoothing procedure. The free-surface boundary condition i s rearranged to y i e l d non-linear solutions under a small time step assumption. The fluid-body i n t e r s e c t i o n s i n g u l a r i t y problem i s handled by an i n t e r p o l a t i o n function that requires information from both surfaces. Numerical r e s u l t s are comparable to l i n e a r theory as well as published experimental work. Limitations are found i n modelling a three dimensional d i f f r a c t i o n problem. Solutions f o r high amplitude waves are r e s t r i c t e d by the use of q u a d r i l a t e r a l elements and coarse g r i d s i z e . Wave forces c a l c u l a t e d on the surface-p i e r c i n g c y l i n d e r are influenced by the chosen tank length. As the simulated wave length becomes comparable to the tank length, the cal c u l a t e d force experienced by the cyli n d e r diverges from the experimental r e s u l t . When the wave length i s short compared to the tank, good r e s u l t s are obtained. 91 REFERENCES Banerjee, P.K. and B u t t e r f i e l d R. 1981. Boundary Element Methods in Engineering Science, London : McGraw H i l l . Brebbia, CA. 1978. The Boundary Element Method for Engineers, London : Pentech Press. C a l i s a l , S.M. and Chan, J.L.K. 1987. Breaking Waves Simulation. Second International Workshop on Water Waves and Floating Bodies, Dept. of Math., U. of B r i s t o l , U.K., Report no. AM-87-06, pp. 9-12. C a l i s a l , S.M. and Chan, J.L.K. 1986. S t a b i l i t y of Low L/B Vessels. Unpublished Report, Dept. of Mech. Engineering, U. of B r i t i s h Columbia, Vancouver, B.C. Chan, J.L.K. 1984. Hydrodynamic C o e f f i c i e n t s of Axisymmetric Bodies at F i n i t e Water Depth. Master Thesis, Dept. of Mech. Engineering, U. of B r i t i s h Columbia, Vancouver, B.C. Dean, R.G. and Dalrymple, R.A. 1984. Water Wave Mechanics for Engineers and S c i e n t i s t s , New Jersey : Prentice H a l l . Fredholm, I. 1903. Sur une class d'equations f o n c t i o n e l l e s , Acta. Math., Vol 27, pp 365-390. Hogben, N. and Stanging, R.G. 1974. Wave Loads on Large Bodies. Proc. Int. Symp. on the Dynamics of Marine Vehicles and Structures in Waves, Univ. College, London, pp 258-277. Hunt, B. 1980. The Mathematical Basis and Numerical P r i n c i p l e s of the Boundary Integral Method for Incompressible P o t n t i a l Flow over 3-D Aerodynamic Configurations. Numerical 92 Methods in Applied Fluid Dynamics, London : Academic Press, pp. 49-136. Jawson, M.A. 1963. Integral Equation Method in Potential Theory, I. P r o c , R. Soc. Lond., Series A 275, pp 23-32. Lin, W. 1984. Non-linear Motion of the Free Surface Near a Moving Body. Ph D Thesis, Dept. of Ocean Engineering, MIT, Cambridge, Mass. Longuet-Higgins, M.S. and Cokelet, E.D. 1976. The Deformation of Steep Surface Waves on Water, I. A Numerical Method of Computation. Proc. R. Soc. Lond., Series A, 350, pp. 1-26. Mogridge, G.R. and Jamieson, W.W. 1975. Wave Forces on C i r c u l a r Caissons : Theory and Experiment. Can. J. Civil Eng., Vol 2, pp 540-548. Newman, J.N. 1977. Marine Hydrodynamics. Cambridge, Mass., MIT Press. O g i l v i e , T.F. 1972. The Wave Generated by a Fine Ship Bow. Dept. of Naval Architecture and Marine Engineering, U. of Michigan, Report no. 127. Sarpkaya, T. and Isaacson, M. 1981. Mechanics of Wave Forces on Offshore Structures. New York : Van Nostrand Reinhold. Vinje, T. and Brevig P. 1980. Non-linear, Two Dimensional Ship Motion. SIS Report, Norwegian Hydrodynamic Laboratory Trondheim. Webster, W.C. 1975. The Flow About A r b i t r a r y , Three Dimensional Smooth Bodies. Journal of Ship Research, Vol. 19, 93 no. 4, pp 206-218. 94 Boundary Condition Figure 1. A Typical Wave Making Problem Figure 2. Numerical Model with Mirror Image at y=-d 95 F i g u r e 3. L i n e a r and N o n - l i n e a r Wave Tank M o d e l s F i g u r e 4. I n t e r p o l a t i o n o f V e l o c i t y Between E l e m e n t s 96 2 . 4 0 e e o I I I h H 1 H t. = 2 . 6 0 ©oo 2. 00 1. 50 1. 00 0. 50 0. 00 -0. 50 -1. 00 •1. 50 -2. 00 0. •t « 2 . 8 0 ©eo ' I I I — J 1 I- I 1 1 H ' • ' ' 3 . 0 0 e e o - t - i I 1 \-00 e. 00 WAVE-MAKER. 16. 00 24. 00 32. 00 40. 00 TANK DIMENSIONS . . 4 0 . X 4 . Wov« Mok*r> Fn*C| . . 0 . 5 Hz Wav* Maken Amp . . . 0 . 3 Figure 5. Wave Generation by a Numerical Wave Tank Model ( dt = 0.05 sec. ) 97 E d = 4 m o d = 6 m 0.20 Frequency ( H z ) Figure 6. Prediction of Wave Length by Linear Numerical Model 10.80,— B.00 8. ZZ 7 . 0 0 _ 6.00 S . 0 0 _ 4.00 a. 0 0 _ d=2m H d = 4 m o d = 6 m B.Bet i i i i I i i i i I i i i i I i i i i I i i i i I i i i i I i i i i I 0.00 0.20 0.40 0.60 0.00 1.00 1.20 1.40 Frequency ( Hz ) Figure 7. Comparison of Phase Velocity Between Solution and Wave Theory A. 00 i-3. 50 - • 3. 00 -E 2. 50 -Shallow Water / 3 <C Aw/Am = kd / - Numerical / Solution \ itude 2. 00 Ampl 1. 50 -Wave 1. 00 -/ ca .-' / Deep Water / y ' Aw/Am = A:d/2 0. 50 / t a - ' 0. 00 i i i i 1 i i i i 1 i i t i 1 i i i i 1 i i i i 1 0. 0 0 0.20 0.40 0.60 0.80 1.00 Hotion Amplitude Am ( m ) Figure 8. Wave Amplitude Predictions by Numerical Solution Figure 9. Wave Profile's Smoothness vs Element Number o o t - 1 8 . 0 0 omo t - 1 8 . 5 0 omo t - 1 9 . 0a mmo - 1 9 . 5 0 mmc 48.(90 04.00 80. 00 TANK DIMENSIONS .. BO. X 4. Vav. Molw Fr^q .. 0.25 Hz Vav* Mah«r Anp ... 0. 10 t - 1 8 . 0 0 e » o - 1 8 . 5 0 w o fc - 1 9 . 0 0 o o o 0. 40 0. 20. -0. 20. -0. 40. 0.00 t - 1 9 . 5 0 e * o 18.00 n—1inoa 32.00 64. 00 r~ w a v « 48. 00 TANK DIMENSIONS .. BB. X 4. Wav. H o l w Fr-«q . . B.ZS H> Wav. N o l w Aap ... 0. 10 80. 00 Figure 10. Comparison Between Linear and Non-linear Wave Tank Results ( dt = 0.05 sec. ) 4. 000 2. 000 /~B. 000 - 2 . 000 _ - 2 . 000 t = 5 sec. t = 6 sec - 4 . i i i i l i i i i I i i i i I i i i i I i i i i I i i i i_ 0 . 0 0 0 2. 000 4. 000 x ( m ) 6. 000 e. 000 10 .000 Figure 11. Non-linear Wave Modelling by Numerical Wave Tank 4. 000 r-2. 000 -/ N e w 0 . 0 0 0 § 0 3 5 ) t=0 t=8 s - 2 . 000 -- 4 . 000 1 1 1 1 1 i i l i i i i I i • • 1 i i • i I i i , . I - 2 . 000 0. 000 2. 000 4. 000 x ( m ) 6. 000 e. 000 10 .000 Figure 12. Motions of Free-Surface's P a r t i c l e s i n Time ( 0 to 8 sec. ) 101 Ay element N+l-i 4-c ' M 1 t; Hr—I 1 t-y = o LINEAR element i Mirror Image Line Figure 13. Single Wave Length Model for Spatially Periodic Wave 102 1.0 r 0.5 _ 0.0 -0.5 -1.0 19.5 t = 20.0 fit = 0.05 sec. H/x H/d 0.10 0.25 "I. o I—1—1—1—1—J—<—•—' i I i i i i-1.0 3.0 5.0 x ( ni ) 7.0 9.0 Figure 14. Simulation of Free-Surface Wave ( Linear Modelling ) 11.0 4.0 2.0 -2.0 U -4.0 I . . . . i • H/X = 0.10 H/d = 1.00 fit = 0.01 sec. = 1.25 sec - J — i — i — i — i — i i i • • t • 0.00 2.00 4.00 6.00 8.00 x ( m ) 10.00 J 12.00 Figure 15. Deformation of High Amplitude Wave ( Non-linear Modelling ) 4.0 2.0 -2.0 HA = 0.125 H/d = 1.25 6t = 0.01 sec. "4 • 0 I I L __ l I J 1 1 1 I I 1 • • ' I • ' • ' I ' • i ' I • i i i I 0.0 2.0 4.0 6.0 8.0 10.0 11.0 x ( m ) Figure 16. Deformation of High Amplitude Wave ( Non-linear Modelling ) 4.0 r -2.0 --2.0 --4.0 H/x = 0.15 H/d = 1.50 6t = 0.01 sec. —1 1 I I I I I I I I I I • I i • 1 ' I • • i ' I ' I I L 0.0 2.0 4.0 6.0 x ( m ) 8.0 10.0 12.0 Figure 17. Deformation of High Amplitude Wave ( Non-linear Modelling ) 104 Figure 18. Waves Generated by a Wedge-shaped Model mirror image line Figure 19. Sectional View of the Considered Domain a.saa r-U = 8 ft/sec a.4M T = 8" 4-> 14-•.aaa a = 7.5° B •.2M / N . non-linear solution •. in linear s o l u t i o n " ^ — - . a. ne . . 1 x ( f t ) 5. am Figure 20. Bow Wave Prediction Calculated Along Wedge-shaped Model 105 Q) » 3 cr Z I 0 o LINEAR A NON-LINEAR • EXPERIMENT • A T = 4" A A • • \ A • * D • O O O A • O O 8 12 U ( f t / s e c ) 16 20 cn CD -C c 6 1 * 0 o LINEAR A NQN-LINEflR • EXPERIMENT - A T = 12" A A D • A D • A • • O o n 6 ° ° d n 8 12 U ( f t / s e c ) 16 20 10 o LINEAR A NON-LINEAR • EXPERIMENT • A T = 8" A a A • A • o o A o • o A D ° . 1 8 12 U ( f t / s e c ) 16 20 ra -C o c 8 & 4 2 0 o LINEAR A NON-LINEAR a EXPERIMENT A T = 16" A A a • A A • • O A D O O 8 D 6 D • 8 12 U ( f t / s e c ) 16 20 Figure 21. Bow Wave Predictions at Different Draft Calculated Along a 15° Wedge ( a = 7.5° ) > o LINEAR A NON-LINEAR • EXPERIMENT T = 4" A A • D O o o o o o o o o 2 3 F x I 2 >-o LINEAR A NON-LINEAR • EXPERIMENT A A T = 12" ^ A A • A _ O * A p ° ° ° o o—D—e e o o 0.0 0.5 1.0 1.5 F 2.0 2.5 3.0 x g 2 I • o LINEAR T = 8" A NON-LINERR A d EXPERIMENT A A A • • • • • " O O O O O O D O -1 2 3 4 5 F x g 2 o LINEAR A NON-LINERR o EXPERIMENT A T = 16" • O -o o—a—o e o o-0 0.5 1.0 1.5 F 2.0 2.5 3.0 Figure 22. Non-dimensional Bow Wave Amplitude Comparisons at Different Drafts ( a = 7.5° ) 2 • A NON-LINEAR D EXPERIMENT • A A A n o ° ° a A O Q A • • • T = 4" • 1 2 F 3 X a E >- 2 • A NON-LINEflR a EXPERIMENT A a a A • D T = 12" • o 0.5 1.5 F 2.0 2.5 3.0 5 A NON-LINEflR a EXPERIMENT 4 3 X 2 >-A A 2 A • D T = 16" • 1 n a 8 0 0.5 1.0 1.5 2.0 2.5 3.0 F Figure 23. Non-dimensional Bow Wave Amplitude Comparisons at Different Drafts ( a = 15 ) (0 0) sz a c X g X 35 30 25 20 15 10 5 0 o LINEAR A NON-LINEAR • EXPERIMENTAL • o A O T = 4" O A • O A • • • o A o A • o A A • • • o a » 8 12 U (ft/sec) 16 20 35 X co ID "§ 201-15 x x 10 5 0 o LINEAR o A A NON-LINEAR T = 12" o A • • • EXPERIMENT A 8 12 U (ft/sec) 16 20 Figure 24. Longitudinal Position ( Xmax (0 0) sz o c X i X o LINEAR A NON-LINEAR • EXPERIMENT A O T = 8" A O A O A D O A a a A O A O a 2 G • 8 12 U (ft/sec) 16 20 35 X (0 at u 20 c 151-x I 10 5 0 o LINEAR o A O A A NON-LINEAR T = 16" o A a A • a EXPERIMENT A 8 12 U (ft/sec) 16 20 0 ) of Wave Peak ( a = 15 ) Linear Solution Non-linear Solution Figure 25. Bow Wave Modelling at Low Model Speed ( T = 12" , a = 7.5°, U = 5 ft/sec. ) 110 Linear Solution Non-linear Solution Figure 26. Bow Wave Modelling at Intermediate Model Speed ( T = 8" , a = 7.5° , U = 8 ft/sec. ) Non-linear Solution Figure 28. Formation of Water Spray at High Model Speed o ( T = 12" , a = 7.5 , U = 14 ft/sec. ) Figure 29. Numerical Modelling of Wave Tank with a Surface-Piercing Cylinder ( 3-D ) Figure 30. Quadrilateral Element Figure 31. Half-cylinder Modelled by Quadrilateral Panels 113 Figure 32. Free-Surface of Numerical Tank Represented by Quadrilateral Elements at Time t=0. a) at Free-Surface b) at Cylinder Wall and Free-Surface Intersection Figure 33. Interpolation of Velocity Components Between Elements ( 3-D ) 114 LINEAR SOLUTION NON-LINEAR SOLUTION Figure 34.a Numerical Simulations of Diffracted Wave in a Tank at Frequency of 0.06105 Hz 115 LINEAR SOLUTION NON-LINEAR SOLUTION Figure 34.b Numerical Simulations of Diffracted Wave in a Tank at Frequency of 0.06105 Hz LINEAR SOLUTION NON-LINEAR SOLUTION Figure 34.c Numerical Simulations of Di f f racted Wave in a Tank at Frequency of 0.06105 Hz 117 Figure 35. Deformation of Free-Surface Elements at t = 36 sec. 118 t = 8 sec t = 12 sec NON-LINEAR SOLUTION Figure 36.a Diffracted Waves are Studied in a Short Tank ( 0.06105 Hz ) 119 t = 24 sec t = 32 sec NON-LINEAR SOLUTION Figure 36.b Diffracted Waves are Studied in a Short Tank ( 0.06105 Hz ) 120 t = 0 sec Figure 37. Deformation of Free-Surface Elements Affected by Shorter Tank Length at Time t = 32 sec. 121 LINEAR SOLUTION NON-LINEAR SOLUTION Figure 38.a Simulations of Diffracted Waves at Higher Frequency ( 0.1221 Hz ) 122 LINEAR SOLUTION NON-LINEAR SOLUTION Figure 38.b Simulations of Diffracted Waves at Higher Frequency ( 0.1221 Hz ) 123 t = 28 sec t = 36 sec LINEAR SOLUTION NON-LINEAR SOLUTION Figure 38.c Simulations of Diffracted Waves at Higher Frequency ( 0.1221 Hz ) 124 I 4 sec 125 Figure 39. Failure of Non-linear Model at High Wave Amplitude 3.0 o NUM. SOLN. A MOGRIDGE 1975 • HOGBEN 1974 2.0 -1.0 0.0 0.0 0.5 .0 k a 1.5 2.0 Figure 40. Non-linear Forces Calculated at the Cylinder at Various ka values APPENDIX I CONSTRUCTION OF A SYSTEM OF LINEAR EQUATIONS FROM THE BOUNDARY INTEGRAL EQUATION. From equation (13) : 0(P> + I <£(Q) G dS = f ^ (Q) G dS (13). s s I f the control boundary, S, i s divided into N elements, and 4> i s assumed constant on each of these elements, equation (13) can be rewritten as : N S 4 + ) <f>. G SS =/~<t>< G 6S. where N 0 = 1 Subscript i and j are element numbers from 1 to N . Since there are N elements, N equations of N unknowns 127 can be written. As a consequence, the following matrix can be obtained. [ A ] <j> = [ c ] <f> n where <h and d> are column vectors. The associated entries of matrix A and C can be computed as : A = I + i j i j + S j G dS i j , n S j + S . C = f G dS i j J i J S j where and SS = 2s j j I = 1 when i=j i j J I = 0 when i * j i j I i s the i d e n t i t y matrix. The above i n t e g r a l terms f o r A and C can be 6 i j i j c a l c u l a t e d by any numerical i n t e g r a t i o n technique. A Gaussian 5 t h order i n t e g r a t i o n method i s applied here. D e t a i l s of 128 these i n t e g r a l s are found i n Brebbia (1978) as well as Banerjee and B u t t e r f i e l d (1981). Although the involved i n t e g r a l equation f o r free-surface problems i s usu a l l y more complicated than equation (13), the p r i n c i p l e of w r i t i n g down the matrix from the i n t e g r a l equation i s very s i m i l a r and st r a i g h t forward. F i n a l l y , i t i s mentioned i n both Brebbia (1978) and Banerjee (1981) that, when carrying out the numerical i n t e g r a l f o r A and C , s p e c i a l consideration should be ij i j made to take care of the s i t u a t i o n when i equals to j . These terms, A and C cannot be cal c u l a t e d d i r e c t l y as the other i i i i J terms. A n a l y t i c solutions are provided by Brebbia and Banerjee, and should be used f o r the diagonal terms of both matrix A and C. 129 APPENDIX II INTERPOLATION OF VELOCITY COMPONENTS BETWEEN ADJACENT ELEMENTS BY EQUATION (24). From equation (24) : = dx + ey +f (24). I f three sucessive elements with number i-1 , i , i+1 are considered, and t h e i r v^ values are designated r e s p e c t i v e l y by V y | ^ i i v y l f ' a n c * V y l i + 1 ' t* i e n> t n e following equations can be written : d X. + 1-1 e y i - i + f = V y i . . d X . + 1 e y i + f ' V y i d X. + 1+1 e y i + i + f = V y x. y i - l 1 " d " V y i - l X. 1 y i 1 e = V y i L 1+1 y i + l 1 - f V L y i + i -Upon sol v i n g t h i s set of l i n e a r equations, equation (24) becomes a function of x and y only, and can be used to 130 interpolate the v e l o c i t y value between these three elements. The same procedure i s applied to interpolate , i n which the column vector on the right-hand-side of the above equation i s replaced by v 131
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical procedure for potential flow problems with...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical procedure for potential flow problems with a free surface Chan, Johnson Lap-Kay 1987
pdf
Page Metadata
Item Metadata
Title | Numerical procedure for potential flow problems with a free surface |
Creator |
Chan, Johnson Lap-Kay |
Publisher | University of British Columbia |
Date Issued | 1987 |
Description | A numerical procedure based upon a boundary integral method for gravity wave making problems is studied in the time domain. The free-surface boundary conditions are combined and expressed in a Lagrangian notation to follow the free-surface particle's motion in time. The corresponding material derivative term is approximated by a finite difference expression, and the velocity terms are extrapolated in time for the completion of the formulations. The fluid-body intersection position at the free surface is predicted by an interpolation function that requires information from both the free surface and the submerged surface conditions. Solutions corresponding to a linear free-surface condition and to a non-linear free-surface condition are obtained at small time increment values. Numerical modelling of surface wave problems is studied in two dimensions and in three dimensions. Comparisons are made to linear analytical solutions as well as to published experimental results. Good agreement between the numerical solutions and measured values is found. For the modelling of a three dimensional wave diffraction problem, results at high wave amplitude are restricted because of the use of quadrilateral elements. The near cylinder region of the free surface is not considered to be well represented because of the coarse element size. Wave forces calculated on the vertical cylinder are found to be affected by the modelled tank length. When the simulated wave length is comparable to the wave tank's dimension, numerical results are found to be less than the experimental measurements. However, when the wave length is shorter than the tank's length, solutions are obtained with very good precision. |
Subject |
Water waves -- Mathematical models Water waves -- Forecasting -- Mathematical models Ocean waves -- Forecasting -- Mathematical models |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-09-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0097958 |
URI | http://hdl.handle.net/2429/28637 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1988_A1 C42.pdf [ 5.66MB ]
- Metadata
- JSON: 831-1.0097958.json
- JSON-LD: 831-1.0097958-ld.json
- RDF/XML (Pretty): 831-1.0097958-rdf.xml
- RDF/JSON: 831-1.0097958-rdf.json
- Turtle: 831-1.0097958-turtle.txt
- N-Triples: 831-1.0097958-rdf-ntriples.txt
- Original Record: 831-1.0097958-source.json
- Full Text
- 831-1.0097958-fulltext.txt
- Citation
- 831-1.0097958.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0097958/manifest