Numerical Study of Landslide-Generated Waves By Gang Yang B.E., Dalian University of Technology, P.R.China, 1992 M.A.Sc, Dalian University of Technology, P.R.China, 1995 M.A.Sc, University of British Columbia, Canada, 1996 A thesis submitted in partial fulfillment of the requirement for the degree of DOCTOR OF PHILOSOPHY In The Faculty of Graduate Studies Department of Civil Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH C O L U M B I A April, 2002 © Gang Yang, 2002 In presenting this thesis in partial fulfillment of the requirements for an advance 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 the 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. Department of Civil Engineering, The University of British Columbia, 2324 Main Mall, Vancouver, B.C., Canada, V6T 1Z4 Date: ppj,„ 2fr. 7^2- n ABSTRACT This thesis describes a time-domain boundary element method to numerically simulate linear waves generated by a horizontal-moving landslide of an arbitrary profile. The approach to setting up the numerical models is based on a Green's function method and a time-stepping procedure. The numerical models include a boundary condition at the landslide surface that accounts for wave generation, a radiation condition in the far field that accounts for open water, and free surface boundary conditions in a time-stepping procedure that account for wave propagation. By solving the first-order boundary value problem in the time domain, waves generation and propagation are simulated numerically. Although the present numerical method can be applied to a variety of landslide configurations, only vertical wall with horizontal movement is considered in the present study. The corresponding numerical models in two dimensions and three dimensions are developed. The numerical models are initially used to simulate the waves generated by a piston-type wavemaker with a periodic wave paddle movement, and the results are validated against a number of analytical solutions as well as available experimental results. Specific two-dimensional landslide problems are also simulated on a horizontalmoving wall movement and comparisons of numerical results with analytical solutions have demonstrated good agreement. A number of numerical simulations on landslide-generated waves have been performed and a selection of numerical results is presented and analyzed. The results exhibit various features of interest that are discussed. Engineering curves have been developed based on the numerical results to estimate order-of-magnitude wave height caused by landslides. A case study has been carried out to demonstrate the application of numerical models. The numerical method can be applied to account for a landslide of an arbitrary shape and movement, and a water body surrounded by surfaces of various geometries. Ill The approach and numerical models can also be extended to account for nonlinear wave effects. Overall, the numerical method and the numerical models are found to be able to provide a reasonably flexible and reliable means of studying and predicting landslidegenerated waves. iv TABLE OF CONTENTS ABSTRACT ii T A B L E OF CONTENTS iv LIST OF T A B L E S vi LIST OF FIGURES vii LIST OF S Y M B O L S ACKNOWLEDGMENTS C H A P T E R 1 INTRODUCTION 1.1 General 1.2 Literature review 1.3 Objectives C H A P T E R 2 THEORETICAL F O R M U L A T I O N x xiii 1 1 2 12 14 2.1 Three-dimensional problem 14 2.1.1 Governing equations 14 2.1.2. Linear problem 16 2.1.3 Green's function method 17 2.2 Two-dimensional problem 18 2.3 Time integration and initial conditions 19 2.4 Landslide movement 21 2.5 Sponge layer damping technique in wavemaker problems 22 CHAPTER 3 NUMERICAL FORMULATION 3.1 Field solution 24 24 3.1.1 Three-dimensional problem 24 3.1.2 Two-dimensional problem 28 3.2 Temporal solution 29 3.2.1 Free surface boundary conditions 29 3.2.2 Radiation condition 30 C H A P T E R 4 V A L I D A T I O N OF N U M E R I C A L M O D E L S 4.1 Two-dimensional wavemaker problem 4.1.1 Theoretical solution 34 34 34 V 4.1.2 Numerical results and comparisons 4.2 Two-dimensional landslide problem 35 36 4.2 1 Theoretical solution 36 4.2.2 Numerical results and comparisons 38 4.3 Three-dimensional wavemaker problem 38 4.3.1 Theoretical solution 38 4.3.2 Experimental Work 39 4.3.3 Numerical results and comparisons 40 4.4 Three-dimensional landslide problem 41 C H A P T E R 5 RESULTS A N D DISCUSSION 42 5.1 Computational considerations 42 5.1.1 Discretization schemes 42 5.1.2 Computing cost 43 5.1.3 Accuracy of matrix solution 44 5.1.4 Stability of time-stepping procedure 45 5.1.5 Numerical accuracy of discretization 46 5.2 Dimensional analysis 50 5.3 Two-dimensional landslide-generated waves 51 5.4 Three-dimensional landslide-generated waves 53 5.5 Limitations of the numerical models 54 C H A P T E R 6 ENGINEERING APPLICATIONS 58 6.1 Design curves for engineering applications 58 6.2 Case study: Skagway Harbour landslide 59 6.3 Potential of landslide-generated waves in Canada 61 C H A P T E R 7 CONCLUSIONS A N D R E C O M M E N D A T I O N S 7.1 Summary 63 63 7.2 Conclusions 65 7.3 Recommendations for further study 66 REFERENCES 68 vi LIST OF T A B L E S Table 1 Numerical simulations of two-dimensional wavemaker problem (d = 1.0 m, a = 0 0.02 m) Table 2 Test cases for three-dimensional wavemaker problems 73 74 Table 3 Accuracy of numerical discretization for two-dimensional wavemaker problem. 75 Table 4 Accuracy of numerical discretization for three-dimensional wavemaker problem. 76 Vll LIST O F FIGURES Figure 1 Sketches of various types of landslide-generated waves 77 Figure 2 Definition sketch of the problem of landslide-generated waves - plan view 78 Figure 3 Definition sketch of the problem of landslide-generated waves - side view 78 Figure 4 Plan view of simplified three-dimensional wave tank (landslide problem) 79 Figure 5 Side view of simplified three-dimensional wave tank (landslide problem) 79 Figure 6 Plan view of three-dimensional wave tank with sponge layer (wavemaker problem) 80 Figure 7 Side view of three-dimensional wavemaker problem and sponge layer technique (wavemaker problem) 80 Figure 8 Discretization of three-dimensional landslide problem 81 Figure 9 Discretization of two-dimensional landslide problem 82 Figure 10 Numerical simulation of two-dimensional wavemaker problem: wave profile at different times (d = 1 m, ao = 0.02 m, T = 2.3 s) 83 Figure 11 Numerical simulation of two-dimensional wavemaker problem: wave elevation at fixed locations (d = 1 m, ao = 0.02 m, T = 2.3 s) 84 Figure 12 Comparisons of dimensionless wave height between theoretical solution and numerical results for steady state two-dimensional wavemaker problem 85 Figure 13 Comparisons of dimensionless wave amplitudes at fixed locations between Noda's theoretical solution and numerical results on two-dimensional landslide problem when simplified landslide is moving horizontally at a constant speed 86 Figure 14 Experimental setup by Takayama (1984).. 87 Figure 15 Numerical simulation of three-dimensional wavemaker: wave profiles at different times (wave elevation in metres) 88 Figure 16 Numerical simulation of three-dimensional wavemaker: wave elevation at fixed points (numbers in brackets showing the coordinates of fixed points in metres). 89 Figure 17 Comparisons between theoretical solutions and numerical results of threedimensional wavemaker problem for b/L = 0.5 (contours showing values of H / a ) . 0 90 Vlll Figure 18 Comparisons between theoretical solutions and numerical results of threedimensional wavemaker problem for b/L =1.0 (contours showing values of H/ao).91 Figure 19 Comparisons between the numerical results and Takayama's theoretical and experimental results in x direction for three-dimensional wavemaker problem (T = 1.0 s) 92 Figure 20 Comparisons between the numerical results and Takayama's theoretical and experimental results in y direction for three-dimensional wavemaker problem (T = 1.0 s).... 93 Figure 21 Comparisons between the numerical results and Takayama's theoretical and experimental results in x direction for three-dimensional wavemaker problem (T = 1.4 s) ! 94 Figure 22 Comparisons between the numerical results and Takayama's theoretical and experimental results in y direction for three-dimensional wavemaker problem (T = 1.4 s) 95 Figure 23 C P U time as a function of the rank of the matrix equation for the twodimensional problem (100 time steps) 96 Figure 24 C P U time as a function of the rank of the matrix equation for the threedimensional problem (100 time steps) 97 Figure 25 Effect of number of segments for two-dimensional wavemaker 98 Figure 26 Wave elevations at the centerline of the wave basin for three discretization schemes (t = 0.5 s) 99 Figure 27 Wave elevation as a function of distances for two-dimensional landslidegenerated waves (d = 0.63 m, a = 0.01 m, T/2 = 0.9 s) 0 100 Figure 28 Wave elevation as a function of time at fixed locations for two-dimensional landslide-generated waves (d = 0.63 m, a = 0.01 m, T/2 = 0.9 s) 0 101 Figure 29 Numerical results of dimensionless wave heights as a function of dimensionless distance for various d/gT 2 values on two-dimensional landslide- generated waves 102 Figure 30 Numerical simulation of three-dimensional landslide-generated waves (d = 0.63 m, a — 0.01 m, T/2 — 0.5 s, dimensions are in metres) 0 103 IX Figure 31 Contours of dimensionless wave height H / a as function of dimensionless 0 location for b/d = 1.0 and various d/gT 2 values from numerical simulations of three-dimensional landslide-generated waves 104 Figure 32 Contours of dimensionless wave height H / a as function of dimensionless 0 location for b/d = 2.0 and various d/gT 2 values from numerical simulations of three-dimensional landslide-generated waves 105 Figure 33 Contours of dimensionless wave height H / a as function of dimensionless 0 location for b/d = 4.0 and various d/gT 2 values from numerical simulations of three-dimensional landslide-generated waves 106 Figure 34 Dimensionless wave heights at the centerline of wave tank (y/d = 0) as function of x/d for various b/d and d/gT 2 values from numerical simulations on three- dimensional landslide-generated waves 107 Figure 35 Fitting curve for dimensionless wave height as function of dimensionless distance for two-dimensional landslide-generated waves when d/gT =0.02 2 108 Figure 36 Design curves of a and P parameters in determining two-dimensional landslide-generated wave heights governed by H / a = a(x/d) for various d/gT 2 values 109 p 0 Figure 37 Design curves of a and p parameters in determining three-dimensional landslide-generated wave heights at the centerline (y/d = 0) governed by H / a =a(x/d) for various b/d and d/gT p 0 2 values 110 Figure 38 Location of Skagway Harbour and the landslide area (Kulikov et al. 1996). I l l Figure 39 Geologic section in the middle of landslide area (Cornforth et al. 1996) 112 LIST O F S Y M B O L S a = = b = Ajj, Bjj == c = C n = displacement of landslide; displacement amplitude of landslide; landslide width; calculated coefficients; wave celerity; wave celerity in the normal direction; C = constant; d = Di, D , D 2 3 E water depth; = discretization scheme; = error in numerical solution; g = acceleration of gravity; G = Green's function; J - k = m = Jacobian matrix; wave number; total number of facets; w = number of facets on the wall surface; m 0 - number of facets on the still water surface m c = number of facets on control surface; m n distance in the direction of normal; N shape function; p,q local coordinate system - p = parameter, see Eq. 3.13 Q = parameter, see Eq. 3.13 w parameter, see Eq. 3.13 parameter, see Eq. 3.13 Pv r S s b radial distance measured from the moving wall; - boundary surface; tank bottom surface, see Fig. 3; xi S = c control surface, see Figs. 2 and 3; Sf = free water surface, see Figs. 2 and 3; S = still water surface, see Figs. 2 and 3; - still wall surface, see Figs. 2 and 3; = landslide or wave paddle surface, see Figs. 2 and 3; 0 S S t w t = t' = time; dimensionless time; T = wave period; Tj = dimensionless time interval; u = dimensionless wave number; V = velocity of wave paddle or landslide; n = normal component of landslide velocity; v' = dimensionless V Vj = w , w m n = x, y, z = dimensionless constant speed; weight functions at Gauss points; Cartesian coordinate system; x' = dimensionless distance in x direction; X = point at (x, y, z); X' = point at (£, \|/, Q ; a, P = fitting parameters of design curves for landslide-generated waves; dd, Pd = tuning factors for damping parameter; P> PP = parameters, see Eqs. 4.12 and 4.13; = Courant number; A p 0 r = an integration of Green's function; r\ = free surface elevation; r| e = exact water surface elevation; r| 0 = wave amplitude; (j) = Q velocity potential; •= landslide volume; u. = damping parameter; xii v = parameter, see Eq. 4.11 a dimensionless wave frequency; = co = \|/, C, - angular wave frequency; distances in the x, y and z directions; Xlll ACKNOWLEDGMENTS The author would like to express his sincere appreciation and gratitude to his supervisor, Dr. Michael Isaacson, for his advice and guidance throughout this research. Thanks are extended to Dr. Sander Calisal for his brilliant teaching and advice on the Boundary Element Method, to Dr. Carl Ollivier-Gooch for his timely comments on the matrix solution and stability of numerical procedures, and to Dr. Ricardo Foschi and Dr. Greg Lawrence for their useful comments and suggestions. The author is also grateful to his friends, including Dr. Sundaralingam Premasiri and Mr. Eric Morris, for their support and encouragement.. The author wants to give special thanks to his father for his love and inspiration to pursue education. Finally the author wants to thanks his wife, Yajun Shi, for her love and support. Chapter 1 Introduction l CHAPTER 1 INTRODUCTION 1.1 General Water waves generated by landslides have been of great interest to ocean, coastal and hydraulic engineers for a number of years. Devastating damages caused by landslidegenerated waves have been recorded all around the world. However, the prediction of landslide-generated waves, especially in three-dimensional cases, is still in the research area. A landslide is often defined as the downward sliding or falling of a relatively dry mass of earth, rock, or mixture of the two. The slope is usually steep and the movement rapid or moderately rapid. Landslides may be categorized into several types according to the rate and character of movement and the kind of material. In coastal or waterfront areas, depending on the geometry of shoreline, landslides are often categorized into vertical rock fall, horizontal wall movement, massive slump, vertical seabed movement due to earthquake, sliding wedge down a slope, or submarine landslide (see Fig. 1). Experimental study, theoretical development, and numerical modeling have been the three main approaches to studying landslide-generated waves. Experimental study is often very time-consuming and costly, thus has mostly been used to validate theoretical solutions and numerical model results. Theoretical development becomes very difficult due to the complicated geometry and moving boundary of landslide, thus is mostly limited to the two-dimensional cases. Numerical modeling, on the other hand, has the advantage of being able to handle complicated boundary geometry and simulate various forms of landslide movement. The drawback of numerical modeling is its limitation on the size of the domain and its high need for computer hardware resources (CPU speed, memory, etc.). With the help of fast-growing computer hardware industry, numerical modeling has been more and more favorable to researchers. Finite Difference Method (FDM), Finite Element Method (FEM), and Boundary Element Method (BEM) are the three main numerical modeling methods. F D M and F E M both digitize the problem domain into a finite number of elements, and solve the Chapter 1 Introduction _2 governing differential equations through some approximation method. B E M on the other hand, digitize only the boundary of the problem domain, and solve the governing differential equations based on the boundary values. In many numerical simulations in coastal engineering, B E M has become favorable for many researchers. The following section gives a brief review of efforts made on the study of landslidegenerated waves and related topics. 1.2 Literature review The problem of landslide-generated waves has been treated as various theoretical and numerical models by different authors. These landslide models include horizontally moving wall, vertically falling box, two-layer flow, and a solid body sliding down a slope. Some experimental work has been performed to verify the models or to simulate real cases. Most of the previous studies have been focusing on two-dimensional cases, while tentative solutions for three-dimensional cases are limited. Previous studies on landslide-generated waves and some related topics are grouped into following categories: • Two-dimensional horizontally moving wall or vertically falling box; • Two-dimensional vertical deformation at the bottom of the seabed or reservoir; • Two-dimensional or three-dimensional solid body sliding down a slope; • Two-dimensional or three-dimensional two-layer flow; • Experimental work and other models; • Wavemaker Problems; • Relevant studies based on Green's function method and time-stepping method. These are briefly reviewed in turn as follows. Chapter 1 Introduction 3 Two-dimensional horizontally moving wall or vertically falling box Wiegel et al. (1970) studied waves generated by vertically falling rockfall (referring to Fig. 1(a)). A two-dimensional linear theoretical approach was presented and some experiments were carried out. The rockfall was simulated as a rectangular box falling down from the above the water surface or from the water surface. The nonlinear effects caused by the dimensions of the rockfall were discussed. Noda (1970) simplified the problem of landslide-generated waves into two distinct types and presented corresponding analytic models. One type of landslides was a vertical landslide, in which case a landslide entered and moved through the water vertically (referring to Fig.l (a)). The model used to approximate vertical landslides was a twodimensional box falling to the bottom at the end of semi-infinite channel. The other type was called horizontal landslide, in which case a landslide entered and moved horizontally through the water (referring to Fig.l (b)). The horizontal landslide was modeled by a two-dimensional moving wall. The analytic solutions were based on the linearized theory of gravity surface waves and the applications of such solutions were determined by the comparison with the experimental work. .For the vertical landslide problem, Noda compared his theoretical results with experiments performed by Wiegel et al. and that by himself, and used his linear theory solution to predict the nonlinear waves generated by Lituya Bay landslide in Alaska. He concluded that the landslide-generated wave surface profile was only dependent on the volumetric inflow rate, and thus the distribution of velocity-depth under a potential landslide could be neglected. While Noda neglected impact effects, experiments carried by Wiegel et al. showed that for reasonably large distances from the box, the wave surface profile was essentially independent of the impact and splash produced by the box. Although a vortex was formed next to the box as the box fell during the experiments, the linearized theory which assumed irrotational motion still predicted very well the wave profile in the oscillatory region. As for horizontal landslide problem, the theoretical model results were not validated due to limited experimental data. Chapter 1 Introduction 4 Gozali and Hunt (1989) used the method of characteristics to obtain numerical solutions for water waves generated in a reservoir by a relatively local landslide. The landslide was modeled by a vertical wall moving a horizontally over a distance of at least ten reservoir depths into the reservoir and following an assumed velocity variation with time (referring to Fig.l (b)). The generated waves were approximated by nonlinear, nondispersive long-wave approximations. The results were shown with dimensionless plots of maximum wave heights as functions of time and distance for both one-dimensional and axis-symmetric solutions. Since the landslide was much simplified, and no experimental results or other source of data were compared, the use of the numerical model was limited to some certain conditions and questions about the validation remain. Two-dimensional vertical deformation at the bottom of the seabed or reservoir Raney and Butler (1976) gave a two-dimensional numerical approach to the landslide-generated wave problem by neglecting the vertical component of velocity in the governing hydrodynamic equations, using a finite difference method. In this model, the landslide was represented by a time-dependent vertical deformation of the bottom of the reservoir plus additional terms to represent the effect of the landslide due to viscous and inertia forces (referring to Fig. 1(c)). The numerical results were compared with the experimental data from a physical model of Libby Dam and Kootenai River in western Montana, U S A , and good agreement was obtained. They concluded that, for small landslide velocities, the viscous drag and pressure drag contribution to wave heights and velocities were small, and the physical displacement of the water by the landslide was the predominate factor. For large landslide velocities, significant contributions to wave heights and propagation velocities were produced by the viscous drag and pressure drag. Hunt (1988) modeled waves generated in a reservoir by a landslide by injecting an instantaneous point source of fluid through the bottom of a reservoir of infinite extent based on linear wave theory (referring to Fig.l (d)). It was assumed that times were large relative to the duration time of the landslide and that distances were large relative to a characteristic horizontal dimension of the volume of water displaced by the slide material. Solutions were obtained for two dimensional and axisymmetric problems. The Chapter 1 Introduction 5 expressions for maximum wave amplitudes, length, period and position were obtained. The calculated results were compared with Noda's (1970) experimental results and good agreement was obtained. The results may be useful for providing boundary conditions for either numerical or experimental models of wave run-up on sloping beaches. Two-dimensional or three-dimensional solid body sliding down a slope Harbitz (1992) developed a three-dimensional mathematical model based on the hydrodynamic shallow water equations for numerical simulation of water waves generated by the submarine Storegga Slides on the Norwegian continental slope. He focused on the total water displacement and the shear stress acting between the slide masses and the fluid. The slide was assumed to be a box with smoothed sides, and it's motion assumed to be a simple function of time and the displacement. The equations were solved numerically by a finite difference technique. The wave run-up heights generated by the Second Storgga Slides was computed on the basis of the numerical model and the results were in remarkably good agreement in some locations. It was concluded that the generated wave heights increased significantly with the acceleration of the slide, and a slide with a smaller volume might easily cause higher waves than a bigger one i f the smaller slide had a sufficiently higher initial acceleration. Also, the analysis of the effects of shear stress at the interface between the water and the slide body proved that this effect must be included, except for run-up zones with gentle beach slopes. As applying the numerical model, a correct estimate of the maximum slide velocity is essential. Harbitz et al. (1993) developed another similar three-dimensional numerical model based on the hydrodynamic shallow water equation for numerical simulation of slidegenerated waves in fjords. The numerical model was used to simulate the slide catastrophe in Tafjord, western Norway, 1934, and showed close agreement. As a result of the numerical calculation, wave reflection and standing wave oscillations were predicted due to the limited area of the fjords. Heinrich (1992) set up a two-dimensional numerical model based on the program Nasa-Vof2D to study the generation, propagation, and run-up on the shore of water Chapter 1 Introduction waves created by landslides. 6 The Nasa-Vof2D, is a nonlinear Eulerian code, which solves the complete incompressible Navier-Stokes equations by a finite difference method. A n experimental study on nonlinear waves generated by a two-dimensional triangular body sliding along a 45° inclined plane was conducted to examine the numerical Nasa-Vof2D extension (referring to Fig. 1 (e)). Two types of landslides, a submarine landslide and an aerial landslide were studied. The computed wave profiles showed very close agreement with the experimental ones, except when free-surface turbulence occurs, which the numerical method cannot simulate. The advantages of this numerical model also include the abilities of dealing with nonlinear wave propagation and vertical deformation of channel bottom. A three-dimensional Nasa-Vof method would be ideal to deal with more realistic problem, taking the attribution of turbulence into account, yet is very complicated and extremely expensive computationally. Two-dimensional or three-dimensional two-layer flow Jiang and LeBlond (1992) studied the submarine landslide and the surface waves that it generates. They presented a formulation of the dynamics of the problem, where the landslide was treated as the laminar flow of an incompressible viscous fluid on a gentle uniform slope and water motion was assumed irrotational (referring Fig. 1(f)). Long-wave approximation was adopted for both water waves and the mudslide and the dispersion of waves was neglected. A finite-difference method was used to solve the resulting differential equations. One-way coupling (bottom deformations affect the free surface) and full coupling (surface pressure gradients react on the mud flow) are considered with the behavior of the mud flow. The calculations showed that three major waves were generated. The first wave was a crest which propagated into deeper water, followed by a trough which was a force wave propagating with the speed of the slide front to deep water. The third wave was a small trough which propagated shoreward. They found that two major parameters determined the interaction between the slide and the waves: the density of sliding material and the depth of initiation of the mudslide. They also concluded that the amplitudes of the waves depended primarily on the density of mud, the initiation depth of the mudslide, the volume of the landslide, and the viscosity of mud. The slide with larger density and volume, and smaller viscosity in shallower Chapter 1 Introduction 7 water generated larger waves. The ratio of the energy transfer from the slide to the waves was found to have a maximum value at the beginning stage and decreased to 2-4% when the waves propagated away from the slide site. This study was the first one that addressed the details of underwater landslide dynamics reactively coupled to water wave dynamics. Later, Jiang and LeBlond (1993) presented another numerical model that simulates the coupling of Bingham plastic mudslide on a gentle uniform slope with the surface waves which it generates. The landslide was treated as an incompressible Bingham plastic flow and the water motion was assumed irrotational. The long-wave approximation was adopted for both water waves and the mudslide and the dispersion of waves and potential turbulent mixing wee not considered. The resulting differential equations were solved by a finite difference method. They found that three parameters dominated the magnitude of the waves (with mud volume fixed): the density of mud, the yield stress of the mud, and the depth of water at the mudslide site. It was concluded that the Bingham plastic behavior of the mud significantly reduced the extent and the speed of the mudslide and also the magnitude of the surface waves generated. The comparison of the solution of the Bingham plastic mudslide model with a snow flow test was presented. It is noted that i f the yield stress of the mud is zero, the Bingham fluid model reduces to the viscous fluid model, which is mentioned above. The evolution of the surface waves generated by Bingham plastic slides is similar to that of the waves generated by a viscous slide. The wave current speed and the wave amplitude, however, are much smaller. It is noted that the underwater slide problem that were treated are two-dimensional and the solution is only valid for a slide and waves on a gentle slope. Jiang and LeBlond (1994) extended their numerical model into three-dimensions. They assumed that the submarine slides behaved as a viscous fluid and applied the longwave approximation for both the water waves and the slide. The slide was assumed to have parabolic surface with a rectangular bottom periphery initially. The numerical results indicated that three-dimensional effects were important in the effectiveness of tsunami generation. It is noted that no experimental data or field data was compared with the numerical results. Chapter 1 Introduction 8 Rubino et al. (1994) studied the generation of tsunamis by mudslides and the near and far-field wave propagation in a three-dimensional framework, based on a two-layer shallow water model. The generation and initial stage of propagation was described through a nonlinear, hydrostatic, turbulent, two-layer, shallow water model, with a local deep layer representing the sliding sediments. The far-field propagation was simulated by a non-hydrostatic, nonlinear, regularized Kadomtsev-Petviashvili model, which is able to describe the evolution of long, nonlinear, weakly dispersive and weakly 3D waves. They concluded that the phase dispersion effects were important in some realistic circumstances, and thus could not be neglected. Experimental work and other models Eie et al. (1971) did some experiments modeling landslides in a reservoir behind the Viddalsvatn rock fill dam, Aurland power plant, and rock slides into the lake Ardalscvatn. The landslide was modeled by a solid mass in position attached to a string which ran over two pulleys and carried an adjustable counterweight at the other end to control the velocity of the slide. The results indicated that the wave heights attain a maximum value for a certain slide velocity, and an increase in thickness of the slide tended to increase the wave heights, even i f the mass of the slide was unchanged. The model results showed that slide form and the slide velocity were significant parameters. Rockslide was modeled by a piece of concrete which was suddenly released and feel by gravity into the water. No enough results were obtained to lead to conclusions. Koutitas (1977) studied the water waves generated by landslides in a long channel, and set up a one-dimensional numerical model using a Finite Element Method. The numerical model was based on the classical St. Venant equations of unidirectional flow, and the inclusion of the landslide effect was achieved by the introduction in the continuity equation of a term describing the temporal variation of the cross section of the conduit. The numerical model was not validated by any experimental work or field data. Chaudhry et al. (1983) conducted a physical model and set up a mathematical model to study the generation and propagation of waves created by hypothetical movement of the 1.4 billion m Downie slide into the Revelstoke reservoir, located on the 3 Chapter 1 Introduction 9 Columbia River, British Columbia, Canada. Two physical models were constructed reproducing about 17 miles of reservoir and simulating the movement of the slide. Because of the narrow width of the canyon, wave propagation becomes approximately one-dimensional once the waves have traveled some distance down the reservoir. Therefore, a mathematical model using an explicit finite difference method to numerically integrate the one-dimensional Saint Venant Equations was used to simulate the wave propagation of the waves in the reservoir about forty miles on either side of the slide. Comparisons of the computed and measured water levels in the model showed good agreement. Watts (1998) studied water waves generated by unsteady motion of a submerged object in two-dimensional cases. Underwater landslide was assumed to be governed only by a characteristic distance and a characteristic time and characteristic wave amplitude was obtained. A non-dimensional wavemaker curve of characteristic water wave amplitude was constructed from the landslide length, the initial landslide submergence, the incline angle measured from the horizontal, the characteristic distance of landslide motion, and the characteristic duration of landslide motion. The wavemaker curve was confirmed by two-dimensional experiments. The underwater landslide experiments were conducted with a P V C solid block that had a horizontal top face and vertical front sliding down a slope in a wave tank. Wavemaker data was obtained from experimental measurements of the maximum near-field wave amplitude as well as solid block initial acceleration and terminal velocity on a 45-degree incline. Wavemaker Problems Wavemaker problems have been studied extensively both theoretically and numerically. Wavemaker problems have similar phenomena with landslide-generated waves in terms of governing equation and boundary conditions. A slope failure beneath a vertical quay wall may introduce a horizontal movement of the quay wall. This movement can be simulated as a horizontal moving panel, which is very similar to a wavemaker problem, except that the panel movement is no longer periodic. Chapter 1 Introduction 10 Dean and Dalrymple (1984) gave close-form solution of a two-dimensional wavemaker problem for both piston and flap type wavemaker motions. Takayama (1984) studied a three-dimensional wavemaker problem in a wave basin. A close-form solution was developed for the wave generation by a single wave panel and the results were verified against experimental results. Relevant studies based on Green's function method and time-stepping method Majority of numerical wave tank employ boundary integral equation method based on Green's theorem to determine the velocity potential of the flow field. The boundary surfaces are discretized by boundary element method. The boundary element method may be classified as constant panel method, linear boundary element method, and higher-order boundary element method such as quadratic boundary element method and cubic boundary element method. The constant panel method is the easiest in formulation, but fails to model the intersections between the Dirichlet and Neumann boundaries accurately. Higher-order boundary element method is more complicated in formulation, but can be used to avoid the singularity caused by the constant panel method at the intersections between the free surface and solid boundaries. Time integration of the free surface conditions was performed by Longuet-Higgins et al. (1976) using fourth-order Adam-Bashforth-Moulton predictor-corrector scheme and fourth-order Runge-Kutta method. Many researchers have used the same technique recently. Dommermuth et al. (1987 and 1988) studied the two methods and obtained Courant condition for the fourth-order Runge-Kutta scheme and found week instability in the fourth-order Adam-Bashforth-Moulton scheme. Open boundary condition is another important issue in numerical wave tank simulation. Radiation condition, damping beach, active absorber and some other techniques have been employed. Some researchers have employed combinations of these techniques. Some recent representative works on numerical wave tank are listed below. Chapter 1 Introduction 11 Isaacson and Cheung (1991) developed a numerical method based on potential flow theory for simulating transient, second-order interactions of ocean waves with large fixed bodies of arbitrary shape in two and three dimensions. The flow potential to second order was defined with respect to a time-independent boundary which included the still water surface, and its solution involved a numerical discretization of an integral equation based on Green's function theorem. The free surface boundary condition and the radiation condition were satisfied to the second order by a numerical integration to account for transient effects. Applications of the method were made to studies of regular wave diffraction around a fully submerged and a semi-piercing circular cylinder in three dimensions. The stability and numerical accuracy of the proposed solution and the treatment of the radiation condition to second order were examined. Comparisons of computed wave forces and runup were made with previous theoretical and experimental results and they indicated favorable agreement. Isaacson et al. (1993) developed a numerical method on the study of transient, nonlinear wave propagation in a laboratory. The nonlinear free surface boundary conditions and the wave generator boundary condition are expanded about the corresponding equilibrium positions by perturbation expansions. The Sommerfeld radiation condition applied on the downwave control surface is modified to incorporate a time-dependent celerity to account for transient effects. The boundary conditions are then satisfied to second order by a method based on Green's theorem. A series of laboratory tests have been performed to verify the numerical model. The case of irregular wave packet was also studied. The experimental results of the free surface elevations obtained at twelve different locations along the flume were compared with the corresponding numerical predictions and favorable agreement was indicated. Boo et al. (1995, 1996) studied nonlinear diffraction of a vertical cylinder problem in a three-dimensional numerical wave tank. The numerical wave tank was modeled using a quadratic element method. A discontinuous element technique was employed at the intersections between water surface and solid boundaries. Open boundary condition was treated by introducing an artificial damping and potential stretching technique. The fourth-order Adam-Bashforth and Adam-Moulton scheme was employed in time Chapter 1 Introduction 12 integration. The open boundary condition was tested by simulating an Airy wave for long duration and proven to be effective. Skourup et al. (1996, 1998) proposed an active absorption in a numerical wave tank. The same techniques by Isaacson and Cheung (1991) were used to develop the second-order wave tank using perturbation series and Taylor expansion. A fourth-order Adam-Bashforth method and a fourth-order Adam-Moulton method were used in time domain integration. Active wave absorption was proposed based on using an array of independently controlled two-dimensional wavemakers along the lateral boundary as active wave absorbers. Some numerical examples were presented. 1.3 Objectives Most of the previous studies on landslide-generated waves have been found to be either algebraically complicated or numerically expensive. Due to the variation of landslide geometry and movement, the problem becomes sophisticated and theoretical solutions become extremely difficult to obtain. On the other hand, numerical methods, such as the Finite Element, Finite Difference, and Boundary Element Methods, have been successfully applied to many complicated boundary value problems. Some researchers have applied a Finite Difference Method (such as the V O F method) and Finite Element Method to solve the Navier-Stokes equations. However, Finite Difference and Finite Element Methods are expensive computationally in solving three-dimensional problems hence are mostly applied to two-dimensional problems. The Boundary Element Method, with discretization applied only to the surfaces, has the advantage of less computational effort when a basic function is available for the governing equations. The primary aim of the present investigation is to study landslide-generated waves using a Boundary Element Method based on a Green's function method in conjunction with a time-stepping procedure, along the lines proposed by Isaacson et al. (1994). The boundary value problem includes a boundary condition at the landslide surface that accounts for wave generation, a radiation condition in the far field that accounts for open water, and free surface boundary conditions that account for wave propagation. The latter are applied to a time-stepping procedure. Linear waves and a constant water depth Chapter 1 Introduction are assumed. 13 Although the present method can be applied to various landslide configurations, only vertical wall with horizontal movement is to be explored in the present study. Thus, although the general numerical method has been developed previously in the context of other hydrodynamic problems, such as wave interactions with fixed and floating structures, the extension of the method to the case of landslidegenerated waves and the associated findings are considered to be new contributions. The numerical models for both two-dimensional and three-dimensional cases are to be set up for the simplified landslide configurations and to be validated against a number of theoretical and experimental results. A number of numerical simulations are to be performed to study the characteristics of landslide-generated waves. Based on the numerical simulation results, parametric engineering curves are to be developed to estimate order-of-magnitude wave heights for the first time. Numerical simulation results are expected to supply wave inputs for other numerical wave models that account for wave diffraction, refraction and shoaling so as to estimate wave propagation away from the landslide site. Finally, the method developed in this thesis may be extended in the future to account for non-linear waves, various landslide movements, water domains with arbitrary boundaries and various other scenarios. Chapter 2 Theoretical Formulation 14 CHAPTER 2 THEORETICAL FORMULATION The problem of landslide-generated waves becomes complicated largely due to the uncertainty of the movement of landslide as well as the difficulty in determination of the instantaneous free surface boundary conditions, especially when the nonlinear component of wave height becomes significant, mostly in shallow water cases. Simplifications are often made to resolve these problems in simulating the movement of landslides and the free surface. In the present study, the problem of landslide-generated waves is simplified into a horizontal moving wall in a water tank. The following sections describe the theoretical formulation of both three-dimensional and two-dimensional cases. 2.1 Three-dimensional problem The landslide is treated as a block of mass moving down a slope into a three-dimensional wave tank of constant depth. A definition sketch of the problem is shown in Figs. 2 and 3. A Cartesian coordinate system (x, y, z) is defined with x measured in the direction of wall movement, y measured perpendicular to x in a horizontal plane, and z measured upwards from the tank bottom. 2.1.1 Governing equations Let t denote time and n the free surface elevation above the still water level. The tank bottom is assumed impermeable and horizontal. With the fluid assumed incompressible and inviscid, and the flow irrotational, the fluid motion can be described by a velocity potential (|)(x, y, z, t) which satisfies the Laplace equation within the fluid domain at any time t: V c|) = 0 at any timet 2 (2.1) The potential § is subject to boundary conditions on the impermeable tank bottom Sb, the moving slide surface S , the still slope surface S surrounding the slide surface, and the w t free surface Sf, which are given respectively as (again, at any time t): Chapter 2 Theoretical Formulation 15 on S (z = 0) (2. 2) on Sf (z = d + n) (2. 3) on Sf (z = d + n) (2. 4) on St(|y|>b/2) (2. 5) on S ( y <b/2) (2. 6) b d§ dn c9<() 3z dt dx dx t 9 n 3<j) dr) =0 dy dy 1 2 + gn + -|V<t)| =0 dx w where d is water depth; b is the width of landslide; g is the acceleration due to gravity; a(t) is the horizontal displacement of the landslide; V (t) is the normal component of the n landslide velocity V(t); and n is the distance in the direction of the unit normal vector which is directed outward from the fluid region. Equations 2.2 and 2.5 correspond respectively to the kinematic boundary conditions on the impermeable tank bottom and the still wall surface where the normal fluid velocity component vanishes. Equation 2.3 is the kinematic free surface boundary condition requiring that the fluid particle velocity normal to the free surface is equal to the velocity of the free surface itself in that direction. Equation 2.4 is the dynamic free surface boundary condition stating that the pressure at the free surface expressed in terms of the Bernoulli equation is equal to the atmospheric pressure with neglecting the effects of surface tension. In the present application, the atmospheric pressure can be set to zero without inducing any loss of generality. Equation 2.6 is the kinematic boundary condition on the moving wall surface where the normal fluid velocity component is equal to the horizontal wall velocity. In addition, the boundary condition in the far field from the landslide has to be defined. Here, the Sommerfeld radiation condition, originally derived on the basis of spatially and temporally periodic conditions of a potential function to account for unsteady wave motions, is applied in the far field. In the time domain, the Sommerfeld radiation condition for the wave potential is given by: Chapter 2 Theoretical Formulation Sty Vr] \dt df dr j =0 16 as r ^ co (2.7) where c is the celerity of the waves and r is radial distance measured from the moving wall. This has been modified by Orlanski (1976) so as to apply on a control surface S c located at a finite distance form the moving wall. The term Vr is no longer necessary, and the condition is applied in the direction normal to the control surface: a<> i a* ^ +c £= 0 dt on n onS c (2.8) where c is the time-dependent celerity of the scattered waves on the control surface S in n c the direction of the normal n. 2.1.2. Linear problem The major difficulties in solving the problem described in the previous section arise primarily from the poorly defined radiation condition involving transient waves and from the two nonlinear free surface boundary conditions applied on the instantaneous free surface which itself is unknown a priori. Considering linear waves only, the nonlinear terms in the boundary conditions on the free water surface are neglected. To further simplify the problem, assuming the wave height is moderate and the water depth is large compared with a typical wave length, the free water surface boundary conditions are applied to the still water surface. Thus Eqs. 2.3 and 2.4 are simplified as: **-*Uo o n S ( z = d) (2.9) ^ + gr) = 0 dt o n S ( z = d) (2.10) dz dt 0 0 Also, the slide is further simplified as a vertical wall, and it is assumed that the displacement of the moving wall is small compared with the typical water depth. The boundary condition on the moving wall is then evaluated on the original still wall surface (see Figs. 4 and 5), thus Eq.2.6 becomes: Chapter 2 Theoretical Formulation ox = V(t) 17 on S (x = 0, w lyl < b/2) (2. 11) 1 1 2.1.3 Green's function method A well-posed boundary value problem at a time t has been set up with Eqs. 2.1, 2.2, 2.5, 2.8, 2.9, 2.10 and 2.11. A boundary integral method involving a Green's function may be used to solve the boundary value problem at the fixed time t. Since the potential (j) is a harmonic function, the second form of Green's theorem may be applied over a closed boundary containing a fluid region, relating the values of < > j and its normal derivative 5<j>/3n along the boundary. The potential (j)(X) at a point X ( x , y , z ) on the surface S is expressed as: • M = i-j[G(X,X')|i(X')-*(X')|^(X,X')kiS 2n (2. 12) s Here X ' represents a point (E,,\y,C) on the surface S over which the integration is performed, and n is measured from the point X'(^,i|/,C). For X located inside the boundary, the factor 1/271 is simply replaced by 1/4TC. G ( X , X ' ) is a Green's function which can be interpreted as the potential at any field point X ( x , y , z ) due to a source of unit strength at X'(^,i|/,C) and which satisfies the Laplace equation in the domain except at the point X(x, y, z) = X'(£,, \\i,C): V G ( X , X ' ) = 8(x-^)5(y-\(/)5(z-c;) inD 2 (2.13) where 8 is a dirac-delta function. The solution to Eq. 2.13 is as follows: G(X,X') = r (2.14) where r is the distance between points X(x,y,z) and X'(t„\\i,C). In the present problem, the surface S would comprise of the still wall surface S , the t moving wall surface S (x = 0), the still water surface S , the control surface S and the w tank bottom Sb, as indicated in Figs. 4 and 5. 0 c However, with the moving wall and Chapter 2 Theoretical Formulation 18 subsequently the wave field symmetric about the x-z plane and the horizontal tank bottom, it is more efficient to simulate one half of the total surface with tank bottom excluded from S and to choose a Green's function that accounts for the double symmetry about the x-z: plane and the tank bottom. Thus the Green's function can be written as: G(X,X')=t- (2.15) k=l k r where rk is the distance between the points X and X' . k The source points X ' are defined k such that XJ = (t,,\\f,C) is a point on S; X ' = (^,-\\),C) is the image of XJ about the x 2 axis; XJ, = (t,,-\\i,-C) is the image of X ' 2 about the tank bottom; and X ' = (^xy-t) is 4 the image of XJ about the tank bottom. Thus, r± is given by: r 1 = [fe-x) r 2 = fe-x) 2 + (v,;-y) (M/ 2 + 2 y) (c;-z) ] 2 + (;-z) J 2 + 2 / 2 1 / 2 + r = [fe - x ) + ( + y) + fe + z) f 2 2 3 r =[fe-x) (H/-y) 2 4 2 2 V + 2 + (c; z) ] 2 + l/2 (2.16) (2.17) (2. 18) (2.19) Thus Eq. 2.12 actually represents an integration of the surfaces in the four quadrants of the doubly symmetric geometry. 2.2 Two-dimensional problem A simplification to the formulation described above may readily be adopted to treat the corresponding two-dimensional problem in the x-z plane with all boundary conditions remain unmodified (see Fig. 5). In adapting the method used here, the surface integral equation deriving from Green's theorem reduces to a line-integral equation in which the influence of the coordinate y is absent. When X is located on a smooth part of the boundary, the integral equation for the potential is then given as: Chapter 2 Theoretical Formulation 19 Hx) = -ij[G(X,xO|i(xO-ct»(xO^(X,X')ldS TT -\_ on on (2.20) Here X ' represents a point (^,C) on the surface S over which the integration is performed, and n is measured from the point X ' . For X located inside the boundary, the factor -1/71 is simply replaced by - 1/2TC . The Green's function for two-dimensional problems is expressed as: (2.21) G ( X , X ' ) = ln(r) In the present context the surface S would comprise of the still body surface S t below the still water level, the still water surface S , the control surface S surrounding 0 the water domain and the tank bottom. c However, because of the assumption of a horizontal tank bottom, the tank bottom can be excluded from S and a Green's function can be chosen to account for the symmetry about the tank bottom. Such a Green's function can be expressed as: G(X X') = i > ( r ) f (2.22) k k=l where rk is the distance between the points X and X ' . The source points X ' are defined such that XJ ={£,,C) is a point on S; X ' =(^,-C) is the image of XJ about the tank 2 bottom. Thus rj is given by: = = fe-x) (C-z)f 2 + [fe-x) fe+ z)2]1/2 2 + (2.23) (2.24) The integral equation 2.20 actually represents an integration of the surface S and its image about the tank bottom. 2.3 Time integration and initial conditions In the preceding formulation, the boundary conditions are defined on the still wall surface, the control surface and the still water surface, and the solution is expressed in Chapter 2 Theoretical Formulation 20 terms of an integral equation relating the velocity potential and its normal derivative on the boundary. For a boundary value problem in this case, i f either the velocity potential or its normal derivative on each portion of the boundary is known, the solution to the integral equation can therefore be obtained by a numerical procedure from which the remaining unknowns can be evaluated. In the present application, the normal derivative of the velocity potential on the wall surface is known and given by Eqs. 2.5 and 2.11. The normal derivative of the velocity potential on the control surface is given in the radiation condition, Eq. 2.8. However, this normal derivative depends on the temporal derivative of velocity potential on the control surface, which is unknown before the solution is obtained. Also, the normal derivative of the velocity potential given in the kinematic free surface boundary condition, Eq. 2.9, depends on a temporal derivative which is unknown a priori. On the other hand, the velocity potential on the control surface and the still wave surface at an advanced time can be evaluated by a time-integration of the temporal terms respectively in the radiation condition, Eq. 2.8, and the dynamic free surface boundary condition, Eq. 2.10. With the flow field know at time t, the velocity potential on the control surface and the still water surface at an advanced time (t + At) is given by: c()(t + At) = (t)(t)+ f - ^ d t on S and S c (2.25) 0 in which d§/dt is evaluated from the corresponding boundary condition on S or S , Eqs. c 0 2.8 and 2.10. According to Eq. 2.10, the time-integration of the wave potential on the still water surface requires knowledge of the free surface elevation, which can be obtained by a time-integration of the kinematic free surface boundary condition, Eq. 2.9. Likewise the free surface elevation at an advanced time (t + At) is given by: t+At p. n(t + At) = ri(t)+ f^t on S 0 a in which dr\/dt is evaluated from the kinematic free surface boundary condition. (2.26) Chapter 2 Theoretical Formulation 21 Initial conditions are needed to carry out the computation. At time t = 0, the flow field is known to be still, i.e., the velocity potential and the first derivative of the potential are zero throughout the field, so as the free surface elevation. For t = At, the moving wall boundary condition, i.e. the derivative of the velocity potential on S , is introduced, and w thus the Green's function method (Eq. 2.12) is applied to solve for (j) on the wall boundary and d§/dn on both the water surface and the control surface. Then the velocity potential on the water surface at time t = At is obtained by the water surface boundary conditions, Eqs. 2.9 and 2.10, and the time-integration of the flow parameters in Eqs. 2.25 and 2.26. At the same time, the velocity potential on the control surface is numerically computed in Eq. 2.8. The velocity potential on the water surface and the control surface, together with the known moving wall boundary condition and the still wall boundary condition, i.e. the normal derivative of the potential (or wall velocity) at time t = At, will provide the boundary conditions at the water surface for the new application of Green's function on the flow field at an advanced time t = 2At. The simulation is continued as the waves generated by the moving wall grow up. 2.4 Landslide movement For the purpose of simplification, the displacement and velocity of the horizontal-moving wall are assumed to follow a half-sinusoidal curve and satisfy the following equations: 0 v(t)= where a 0 a cosin(ot) 0 2 0 <(t °>:T/2) i7i / 2 ) (2.27) (0<t<T/2) (t*T/2) ( 2 - 2 8 ) is the landslide displacement amplitude, © = 27t/T is the angular wave frequency, and T/2 is the duration of the landslide. By changing the values of a and T, 0 various combinations of displacement and duration of landslides can be simulated. Also, piston-type wavemaker problems can be simulated by implementing a periodic vertical wall movement with a wavemaker amplitude a /2 (stroke of a ) and a wave period T. 0 0 Chapter 2 Theoretical Formulation 2.5 22 Sponge layer damping technique in wavemaker problems Problems of landslide-generated waves are very close to wavemaker problems. Wavemaker problems are studied in detail in this thesis for verification purposes. Numerical simulation of wavemaker problems, often referenced as Numerical Wave Tank (NWT), has been widely studied in recent years. Some key issues in the numerical simulation, such as discretization schemes, boundary element method, time-stepping integration, radiation (open) boundary, absorbing beach or sponge layer technique, etc. have been comprehensively studied and discussed (Kim, 1995 and K i m et al. 1999). The treatment of the open boundary has been one of the great interests. Many N W T studies indicate that the use of the Sommerfeld / Orlanski scheme at the open boundary induces wave reflection. This is mainly due to the numerical error in numerical evaluation of the wave celerity. Apart from numerical Sommerfeld / Orlanski scheme on the radiation condition, various numerical absorbing beach techniques have been proposed and studied. These numerical beach techniques include sponge layer technique (Cointe et al., 1990, Boo and Kim, 1995), active wave absorber technique (Skourup, 1996). In this study, the sponge layer damping technique combined with Sommerfeld / Orlanski scheme are implemented in the wavemaker numerical models. The concept and theoretical development of the sponge layer damping technique are described below. As shown in Figs. 6 and 7, a sponge layer is placed in the longitudinal end as well as the lateral end of the wave tank in order to partially absorb the waves. The linearized dynamic free surface condition on the sponge layer can be written as: ^ at + gn + M .(X> = 0 X <X<X +L 0 0 d (2. 29) d (2. 30) where p.(X) is the damping parameter, expressed as: u.(x) = <xa> — -p d d X <X<X +L 0 J 0 Chapter 2 Theoretical Formulation 23 where a d and Pd are tuning factors, co = 27t/Tis the circular wave frequency, L is the wave length, and L = P L is the sponge layer length denoted as L and L in x and y d d x y directions respectively, and Xo = (X , Yb) is the coordinate of inner boundary of the D sponge layer on the z = d plane. The principle of the sponge layer is to absorb the incident wave energy before it can reach the open-end control surface. As discussed by Kim et al (1999), i f the absorption is too weak, part of the incident energy will reach the control surface and be reflected. Inversely, i f the absorption is too strong, the sponge layer itself will reflect part of the energy. Some N W T studies (Boo and Kim, 1994, 1995, 1996) show that a = 1 and p = 1 in Eq. 2.30 with a sponge layer length of one wave length absorbs reasonably well over a large range of frequencies. In this study, a = 1, p = 1 and L = L are adopted. d Eq. 2.29 can be re-written as follows: |L-gn-u.(x)|> X <X<X +L d (2.31) X <X<X +L d (2.32) 0 0 where p.(X) is simplified as: U(X)=<D 0 J 0 0 This sponge layer damping technique is implemented in both two-dimensional and three-dimensional numerical models of wavemaker problems. The numerical results and the comparison with theoretical solutions are described in the following sections. In simulating the problems of landslide-generated waves, as the numerical error induced by the Sommerfeld / Orlanski radiation condition due to wave reflection at the control surface occurs only when wave crest/trough passes the control surface, thus the error can be avoided by limiting the simulating time so that the leading wave crest stays within the computation domain. Without a sponge layer damping technique, it is possible to maintain an as-large-as-possible wave domain, as computation in three-dimension is extremely expensive on computer resources and very time-consuming. Chapter 3 Numerical Formulation 24 CHAPTER 3 NUMERICAL FORMULATION In Chapter 2, the field solution at a time instant is expressed in term of an integral equation relating the potential and its normal derivative on the boundary. In practical engineering applications, the direct solution to the integral over the boundary is difficult to obtain. Instead, a numerical method can be applied. In the numerical method, the boundary is discretized into a series of facets over which the potential and its normal derivative are assumed to vary according to an interpolation function. The integrals over each facet are then evaluated by a numerical integration in terms of the potential and its normal derivative at a number of nodes where the boundary condition are applied. The application of the discretized equation to all the nodes in turn gives rise to a system of linear algebraic equations from which a numerical solution can be obtained. Then a timedifference method can be used to solve the integral equation involving the time-integral of boundary conditions. 3.1 Field solution 3.1.1 Three-dimensional problem The surface S in Eq. 2.12 over which the integration is performed comprises of the slope surface S , the moving wall surface S , the still water surface S , and the control surface t w 0 S . A l l these surfaces are discretized into finite numbers of planar quadrilateral facets as c shown in Fig. 8. The corresponding values of § and dfy/dnare then taken as constants over each facet and applied at the facet centroid. Eq. 2.12 becomes: J - 1 ASj ASj (i=l,2,...,N) (3. 1) in which i and j are the indices associated with X and X ' respectively; m is the total number of facets in the formulation; and ASj is the area of the j-th facet. Chapter 3 Numerical Formulation 25 As G ( X X j ) is expressed in Eq. 2.15, to solve Eq. 3.1, the relationship between r i 5 and dr/dn needs to be studied. This can be expressed as: where (Xj - X j ) denotes the distance between the points X j and XJ; and rij denotes the normal vector of the boundary at the point X •. Therefore, £ ( dn X , X l i ) - - i & ^ r, k=T " (3.3) k Assume the number of facets on the wall surface, the still water surface, and the control surface is m , m , and m respectively. w 0 It's understood that the integration c performed on the surfaces is in the order of S , S , S and S as the indices i and j increase t w 0 c from 1 to m. Substitute Eqs. 3.2 and 3.3 into Eq. 3.1, re-arrange terms so that unknown values of § and d$/dn are put at the left-hand-side, the following equation can be obtained for rectangular facets: t>,*fc) j = l B,f 1 + j = m „ + l fc)=-£B,ftfc)- On j = On 1 EA, (X;) J + j = m w (3.4) + l in which the coefficients By and Ay respectively correspond to integrals of the Green's function and its normal derivative over area of the j-th facet as in Eq. 3.1. When i = j , the coefficients are evaluated by a closed form integration and when i * j , the integration may be preformed by a numerical integration whereas in this study, the commonly used Gaussian integration is implemented. The coefficients Ay and By are thereby approximated by: £IEdet(j)w w 4 M M m k = l A iJ = m = l n fx' l X - X Vn *™ n = l r (i * j) k m n (3.5) - 2 XZZdet(j)w w„fe=^il^* M M k=2 m = l n = l 4 I+ m r^ (i = j) Chapter 3 Numerical Formulation 4 26 M M EI2>M W mW „n 4 (i*j) m k=l m=l n=l kmn (3.6) M M w 'W in k=2 m=l n=l n (i = j) kmn where m and n are the indices of Gauss points in the k-th image of the j-th facet; M is the number of Gauss points to be used in the numerical integration in one direction; and w m and w are the weight functions at each of the Gauss points. The term J is known as the n Jacobian matrix which relates the global and local coordinate systems of a facet in the numerical integration. The determinant of the Jacobian matrix, det(J), is evaluated at the location of each Gauss points. For a quadrilateral facet, the Jacobian matrix is given by: J= M dp 1^ dp £f4.. tt dq (3.7) ± ^ dq in which Nj is the shape function relating the global and local coordinate system of a facet and is defined as: N,=^(l-pXl-q) (3.8) N = i ( l + p)(l-q) (3.9) 2 N =i(l + pXl + q) (3. 10) N =l(l-pXl + q) (3.11) 3 4 where p and q are local coordinates of a square with the origin at the centroid and a lateral length of unit 2. For a rectangular facet, det(J) is constant and is equal to AS/4. The term T in Eq. 3.6 represents the integration of the Green's function over the i area of the i-th facet when i = j , and is expressed as: Chapter 3 Numerical Formulation I> 27 J-dS (i=l,...,m) (3.12) AS: In principle, a Gaussian integration can consistently be applied to evaluate Eq. 3.12. However, the numerical integration is inaccurate near the singular point r = 0 and a substantial number of Gauss points are required to yield satisfactory results. On the other hand, a closed form solution of the integral for any polygonal facet shape has been given by Hogben & Standing (1974) and is adopted here. For a quadrilateral facet with vertices (p4, q4) in a local coordinate system (p, q) and each making angles 0 i , at (pi, qi), 05 (0i = 0 and 05 = 2TC) in turn with the p-axis at the origin, T is given as: JL, w r= ^¥ k=i in which P = p Pcos(0/2) + Q sin(6/2) - R cos(0/2)| In n8k+ ' (3- 13) P cos(0/2) + Q sin(9/2) + R cos(0/2)| K k+1 -p , Q=q k k+1 -q , W=p k k + 1 q -p q k k and R = (p + Q ) ' . 2 k + 1 2 /2 When this boundary value problem is being solved for Laplace equation, each boundary facet has two variables, potential < > | and its normal derivative d§/dn, and one of these two variables is always given with the other unknown. In our case, d§/dn on the wall surface is known as either 0 or V(t); § on the still water surface and the control surface is known or can be solved by time integration. To solve for (j) on the wall surface and dfy/dn. on the still water surface and control surface, assuming the total facet number is m, m equations can be written in the format of Eq. 3.4 with each corresponding one facet. Thus m variables of and 5<j)/<3n can be solved. It is noted that, with the linear wave assumption, the free surface boundary conditions are applied at the still water surface where the discretization is made. Since the matrix coefficients Ay and By in Eq. 3.4 are dependent only on the facet coordinates, as shown in Eqs. 3.5 and 3.6, Ay and By remain constant throughout the time-stepping procedure. Thus L U factorization can be performed initially on the matrix, and the resulting L U coefficient matrix can be stored and used to solve the linear system, which Chapter 3 Numerical Formulation 28 implies that much less effort is needed in solving the matrix equation in the time-stepping procedure. 3.1.2 Two-dimensional problem The integral equation (2.20) for the two-dimensional problem is solved by a similar numerical procedure in which the boundaries S , S , and S are discretized into finite w c 0 numbers of facets made up of straight-line segments, as shown in Fig. 9. The corresponding values of § and cty/dn are then taken as constants over each facet and applied at the facet center. The numerical discretization gives rise to the same system of equations as in Eq. 3.4 but the coefficients Ay and By are evaluated by different expressions on account of the different integral equation and Green's function that are now applicable. These are now given as: 2 M jk w„ k=l m=l A, 2 M 7 t -k=l X £m=ld e t ( J ) 2 (i*j) 'km ( jkm X w„ (3. 14) - X i)"fijk 0 = j) 'km M ££det(J)w ln(rJ m (3- 15) k=l m=l M m=l (i = j) where m is therindex of the Gauss points in the k-th image of the j-th segment; M is the i+Xdet(J)w mln(r 2m) number of Gauss points used in the numerical integration; and w is the weight function m at each of the Gauss points. The determinant of the Jacobian, det(J), is simply ASj/2 for facets make up of straight-line segments. The term J^ in Eq. 3.15 represents the integration of the Green's function over the length of the i-th segment when i = j and k = 1, and is given by a closed form integration as: Chapter 3 Numerical Formulation ASi/2 r f = Jln|r|dr = A S j In 29 AS, (i = l , . . . , N ) (3- 16) -ASi/2 3.2 Temporal solution As described in Chapter 2, a time integration of boundary conditions is performed on the still water surface and the control surface to provide proper boundary conditions for the application of Green's function on the flow field at an advanced time. The numerical procedures are discussed in details in below. 3.2.1 Free surface boundary conditions The free surface boundary conditions are used in a time-stepping procedure to provide free surface elevation and velocity potential on the still water surface at each time step. Various time-stepping equations are available. One of the procedures, the Adam- Bashforth-Moulton method, has been successfully implemented by Isaacson and Cheung (1991) in their study on time-domain solution for second-order wave diffraction problems. As a simple and reliable procedure, the Adam-Bashforth-Moulton method is adopted in this thesis. The method also has the advantage of providing an optional iteration scheme to account more accurately for the rapid variation of the solution in time. Initially, the free surface elevation at new time (t + At) is first evaluated in terms of the know solution up to time (t - At) by the second-order Adam-Bashforth equation as: r,(t + At) = n(t) + At vdty, V (3. 17) dt y,_ A t where dr\/dt can be obtained by the kinematic free surface boundary condition Eq. 2.9. The dynamic free surface boundary condition, Eq. 2.10 for landslide problems or Eq. 2.31 for wavemaker problems with a sponge layer, can now be used to provide d$/dt at (t + At). In turn the potential <> ( may then be obtained by the first-order Adam-Moulton equation: 4>(t + At) = 4>(t)+ At v dty + (3. 18) vdt7t At + Chapter 3 Numerical Formulation 30 With the normal derivative of the potential on S and S known from the boundary t w conditions on the wall surface Eqs. 2.5 and 2.11, and the potential on S and S evaluated c 0 from the time-stepping equations, the flow in the domain at the advanced time (t + At) can now be solved from the discretized boundary integral equation 3.17. The output vector provides 5<j)/5n on S at (t + At). The kinematic free surface boundary condition, 0 Eq. 2.9, can now be used to provide the corresponding values of dn/St at the advanced time (t + At). This completes an initial calculation of all the boundary parameters at time (t + At) and the computation can then proceed to the next time step. 3.2.2 Radiation condition In order to simulate a sufficiently long duration in a reasonably sized domain, the application of the Sommerfeld radiation condition on the control surface requires particular consideration. When the control surface is located at a sufficiently far distance from the wavemaker, the Sommerfeld radiation condition (Eq. 2.8) can be rewritten in term of the normal derivative of the potential and applied on the control surface. A numerical scheme of Sommerfeld radiation condition proposed by Orlanski (1976) has been widely adopted and is described as followings. In order to develop a numerical solution, Eq. 2.8 is expanded by a central difference in time and a leap-frog difference in space: Hx,t + A t ) ^ (3. 19) where X is a point on the control surface; At is a time step size; An is a small distance of the order of a characteristic facet length; and n is the normal vector at X . From Eq. 3.19, the potential on the control surface at an advance time (t + At) can be expressed in terms of the solutions obtained at previous time-steps: Chapter 3 Numerical Formulation 4>(X, t + At) = 1 - A <|>(X, t - At) + 31 Po <|)(X - nAn, t) (3. 20) in which 'AO v (3.21) Any where p is the Courant number on the control surface which should be maintained at less o than one for stability reasons. Literally, the boundary values of (j) is extrapolated in time and space from the values of § near the boundary inside the domain at previous time steps through the knowledge of the celerity c , which is evaluated on the control surface n at time t. The celerity c on the control surfaces at a time t is approximated by its value at the n previous time step at locations (X - nAn) which are just within the control surfaces. This is done on the basis of: ai>/dt c„ = at (X - nAn) and (t - At) (3. 22) The numerical value of c at (X - nAn) and (t - At) is given by: n c„ = 'An^ v <KX-nAn,t)-<|)(X " "An, t - 2 At) A t y (|>(X - nAn, t) + <(>(X - nAn, t - 2At) - 2<|>(X - 2nAn, t - At)_ (3. 23) In evaluation of celerity c„, potentials within the boundary <|)(X-nAn) and <()(X - 2nAn) are needed. For three-dimensional problems, these can be obtained by the following boundary integral equations after the boundary value problem is solved by Eq. 3.1: •(X) = ^ | G ( X , X 0 £ ( X 0 - * ( X ' ) ^ ( X , X ' ) (3. 24) Chapter 3 Numerical Formulation 32 As all the surfaces are digitized into finite numbers of planar quadrilateral facets and the corresponding values of are taken as constants over each facet and applied at the facet centroid, Eq. 3.24 can be numerically represented as: (3. 25) where Ay and By are defined in Eq. 3.5 and Eq. 3.6 (i * j). Note in Eq. 3.25, boundary values <))(Xj) and ^ " ( ^ j ) a r e either defined by boundary conditions or obtained by solving Eq.3.2. For two-dimensional problems, the boundary integral equation and the numerical representation are as follows, •(x) = - i - J G ( X , X ' ) § ( X ' ) - < K X ' ) § ( X , X 0 dS (3. 26) (3. 27) where Ay and By are defined in Eq. 3.14 and Eq. 3.15 (i * j ) . There is a numerical difficulty in using Eq. 3.23 in that both dfy/dt and d§/dx are zero at locations where r) = 0 and Eq. 3.23 becomes undefined. Even at locations near zero-crossing points the numerical values of c calculated by Eq. 3.23 are inaccurate. In n order to overcome these numerical difficulties, the following numerical scheme was proposed by Isaacson and Cheung (1991) and is adopted. When the absolute values of 8§/dx are less than a certain prescribed value related to the chosen computational precision, c = 0 is used because <j) changes very little in one time step on account of n d§/dt being small. The Courant number on the control surface given by Eq. 3.21 is maintained at less than one for stability reasons. The maximum value of the celerity calculated by Eq. 3.23 is then limited to have a maximum value of A x / A t . Thus, the numerical values of c are given by: n Chapter 3 Numerical Formulation c =0 33 for c' < 0, or when dfdX and dfy/dx are small (3. 28) c = c' for 0 < c' < Ax/At (3.29) c = Ax/At for c'„ > Ax/At (3.30) n n n n n n where c' is the initially calculated value of the celerity from Eq. 3.23. n Although the above numerical scheme was originally developed for hyperbolic flows, its application in parabolic problems has been shown to be effective by Han et al. (1983) and Isaacson and Cheung (1991). Chapter 4 Validation of Numerical Models 34 CHAPTER 4 VALIDATION OF NUMERICAL MODELS Numerical models were setup based on the theoretical development and the numerical formulation described in Chapters 2 and 3. In order to validate the numerical models, numerical simulations were carried out on piston-type wavemaker problems as well as landslide problems. In simulating piston-type wavemaker problems, sinusoidal periodic horizontal wall movement in a rectangular water tank was introduced. A sponge layer damping technique is applied to the water surface in areas adjacent to the control surface to absorb wave energy and reduce reflections. The numerical results in two-dimensional cases are compared with the steady state theoretical solutions, and the numerical results in three-dimensional cases are compared with theoretical solutions as well as experimental results developed and obtained by Takayama (1984). Two-dimensional landslide problems were also simulated using the numerical models and the results were compared with the theoretical solutions derived by Noda (1970). In all the validation cases, numerical results agree well with theoretical/experimental results. These are described and discussed in the following sections. 4.1 Two-dimensional wavemaker problem Two-dimensional piston-type wavemaker problem was studied. The theoretical solution presented by Dean (1985), the numerical results in testing runs using the numerical models developed in this thesis, and their comparisons are described and discussed in the following sections. 4.1.1 Theoretical solution The theoretical solution for a piston-type two-dimensional wavemaker problem was obtained by Havelock in solving a boundary value problem. The definition sketch of the problem is the same as shown in Fig. 7 in two-dimensional scenario. This solution was presented by Dean (1985) and is summarized below. The wave elevation generated by a two-dimensional piston-type wavemaker far from the wave paddle is given by: Chapter 4 Validation of Numerical Models 35 a (cosh2kd-l) t| = — — -cos(kx-cot) sinh2kd + 2kd ' n A x (4- ) v where n is wave elevation; a i s wave paddle displacement magnitude (stroke); 0 co = 2TC / T is circular frequency; T is wave period; k = 2TC / L is wave number; L is wave length; and d is water depth. The relative wave height, defined as the ratio of wave height and wave paddle displacement magnitude H / a , is given by: 0 H _.2(cosh2kd-l) a 0 sinh2kd + 2kd 4.1.2 Numerical results and comparisons The numerical formulation described in Chapter 3 was coded and utilized to simulate piston-type wavemaker problems by applying a periodic movement to the wave paddle. Similarly, the wave paddle boundary, the free surface boundary including a sponge layer, and the control surface boundary were digitized into small segments, as shown in Fig. 9. The periodic wave paddle displacement a(t) and velocity v(t) are governed by the following equations: a(t) = ^[l-cos((ot)] (4.3) v(t) = -^-cosin(cot) (4.4) where co = 2TC/T is the circular wave frequency. Eq. 4.4 provides the boundary condition for the wave paddle. Numerical simulations on two-dimensional wavemaker problems were carried out and one set of results are presented in Figs. 10 and 11. In this testing case, the dimension of the tank was 40 m in length and 1 m in depth. The wave paddle with a motion magnitude of 0.02 m was moving at a period of 2.3 seconds. The sponge layer length was 6.28 m at the downstream end of the tank. The wave paddle surface boundary and the control surface boundary were discretized into 10 segments, and the still water surface was discretized into 200 segments. Wave elevations were computed at a 0.01 second time increment. Figure 10 shows wave profiles along the wave tank at varying Chapter 4 Validation of Numerical Models 36 times, while Fig. 11 shows wave elevation changes as a function of time at three different fixed points. It is seen that the wave elevations obtained are stable. More numerical simulations were carried out for varying wave periods and the results were compared with the steady state theoretical solutions presented by Dean (1984). In all the testing cases, the wave tank depth was 1 m, and wave paddle motion magnitude was 0.02 m. The wave periods, the computed wave lengths and wave elevations, dimensionless wave frequency kd and the dimensionless wave height H/ao of these cases are listed in Table 1. Comparisons between these numerical results and the theoretical solutions on relative wave height H/ao versus dimensionless wave frequency kd are shown in Fig. 12. Very good agreement was obtained between the numerical results and theoretical solutions. 4.2 Two-dimensional landslide problem To further validate the two-dimensional numerical model, test cases were carried out using the two-dimensional numerical model to simulate a horizontal wall movement with constant speed, and the results were compared with the theoretical solution derived by Noda (1970). These are described and discussed in the following sections. 4.2 1 Theoretical solution The theoretical solution for a two-dimensional horizontal landslide problem was derived by Noda (1970) using a direct-integral method by treating landslide as a horizontal moving wall, the velocity of which was approximated by N straight lines as governed by the following expression: V 0 < t' < T, T, < t' < T 2 (4. 5) T A N-1 < t' < TN — 1 — t' > T 1 ] Chapter 4 Validation of Numerical Models 37 in which, v'(t') is dimensionless wall speed; V i , V2, ... V N are dimensionless constant speeds; t' is dimensionless time defined as t' = t^/g/d . The derived wave elevation by a direct-integral solution is expressed as: •n(x't') J k T - i d ^ [2 r « V tanhucosux' . n = Sl-Jo ^[TI r /, „ \ . / , ^ \i. ] [sina(t'-T J-sina(t'-T )]du — n (4.6) n au J O „ J in which x' = x / d , u and a are dimensionless wave number and dimensionless wave frequency, expressed as: a = Vgktanh(kd)d7g ^ / v u = kd ( 4 8 ) In Eqs. 4.7 and 4.8, k = 2u/L is wave number, L is wave length. Considering when u > 9, tanh u « 1, thus a ri(x',t') _ Z J I T d ^[2 , Eq. 4.6 can be further re-written as: r V tanhucosux' . 5111 / , t ^T = Z - Jo — " + \ 9 n r [71 au 2 r» V„ tanhucosux' . 8111 /, CT I °( - "-. ) - _ TY I , t - n ft " 1 J 0 -nJ -. u 3 / 2 (4. 9) sinu>(t'-T _,)-sinu2(t'-T ) du n n The first integral in Eq. 4.9 can be solved numerically, while the second integral can be integrated by parts when x' > 0. Note this analytical solution was based on a linearizing assumption that the displacement of the wall is small in comparison to the water depth, so that the integrals can be performed on the still wall surface instead of the moving wall surface. This assumption is the same assumption utilized to derive the numerical solutions in this thesis. In a strict sense the derived analytical solution are for a water inflow problem at x = 0, where v(t') is the velocity of water inflow. The analytical solution, Eq. 4.9, was utilized to validate the two-dimensional numerical model. The theoretical solutions, numerical results and the comparisons are described and discussed in the following section. Chapter 4 Validation of Numerical Models 38 4.2.2 Numerical results and comparisons In the numerical test case, the wave tank was 1 m deep and 20 m long; the water surface was digitized into 2,400 segments; the time step was 0.01 second and the simulation duration is 3.2 seconds; the wall was moving at a constant speed of 3.13 m/s (Vi = 10). Two probes were placed atx = 2 m ( x ' = 2), and x = 5 m ( x ' = 5)to record time history of wave elevation. The recorded dimensionless wave elevations as a function of dimensionless time in comparison with the results from the theoretical solution by Noda (1970) are presented in Fig. 13. It is seen that very good agreement was obtained. In summary, the validation test runs on wavemaker problems and landslide problems in two-dimensional cases successfully validated the two-dimensional numerical model. 4.3 Three-dimensional wavemaker problem A three-dimensional piston-type wavemaker problem was studied. The theoretical solution and the experimental work obtained by Takayama (1984), the numerical results in testing runs using the numerical models developed in this thesis, and their comparison are described and discussed in the following sections. 4.3.1 Theoretical solution The theoretical solution for a three-dimensional piston-type wavemaker problem was derived by Takayama (1984). The definition sketch of the wavemaker problem in Takayama's derivation is the same as shown in Figs. 6 and 7, excluding the sponge layers. The derived elevation of waves generated by a piston-type single wave paddle are given by: fa l J_^ No(V(kx-q) +(ky) )dq 2 p + 2 2 2 (4. 10) Chapter 4 Validation of Numerical Models 39 where J (x), N (x), and K (x) are Bessel, Neumann and modified Bessel functions with 0 0 0 index zero, respectively; values of v are the positive roots of the following equation: a> = - g v tanh(vd) 2 (4.11) ct and f3 are defined by the following equations: p p sinh kd a p = (4. 12) r sinh 2kd^| kd 1 + 2kd J sinh vd 2 (4. 13) sinh 2vd vd 1 + 2vd The relative wave height as a function of location in the x-y plane, defined as the ratio of wave height and wave paddle displacement magnitude H/ao, is given by: + ZPpJ*^(V(vx-q) (v))d' 2 0 + ! y q (4. 14) 1/2 +a. ^ J (V(kx-q) (ky) ]dq 2 2 2 0 2 + The theoretical solution presented in Eq. 4.14, together with some experimental results by Takayama (1984), which were used to validate the three-dimensional numerical model, are described in the following section. 4.3.2 E x p e r i m e n t a l W o r k The experimental work on waves induced by a single wave paddle was carried out in a wave basin and the results were presented by Takayama (1984). The experimental setup by Takayama is illustrated in Fig. 14(a). The experimental work was carried out in a wave basin of 25 m long, 15 m wide, and 1 m deep. Rubble mound and wave absorbers were placed along the tank walls to reduce wave reflection. Three 5-metre long piston type wave paddles were placed on one side of the tank. Only the middle wave paddle Chapter 4 Validation of Numerical Models 40 was operated in the experiments. Wave elevations were measured at 140 points in the wave tank, as shown in Fig. 14(b). The water depth was 0.63 m for all experiments. Two experiments corresponding to a wave period of 1.0 second and 1.4 seconds were carried out and the results were presented by Takayama. These experiments were numerically simulated in two test cases and comparisons are presented and discussed in the following section. 4.3.3 Numerical results and comparisons In the three-dimensional numerical model, the wavemaker was simulated as a vertical wave paddle moving periodically in a rectangular wave basin (see Fig. 7). The wall surface, the wave paddle surface and the still water surfaces including the sponge layer were digitized into a number of facets (see Fig. 8). The wave paddle's displacement and velocity are governed by Eqs. 4.3 and 4.4 respectively. Numerical results of four test cases are presented in this thesis, corresponding to the experiments performed by Takayama (1984). In all these test cases, water depth was 0.63 m, wave paddle motion magnitude was 0.02 m, duration of simulation was 10 seconds, the water surface was digitized into a 48 x 48 mesh, water depth was digitized into 4 segments, and the total number of facets was 2,880. Wave basin dimension, wave paddle length, and wave period for each test case are listed in Table 2. Numerical simulation test cases A and B correspond to theoretical results presented by Takayama, with relative wave paddle length b/L = 0.5 and 1.0 respectively (L is the wave length). Wave profiles generated in test case A are shown in Figs. 15 and 16. Figure 15 shows the three-dimensional wave profile as a function of location at varying time instances. Figure 16 shows the simulated water surfaces as a function of time at six fixed points where numerical probes were located (locations of the numerical probes are shown in Fig. 14(a)). It is seen from Figs. 15 and 16 that the generated waves are stable. The relative wave height H/ao from the numerical simulation test case A , in comparison with that of the theoretical solutions derived by Takayama (1984), are presented in Fig. 17. The relative wave elevation from the numerical model test case B , in comparison with that of the theoretical solutions and the results, are shown in Fig. 18. Note that wave Chapter 4 Validation of Numerical Models 41 heights at only half of the wave basin are presented due to symmetry. It is seen from Figs. 17 and 18 that for both b/L = 0.5 and b/L = 1.0, the numerical simulations generally agree with the theoretical solution. The numerical test cases C and D correspond to the experimental work performed by Takayama. From test case C (T = 1.0 second), the simulated relative wave height H/ao as function of relative distance x/L from the wave paddle at varying y values, in comparison with the theoretical and experimental results, are shown in Fig. 19. Similarly, the simulated relative wave heights H/ao as a function of relative distance y/L from the centerline of the wave basin at varying x values, in comparison with the theoretical and experimental results, are shown in Fig. 20. From test case D (T = 1.4 second), the simulated relative wave heights H/ao as a function of relative distance y/L from the wave paddle at varying x values, in comparison with the theoretical and experimental results, are shown in Fig. 21. And the simulated relative wave heights H/ao as a function of relative distance x/L from the centerline of the tank at varying y values, in comparison with the theoretical and experimental results, are shown in Fig. 22. It is seen from Figs. 1 9 - 2 2 that the numerical results agree with the theoretical and experimental results. 4.4 Three-dimensional landslide problem As described in Chapter 1, the studies on problems of the three-dimensional landslidegenerated waves are very limited at the present time. Thus the validation of the numerical model on the three-dimensional cases was carried out mostly by visual observations. The numerical modeling results, however, were compared with the visual observation and the field measurement in a case study, which are described and discussed in Chapter 6. Chapter 5 Results and Discussion 42 C H A P T E R 5 R E S U L T S A N D DISCUSSION Further to the validation of the numerical models, some test runs were performed to study the computational characteristics of the models. Dimensional analysis and a number of numerical simulations on two-dimensional and three-dimensional landslide-generated waves were then carried out and a selection of results is presented. These are discussed in the following sections. 5.1 Computational considerations The computational considerations include discretization schemes, computing cost, accuracy and stability of matrix solution, stability of time-stepping procedure, and numerical accuracy of discretization. 5.1.1 Discretization schemes The three-dimensional numerical model corresponds to that of a wave basin with a part of moving boundary located at one end. The moving wall surface, the still wall surface, the still water surface and the control surface are discretized into finite number of planar quadrilateral facets (see Fig. 8). The wave basin bottom is excluded due to the choice of the Green's function that accounts for the symmetry. A F O R T R A N code L G W 3 D has been developed to perform the surface discretization, generate the matrix coefficients, solve the matrix equation and time-step the development of waves. A Gaussian integration is used to compute the matrix coefficients although other numerical integration methods may be feasible. For the purpose of high order of accuracy, a 64point Gaussian integration and double precision is used throughout the computation. Typically 2,000 to 3,000 facets are used in the discretization of the surfaces. The two-dimensional numerical model corresponds to that of a narrow wave flume with a moving boundary on one end. The moving wall surface, the control surface and the still water surface are discretized into numbers of straight-line segments (see Fig. 9). The tank bottom surface is excluded due to the choice of the Green's function that accounts for the symmetry. A F O R T R A N code LGW2D has been developed to perform Chapter 5 Results and Discussion 43 the surface discretization, generate the matrix coefficients, solve the matrix equation and time-step the development of waves. A n 8-point Gaussian integration and double precision is used throughout the computation. About 200 to 400 facets are generally used in the discretization of the surfaces. Both L G W 3 D and LGW2D models can be used to either simulate landslidegenerated waves with a radiation condition or piston-type wavemaker problems with the radiation condition in combination with sponge layer technique by changing a control parameter in the input file. 5.1.2 Computing cost When the rank of a matrix equation is large, computing cost might pose an important issue affecting the applicability of the method. Typically, with the rank of the matrix equation denoted by N , the number of floating-point operations for the generation of matrix coefficients is proportional to N , while the number of floating-point operations for solving the matrix equation is proportional to N . For LGW2D and L G W 3 D models, since the elements of the left-hand side matrix are functions only of facet locations, orientations and relative distances between the facets, the matrix remains constant throughout the time-stepping procedure. Thus L U factorization can be performed initially on the matrix, and the resulting L U coefficient matrix can be stored and used to solve the linear system, which can significantly speed up the procedure of solving the matrix equation. The programs LGW2D and LGW3D were coded and compiled using the Microsoft Fortran Power Station (Version 4.0) and can be run on personal computers running Microsoft Windows operating systems. As the computing cost depends on a number of factors such as the C P U speed, the R A M size and speed, and sometimes the hard drive speed when temporary storage on the hard drive is needed, the C P U time for a specific problem may vary slightly from time to time depending on the computing environment. Thus, the results presented in this section can only be used for reference purposes. Figures 23 and 24 show the C P U time as a function of the rank of the matrix equation for the program L G W 2 D and L G W 3 D running on a P C equipped with Intel PIII-600 C P U and 256 Mb R A M . Note the C P U time does not include the time spent on Chapter 5 Results and Discussion 44 the generation of matrix coefficients and the L U factorization, since this time duration occurred only once throughout the simulation period. It is noted that for both twodimensional and three-dimensional problems, the C P U time is proportional to approximately N ' , where N is the total number of segments used. 5.1.3 Accuracy of matrix solution When the rank of a system of simultaneous equations is large, the round-off error in a matrix solution requires particular attention before it can be applied more generally to the stepping procedure. A number of different techniques can be used to indicate the magnitude of the round-off error in a matrix solution. One of the most frequently used is the condition number which has the property that the higher the condition number the less accurate is the matrix solution. In this thesis, condition numbers of matrices were evaluated using routine D L F C R G in the IMSL Math Library included in the Microsoft Fortran Power Station (version 4.0) package. Routine D L F C R G performs an L U factorization of a real general coefficient matrix and to estimate the condition number of the matrix. The solution of the linear system is then found using the iterative refinement routine DLFSRG. Routine D L F C R G fails i f U , the upper triangular part of the factorization, has a zero diagonal element or i f the iterative refinement algorithm fails to converge. These errors occur only i f the matrix is singular or very close to a singular matrix. If the estimated condition number is greater than 1/s (where 8 is machine precision), a warning message is issued, which indicates that very small changes in the matrix can cause very large changes in the solution, thus the matrix is ill conditioned. If no warning is issued, the solution is considered accurate to at least the machine precision. In this thesis, warning messages from routine D L F C R G for cases of both twodimensional and three-dimensional model simulations were monitored. For all the test cases, no warning message has ever been encountered. Thus the matrix solutions for all the test cases presented in this thesis are considered to be accurate. It is seen from Eqs. 3.4, 3.5, 3.6, 3.14, and 3.15, that the left-hand side matrix in the matrix equation consists of elements that are functions of only facet locations, Chapter 5 Results and Discussion 45 orientations and relative distances between the facets. Thus the geometry and discretization scheme of a specific problem domain are dominant factors for the accuracy of matrix solutions. Since the geometry and discretization scheme can be chosen in any form theoretically, there is no guarantee that the formed matrix of a specific problem is well conditioned. Thus condition numbers and error messages from routine D L F C R G should always be monitored in any numerical simulations to guarantee the accuracy of a matrix solution. 5.1.4 Stability of time-stepping procedure The stability of the proposed solution requires that the time-step size At is sufficiently small in relation to a characteristic facet diameter AS. This condition corresponds to a minimum value of T/At, for any given value of L/AS, where T and L denote a representative wave period and length. Isaacson and Cheung (1991) studied the stability of the time-stepping procedure, which is adopted in this thesis. It was found that in two-dimensional cases, stable solutions could be obtained for L/AS and T/At as small as 10 and 7 respectively. Isaacson and Cheung (1991) also studied the Courant stability condition which specifies that the Courant number given by cAt/AS = (L/AS)/(T/At) is less than one, and found that this appeared to apply to the cases when L/AS is approximately less than 30. For larger values of L/AS, a stable solution could be obtained for a time-step size which was significant larger than that required by the Courant stability condition. In the context of this thesis, the Courant stability condition has been enforced for all numerical simulations of wavemaker problems. In cases of problems of landslide- generated waves, the wave length of each wake wave following the leading wave is smaller than its precedent wake wave. Because of the decreasing wave length, the Courant number cannot be defined explicitly. Nevertheless, in engineering applications, the first several leading waves generally have the largest impact on coastal structures and are the focus of this thesis. The Courant stability condition provides a good guideline in the numerical simulations. In this thesis, time step is always kept small and visual observations have been performed in each test case to ensure satisfying stability. Chapter 5 Results and Discussion ; 46 5.1.5 Numerical accuracy of discretization Apart from the issue of stability, the numerical accuracy of discretization must also be assessed in terms of the chosen segment or facet size. The computed wave elevations were found to depend on the number of segments or facets used to represent the water surface. To illustrate the effects of the number of segments or facets used to represent the water surface, convergence of free surface wave elevations is examined for wavemaker problems. Two-dimensional problem Let r) denote the exact wave elevation, and n ^ denote the numerically simulated e wave elevation on a discretization scheme with a spacing of Ax . The error Ej in the numerical solution at point i is: E i.Ax = ^ l A x ( i ) - l e ( x ) x T i (5.1) To summarize the error as a single number, a commonly used norm, L i norm, defined by the following equation (5.2) is used: n IN Ml n Suppose numerical solutions on three discretization schemes D i , D2, and D3 are available, where D2 has triple as many segments as D i , and D3 has triple as many segments as D2. The use of the tripling schemes allows the differences among the schemes on specific locations to be computed without additional approximation error, as wave elevations are represented by the values at the center of each segment. Assuming that the error in each solution is proportional to its segment spacing to some power k, the solutions can then be written as: n| Di =n +CAx (5.3) k e 'Ax le+Cj^-l N k (5.4) Chapter 5 Results and Discussion f V = T le+C 47 AxV (5.5) where C is a constant parameter in the above equations. Taking the norm of the difference of the solutions between schemes D i and D2, and that between D2 and D 3 , we have: > < -< < -< = CAx J_ l 3 3k (5.6) j (5.7) gk The ratio of Eq. 5.7 and Eq. 5.6 can be expressed as: < < -< 3 (5.8) k Eq. 5.8 indicates that the convergence order k can be evaluated without knowing the exact solution. Moreover, the norm error of the finest scheme D3, | | E II = C A x X-, can k n II 3 II (jk ' be evaluated by taking the ratio of this norm to that in Eq. 5.7, i.e. E o 1 3 < -< (5.9) 3 -l k Eq. 5.9 can then be used to compute the norm of scheme D 3 after the norm and the convergence order k are obtained. To further reduce the numerical error caused by the time-stepping procedure, the time step for the finer mesh is reduced to 1/3 of the time step for the corresponding 2 coarser scheme, to ensure that the order of the accumulation of this numerical error is no larger than the order of errors caused by the discretization schemes. Chapter 5 Results and Discussion 48 The above methodology is applied in a series of numerical simulations to study the accuracy of discretization schemes of the numerical models. These simulations correspond to a wave tank of 1.0 m water depth, 30 m in length with a 5 m long sponge layer at the downstream end of the tank. The wave paddle is moving at a period of 2.3 seconds, with magnitude of 0.02 m. The water surface is discretized into 60, 180, and 540 segments respectively. The time step for each discretization scheme is set at 0.081, 0.009, and 0.001 second, respectively. The wave elevations corresponding to these discretization schemes at time t = 8.1 second are shown in Fig. 25. The norms of differences in wave elevations between the coarser and finer discretization schemes, the convergence. order k, and the numerical error of the finest discretization scheme at varying time instances are shown in Table 3. It is seen from the results that the estimated order of convergence varies from 1.79 to 2.03, which indicates that the numerical scheme converges nicely. The error norm of the finest scheme is in the order of 10" m to 10" m. 6 5 As the wave height generated in the wave tank is in the order of 10" m, the error caused by discretization schemes is negligible. Three-dimensional problem The accuracy of discretization was studied following the similar method described in the previous section. Again, assuming three discretization schemes, D i , D2 and D3 are available, where D2 is obtained by dividing each facet of the D i scheme into nine equal facets, and D3 is obtained by dividing each facet of the D scheme into nine equal facets. 2 This facet refining method allows the accurate comparison among different discretization schemes be carried out at fixed specified locations, since wave elevations are represented by values at the centroid of facets. Thus, on the water surface, the wave elevation from numerical solutions of three discretization schemes can be expressed as: ru+cCAxAy) 1 AxAy tle+C 9 (5.11) J AxAy Tle+C V (5. 10) 81 J (5. 12) Chapter 5 Results and Discussion 49 Following the similar method described in the previous section for two-dimensional problems, the ratio of the norm of the differences in the solution between the finer schemes and that between the coarser schemes can be obtained as: < -< < -< (5. 13) and the norm error of the finest scheme D E "ID 3 j|E 1| = C(AxAy) D3 k , can be evaluated as: 1 3 -Tl D3, L 9 k -1 (5. 14) Similar to the two-dimensional case, the time step for the finer mesh is reduced to 1/9 (first order) of the time step for the corresponding coarser scheme, to further reduce the numerical error caused by the time-stepping procedure. Strictly speaking, to ensure the order of the accumulation of this numerical error is no larger than the order of errors caused by discretization schemes, the time step for the finer mesh should be reduced to 1/9 (second order) of the time step for the corresponding coarser scheme. However, 2 such a time step reduction causes significant increase in simulation time and makes the numerical simulation for the finest mesh infeasible. Nevertheless, the first order time step reduction is considered sufficient for the present purpose, and thus is adopted. To study the accuracy of the three-dimensional numerical model, sample test cases on the three-dimensional wavemaker problem were simulated. The test cases correspond to a wave basin of 1.54 m long, 3.08 m wide, with a water depth of 0.63 m. A 0.77 m wide sponge layer was implemented at the open ends of the tank. Wave paddle motion magnitude was 0.02 m, and the period was 1.0 second. The water surface (1.54 mx 1.54 m) was discretized into 6x6, 18x18, and 54x54 facet meshes for the three schemes respectively. The time step in the simulations for each discretization scheme is 0.081, 0.009, and 0.001 second, respectively. The simulation duration was controlled to be within 1.2 second to make sure the numerical results were not affected by the reflection from the open boundary and the sponge layer as these boundary conditions could not Chapter 5 Results and Discussion 50 provide perfect wave passing, as discussed in Chapter 3. The wave elevations at the centerline of the wave basin (perpendicular to the wave paddle) at time t = 0.6 second are shown in Fig. 26. The norms of differences in wave elevations between the coarser and finer discretization schemes, the convergence order k, and the numerical error of the finest discretization scheme at varying time instances are shown in Table 4. It is seen that the estimated order of convergence varies from 0.85 to 0.98, which indicate the numerical scheme converges nicely. The error norm of the finest scheme is in the order of 10" m. This compares with the approximate 10" m generated wave height, and indicates that the error caused by discretization schemes is negligible. Thus the threedimensional numerical model converges nicely as the number of facets representing the water surface increases. This is consistent with the accuracy study results on the twodimensional numerical models. In summary, the studies on the discretization schemes of both two-dimensional and three-dimensional numerical models indicate that sufficient accuracy can be obtained from the discretization schemes. 5.2 Dimensional analysis The factors that affect the amplitude of landslide-generated waves at a certain location include the water depth, the acceleration of gravity, the distance from the landslide source, and the landslide motion such as volume, width, displacement, speed and duration. Among these factors, the landslide volume Q , width b, displacement a(t), speed V(t) and duration T/2 are correlated. Assuming that a simplified landslide is treated as a vertical wall moving horizontally following a half sinusoidal motion, the slide displacement can be defined by: (0 < t < T/2) (t > T/2) Thus the landslide speed V(t), and the landslide volume can be expressed as: (5. 15) Chapter 5 Results and Discussion v(t) = a 7T 0 '2TI ^ sin — t T j 0 51 (0 < t < T/2) (5. 16) (t > T/2) Q = a bd (5. 17) 0 It is seen from Eqs. 5.15 - 5.17, the landslide motion can be represented by the landslide displacement a , the width b, and the duration T/2, although other parameters 0 such as volume and speed are also representative. Therefore, the amplitude of landslidegenerated waves can be identified as a function of location, water depth d, landslide duration T/2, landslide width b, landslide displacement amplitude a , and acceleration of 0 gravity g. Hence the wave height generated by the simplified landslide can be described as: H = f(x,d,a ,T,g) two-dimensional landslide (5. 18) H = f(x,y,d,a ,b,T,g) three-dimensional landslide (5. 19) 0 0 where H is the wave height generated by landslide, and can be obtained by means of evaluating the free surface elevation. As wave height is linearly proportional to landslide amplitude, by choosing water depth d, acceleration of gravity g, and double duration T as basic parameters, dimensional analysis gives the following relationship: x dgT a„ H d =f 2 x y b d two-dimensional landslide (5. 20) three-dimensional landslide (5.21) d'd'd' r g For different categories of landslide, Eqs. 5.15 - 5.17 can be slightly modified to fit certain types of landslide. 5.3 Two-dimensional landslide-generated waves A number of numerical simulations were performed to study the characteristics of landslide-generated waves in two-dimensional cases. The focus of the study is the Chapter 5 Results and Discussion 52 maximum wave amplitude generated by landslide, which has the most significant impact on coastal structures. The results from the numerical simulations indicate that the highest wave in the landslide-generated wave train is always the leading wave. Also, a number of smaller "wake waves" are generated following the leading wave, and the wave height decreases as the waves propagate. Figs. 27 and 28 demonstrate these phenomena. In this case, the simulated water tank was 40 m long, and the water depth was 0.63 m. The landslide displacement was 0.01 m, and the landslide duration was 0.9 seconds. The water surface was modeled using 200 straight-line segments. Time step in the simulation was 0.01 second. It is noted that the wave periods and the wave lengths of the landslide-generated waves tend to increase as they propagate. Thus the dispersion relation in the linear wave theory cannot be used to evaluate the wave speed. However, wave speed can be numerically determined by monitoring the crest positions of the waves, although this was not attempted in this study. As discussed in the previous section, in two-dimensional cases, the dimensionless wave height H / a of landslide-generated waves depends on the relative location x/d, and 0 dimensionless landslide duration d/gT (Eq. 5.20), where T/2 denotes landslide duration. 2 To study this dependency, numerical simulations were carried out based on varying landslide duration, and the resulting wave heights at various locations were analyzed. Figure 29 shows the maximum dimensionless wave height as a function of the relative location x/d for various d/gT values. It is noted that with the same landslide amplitude, 2 the shorter the landslide duration, the faster the landslide is moving. It is seen in Fig. 29 that landslides with higher speed (or shorter landslide duration) generate higher waves, especially in the vicinity of landslides. Also, landslide-generated wave height decreases as the distance from the landslide increases. Far away from the landslide, for the d/gT 2 values within the range of 0.01 to 0.06, wave height converges to a constant value of approximately 0.4a regardless of landslide duration or speed. In the extreme case of 0 very low d/gT 2 values (i.e. very slow slide moment), a zero wave height is expected in the far field. Therefore, a convergence of 0.4a is valid only for d/gT 0 2 above a certain Chapter 5 Results and Discussion 53 value. And indeed, as expected, the lowest curve in Fig. 29 shows that the wave height for d/gT =0.005 is smaller than those for larger d/gT 2 5.4 2 values. Three-dimensional landslide-generated waves A number of simulations were carried out to study the characteristics of landslidegenerated waves in three-dimensional cases. Typical generation and propagation of landslide-generated waves in three-dimensional cases are demonstrated in Fig. 30. This simulation was carried out in a half wave basin of 5 m long, 5 m wide, and 0.63 m deep in water depth. The landslide was treated as a 0.63 m wide wall block moving into the wave basin with a displacement of 0.01 m and a duration of 0.5 second. The water surface of the half wave basin was digitized into 36x36 facet mesh, and the simulation time step was 0.02 second. From the dimensional analysis, the dimensionless wave height H / a 0 is a function of relative location (x/d, y/d), relative landslide width b/d, and dimensionless landslide duration d/gT , as defined in Eq. 5.21. To study the effect of these factors, numerical 2 simulations were carried out for three landslide widths (b/d = 1.0, 2.0 and 4.0) and three landslide durations (d/gT = 0.02, 0.04, 0.06) respectively. 2 A l l simulations were performed in a half wave basin of 5 m long, 5 m wide and 0.63 m deep in water depth. The water surface of the half wave basin was digitized into a 36x36 facet mesh, and the simulation time step was 0.02 second. The results of the dimensionless wave height for varying landslide dimensions and durations are shown in Figs. 31 - 33. The results of the dimensionless wave height along the centerline of the landslide for these cases are shown in Fig. 34. It is seen from Figs. 3 1 - 3 4 that for a fixed landslide width, the shorter the landslide duration or the higher the landslide speed is (higher d/gT the waves are generated. 2 values), the higher For b/d = 1.0, in the vicinity of landslide, the wave height changes from approximately l.Od to 0.4d as d/gT value decreases from 0.06 to 0.02. A 2 similar trend is noted for cases with b/d = 2.0 and 4.0. Chapter 5 Results and Discussion 54 When the landslide duration is relatively short or the landslide speed is relatively high ( d / g T = 0.06 and 0.04), the waves along the centerline of the wave tank are higher 2 than the waves at other locations. On the other hand, when the landslide duration is relatively long or the landslide speed is relatively low, the waves propagate and scatter almost uniformly in the wave basin. It is also seen from the figures that the waves decay fairly rapidly as they propagate in the wave basin. This agrees with the fact that, in open water, the amplitudes of the scattered waves decay like l / Vr due to the directional spread of wave energy, where r is radian distance measured from a point of interest (e.g. Sarpkaya and Isaacson, 1981). The effect of the landslide dimension (in this case, the width) can be seen by comparisons of Figs. 3 1 - 3 4 for fixed d/gT values. As expected, waves generated by 2 larger landslides (b/d = 4.0) are greater than waves generated by smaller landslides (b/d = 1.0). 5.5 Limitations of the numerical models The numerical models have been validated and the numerical simulations show some results of interest. However, in application of the theoretical formulation and utilizing the numerical models to study landslide-generated waves, some limitations need to be considered. These limitations are as follows: • Constant water depth The numerical models are set up based on the assumption of constant water depth. This assumption provides symmetry with respect to the domain bottom, which allows for a lower number of discretized facets on the surfaces and enhances the efficiency of the numerical model. However, it is noted that the numerical method is in fact not restricted by the assumption of constant water depth. When the water depth is not constant, the bottom boundary has to be included in the discretization scheme. The Green's function has to be revised accordingly, and a larger number of facets will be needed in the numerical model. In any event, it is emphasized that the present Chapter 5 Results and Discussion 55 problem is intended to relate to the more immediate vicinity of the landslide, and not to account for diffraction, refraction and shoaling so as to predict wave propagation away from the landslide site. • Linear waves The current numerical method and numerical models are based on the assumption that landslide-generated wave amplitudes are relatively small compared to water depth or wave length. This assumption may no longer be valid for more extreme events when the speed of the landslide is so large that breaking waves are generated right in front of the landslide. The current numerical methods and models do not take into account breaking wave effects. Thus, the breaking wave limit, which is normally assessed by breaking water depth or wave steepness (Sarpkaya and Isaacson, 1981), should always be verified during the numerical simulations. Attention also needs to be paid to the cases where nonlinear wave effects may be significant. To assess this in the context of shallow water waves, reference may be made to the Ursell parameter, defined by U = H L / d . 2 r 3 The Ursell parameter is of order unity or is moderately large for high waves in shallow water, or when nonlinear wave effect becomes significant (Sarpkaya and Isaacson, 1981). For small values of Ursell parameter, linear wave theory is generally applicable. Nevertheless, it is worthwhile pointing out that linear wave theory has been demonstrated to be robust and sufficient in many engineering applications, even when waves are moderately steep, and has been adopted as the fundamental wave theory in many coastal engineering design and reference manuals, such as the Shore Protection Manual (1984). • Horizontal wall movement While the present numerical method is applicable to a landslide of arbitrary profile, the numerical results presented in this thesis are based on the further assumption that the landslide can be simplified into a vertical wall and the motion follows a half sinusoidal curve. For a slump-type landslide with a dominant horizontal movement, it is expected that the slide may be simplified into a horizontal moving wall by equalizing the slide volume with the wall motion amplitude and duration, so that the Chapter 5 Results and Discussion 56 present numerical method may be applicable. A case study on Skagway Harbour landslide described in Chapter 6 indicates that such a simplification appears feasible. Indeed, the validity of this assumption may be examined numerically using the present models, although this is not attempted in this study. Other types of landslide, such as fockfall, submarine landslide, and vertical seabed movement, can also be simulated with slight modification to the numerical models. More realistic landslide motion can also be simulated provided the time-history of the slide displacement is available. • Geometry of boundaries Theoretically the numerical methods can be applied to a water domain of arbitrary boundaries. However, attention has to be paid on the overall geometry, because the accuracy of matrix solution is dependent on the condition number of the left-handside matrix, which is function of boundary surface facet coordinates. It is possible that in some particular situations the boundary geometry gives rise to an i l l conditioned left-hand-matrix which causes non-negligible errors in the matrix solution. Warning messages from the numerical models will be issued i f such cases occur. It is suggested that warning messages from the numerical models be monitored to guarantee the accuracy of the simulation results. • Discretization scheme In this thesis, straight line segments and rectangular facets are used to digitize the boundaries and the water surface of the problem domain. The velocity potential at the center of each line segment or at the centroid of each rectangular facet is obtained through a Gaussian numerical integration. Water surface elevations are also represented by values at these centroid locations. Thus, to obtain a reasonably accurate wave profile, it is essential to ensure that the digitized water surface mesh is fine enough to represent waves. In most cases, this conflicts with the efficiency of the numerical model and the computer hardware capacity and needs to be resolved with caution. Other discretization schemes, such as higher order element representation, are applicable with some modification to the numerical models. Chapter 5 Results and Discussion 57 Time step and stability In using the numerical models to simulate landslide-generated waves, time step needs to be sufficiently small to ensure the stability of the time-stepping procedure. As described in Section 5.1.4, the Courant stability condition provides a good guideline in the numerical simulations. • Wave reflection at open boundary It has been noted that the numerical implementation of the Sommerfeld / Orlanski radiation condition at the open boundary may cause partial wave reflection when wave trough or crest passes thought the open boundary due to the numerical error in calculating the wave celerity. Since the numerical simulation results indicate that the leading wave has the maximum amplitude, the wave reflection can be avoided by adjusting the simulation duration to ensure the leading wave remains within the water domain. This can be achieved from visual observation on the wave elevation from the numerical simulation results. In summary, the limitations of the current numerical method and models are closely related to the complexity of the problem domain geometry, the movement of the landslides, and limitations caused by computer resources. Some of the limitations may be overcome by simplifying the problem domain, adjusting simulation duration and modifying discretization scheme. Other limitations, such as nonlinear wave effects and inconstant water depth, require some modification and extension to the numerical method. Within these limitations, the current numerical method and models provide a reliable means to study landslide-generated waves. Chapter 6 Engineering Applications 58 CHAPTER 6 ENGINEERING APPLICATIONS Based on the numerical results presented and discussed in the previous chapter, tentative solutions for real engineering applications have been sought in an attempt to present coastal engineers a quick way to obtain order-of-magnitude estimates of landslidegenerated waves. Design curves of such solutions have been developed, and the results of a case study are discussed in the following sections. 6.1 Design curves for engineering applications It is seen from Figs. 29 and 34 that the dimensionless wave height H / a 0 decays rapidly as the distance from the landslide source increases. A simple exponential relationship simulating such a phenomenon is as follows: H =a ix (6. 1) where a and p are two parameters that are dependent on the dimensionless parameter d/gT 2 in two-dimensional cases, and on the dimensionless parameters d/gT , b / d , and 2 y / d in three-dimensional cases. Equation 6.1 can be converted as: In v oy a = lncc + p l n vdJ (6.2) The above relationship can be readily applied to find a and P parameters using the leastsquares method. Figure 35 illustrates such a fit for the two-dimensional landslide- generated wave problem when d/gT = 0.02. 2 For various d/gT 2 values, more fits were performed on the dimensionless wave height for both two-dimensional and three-dimensional cases and design curves for a and P parameters were developed. Figure 36 shows the design curves for a and P parameters as a function of d/gT 2 in two-dimensional case. Figure 37 shows the design curves for Chapter 6 Engineering Applications a and p parameters as a function of d/gT 59 2 for waves at the centerline (y = 0) for different b/d values in three-dimensional cases. Note that Eq. 6.1 and the engineering curves presented in Figs. 36 and 37 for a and P parameters present only order-of-magnitude estimates for landslide-generated waves in a rectangular wave basin within the context of this thesis. Also, Eq. 6.1 produces an infinitely large H / a 0 value when x/d becomes infinitely small, which in reality is not true. Thus, Eq. 6.1 and Figs. 36 and 37 may be used cautiously to predict waves only at locations that are a certain distance away from the landslide source. This restricting distance may vary in different circumstances. A distance that is equal to one water depth may be used as a reference restricting distance. 6.2 Case study: Skagway Harbour landslide Landslides in coastal areas and reservoirs have been recorded in many countries during the past a few decades. However, the detailed information on landslide formation and movement as well as magnitude of landslide-generated waves is rare. In this section, one of the most detailed cases on landslide-generated waves, the Skagway Harbour landslide, is presented and studied with reference to the numerical modeling results. At approximately 7:10 p.m. Alaska standard time (AST) on November 3, 1994, a landslide occurred in Skagway Harbour (59°26'45"N, 135°19'30"W, see Fig. 38), located at the north end of Taiya Inlet, in southeast Alaska. The landslide was triggered by a 260 m length of the Pacific and Arctic Railway and Navigation Company (PARN) dock rapidly sinking into the water. project. The dock was in the early stage of a reconstruction The failure occurred coincidentally with an extreme low tide. The landslide caused a series of large-amplitude sea level oscillations estimated by eyewitnesses to be 5 - 6 m in height in the inlet, and 9 - 11 m at the shoreline. A National Oceanic and Atmospheric Administration (NOAA) analog tide gauge located on the west side of the harbour, opposite the landslide, recorded waves with periods of roughly 3 minutes and maximum trough-to-crest wave heights of about 2 m (Kulikov et al., 1996). The landslide caused one fatality, and the generated tsunami caused approximately Chapter 6 Engineering Applications 60 $1,000,000 in damage to the Ferry Terminal and about $100,000 in damage to the Skagway Boat Harbour. According to Cornforth et al. (1996), the P A R N dock on the east side of the harbour below steep cliffs was undergoing a major renovation to the northern half of the 400 m long structure (see Figs. 38 and 39). The timber deck and the supporting timber piles were being dismantled and replaced with sheet pile cellular bulkheads. The bulkheads were tied back to the shoreline by tail walls approximately 23 m long. The cells were to be backfilled with compacted dredged sand. Construction work began in late September 1994. The dock owner brought fill and granite riprap to the site, built a platform of approximately 90 m to 100 m long, and by October 15, 1994, stockpiled riprap on top of the platform to a maximum height of about 5.5 m (see Fig. 39). The volume of the platform fill was approximately 5,120 m , and the volume of the riprap 3 stockpile was approximately 4,240 m . On November 3, an extreme low tide of elevation 3 -1.3 m occurred at 6:47 p.m. At about 7:10 p.m., the entire south end of the existing dock and the new construction slid into the water causing the death of a member of the working crew. Eyewitness described the sheet pile wall snapping with the ground sliding out from under the cell, and the timber piling disappearing as a huge wall of water and debris came crashing towards the shoreline. The job superintendent estimated the total time from the beginning to total loss of ground below water as 15 to 20 seconds. The background information on the Skagway Harbour landslide provides some data that can be used in numerical studies. The initial stage of the landslide can be treated as a slump type landslide, which can be approximately simplified as a horizontal moving wall. The water depth in front of the dock was approximately 60 m. The duration of the landslide is considered to be 15 to 20 seconds. As an approximation, the landslide thickness is treated as the displacement amplitude. Measured from Fig. 39, this displacement amplitude is approximately 10 m. The length of the landslide was 260 m. Based on the above data, the Skagway Harbour landslide corresponds to one of the test cases presented in Chapter 5, which has d/(gT ) = 0.02, b/d = 4. Figure 33 presents 2 the contours of dimensionless wave height that is generated by such landslides and Fig. Chapter 6 Engineering Applications 61 34 shows the dimensionless wave height along the centerline of the landslide (perpendicular to the dock face). As seen from Figs. 33 and 34, along the centerline of the landslide and adjacent to the slide source, the maximum wave height can be as high as 0.6ao ~ 0.8ao, or 6 m to 8 m. This agrees with the 5 m ~ 6 m wave height observed by the eyewitnesses (Kulikov et al., 1996). The tide gauge, which is located at the west side of the harbour, is approximately 450 m away from the center of the landslide. Taking 400 m and 150 m as the distance along the centerline of the slide (x direction) and along the dock face (y direction) respectively, results in x/d = 6 and y / d = 2.5. Referring to Fig. 33, by extending the contour lines in x direction, the numerical simulation results show that the wave height at the location of the tide gauge is approximately 2.0 m. This also agrees well with the tide gauge measurement (Kulikov et al., 1996). A more accurate numerical simulation taking into account the slope effect is possible following the numerical methods developed in this thesis, but will be very timeconsuming. On the other hand, it is demonstrated that by approximating the slump-type landslide as a horizontal moving body, the numerical modeling results presented in this thesis can be used to estimate the magnitude of waves to a certain degree of accuracy. 6.3 Potential of landslide-generated waves in Canada There has been a record of landslide-generated waves in various locations in Canada. In many cases, waves generated by landslides entering rivers, lakes or other bodies of water have caused substantial damage. To mitigate such damage, the potential of landslidegenerated waves needs to be evaluated. Landslide entering rivers may be devastating. In 1908, at Notre Dame de la Salette, Quebec, a landslide into the frozen Lievre River generated a wave of water containing ice blocks that destroyed 12 houses and claimed 33 lives. It was eastern Canada's largest landslide disaster. Underwater landslides occur on the steep underwater slopes of deltas in lakes and the sea as well as on the slopes of the continental shelf. Damage may be caused to structures sited on the edge of deltas when underwater landslides involve the delta Chapter 6 Engineering Applications 62 surface. At Woodfibre, British Columbia, such a landslide caused damage to wharves and warehouses in 1955. In 1975, at Kitimat, British Columbia, a large failure on the underwater slope of the Kitimat delta generated an 8 m wave that swept around the shores of Kitimat Arm and damaged wharf facilities. Destructive potential of landslide-generated waves has also been recognized at many other locations in the province of British Columbia. Such locations include the Fraser River, Sturgeon Bank, Harrison Lake, Okanagan Lake, Columbia River, and Kootenay Lake. Some experimental work and studies of such potential have been carried out by Northwest Hydraulic Consultants Ltd. (1970, 1976, 1980 and 1990) as well as other authors (Chaudhry et al 1983). Floating bridges, dams, wharfs and marinas, and residents at waterfront locations are all under threaten i f catastrophic events such as landslide-generated tsunamis occur. To mitigate damages caused by landslide-generated waves, comprehensive studies involving cooperative works by geologists, geotechnical experts and coastal engineers are in great need. The theoretical method and numerical models presented in this thesis provide a relative easy-to-use, flexible and powerful means for such tasks. Chapter 7 Conclusions and Recommendations 63 C H A P T E R 7 CONCLUSIONS AND RECOMMENDATIONS 7.1 Summary A time-domain boundary element method has been developed to treat linear waves generated by a landslide of arbitrary shape into a water body of constant depth. B y solving the first-order boundary value problem in the time domain, the present method has greater versatility, less complicated algebra and requires less computer resources compared to other numerical methods. Numerical models to simulate waves generated by either wave paddles or horizontal-moving landslides have been set up and validated by comparisons with a number of available analytical, numerical and experimental results. The major difficulties associated with the full nonlinear problem on landslidegenerated waves arise primarily from the two nonlinear free surface boundary conditions applied on the instantaneous free surface, which itself is unknown a priori, and from the poorly defined radiation condition in the time-domain. Considering linear waves only and based on the assumption that the wave amplitude is moderate, the nonlinear terms in the free surface boundary conditions are neglected and the simplified boundary conditions are applied to the still water surface. The field solutions on a time- independent boundary which includes the still water surface are then obtained by the application of a boundary integral equation. Since the boundary is invariant in time, the solution to the system of simultaneous equations obtained through the discretization of the integral equation is required only once, rather than at each time step. The flow potential at each time step is then obtained through solving the integral equation. The Sommerfeld radiation condition applied to the open boundary is modified to include a time-dependent celerity that is numerically assessed. For wavemaker problems, a combined Sommerfeld radiation condition and sponge layer damping technique is applied to the open boundary to achieve steady wave profiles. The initial conditions correspond to a simplified horizontal landslide movement following a sinusoidal speed profile, and the velocity potential is developed in time and space through the gradual imposition of Chapter 7 Conclusions and Recommendations body surface boundary conditions. 64 The free surface boundary conditions and the radiation condition are then satisfied by a numerical integration in time. Although the computer program for the three-dimensional problem which generates and solves the matrix equation is computationally intensive, the left-hand-side matrix in the matrix equation involved in the formulation is a function of geometry and discretization only. Results of the left-hand-side matrix generated by the program can be stored in the disk memory, and used subsequently in the time-stepping program for different landslide scenarios and varying simulation durations. Because of this, the application of the present method with a larger number of facets is economically feasible. The numerical models are coded in FORTRAN90 programming language, and utilize numerical routines in the IMSL Math Library included in the Microsoft Fortran Power Station (version 4.0) package. These numerical routines assess the accuracy of the matrix solution by evaluating condition numbers and issue a warning message when the matrix is ill-conditioned. If no warning message is issued, the solution is considered to be accurate to at least the machine precision. Features of dynamic allocated arrays are implemented in the programs so that the rank of the matrix equation can be specified by users at runtime. The numerical models are validated by comparisons of numerical results with a number of analytical, numerical and experimental results from previous studies for both the wavemaker problems and the landslide problems. Good agreement has been obtained in all cases; Dimensional analysis and a number of numerical simulations were carried out to study the characteristics of waves generated by a horizontal landslide movement. Engineering curves were also developed based on the numerical simulation results to provide order-of-magnitude estimate on wave heights generated by such slides. Finally, a case study was carried out to demonstrate the application of the numerical results. Chapter 7 Conclusions and Recommendations 7.2 65 Conclusions The following conclusions have been obtained through the numerical study of the landslide-generated waves: • The numerical results indicate that, for horizontal landslide movement, the generated leading wave has the maximum amplitude. A number of smaller "wake waves" are also generated following the leading wave. • For a specified landslide dimension and amplitude, the shorter the landslide duration or the higher the landslide speed, the higher are the waves generated. • For a specified landslide duration and amplitude, the larger the width (dimension) of the landslide, the higher are the waves generated. • Landslide-generated waves decay as they propagate. In two-dimensional cases, wave heights tend to converge to 0.4a 0 far away from the landslide unless d/gT 2 is unduly small (say less than 0.01), where a is the landslide movement amplitude. In 0 three-dimensional cases, waves decay much faster. • In three-dimensional cases, when the landslide duration is relatively short or the landslide speed is relatively high, waves along the centerline of the landslide are generally higher than those in other directions. When the landslide duration is relatively long or the landslide speed is relatively low, the generated waves propagate evenly on the water surface. • For three-dimensional landslide-generated waves, contours of dimensionless wave height as a function of location for various b/d and d/gT values have been 2 developed. These contours can be used to predict landslide-generated waves in engineering applications. • The wave heights in two-dimensional cases or those at the centerline of the slide in three-dimensional cases as a function of distance from the slide source have been Chapter 7 Conclusions and Recommendations 66 found to follow an approximate exponential relationship. Parametric engineering curves of such relationships have been developed and presented to estimate order-ofmagnitude landslide-generated waves for the purpose of engineering applications. With these curves, engineers are able to quickly estimate wave amplitudes at locations that are a certain distance away from the slide source given the slide duration and displacement. • It has been demonstrated that some slump-type landslides can be approximated to a horizontal moving body, and hence be simulated using the numerical models developed in this thesis. A case study on the Skagway Harbour landslide was carried out in which the slide was treated as a slump-type slide. Numerical simulation results on wave heights show good agreement with the field measurements. 7.3 Recommendations for further study In this thesis, a time-domain boundary element method has been developed to study the landslide-generated waves. The numerical models have been set up to simulate waves generated by a horizontal landslide movement, and have been validated against a number of available analytical, numerical and experimental results on both wavemaker and landslide problems. The numerical simulation results demonstrate a number of features of interest, and may be used to supply detailed wave information for other wave models that account for diffraction, refraction and shoaling effects as waves propagate away from the landslide site. The work carried out and the study performed of the related literature indicates that there are several extensions of the present research, which may improve the overall understanding and prediction of landslide-generated waves. Prospective extensions of the present study are as follows: • A n experimental study on waves generated by horizontal landslide movement would further validate the numerical models. Additional experimental studies on other types of landslide, such as slump-type landslide, would be useful to assess the extent to which the numerical simulation results for vertical walls and horizontal landslide movement may be applied to slump-type landslide movement. Chapter 7 Conclusions and Recommendations • ; 67 With modifications to the discretization scheme, the numerical method can be extended to simulate waves generated by landslides of arbitrary shapes, submarine landslides, earthquake-induced seabed displacements, and other scenarios. • The landslide motion has been assumed to follow a half sinusoidal curve in the present study. However, various real-time landslide motions can also be simulated provided the time-histories of such motions are available. Examination of the influence of landslide motions can then be assessed through numerical simulations. • In the present study, only linear waves are considered. For nonlinear wave effects, the present method can be extended based on the Stokes second-order expansion procedure. The nonlinear free surface boundary conditions defined on the instantaneous free surface can be expanded at about the still water level by a Taylor series expansion and terms up to the second order be retained through a perturbation procedure (Isaacson, 1991). . 68 REFERENCES Bleistein, N . (1984). "Mathematical methods for wave phenomena." Academic Press, Inc., Orlando, Florida, 341 pp. Boo, S.Y., and K i m , C H . (1994). "Simulation of fully nonlinear irregular waves in a 3D numerical wave tank." Proceedings of the Fourth International Offshore and Polar Engineering Conference, Osaka, Japan, 1994, Vol. Ill, 17-24. Boo, S.Y., and Kim, C H . (1995). "Weakly nonlinear diffraction due to vertical cylinders in a 3-D numerical wave tank." Proceedings of the Fifth International Offshore and Polar Engineering Conference, Hague, Netherlands, 1995, Vol. Ill, 19-25. Boo, S.Y., and Kim, C H . (1996). "Fully nonlinear diffraction due to a vertical circular cylinder in a 3-D H O B E M numerical wave tank." Proceedings of the Sixth International Offshore and Polar Engineering Conference, Los Angeles, U S A , 1996, Vol. Ill, 23-30. Chaudhry, M . H . , Mercer, A.G., and Cass, D. (1983). "Modeling of slide-generated waves." Journal of Hydraulic Engineering, ASCE, 109(11), 1,505-1,520. Chen, S., and Sharma, S.D. (1996). "Far field condition for nonlinear scattered waves." Proceedings of the Sixth International Offshore and Polar Engineering Conference, Los Angeles, U S A , 1996, Vol. Ill, 50-56. Cheung, K.F. (1991). "Time-domain solution for second-order wave diffraction." PhD thesis, Department of Civil Engineering, University of British Columbia. Contento, G., and Casole, S. (1995). "On the generation and propagation of waves in 2D numerical wave tanks." Proceedings of the Fifth International Offshore and Polar Engineering Conference, Hague, Netherlands, 1995, Vol. Ill, 10-18. Cornforth, D.H. and Lowell, J.A. (1996). "The 1994 submarine slope failure at Skagway, Alaska." Proceedings of the Seventh International Symposium on Landslides, Trondheim, Norway, 1996, Vol.1, 527-532. Davidson, D.D. and Whalin, R.W. (1974). "Potential landslide-generated water waves, Libby Dam and Lake Koocanusa, Montana." Technical Report H-74-15, U.S. Army Engineer Waterways Experimental Station. Dean, R.G. (1982). "Water Wave Mechanics for Engineers and Scientists." Prentice-Hall Inc., Englewood Cliffs, New Jersey 07632. DeSanto, J.A. (1992). "Scalar wave theory." Springer-Verlag Berlin Heidelberg New York, 189 pp. 69 Dommermuth, D.G., and Yue, D.K.P. (1987). " A high-order spectral method for the study of nonlinear gravity waves." Journal of fluid Mechanics, Vol. 184, 276-288. Eie, J., Solberg, G., Tvinnereim, K . , and T<j)rum, A . (1971). "Wave generated by landslides." Proceedings of the 1st International Conference on Port and Ocean Engineering under Arctic Conditions, Trondheim, Norway, 1, 489-513. Gozali, S., and Hunt, B. (1989). "Water waves generated by close landslides." Journal of Hydraulic Research, 27(1), 49-60. Han, T.Y., Meng. J.C.S., and Innis, G.E. (1983). " A n open boundary condition for incompressible stratified flows." Journal of Computational Physics, 49, 276-297. Hammack, J.L. (1980). "Long waves generated by complex bottom motions." Coastal Engineering, 1980, 639-651. Harbitz, C.B. (1992). "Model simulations of tsunamis generated by the Storegga Slides." Marine Geology, 105, 1-21. Harbitz, C.B., Pedersen, G., and Gjevik, B. (1993). "Numerical simulations of large water waves due to landslides." Journal of Hydraulic Engineering, A S C E , 119(12), 1,325-1,342. Heinrich, P. (1992). "Nonlinear water waves generated by submarine and aerial landslides." Journal of Waterways, Port, Coastal and Ocean Engineering, A S C E , 118(3), 249-266. Huang, C.J., Zhang, E.C, and Lee J.F. (1995). "Nonlinear viscous wave fields generated by a piston-type wavemaker." Proceedings of the Fifth International Offshore and Polar Engineering Conference, Hague, Netherlands, 1995, Vol. Ill, 34-41. Hunt, B . (1988). "Water waves generated by distant landslides. " Journal of Hydraulic Research, 26, 307-322. Isaacson, M . , and Cheung, K.F. (1991). "Second order wave diffraction around twodimensional bodies by time-domain method." Applied Ocean Research, V o l . 13, 175-186. Isaacson, M . , Cheung K.F., Mansard E., and Miles. M . D . (1993). "Transient wave propagation in a laboratory flume propagation." Journal of Hydraulic Research, 31, 665-680. Jiang, L., and LeBlond, P.H. (1992). "The coupling of a submarine slide and the surface waves which it generates." Journal of Geophysical Research, 97(C8), 12,731-12,744. Jiang, L . , and LeBlond, P.H. (1993). "Numerical modeling of an underwater Bingham plastic mudslide and the waves which it generates." Journal of Geophysical Research, 98(C6), 10,303-10,317. 70 Jiang, L . , and LeBlond, P.H. (1994). "Three-dimensional modeling of tsunami generation due to a submarine mudslide." Journal of Physical Oceanography, 24, 559-572. K i m , C H . (1995). "Recent progress in numerical wave tank research: a review." Proceedings of the Fifth International Offshore and Polar Engineering Conference, Hague, Netherlands, 1995, Vol. Ill, 1-9. K i m , C H . (1999). "Recent research and development of numerical wave tanks - a review." International Journal of Offshore and Polar Engineering, ISOPE , V o l . 9, No. 4, December, 1999, 241-256. K i m , Y.J., K i m D.J., and Hwang, J.H. (1994). "Calculation of nonlinear free-surface flows using two-dimensional numerical wave tank." Proceedings of the Fourth International Offshore and Polar Engineering Conference, Osaka, Japan, 1994, Vol. Ill, 25-31. Koutitas, C G . (1977). "Finite element approach to waves due to landslides." Journal of the Hydraulics Division, A S C E , 103, 1,021-1,029. Kulikov, E . A . , Rabinovich, A . B . , Thomson, R.E. and Bornhold, B.D. (1996). "The landslide tsunami of November 3, 1994, Skagway Harbor, Alaska." Journal of Geophysical Research, Vol. 101, No. C3, 6,609-6,615. Liu, Y . , Dommermuth, D.G., and Yue, D.K.P. (1992). " A high-order spectral method for nonlinear wave-body interactions." Journal of fluid Mechanics, Vol. 245, 115-136. Longuet-Higgins, M.S.and Cokelet, E.D. (1976). "The deformation of steep surface waves on water. I. A numerical method of computation." Proceedings of Royal Society of London, A350, 1-26. Naito, S. (1994). "Research on element wavemakers and wave field generated by their combination." Proceedings of the Fourth International Offshore and Polar Engineering Conference, Osaka, Japan, 1994, Vol. Ill, 8-16. Noda, E. (1970). "Water waves generated by landslides." Journal of the Waterways, Harbors and Coastal Engineering Division, A S C E , 96(4), 835-853. Northwest Hydraulic Consultants Ltd. (1970). "Hydraulic model studies: wave action generated by slides into Mica Reservoir." Technical report. Northwest Hydraulic Consultants Ltd. (1976). "Hydraulic model study of waves from Downie slide." Technical report. Northwest Hydraulic Consultants Ltd. (1980). "Review of landslide generated waves." Technical report. 71 Northwest Hydraulic Consultants Ltd. (1990). " M T . Breakenridge landslide study: landslide generated wave analysis." Technical report. Orlandski, I. (1976). " A simple boundary condition for unbounded hyperbolic flows." Journal of Computational Physics, 21,251 -269. Raney, D.C., and Butler H.L. (1976). "Landslide generated water wave model." Journal of the Hydraulics Division, A S C E , 102, 1,269-1,282. Roach, G.F. (1982). "Green's functions." Cambridge University Press, Cabridge, 325 PPRubino, A., Backhaus, J.O. and Pierini, S. (1994). "Tsunamis generated by mud-slides." International Sumposium: Wave-Physical and Numerical Modelling, Vancouver, Canada, 1994, 466-475. Sarpkaya, T., and Isaacson, M . (1981). "Mechanics of wave forces on offshore structures." Van Nostrand Reinhold, New York, 651 pp. Skourup, J. (1996). "Active absorption in a numerical wave tank." Proceedings of the Sixth International Offshore and Polar Engineering Conference, Los Angeles, U S A , 1996, Vol. Ill, 31-38. Skourup, J., and Schaffer, H.A. (1998). "Simulations with a 3D active absorption method in a numerical wave tank." Proceedings of the Eighth International Offshore and Polar Engineering Conference, Montreal, Canada, 1998, Vol. Ill, 248-255. Stassen, Y . , Boulluec, M . L , and Molin, B. (1998). " A high order boundary element model for 2D wave tank simulation." Proceedings of the Eighth International Offshore and Polar Engineering Conference, Montreal, Canada, 1998, Vol. Ill, 348355. Stoker, J.J. (1957). "Water waves." Interscience Publishers, Inc., New York, 567 pp. Takayama, T. (1984). "Theory of oblique waves generated by serpent-type wave-maker." Coastal Engineering in Japan, Vol. 27, 1-19. U.S. Army Corps of Engineers (1984). "Shore Protection Manual", Coastal Engineering Research Center, Waterways Experiment Station, Vicksburg, Mississippi. Wang, Y . (1994). " A study on periodic wave breaking by absorbed numerical wave channel." Proceedings of the Fourth International Offshore and Polar Engineering Conference, Osaka, Japan, 1994, Vol. Ill, 32-36. Watts, P. (1998). "Wavemaker curves for tsunamis generated by underwater landslides." Journal of Waterway, Port, Coastal and Ocean Engineering, A S C E , 127-137. 72 Wiegel, R.L., Noda, E.K., Kuba, E . M . , Gee., D . M . , and Tornberg, G.F. (1970). "Water waves generated by landslides in reservoirs." Journal of the Waterways and Harbors Division, A S C E , 96, 307-333. Williams, A . N . , and Crull, W.W. (1996). "Theoretical analysis of directional waves in a laboratory basin." Proceedings of the Sixth International Offshore and Polar Engineering Conference, Los Angeles, USA, 1996, Vol. Ill, 65-72. 73 No. Wave Period (s) Computed Wave Length (m) Computed Wave Amplitude no (m) kd H/a 1 20.0 64.6 0.001 0.10 0.1 2 10.0 32.0 0.002 0.20 0.2 3 5.15 15.9 0.004 0.40 0.4 4 3.10 9.1 0.007 0.69 0.7 5 2.30 6.2 0.010 1.01 1.0 6 1.72 4.2 0.014 1.50 1.4 7 1.45 3.1 0.017 2.03 1.7 8 1.28 2.5 0.019 2.51 1.9 9 1.16 2.1 . 0.020 2.99 2.0 10 1.08 1.8 0.021 3.49 2.0 11 1.00 1.6 0.021 3.93 2.0 0 Table 1 Numerical simulations of two-dimensional wavemaker problem (d = 1.0 m, a = 0.02 m). 0 74 Test Case A B C D Wave Basin Dimension 6.16 m x 12.32 m 8.22 m x 16.44 m 6.16 m x 12.32 m 10.96 m x 21.92 m Wave Paddle Width 0.77 m 2.74 m 5m 5m Wave Paddle Motion Period 1.0 s 1.4 s 1.0 s 1.4 s Note in all test cases: • Water depth was 0.63 m; • Wave paddle motion magnitude was 0.02 m; • Simulation duration was 10 s; • Water surface was digitized into a 48 x 48 mesh; • Water depth was digitized into 4 segments; • The total number of facet was 2,880. Table 2 Test cases for three-dimensional wavemaker problems. 75 Time . < 1 3 k 1 4.1 s 6.5 s 8.1 s 9.7 s 3.04E-4 m 3.99E-5 m 5.01E-4m 6.86E-4 m 3.93E-5 m 5.58E-5 m 6.64E-5 m 7.34E-5 m 0.129 0.140 0.133 0.107 1.86 1.79 1.84 2.03 5.84E-6 m 9.08E-6 m 1.02E-5m 8.81E-6m Convergence Order k J 3 111 3 - l ID, ""ID, k Note: For all cases, wave tank length = 30 m, sponge layer length = 5 m, water depth = 1.0 m, wave paddle motion magnitude = 0.02 m, wave period = 2.3 s. • Scheme D i : 60 segments on water surface, At = 0.081 s. • Scheme D2: 180 segments on water surface, At = 0.009 s; • Scheme D 3 : 540 segments on water surface, At = 0.001 s. Table 3 Accuracy of numerical discretization for two-dimensional wavemaker problem. 76 Time < -< , < -< 1 9 k 0.4 s 0.6 s 0.9 s 1.1 s 9.77E-4 m 1.70E-3m 1.31E-3 2.66E-3 m 1.39E-4m 1.99E-4m 2.94E-4 4.09E-4 m 0.142 0.117 0.146 0.154 0.89 0.98 0.88 0.85 2.31E-5m 2.65E-5 m 5.02E-5 7.42E-5 m i Convergence Order k J 9 k 3 111 - 1 Fin. "I, Note: For all cases, wave basin length = 1.54 m, width = 3.08 m, sponge layer length = 0.77 m, water depth = 0.63 m, wave paddle motion magnitude = 0.02 m, wave period = 1.0 s. • Scheme D i : 6x 6 mesh on water surface, At = 0.081 s; • Scheme D2: 18x18 mesh on water surface, At = 0.009 s; • Scheme D3: 54 x 54 mesh on water surface, At = 0.001 s. Table 4 Accuracy of numerical discretization for three-dimensional wavemaker problem. •//////////////////////////////////. (a) vertical rock fall 1 WW//////////////////////////////. (b) horizontal wall movement '////////////„ ,W9m,, (c) slump W77s777777777777777777777?. (d) vertical seabed movement 777777m V777777777, (e) wedge sliding down a slope (f) submarine landslide Figure 1 Sketches of various types of landslide-generated waves. 78 Figure 2 Definition sketch of the problem of landslide-generated waves - plan view. Figure 3 Definition sketch of the problem of landslide-generated waves - side view. 79 Figure 4 Plan view of simplified three-dimensional wave tank (landslide problem). Landslide So S f • V(t) JS W Sb *7SSSSSSSSSSSSSSSSSSSSSS>SSSSSSSSSSJ Figure 5 Side view of simplified three-dimensional wave tank (landslide problem). 80 Sponge Layer Sponge Layer Wave Paddle ^ S, -> x X b Sponge Layer Figure 6 Plan view of three-dimensional wave tank with sponge layer (wavemaker problem). Wave Paddle t Sponge Layer X S b So L f x *• S7SSSSSSSSSSSSSSS///////S//////SSS/ ,—• 1 Figure 7 Side view of three-dimensional wavemaker problem and sponge layer technique (wavemaker problem). X 81 x Figure 8 Discretization of three-dimensional landslide problem. 82 mirror image line Figure 9 Discretization of two-dimensional landslide problem. 83 ^ Sponge Layer^ 0.02 r|(m) o.oo -0.02 0.02 t = 8 sec r| ( m ) o.oo -0.02 0.02 r| ( m ) o.oo t = 10 sec \ \ -0.02 0.02 t = 12 sec TI ( m ) o.oo -0.02 0.02 t = 14 sec r) ( m ) o.oo -0.02 0.02 t = 16 sec /\ r | ( m ) o.oo -0.02 10 20 x(m) 30 Figure 10 Numerical simulation of two-dimensional wavemaker problem: wave profile at different times (d = 1 m, ao = 0.02 m, T = 2.3 s). 40 84 0.02 n(m) o.oo -0.02 1 1 0 1 1 4 1 1 8 I 1 I 1 I 12 16 20 12 16 20 t(sec) 0.02 x= 20.0 m U(m) 0.00 0 4 8 t(sec) Figure 11 Numerical simulation of two-dimensional wavemaker problem: wave elevation at fixed locations (d = 1 m, ao = 0.02 m, T = 2.3 s). 85 Figure 12 Comparisons of dimensionless wave height between theoretical solution and numerical results for steady state two-dimensional wavemaker problem. 86 Figure 13 Comparisons of dimensionless wave amplitudes at fixed locations between Noda's theoretical solution and numerical results on two-dimensional landslide problem when simplified landslide is moving horizontally at a constant speed. 87 Rubble Mound Numerical probes rR > ' ^ ~ * Rubble Mound A ,,,, 1 ' & (a) Plan view of the wave basin (A, B, C, D, E, F are the numerical probes in the numerical simulation). Wove measuring point, l00B4«aQ0' CUnil':'cmJ 'v lOOfl-4 '400 (b) Experimental wave measuring points. Figure 14 Experimental setup by Takayama (1984). Figure 15 Numerical simulation of three-dimensional wavemaker: wave profiles at different times (wave elevation in metres). -0.01 10 0.01 r\ (m) 12 14 B (3.08, 0.0) o.OO -0.01 0.01 tl (m) 0.00 -0.01 0.01 r| (m) o.OO -0.01 0.01 n (m) o.oo -0.01 0.01 r| (m) o.OO -0.01 time (sec) Figure 16 Numerical simulation of three-dimensional wavemaker: wave elevation at fixed points (numbers in brackets showing the coordinates of fixed points in metres). 90 Theoretical solution Numerical results Figure 17 Comparisons between theoretical solutions and numerical results of threedimensional wavemaker problem for b/L = 0.5 (contours showing values of H/a ). 0 91 Theoretical solution Numerical results Figure 18 Comparisons between theoretical solutions and numerical results of threedimensional wavemaker problem for b/L =1.0 (contours showing values of H/ao). 92 Figure 19 Comparisons between the numerical results and Takayama's theoretical and experimental results in x direction for three-dimensional wavemaker problem (T= 1.0 s). Figure 20 Comparisons between the numerical results and Takayama's theoretical and experimental results in y direction for three-dimensional wavemaker problem (T= 1.0 s). 94 Figure 21 Comparisons between the numerical results and Takayama's theoretical and experimental results in x direction for three-dimensional wavemaker problem (T=1.4 s). Figure 22 Comparisons between the numerical results and Takayama's theoretical and experimental results in y direction for three-dimensional wavemaker problem (T= 1.4 s). 96 100 CD E i- 10 CL Fit : o y = 8.094 x 10 _ 5 x' 8 7 J _L 100 3 4 5 6 I I L 7 8 9 Rank of Matrix Equation Figure 23 C P U time as a function of the rank of the matrix equation for the twodimensional problem (100 time steps). 1000 97 1000 E i- 100 Fit: 0_ O 10 100 y = 1.215x!0^x 1.87 _L 4 J 5 I I I I 6 7 8 9 1000 Rank of Matrix Equation _L 2 J 5 I 6 I I I 7 8 9 Figure 24 C P U time as a function of the rank of the matrix equation for the threedimensional problem (100 time steps). 10000 98 99 Figure 26 Wave elevations at the centerline of the wave basin for three discretization schemes (t = 0.5 s). 100 o.oi t = 1.0 sec r| (m) 0.00 -0.01 0.01 t = 2.0 sec Ti(m) 0.00 -0.01 0.01 t = 3.0 sec r|(m)0.00 -0.01 0.01 t = 4.0 sec Tl(m) 0.00 -0.01 0.01 t = 5.0 sec r|(m) 0.00 -0.01 0.01 t = 6.0 sec T|(m) 0.00 -0.01 8 12 16 x(m) Figure 27 Wave elevation as a function of distances for two-dimensional landslidegenerated waves (d = 0.63 m, a = 0.01 m, T/2 = 0.9 s). 0 20 101 Figure 28 Wave elevation as a function of time at fixed locations for two-dimensional landslide-generated waves (d = 0.63 m, a = 0.01 m, T/2 = 0.9 s). 0 102 Figure 29 Numerical results of dimensionless wave heights as a function of dimensionless distance for various d/gT values on two-dimensional landslide-generated waves. 2 103 Figure 30 Numerical simulation of three-dimensional landslide-generated waves (d = 0.63 m, a = 0.01 m, T/2 = 0.5 s, dimensions are in metres) 0 104 Figure 31 Contours of dimensionless wave height H / a as function of dimensionless 0 location for b/d =1.0 and various d/gT values from numerical simulations of three-dimensional landslide-generated waves. 2 105 Figure 32 Contours of dimensionless wave height H / a as function of dimensionless 0 location for b/d = 2.0 and various d/gT values from numerical simulations of three-dimensional landslide-generated waves. 2 Figure 33 Contours of dimensionless wave height H / a as function of dimensionless 0 location for b/d = 4.0 and various d/gT values from numerical simulations of three-dimensional landslide-generated waves. 2 107 b/d = 1.0 d/(gT 2) = 0.06 A d/(gT 2) = 0.04 A d/(aT 2) = 0.02 A H/a 1.0 0 • -. \ . ~~ •" • ;1- . . . • — . b/d = 2.0 d/(gT 2) = 0.06 A d/(gT 2) = 0.04 A d/(aT 2) = 0.02 A H/ao 1.0 b/d = 4.0 d/(gT 2) = 0.06 A d/(gT 2) = 0.04 A d/(gT 2) = 0.02 A 1.5 H/a 10 0 0.5 " - 0 1 2 3 - 4 x/d Figure 34 Dimensionless wave heights at the centerline of wave tank (y/d = 0) as function of x/d for various b/d and d/gT values from numerical simulations on threedimensional landslide-generated waves. 2 108 Figure 35 Fitting curve for dimensionless wave height as function of dimensionless distance for two-dimensional landslide-generated waves when d/gT = 0.02. 2 109 a - - 0 -P 0.01 0.02 0.03 d/gT 0.04 0.05 0.06 2 Figure 36 Design curves of a and P parameters in determining two-dimensional landslide-generated, wave heights governed by H / a =ct(x/d) for various p 0 d/gT 2 values. 110 o - » - b / d = eo b/d = 4.0 - • - b / d = 2.0 - • - b / d = 1.0 P -0.3 • 0 0.01 0.02 0.03 A " 0.04 0.05 ^ 0.06 0.07 d/(gT ) 2 Figure 37 Design curves of a and P parameters in determining three-dimensional landslide-generated wave heights at the centerline (y/d = 0) governed by H / a = a(x/d) for various b/d and d/gT values. p 0 2 Ill Figure 38 Location of Skagway Harbour and the landslide area (Kulikov et al. 1996). 112 I LU: t_ I •iSCALE;H METRES Figure 39 Geologic section in the middle of landslide area (Cornforth et al. 1996).
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical study of landslide-generated waves
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical study of landslide-generated waves Yang, Gang 2002
pdf
Page Metadata
Item Metadata
Title | Numerical study of landslide-generated waves |
Creator |
Yang, Gang |
Date Issued | 2002 |
Description | This thesis describes a time-domain boundary element method to numerically simulate linear waves generated by a horizontal-moving landslide of an arbitrary profile. The approach to setting up the numerical models is based on a Green's function method and a time-stepping procedure. The numerical models include a boundary condition at the landslide surface that accounts for wave generation, a radiation condition in the far field that accounts for open water, and free surface boundary conditions in a time-stepping procedure that account for wave propagation. By solving the first-order boundary value problem in the time domain, waves generation and propagation are simulated numerically. Although the present numerical method can be applied to a variety of landslide configurations, only vertical wall with horizontal movement is considered in the present study. The corresponding numerical models in two dimensions and three dimensions are developed. The numerical models are initially used to simulate the waves generated by a piston-type wavemaker with a periodic wave paddle movement, and the results are validated against a number of analytical solutions as well as available experimental results. Specific two-dimensional landslide problems are also simulated on a horizontalmoving wall movement and comparisons of numerical results with analytical solutions have demonstrated good agreement. A number of numerical simulations on landslide-generated waves have been performed and a selection of numerical results is presented and analyzed. The results exhibit various features of interest that are discussed. Engineering curves have been developed based on the numerical results to estimate order-of-magnitude wave height caused by landslides. A case study has been carried out to demonstrate the application of numerical models. The numerical method can be applied to account for a landslide of an arbitrary shape and movement, and a water body surrounded by surfaces of various geometries. The approach and numerical models can also be extended to account for nonlinear wave effects. Overall, the numerical method and the numerical models are found to be able to provide a reasonably flexible and reliable means of studying and predicting landslidegenerated waves. |
Extent | 5612983 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-09-29 |
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.0063674 |
URI | http://hdl.handle.net/2429/13309 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2002-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2002-732606.pdf [ 5.35MB ]
- Metadata
- JSON: 831-1.0063674.json
- JSON-LD: 831-1.0063674-ld.json
- RDF/XML (Pretty): 831-1.0063674-rdf.xml
- RDF/JSON: 831-1.0063674-rdf.json
- Turtle: 831-1.0063674-turtle.txt
- N-Triples: 831-1.0063674-rdf-ntriples.txt
- Original Record: 831-1.0063674-source.json
- Full Text
- 831-1.0063674-fulltext.txt
- Citation
- 831-1.0063674.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-0063674/manifest