TRANSFORMATION OF WAVES AROUND ARTIFICIAL ISLANDS by KUSHAL K. TALUKDAR B. Tech., Indian Institute of Technology, Kharagpur, 1980 M. Tech., Indian Institute of Technology, Madras, 1982 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in FACULTY OF GRADUATE STUDIES Department of Civil Engineering We accept this thesis as conforming to the required standard UNIVERSITY OF BRITISH COLUMBIA May 1986 © KUSHAL K. TALUKDAR, 1986 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the THE UNIVERSITY OF BRITISH COLUMBIA, I agree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r reference and study. I f u r t h e r agree that permission f o r extensive copying of t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the Head of my Department or by h i s or her r e p r e s e n t a t i v e s . I t i s understood t h a t copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my w r i t t e n p e r m i s s i o n . Department of C i v i l Engineering THE UNIVERSITY OF BRITISH COLUMBIA 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: May 1986 ii ABSTRACT A numerical method based upon the finite element technique is described to determine the transformed wave field around artificial islands for regular unidirectional waves. The present work takes into account various aspects of wave transformation, including wave refraction, reflection, diffraction, shoaling and wave breaking, in a computer model which finally predicts the transformed wave heights and elevations at different locations. A modified version of the original two-dimensional mild slope equation, to include the energy dissipation from wave breaking has been used for this purpose along with criteria for the onset and cessation of wave breaking. Results are presented first for the cases where analytical solutions are available in order to check for accuracy and convergence of the present computer model. Since no previous results are available for the proper island geometries, the computer predictions are compared with measurements of the wave field around a model island in a laboratory' basin. In particular, the wave heights at specified points and the runup around the island are compared. The study indicates the need for including wave breaking in such an analysis and also suggests possible ways to calibrate certain features of the computer model for a better numerical prediction. iii TABLE OF CONTENTS ABSTRACT ii LIST OF FIGURES v ACKNOWLEDGEMENTS vi CHAPTER - I INTRODUCTION 1 1.1 Development of Artificial Islands 1 1.2 Types of Artificial Islands 2 1.3 The Aim of Present Work 3 CHAPTER - II LITERATURE REViEW 5 2.1 Pure Refraction 5 2.2 Refraction and Diffraction 6 2.3 Refraction and Forward Diffraction 8 2.4 Wave Energy Dissipation 10 CHAPTER - i n THEORETICAL FORMULATION 11 3.1 Problem Definition 11 3.2 Simplification of the Problem 16 3.3 Finite Element Formulation 18 3.4 Determination of the Energy Dissipation Factor 23 3.5 Wave Heights of the Transformed Field 24 CHAPTER - TV NUMERICAL PROCEDURE 25 4.1 The Element 25 4.2 Integration Technique 27 4.3 The Range of Wave Breaking 28 4.4 Arrangement of Global Matrix 29 4.5 The Solution Routine 30 4.6 Organization of the Program 32 iv CHAPTER - V RESULTS AND DISCUSSION 35 5.1 Verification Against Known Solutions 36 5.2 Convergence Tests 40 5.3 Results for a Typical Island 42 5.4 Experimental Set Up 42 5.5 Discussion 52 CHAPTER - VI CONCLUSION AND RECOMENDATIONS 56 REFERENCES 58 APPENDIX - I USER'S MANUAL 61 A-1.1 Description of Input Data File 61 A-1.2 Execution of Data Generating Program 63 A-1.3 Execution of Finite Element Program 63 A-1.4 Transfer of Graphics Files to VAX/VMS 64 A-1.5 Execution of the Graphics Program 65 A-1.6 Limitation of the Finite Element Program 65 V LIST OF FIGURES FIGURE PAGE 3.1 Sketch of the problem 12 3.2 Sketch of simplification 16 4.1 The mapping of the element 26 4.2 The local and global node numbering sequence 30 4.3 Flow chart of the finite element program 33 5.1 Geometry of the island with parabolic shoal 37 5.2 The finite element mesh for parabolic island 37 5.3 Relative amplitude and phase lag around parabolic island 38 5.4 Sketch of a circular cylinder on a horizontal seabed 38 5.5 Relative amplitude and phase lag for a vertical cylinder, «a = 0.6 39 5.6 Runup and phase lag for a vertical cylinder, «a = 1.0 39 5.7 Accuracy of runup with mesh size for a vertical cylinder 41 5.8 Accuracy of runup with boundary distance for a vertical cylinder 41 5.9 The geometry of a typical island 43 5.10 The runup profile around island, wave period = 10 seconds 43 5.11 The runup profile around island, wave period = 20 seconds 44 5.12 Sketch of the experimental model 44 5.13 Comparison of runups around model, wave period = 0.7 seconds 46 5.14 Comparison of runups around model, wave period = 1.4 seconds 46 5.15 Comparison of relative wave amplitudes inside the model berm 47 5.16 Comparison of relative wave amplitudes outside the model berm 47 5.17 Wave height distribution around model (with breaking) 48 5.18 Wave height distribution around model (without breaking) 49 5.19 Wave elevations around model island (with breaking) 50 5.20 Wave elevations around model island (without breaking) 51 vi ACKNOWLEDGEMENTS In presenting this thesis, I wish to express my heartfelt gratitude to Dr. M. de St Q. Isaacson for his invaluable guidance and suggestions during the entire course of this wort 1 would also like to thank Dr. N. D. Nathan sincerely for his critical review of the theoretical formulation during the early stages. Discussions with P. S. Han, P. Pattani and S. Abayakoon were of great help to me and I can offer here only an inadequate acknowledgement of my appreciation to them. I am also indebted to C. Padmore for his patient help with the experiments and to A. Martin for modifying the Civil Engineering Graphics Package to suit this particular problem. Finally, financial support in the form of a research assistantship from the Natural Sciences and Engineering Research Council, Canada is very much appreciated. 1 CHAPTER - I INTRODUCTION Artificial islands have been developed and used as a result of the exploration for hydrocarbons in the Arctic sea. From both economic and technical points of view, they have been considered to be most suitable for the search environment and a number of them have been built since 1972. Still to this day a thorough understanding of the interaction between the islands and the surrounding sea is not complete. The present work is an attempt to analyze the effects of such islands on the incident wave field by using a mathematical model which in turn is expected to lead to better design criteria. 1.1 DEVELOPMENT OF ARTIFICIAL ISLANDS The search for and exploitation of hydrocarbons from the seabed dates back to 1940's when the first offshore structures in the form of steel jackets were introduced in the Gulf of Mexico. Since then, the development of offshore resources has involved deeper waters and implementation of new structural designs to meet the challenges of different, more hostile environments. As the areas of search widened, new areas with potential for finding giant petroleum reservoirs became a distinct economic possibility, which were not taken into consideration before. One such area is the Mackenzie Bay area of the Beaufort sea, at the northern boundary of Canada. The unique nature of the environment there produced the design of 'artificial islands', the Canadian 2 contribution to the world list of offshore structures. A number of other drilling platform alternatives such as monopods, cones and several types of floating rigs have also been studied, but the reasons which favoured artificial islands can be summarized (Ref. Jong, 1976) as follows : i) Shallow depth of water, extending to about the 18.0 meter depth contour. ii) Short (two to three months) working season results in high standby cost of floating rigs in winter. iii) Islands are considered to be the safest option to counter environmental risks. iv) The initial capital investment is very high for other structures compared to artificial islands. v) The technology for building artificial islands is well known and construction work requires only standard construction equipment A study carried out by U.S Minerals Management Service for the U.S.A Alaskan Beaufort sea as reported in Offshore Engineer (April.1985) also confirms the suitability of artificial islands to a water depth of 20 meters. 1.2 TYPES OF ARTIFICIAL ISLANDS Artificial islands have two main sections. Firstly, the body of the island provides a base for drilling operations. Secondly, slope protection materials protect the slopes of the island against wave forces in the summer season and resist ice forces in winter. Depending upon their design, artificial islands may be divided into the 3 following categories. Sacrificial Beach Island : This has an unprotected beach profile with gentle slopes. The beach suffers considerable erosion during heavy storms and a large buffer zone of fill above water is required to prevent erosion of the actual island. Protected Island : The beaches of these islands are protected by various types of materials such as sand bags, gabions or wire netting cages filled with sand bags or rocks etc. Caisson Retained Island : Concrete or steel caisson (such as with Molikpaq described in Ocean Industry June,1985) and also hybrid caissons are ballasted on subsea berms to make the islands. As mentioned in Offshore (Feb,1985) they eliminate the need for elaborate shore proctection and also reduce costs by reducing fill requirements. 1.3 THE AIM OF PRESENT WORK From the above description it is clear that an estimation of the erosion of the beach or the berm is important from the design point of view. A knowledge of this requires an estimate of the longshore current obtained from the transformed wave field around the island. This can be a complicated process if every aspect of the transformed wave field is taken into account Such effects include the randomness of the waves, the effect of wave reflection and diffraction from and around the islands, shoaling on the slope, the presence of a predominant current flow if any, wave energy dissipation and the change in bathymetry. Various simplified models have been proposed taking some of these effects into account and ignoring others. Although most 4 of the mathematical models proposed so far take account of wave reflection, refraction, diffraction and shoaling for regular periodic waves, they do not generally take account of wave breaking due to the reduction in water depth. This results in a loss of energy which has a significant effect on the transformed wave field. The purpose of the present work is to include the effects of wave refraction, diffraction and energy dissipation due to wave breaking simultaneously in a mathematical model. Such a model is used to obtain a numerical solution using finite element techniques in order to obtain a more realistic transformed wave field around the island. 5 CHAPTER - H LITERATURE REVIEW In the first analytical models to determine changes in the wave field, only the effects of shoaling and refraction were taken into account Before this analytical procedure, changes in wave heights were determined by a graphical technique using refraction diagrams. However, if the size of the structure is comparable to the wavelength or if a sudden change in bottom slope occurs, then wave diffraction becomes important and pure refraction theory is no longer adequate. Thus the next analytical development also accounted for wave diffraction. Eventually the breaking of waves and the associated energy dissipation have been considered in order to produce a better and closer approximation to actual conditions. In this chapter various approaches to a description of the transformed wave field as developed over the years are described and their implications discussed. 2.1 PURE REFRACTION When a wave train approaches a sloping seabed at an angle, the section of the wave train in shallower water travels at a slower rate than the portion in the deeper water as indicated by the dispersion relation. Therefore, a wave crest tends to bend and become more aligned with the bottom contours. This is known as wave refraction, with the wave velocity obeying Snell's law _L _ (2.1) Co ainVo 6 where c is the wave velocity, \p is the angle the wave crest makes with the bottom contour and the subscript V denotes deep water conditions taken as a reference. This is the basis for various numerical or graphical schemes for tracing the path of wave orthogonals (lines orthogonal to wave crests) from deep to shallow water over a region where the depth contours are known. Arthur (1946) showed that Eq. (2.1) can be used to obtain a relatively simple analytical solution for circular islands having axisymmetric bathymetry. The paths of orthogonals should satisfy the first order differential equation where r is the radial distance of any point on an orthogonal from the centre of the island, x is the angle between the direction of wave propagation and the direction of r. r Q and x Q are respectively the radial distance and the angle subtended by the point on orthogonal at the beginning of the sloping portion of the seabed. This theory is valid only for very gradual slopes since it cannot take into account the scattering of waves by a vertical wall or the effects of a sudden change in slope. 2.2 REFRACTION AND DIFFRACTION A pioneering study which accounted for refraction and diffraction simultaneously was carried out by Homma (1950), who investigated the effects of long seismic sea waves (tsunamis) around an island with a parabolic shoal. He showed that when wave lengths are very long compared to the size of the island, the increase in wave height is due more to resonance of the water around the island rather than to changes in wave height due to shoaling. The governing differential equation for long waves is of (2.2) 7 elliptic type and if suitable conditions are prescribed on every section of the boundary, a solution to the differential equation exists and is unique in the closed domain. However, the problem of wave scattering by an island involves an open region and for the correct solution of the differential equation the open region must be terminated at some finite value where the correct waveform is known. The only restriction is that the solution should obey the radiation condition, i.e. scattered wave energy must propagate away from the island in a suitable way. Far off from the island the solution may be expressed by an infinite Fourier-Bessel series with unknown coefficients. Hankel functions describing the outgoing waves satisfy the required radiation condition but the precise form of the wave at any radial distance is not known beforehand. As a first approximation, Vastano and Reid (1967) replaced the Hankel functions in the expansion of the far field by their asymptotic forms and the numerical long wave solutions in the neighbourhood of a paraboloid island was obtained by using a finite difference technique. Lautenbacher (1970) also obtained solutions for long waves in terms of an integral derived with the aid of Green's theorem and appropriate Green's function. To overcome the restriction of long waves, a general linear theory of combined refraction and diffraction was first presented by Berkhoff (1972) under the assumptions that waves are of small amplitude (linear approximation) and the sea bed slope is gently varying, so that linear wave theory for constant water depth can be applied at any one location. This gives rise to the mild slope equation V(cc j ,V^) + — u24> = 0 (2.3) where c is the wave velocity, c is the group velocity, CJ is the angular frequency of the waves and V is the horizontal gradient operator. Also 4> is the two dimensional velocity potential function at the still water level. 8 To obtain a solution to the mild slope equation, Berkhoff (1972) employed a finite element technique for the region around the island and at a boundary some distance from the island matched the finite element solution with the far-field radiation condition satisfied by adopting a source distribution for the outer region. Unfortunately the solution destroys the symmetry of the equations and results in large requirements for computer time and storage. Chen and Mei (1974) adopted a special variational statement which includes the far-field series solution and results in a symmetrical set of equations, but they treat only the special case of the shallow water wave equation. Betts and Zienkiewicz (1977) stirnmarized different possible methods of coupling the exterior solution with the interior solution, viz; boundary dampers, exterior analytical or series solution, exterior boundary integral formulation and introduced the 'infinite elements' to take the far-field boundary conditions into account Houston (1981) adopted a hybrid solution technique, i.e. combining a finite element solution with an analytical solution for the far-field and showed that extremely close results can be achieved by using a refined mesh for waves having very short periods. Tsay and Liu (1983) calculated the wave forces and moments acting on docks based upon a finite element solution. 2.3 REFRACTION AND FORWARD DIFFRACTION The mild slope equation (2.3) is elliptical in form and needs to be solved over the whole domain simultaneously. One of the common requirements for the finite element solution of the mild slope equation is that the water depth should become constant at a certain distance away from the scattering object This cannot be satisfied easily in coastal regions with non-zero slopes. The parabolic approximation method developed by Radder (1979) circumvents this difficulty. In this method the wave velocity potential is split into two components, one for the transmitted waves moving 9 forward and the other one for reflected waves primarily propagating in the opposite direction. For most coastal engineering problems the transmitted wave is the important part of the wave, as the reflected waves move in the offshore direction. So only the forward scattered waves are treated in the analysis and this enables a parabolic approximation to Eq. (2.3) to be adopted. This converts the boundary value problem to an initial value problem such that the analyst can calculate the wave amplitude along a phase line, starting at an upwave boundary and marching along the predominant wave direction. An analytical solution for this kind of combined linear refraction and forward diffraction wave field in the neighbourhood of a thin breakwater was reported by Lozano and Liu (1980). The bottom topography was, however, restricted to a uniform plane beach and the breakwater was considered perpendicular to the shoreline. A more general numerical solution, where the orientation of the breakwater can be arbitrary, was obtained by Tsai and Liu (1982) and the accuracy of the method was demonstrated by comparing the numerical solution with laboratory data. The splitting of the velocity potential can be accomplished in a number of ways and accordingly alternative forms'of'the governing equation have been obtained. The linear equation given by Kirby and Dalrymple (1983) is 2iKccg — + 2K{K - icQ)ccg A* + i—{Kccg)A* + — \ccg — J = 0 (2.4) where A* is the amplitude of the velocity potential defined by Also g is the gravitational constant, K is the wave number and K q is the wave number in deep water. Using a suitable perturbation procedure the wave non-linearities can be taken into account and Kirby and Dalrymple (1983) showed that the effect of non-linearities 10 for Stokes waves is governed by an additional term proportional to |A |A. The governing equation can be further modified to take account of the effects of currents. Studies by Booij (1981) and Liu (1983) led to two different versions of the governing equation which includes the effect of currents. Kirby (1984) found the reasons for the discrepency and obtained the correct equation governing the propagation of linearized waves in regions with varying depth and current. 2.4 WAVE ENERGY DISSIPATION None of the models mentioned so far take the dissipation of energy into account in describing the wave transformation in the surf zone. The energy loss may be due to various possible effects, including seabed friction, porous bottom, bottom boundary layers etc., but the primary loss of energy is generally associated with the breaking of waves. Booij (1981) has shown that the inclusion of the term kjW0 in the mild slope equation, in which u is the angular wave frequency and W is an unspecified damping factor, incorporates the effects of wave damping.. The modified mild slope equation is thus: Dalrymple, Kirby and Hwang (1984) solved this modified version using a parabolic approximation and determining the damping factor for cases of localized areas of energy dissipation, such as due to pockets of mud in a sandy bottom, pile clusters etc. For straight shores with mild slopes Dally, Dean and Dalrymple (1984) developed a model for wave breaking and decay, analogous to a hydraulic jump, when waves enter into the surf zone. Isaacson (1985) linked this model of wave breaking with the damping factor and showed how the energy dissipation due to wave breaking can be taken into account. (2.5) 11 CHAPTER - HI THEORETICAL FORMULATION In this chapter, the mathematical problem of determining the transformed wave field is first outlined. Under the assumptions of linear wave theory and a gentle seabed slope the governing equation for the wave motion around objects rising from seabed or deep water may be approximated as a two -dimensional partial differential equation involving a potential function. This equation named the "mild slope' equation is similar to the Helmholtz equation. The solution of the mild slope equation provides wave heights at different locations which is the primary aim of the present work. Since no general analytical technique is available to solve the potential function appearing in the mild slope equation, a finite element discretization procedure is adopted. Thus the next step is to develop the finite element equations from the governing differential equation. The method of weighted residuals is the most general method for obtaining the approximate solutions to partial differential equations. One of the most widely used, Galerkin's method, will be used for this purpose. 3.1 PROBLEM DEFINITION The geometry describing the problem is as sketched in Figure 3.1. Here x is the horizontal coordinate measured in the direction of wave propagation, z is the vertical coordinate measured upwards from the still water level and the y axis lies on the undisturbed water surface, its direction being fixed by the right-hand rule. The 12 The incident wave train is assumed to maintain a permanent sinusoidal form and the further assumptions are that there is no underlying current, the free surface is uncontaminated, the fluid (water) is incompressible and inviscid and the flow is irrotational. The incident wave train has a height HQ, a wavelength X, a period T and a speed c. The angular frequency of waves is defined as CJ = 2ir/T and the wave number as K = 2n7X. The fluid motion may be described by a velocity potential $ which must satisfy the Laplace equation. In the three-dimensional Cartesian coordinate system, this is a 2 * a 2 * a 2 * n For an obstacle arising from deep water or seabed the total wave potential can be written in complex form as (Ref. Sarpkaya and Isaacson, 198 La) 13 where A = -igH /2u o cosh(/c(z + tfl) ^ A cosh*<f * ( * ' y ) c (3.2) i = y M d = Still water depth t = Time <t>(w) = Total potential function. Eq. (3.2) satisfies the boundary condition on the seabed which imposes a zero vertical component on the fluid particle velocity there ? * = 0 at z = -d (3.3) dz The linearized kinematic free surface boundary condition, which requires that fluid velocity normal to the free surface 17 is equal to the velocity of the free surface itself in that direction is: at s = 0 (3.4) dz dt The linearized dynamic free surface boundary condition, which requires that the pressure at the free surface, expressed in terms of Bernoulli's equation is constant, is: £ £ + y , = 0 at z = 0 (3.5) at v' Eqs. (3.4) and (3.5) may be combined to involve # and exclude 17: £ ! * + f f ^ i = 0 at * = 0 (3-6) dt2 vdz Eq. (3.2) also satisfies this combined boundary condition. Since the Laplace equation is a linear homogeneous partial differential equation, the solution for the total wave potential <t> can be obtained as a sum of two potentials, 14 $ = + $ s (3.7) where the incident wave potential $j is given by = A w ^ « » + d)) e[.(«-w«)l (3.8) 1 cosh(/ed) and the scattered potential $g as There are two remaining boundary conditions to be satisfied for a unique solution of the problem. The first one requires that there is no flow across the solid body boundary, and the second one is that the scattered potential should decay far away from the body in an appropriate way. Represented mathematically, the first one is — = 0 on body surface (3.10a) dn where n is the normal vector to the body surface and the second one, also known as Sommerfeld radiation condition, is Bm A (?p- - iK$s) = 0 (3.106) r—oo \ Or J where r is the outward radial direction. Substituting for the total and scattered potentials from Eqs. (3.8) and (3.9) the boundary conditions become — = 0 on body surface (3.11a) dn \lmri(?p--iKts)=0 (3.116) r—oo \ Or ) 15 For $ to satisfy Laplace's equation, <(> should satisfy the original mild slope equation as suggested by Berkhoff (1972) If energy dissipation is to be taken into account then <p should satisfy a modified version of the mild slope equation; which will be the basis of the present work. For a propagating wave train, the average rates of energy transfer across vertical planes a short distance dx apart and parallel to the wave crests can be equated to the energy dissipation occuring between those planes. Thus: where P is the energy flux, which is the average rate of energy transfer across a plane of constant x. D is the average rate of energy dissipation per unit area, E is the energy density, which is the average energy per unit horizontal area of the wave train and W is the energy dissipation factor. Booij (1981) showed that the additional term kjW<£ in the mild slope equation accounts for the wave energy dissipation and the mild slope equation becomes (3.12) dp dx = -D = -WE (3.13) (3.14) To summarize, the problem to be solved is the governing equation (3.15) subject to the boundary conditions on body surface (3.16a) (3.166) 16 3.2 SIMPLIFICATION OF THE PROBLEM The problem as stated applies over an infinite domain, but in order to obtain numerical solution the outer boundary is considered to be located at a finite distance sufficiently far away from the body where the radiation condition is approximately satisfied. Since any finite element problem calls for the solution of a large system of equations, it will be helpful if some symmetry can be used to reduce the size of the problem. y * X Fig. 3.2 Sketch of simplification. Bearing these two factors in mind, the problem is simplified in the following manner. Since the outer boundary can be of any shape a circular boundary at a finite distance has been chosen. Also, only geometries which are symmetrical about the x-z plane will be considered so that only the upper half of the problem has to be solved as shown in Figure 3.2. Only the upper half of the problem can be considered because in a potential irrotational flow the velocity across the line of flow symmetry must be zero. This makes the x axis part of the boundary along which the flow velocity in the y 17 direction is zero. Also for a semicircular outer boundary the normal vector will coincide with the outward radial direction. With reference to Figure 3.2, the simplified problem can be stated as follows. The governing equation is V(ec,V0) + (^u2 + iww) $ ss 0 in domain O (3.17) subject to the boundary conditions dn where is the body boundary, on Si d$ W o c — = — = 0 on S2 8y dn 2 (3.18a) (3.186) where S2 is the boundary formed by the x axis anddf/dy is the velocity across it, and d<f>s dr d$§ — \K$S = — tK<f>s = 0 on S3 (3.18c) where is the semicircular outer boundary. If attention is further restricted to geometries which are axisymmetric about the z axis, i.e. if depth is only a function of the radial distance from the centre of the island, a drastic reduction in the computational effort may be achieved. This may be brought about by applying a discretization in the form of ring elements followed by Fourier series expansions of variables with respect to the angular direction in a cylindrical coordinate system. The incident wave field is not axisymmtric and therefore has to be represented in terms of an expansion series of orthogonal functions. However, use of this technique will restrict the solution only to bodies of revolution and will not be able to handle real life situations where the depth may not be a function of the radial distance alone. So in order to retain the generality of the problem this approach is not pursued in the present work. 18 3.3 FINITE ELEMENT FORMULATION The method of weighted residuals involves two steps. In the first step the general functional behaviour is in some way approximated so that the governing differential equations and the boundary conditions are roughly satisfied. This approximation will therefore result in some error called a residual. This residual should vanish in some average sense over the solution domain in order to obtain a correct solution. The second step is to solve the equation or system of equations which arise from minimizing the residual, and thereby obtaining the solution for the unknown function (Huebner, 1975). Consider a differential equation expressed in a general form as Ifo) - / = 0 (3-19) in a domain and with a boundary S, where L is the differential operator and f is a known function of the independent variables. First the exact solution <t> is approximated by 0 so that m 7 ( 3 - 2 ° ) where Nj = Assumed functions $ I = The unknown variables m = Some positive integer. Substitution of 0 into Eq. (3.19) results in a residual, R: £(?)-/ = * (3-21) In the method of weighted residuals, the m unknowns 0^ are determined such that the residual or error R over the whole solution domain is small. This is achieved by producing a weighted average of the error and making this weighted average vanish 19 over the entire domain. If m linearly independent weighting functions (a) are chosen, then Jn n There are many ways to select suitable weighting functions, which give rise to the variety of weighted residual techniques. In the Galerkin method the weighting functions are taken to be the same as the approximating functions used to represent <f>, i.e. ai = Ni for Z=l,...,m (3.23) Substituting Eq. (3.23) into Eq. (3.22) gives rise to the Galerkin formulation [ [xtf) -/] JVn = o 1=1 m (3-24) The above equation is true for the whole domain; but Eq. (3.19) holds for any point or collection of points in the solution domain defining any arbitrary subdomain or element of the whole domain. For this reason a local approximation for Eq. (3.24) can be defined and it. will be valid for one element at a time. Hence for the Galerkin method the equations governing the behaviour of an element can be written as f [ l ( ^ e ) - / e ] ^ e d n = 0 (3-25) where the superscript V restricts the range to one element For the problem under consideration we have from Eq. (3.17) L = V{ccaV) + (%2 + iuW) / = 0 Substitution in Eq. (3.25) and dropping of the superscript V for convenience results in 20 «o *• c j Now V where i , j are unit vectors in the x, y directions respectively. Therefore so that and vKv»-(iJ +ij)(^ } + «,g,) Hence the Galerkin formulation can be written as Integration of the first term in the equation by parts gives where n x is the x component of the unit vector normal to the boundary and dS is the differential arc length along the boundary. Similarly /X(<4)*W.«*g***-/ 0*8$** (3-28) 21 Substitution of Eqs. (3.27) and (3.28) in Eq. (3.26) and removing the negative sign results in Now in the line integral of Eq. (3.29) dS d<b (dS*. d<b~\ ( » *.\ For elements entirely within the whole domain the line integral will be zero. For the elements on the boundary S, where S = Si + 5^ + S$ the line integral is also equal to zero along and S2 from Eqs.(3.18a) and (3.18b) respectively. So the line integral in Eq. (3.29) reduces to 'Si But from Eq. (3.7) dn dn dn and substitution of ^^s/gn from Eq. (3.18c) gives Substituing in Eq. (3.29) the total integral is obtained as - J #1 (^ iK* ~ d3 = 0 ( 33°) From Eq. (3.20) 22 m where the symbol { } denotes a column vector, it is evident that £0 So substitution in Eq. (3.30) results in [ / 0 h ( { ^ r § + { f } T f ) - M ^ + - ) H dxdy 9t) - I [jf cc„ JV, K {Njf rfaj = j f « , JV, (|£ - (fa (3.31) or in a brief matrix form the finite element equation for an element can be written as [K\eMe = {qi}e (3.32) This is similar to the structural mechanics stiffness formulation where [ K] is like the stiffness matrix and {q^ } is the load vector. In the present case 1*1 = M + [Ka] + [K3] where \K 1 f T ( d N i dNi dNi 9N,\ ca , 1 [Klb=Ja htei; +i?i?)--y Ni N'\ d x d* ^ [K2\ji = -* / ccg KNJ Nt da (Jfjlj, = - » f uWNj Nt dx dy Jn (3.336) (3.33c) (3.33d) 23 These element matricies and vectors can be assembled into a global system by considering each and every element covering the whole domain and can be solved by using any standard technique. 3.4 DETERMINATION OF THE ENERGY DISSIPATION FACTOR Using the study of energy dissipation and wave height decay of shoaling waves breaking on a plane beach carried out by Dally, Dean and Dalrymple (1984), Isaacson (1985) showed that the energy dissipation factor for wave breaking can be obtained from the relation where T, stable wave factor and K', wave decay factor are empirical coefficients. From a study of experimental data Dally et al. found K' to lie between 0.10 and 0.28 for beach slopes within the range 1/80 to 1/30, whereas the range of T is from 0.35 to 0.48. The application of Eq. (3.34) is also dependent upon the conditions for the onset and cessation of wave breaking. As a simple approach the onset of wave breaking is considered to start when the H/d ratio increases over a specified value of 0.78 and cessation of breaking takes place when H/d becomes less than T. The numerical values of the coefficients are for gentle slopes and plane beaches where reflection and scattering are not present Thus these values may not be applicable to artificial islands which involve significant wave reflection and scattering. Different geometries, steeper slopes and the presence of reflection and diffraction can all contribute to different extents of energy dissipation and consequently to changes in the empirical coefficients and the wave breaking condition. Since no experimental or theoretical work is available at present to indicate the way these changes may take place, the original numerical coefficients for the plane beach are used unaltered in the (3.34) 24 present work. 3.5 WAVE HEIGHTS OF THE TRANSFORMED FIELD Once the numerical values of the potential function tf>(x,y) at different locations are known from the solution of the assembly of equations (3.33a) to (3.33d), the total potential can be found by using Eq. (3.2). The wave elevation r\ can be obtained from the linearized dynamic free surface boundary condition given by Eq. (3.5) n = JJA at z = 0 (3-35) * gdt In the present work, all results are obtained at a time when t equals zero for comparison with other available results. At that time the wave elevation becomes •7 = R e [ ^ ( x , y ) ] (3-36) The nondimensional runup R around the island is defined by J * = M = jAb8tf(x > y )] (3-37) •Ho and the relative amplitude A r is defined as The phase angle is calculated from the ratio of the real and imaginary parts of the potential function. 25 CHAPTER - IV NUMERICAL PROCEDURE The development of a numerical procedure for the finite element formulation involves several steps. This chapter deals with various aspects of the computation. The type of element used, the numerical integration technique employed within the element domain, the approach used to determine the height of breaking waves, the arrangement and decomposition of the global assembled matrix to fit available solution routines at UBC, the type of solution routine chosen and finally the organization of the program are discussed in the following sections. 4.1 THE ELEMENT In order to model the curved boundaries effectively with fewer elements, the eight noded quadratic isoparametric element has been adopted. An isoparametric element implies that the curved elements in the (x,y) coordinate system can be mapped into a rectangular element in the (s,t) coordinate system as shown in Figure 4.1 and that parametrically the elements in both systems are equivalent Details of isoparametric elements can be found in Zienkiewicz (1981). The coordinate transformation is achieved from the relations and (4.1) where / takes value from 1 to 8 and the shape functions N. are given by 26 * 8 \ 1 (-1.1) (1.D 4 7 8 3 6 i 1 [5 2 (-1.-1) (1.-D //g 4.1 The mapping of the element. Nt = ]{1 + S3l) (1 + ttt) (S3, + « , - 1) 4 (4.2) The variable 0 is also obtained from its nodal values by the same interpolating functions 4> = £ / V ^ j (4.3) Since the shape function derivatives and the integrations in Eqs. (3.33a) to (3.33d) refer to the coordinate system (x,y) they have to be converted initially to the (s,t) coordinate system in order to carry out the numerical scheme. From the rules of partial differentiation it can be shown that (Segerlind, 1976) dN, dN, = [J\ i t ) (4.4) It dt where [J] is the Jacobian matrix defined as [J] = and for a two-dimensional surface element dxdy = dei[J]dsdt dx a ~5t 9 (4.5) (4.6) 27 4.2 INTEGRATION TECHNIQUE Since the Jacobian matrix [ J ] is generally a function of s and t, it cannot be readily evaluated because the coefficients are polynomials. So [J]~* is never determined explicitly. Instead the element matrices are evaluated using numerical integration techniques. The numerical evaluation of an integral can proceed in one of two basic ways. In the first, the function (F) is evaluated at equally spaced points, which gives rise to the trapezoidal and Simpson's rules. The second approach is to locate the sampling points in such a way as to achieve the best accuracy. This process is known as Gaussian quadrature and is usually employed in finite element analysis because it requires fewer sampling points than the first method for the same degree of accuracy. In Gaussian quadrature with n sampling points, a polynomial of degree (2n-l) can be constructed and integrated exactly. If this is satisfied then the rate of convergence which can be attained by exact integration will be preserved (Zienkiewicz, 1981). Since in Eqs. (3.28a) to (3.28d) the order of the function (F) for the quadratic 4 4 polynomials can be upto (s ) or (t ) because of the term N.N, is the minimum requirement Since n must be an integer, 3 sampling points in each direction are selected in order to avoid loss in convergence. So the integration in the range (-1,1) in each direction in the numerical form will be where Q is the weighting function at the quadrature point The locations and weighting functions for the three points are obtained from Cook (1981) and are as follows: 2n - 1 = 4 or n = 2.5 (4.7) 28 LOCATION WEIGHT ± 0.77459 66692 41483 0.55555 55555 55556 0.00000 00000 00000 0.88888 88888 88889 4.3 THE RANGE OF WAVE BREAKING As mentioned earlier, the incipient wave breaking criterion required in the calculation of W has been taken as H/d = 0.78. However, the wave height H is obtained as part of the solution and is not known beforehand in order to establish the onset of wave breaking. As an initial solution procedure, the increase in wave height over the berm, used only to obtain W, is assumed to be only due to shoaling. From this increased height the range over which the waves break is estimated and then the contribution to the system matrix due to energy dissipation is calculated. Actually the solution technique should be iterative, with the initial solution for wave heights at different locations tested for wave breaking, and a refined contribution from the dissipation factor then calculated. This process should be repeated until the solution converges. However, as an initial procedure, such an iterative approach is not pursued in the present work and the calculations are stopped after obtaining the first set of results. In order to obtain wave heights for use in the wave breaking criterion, the wave height H over the berm up to the point of breaking is taken as due to shoaling only and is given by: (e.g. Sarpkaya and Isaacson, 1981b). Once the waves start breaking, the heights of the decaying waves are calculated by using the expression given by Dally, Dean and Dalrymple (1984) for a plane beach of H f 2cosh 2 («fl V (4.8) H0 2nd + 8inh(2<cd) 29 uniform slope: • H (4.9) where a (4.10) u = The berm slope H^ , = Incipient breaking height and depth respectively. It is to be noted that for K' /m = 5/2 the expression is not valid. For this special case the decay in wave height is given instead Thus the height of the breaking waves can be determined from the Eqs. (4.8) to' (4.12) and the energy dissipation factor can be obtained from Eq. (3.29) within the range specified in Section 3.4. 4.4 ARRANGEMENT OF GLOBAL MATRIX Each of the element matricies occupies a unique position in the global matrix. The way to locate the element matricies in the global system is to identify the global node numbers to the corresponding internal node numbers. For example Figure 4.2 shows the numbering sequence for a typical element, with the internal node numbers written inside and the corresponding global node numbers written outside the element (4.11) where (4.12) Fig. 4.2 The local and global node numbering sequence. Then writing the global node numbers against the element node numbers for the local matrix addresses, it can be seen 74 72 55 57 73 61 56 62 1 2 3 4 5 6 7 8 74 1 "#n # 1 2 #13 #H #16 #17 #18 72 2 K2l #22 #23 — — — — #28 62 8 _#8i #82 #83 — — — — #28 that the element K^g corresponds to a global location of (74,62); element corresponds to the location (72,72) and so on. Repeating this for every element in the domain the final global matrix and the column vector may be set up. 4.5 THE SOLUTION ROUTINE The complex global matrix [ K] given by Eq. (3.32) contains many zeros and it is therefore inefficient to solve the matrix equation using a general solution routine for the whole system. Since available solution routines at UBC do not include any to 31 solve a sparse complex matrix system, the global complex matrix of size (n X n) is decomposed into real and imaginary parts, and the resulting real matrix system has a size (2n X 2n) where n is the total number of nodes. Consider a complex matrix equation W + iV\ ($r + ifi) = {C + iD} (4.13) where [U + iV\ = Complex Global Matrix [C + \D\ = Compkx Column Vector ^ r = Real Part of Solution = Imaginary Part of Solution Eq. (4.13) can be written as: \U] ffr}+i[V] {?,} + »•[(/] {?.•} - [V] {?,•} = {C} + i{D} Equating the real and imaginary parts, one obtains W\ { W - M ft) = {c} or written in the matrix form i v - ; ] { £ } - « } This increases the matrix size and makes it unsymmetric, but on the other hand a very convenient and efficient routine UBC:SPARSPAK is available for solving large linear sparse system of equations. This routine exploits all the zeros in the global matrix system, with only the non-zero elements of the matrix specified for ordering before the actual solution starts. This routine has the option of using a nested dissection method for an unsymmetric matrix, specifically designed for finite element problems. For large problems the reduced factorization time more than compensates for the time required for the solution of a full (n X n) complex matrix (UBC SPARSPAK, 1982). Once the real and imaginary parts have been obtained they can be 32 put back in the complex form to get the complete nodal solution. 4.6 ORGANIZATION OF THE PROGRAM The finite element program developed to solve the present problem uses a number of subroutines and the manner in which they are arranged is shown in the flowchart in Figure 4.3. The functions of the various subroutines are described below. LAYOUT : This routine reads in all the data required for the problem; viz the wave height and period, island geometry, the nodal coordinates and element information. SPSSLV : This is an initialization routine for the sparse matrix solver; UBC:SPARSPAK. ORDRA6 : This is another routine for UBC:SPARSPAK which orders and allocates the storage for the unsymmetric global matrix. STIFF : This subroutine calculates the real part of the element matrix as given in . Eq. (3.33a). BREAK : This subroutine finds the range of wave breaking and calculates the imaginary part of the element matrix due to wave breaking as in Eq. (3.33c). BOUND2 : This subroutine calculates the imaginary part of the element matrix as given in Eq. (3.33b) and the second part of the column vector of Eq. (3.33d). BOUND1 : This subroutine calculates the first part of the column vector as in Eq. (3.28d). SOLV6 : It is a routine of UBC:SPARSPAK that solves an unsymmetric sparse system of linear equations' by nested dissection method. EXPAND : This routine puts together the real and imaginary part of the solution of LAYOUT r SPSSLV t ORDRA6 E-I BREAK REPEAT FOR r ALL ELEMENTS I V BOUND2 BOUND1 _y SOLV6 EXPAND GRAPH > STOP 34 the potential function and obtains the wave height, elevation and phase at every node. GRAPH : This subroutine arranges the nodal data in a form so that a three dimensional plot of wave elevation and height can be obtained when used in conjunction with the plot routines of UBC Civil Engg. Graphics Lab. WAVE : This subroutine is called at various stages to find out the phase and group velocities, the wave number and the incident velocity potential at different locations. 35 CHAPTER - V RESULTS AND DISCUSSION Before attempting to obtain results for a specific problem, any finite element program should first be verified by running it for problems for which previous results are known. In addition the program should also be tested for convergence. For the present program, two cases of independent analytical solutions are known where the effect of wave breaking is not taken into account. One of these corresponds to a circular island with parabolic shoal and the other one corresponds to a vertical cylinder on a horizontal seabed. The sections in this chapter present results for these known cases and also the results of convergence. Results are then presented for a typical artificial island with realistic geometry and taking account of the effects of wave breaking. Since no results are available in the existing literature for comparison with results when wave breaking is taken into account, a few preliminary experiments have been performed with a scale model of an island. The experimental results are compared with the numerical predictions, both with wave breaking taken into account as well as neglected. One section of this chapter describes the experimental set up. Various aspects of the experimental results and the limitations of the theory and the program are discussed in the last section. 36 5.1 VERIFICATION AGAINST KNOWN SOLUTIONS For a circular island with parabolic shoal as shown in Figure 5.1, relative wave amplitudes have been calculated by Houston (1981) on the basis of an analytical solution developed by Homma (1950). For verification, the program developed here has been tested for a long wave of period 720 seconds. The corresponding wavelength is 142600 meters. The radiation boundary is taken to be located at a distance 95000 meters away from the centre of the island. A total of 135 elements with 454 nodes have been used for the solution. The element mesh is shown in Figure 5.2 where the angular interval is 12°. The results of Houston are compared with those obtained from the present program in Figure 5.3. Vastano and Reid (1967) obtained the phase lag for the same problem using a Finite difference scheme. Their results are also compared with the phase lag obtained from the present program in Figure 5.3. The comparisons indicate that the present program predicts the relative amplitude and phase for this case remarkably well. The other analytical solution, known as the MacCamy Fuchs solution, is for a vertical circular cylinder rising from a horizontal seabed as shown in Figure 5.4. The nondimensional analytical solutions for runup and relative wave amplitude are obtained from Sarpkaya and Isaacson (1981) for two cases »ca = 0.6 and /ca = 1.0, where a is the radius of the cylinder. For the dimensions shown in Figure 5.4 these conditions correspond to wave periods of 26.44 seconds and 16.38 seconds respectively. In Figure 5.5 the analytical solution for relative amplitude and phase lag is compared with the finite element solution at a radial distance 75 meters away from the centre of the cylinder for the case Ka = 0.6. The number of elements used in the finite element solution is 72 and the radiation boundary is taken at a distance 295 meters away - from the centre of the cylinder. Results for the shorter wave period of 16.38 seconds 37 y Fig. 5.1 Geometry of the island with parabolic shoal. Fig. 5.2 Finite element mesh for the parabolic island. R B L A T T V E A M P L I T U D E 0 - 9 - C PR8SKNT SOLUTION AJ1ALTTICAL SOLUTION P H A S E L A C * ' t PRESENT SOLUTION X - x - K FINITE DIF7. SOLUTION — r -120 00 00 ANGLE IN DEGREES 180 130 30 Fig. 5.3 Relative amplitude and phase lag around parabolic island. y INCIDENT WAVE • • o J ' 50m-i 2 -y yyy. s / y y yA — r 10m y y^y y s Fig. 5.4 Sketch of a circular cylinder on horizontal seabed. 39 120 SO 00 ANGLE IN DEGREES Fig. 5.5 Relative amplitude and phase lag for a vertical cylinder, KO .= 0.6. " S 3 S3 RUNUP O - e - 0 PRH3BHT SOLUTION I MacCAUY-PUCHS 30LH PHASE L A G PRI3BNT SOLUTION Zl X - * - K M . c C A M Y - P U C H S S O U w a 3C 1 ISO 130 120 BO 00 ANGLE IN DEGREES 30 Fig. 5.6 Runup and phase lag for a vertical cylinder, ica = 1.0. 40 are plotted in Figure 5.6. The runup is shown at a radial distance 25 meters from the centre of the cylinder, i.e. on the cylinder wall. In this case the number of elements is 100 and the radiation boundary is located 225 meters away from the centre. Since no wave breaking is taken'into account, the solutions will not change for different wave heights. As in the case of a cylinder with a parabolic shoal, the present results agree quite well with the established results. 5.2 CONVERGENCE TESTS Since the continuity of the variable 0 is maintained at the element boundaries and also inside the eight noded isoparametric elements and since completeness is assured when the derivatives can achieve arbitrary constant values, convergence to the accurate solution follows. The accuracy will depend upon parameters such as the incident wave height and period, the location of the outer boundary and the degree of refinement in the mesh size. However, an elaborate mathematical analysis to determine the order of convergence and error bound is omitted because the general opinion among finite element users is that they are given in terms unknown *a prion" (Zienkiewicz, 1981) and can also be quite complicated. Therefore a qualitative study of the solution accuracy as a function of the mesh size and the location of the outer boundary has been carried out for a particular geometry and wave period. The vertical circular cylinder of Figure 5.4 is chosen as the object and the wave period is taken to be 26.44 seconds for which the analytical solution for runup is known. Figure 5.7 shows the effect of varying mesh size for a constant outer boundary and Figure 5.8 displays the effect produced by shifting the outer boundary for a constant number of elements. It shows definite improvement by placing the outer boundary further away and by refining the mesh, but no definite mathematical relationship can be developed because the effect of 41 180 130 1Z0 SO 60 30 0 ANGLE IN DEGREES Fig. 5.7 Accuracy of runup with mesh size for a vertical cylinder, tea = 0.6, outer boundary — 295m. CD 180 130 120 90 00 30 0 A N G L E IN DEGREES Fig. 5.8 Accuracy of runup with boundary distance for a vertical cylinder, tea = 0.6, elements = 45. 42 moving the boundary further out improves the solution in itself, but has the effect of making the elements larger, which has a detrimental effect. Thus as a general rule it is desireable to place the outer boundary as far away as possible and to incorporate as many elements as possible within the practical limits of the computer capacity and scope of the problem. 5.3 RESULTS FOR A TYPICAL ISLAND A sketch of a typical artificial island with vertical sides and sloping berm is shown in Figure 5.9. Results of the present method for this case are shown in Figure 5.10 and 5.11 for wave periods of 10 seconds and 20 seconds respectively. Different wave heights have been considered and the results show the effect of wave heights when wave breaking is taking place along with the case when wave breaking is not considered. It is evident that the breaking of waves has a damping effect on the runup around the island and the effect increases with increasing wave heights- as expected. However, this does not indicate how realistic the breaking model is because, as mentioned before, no data are available in the existing literature to compare with the results. Thus, a few sets of preliminary experiments have been conducted to check the predictions of the computer model. 5.4 EXPERIMENTAL SET UP In order to provide a comparison of the numericalcal predictions, a few sets of experiments have been carried out on a geometrically scaled-down model of the island shown in Figure 5.12. The model is made of aluminium sheet and has a length scale factor of 1:100. By dimensional analysis it can be shown that the wave periods will be reduced by a factor 1:10. Apart from the runup, wave heights were measured at 43 — 1 1 1 1 1 1 180 150 120 80 80 30 0 ANGLE IN DEGREES Fig. 5.10 The runup profile around island, wave period = 10 seconds. RUNUP G-e-G WITHOUT BREAKING A A A BREAKING, WAVE HT-4.37m < ' ' BREAKING, WAVB HT-7.82n» 1 % •to J \ * X * ( I » % * % « 1 % i * # # i « » # —Q— « t t ^ \ > \ i 180 190 120 00 60 30 ANGLE IN DEGREES Fig. 5.11 The runup profile around island, wave period = 20 seconds. Jf 1.52m INCIDENT WAVE I " > SLOPE 1:3.6 7* Fig. 5.12 Sketch of the experimental model. 45 seven different locations along radials 30° apart and both 0.15 meters inside and outside the berm as shown in Figure 5.12. The incident wave heights were measured at a point sufficiently far away from the island such that the measurement could be completed before diffracted waves arrived at the measurement location. The experiments were performed in a wave basin which is 13.72 meters long and 4.45 meters wide inside. Waves were generated at one end by a paddle hinged at the channel bottom and driven by a variable speed control through a crank of variable stroke. In order to obtain the required scaled down water depths, the basin was fitted with a raised floor except in the region near the paddle where the floor was inclined to meet the true bottom. A wave filter was located near the paddle and ensured reasonably smooth incident waves. In order to minimize wave reflections, the side walls were covered and the far end of the basin was fitted with a sloping beach made of artificial hair matting. For each test the wave heights at the marked locations were measured using a hook and pointer gage. For runup measurement, half the vertical side of the model was wrapped with a paper sleeve so that the wetted region could be easily distinguished and the line defining this could be marked in. However, no multichannel electronic wave recorder to measure the wave elevations at different points simultaneously was available and therefore information about wave elevations and phase angles could not be obtained. Despite the limitations, the experiments were considered justified in order to provide the basis for comparing the theoretical runup and relative wave amplitude predictions. Profiles of runup and wave heights at marked points were obtained for various conditions and selected results are reproduced together with the corresponding predictions for cases of breaking both taken into account and neglected. Figures 5.13 and 5.14 show the results for runup of waves of period 0.7 seconds and 1.4 seconds 180 ISO 120 BO 60 ANGLE IN DEGREES 30 Fig. 5.13 Comparison of runups around model for wave period = 0.7 seconds. 180 ?1 r f % y i i 1? f T f i RUNUP f % 0 - 0 - 0 WITHOUT BREAXINC t I I BREAXINC. W A V E HT=»0.07m A A d EXPERIMENTAL RESULTS i \ ; \ i V i f \ j t % V ' i i t § i i ! \ * 1 1 ^ 1 % /<\ 5 < i * « f t % ' \ • 9 " i 1 •' » i 1 # « * 9 t 1 • i V » $ 1 I \ % t I T / i t — * — H — V r i r » 130 120 BO ANGLE IN DEGREES 60 30 Fig. 5.14 Comparison of runups around model for wave period = 1.4 seconds. 47 3* U o I k t 1 1 1 1 A i f t 1 1 i » .* 1 » 1 1 f 4 /'. ; • t f I t 1 V J ', T I i f 1 1 1 \ • i • i 1 1 ; t $ 1 Y I i \ i i • • i • ! i i • V, i i .' « i i - * * * % • « • i • • • i ! > \ V • < » ; • \ i • t t < t '*>! 9 i i A TTV IDT n n i n t »' J 0 - 9 - 0 1 +—f ! mTHOUT BREAKING JRKAKING, WAVE HT-0.03 ECP8R1MBNTAL R83ULTS A—A ] A r—7 1 1 1 1 1 1 1 180 130 120 00 00 ANGLE IN DEGREES 30 Fig. 5.15 Comparison of relative wave amplitudes at 0.15m inside the model berm, wave period = 0.7 seconds. 120 00 00 ANGLE IN DEGREES Fig. 5.16 Comparison of relative wave amplitudes at 0.15m outside the model berm, wave period = 1.4 seconds. Incident wave ht = 0.05m, Period = 0.7sec, Elements = 396 Fig. 5.1? Wave height distribution around mode] (irith breaking/. Oo Incident wave ht = 0.05m, Period = 0.7sec, Elements = 396 Incident wave ht = 0.05m, Period = 0.7sec, Elements = 396 52 respectively. The incident wave height in the first case is 0.050 meters and the second case 0.045 meters. The uniform water depths for both cases were maintained at 0.10 meters. Figure 5.15 shows the relative wave amplitudes at 0.15 meters inside the berm for the lower period and Figure 5.16 shows the same at 0.15 meters outside the berm for the higher period. In both the cases a total of 396 elements were used with 1269 nodes for the computer prediction and the radiation boundary was taken to be 2.5 meters away from the centre of the model. Figures 5.17 and 5.18 show the overall view of the wave height distributions around the model with and without wave breaking being taken into account. The wave period is 0.7 seconds. Figures 5.19 and 5.20 show the overall wave elevations for both the cases at time t equal to zero. The waves tend to bend around the island and similar patterns were observed during the model testing. Though the agreement between the theoretical and experimental results are not excellent, it certainly shows that the results when wave breaking is considered are more realistic. There are quite a few reasons which may lead to the difference between the observed and predicted values and some of these are discussed in the next section. 5.5 DISCUSSION i) The wave lengths are dependent upon wave periods and for lesser periods the wave lengths get shorter. It can also be seen from the results that the variation in the runup profile is greater for waves of shorter periods than that of longer ones. Hence, to maintain the accuracy the element grid should be finer for shorter wave periods so that the rapid variation in the wave field with position can be correctly produced. ii) Since the wave properties change with location, it is desirable to have sufficient number of elements within one wave length to take into account the changes 53 correctly. For the same island geometry, more elements are required when the analysis is being carried out for shorter wave lengths. The wave lengths also decrease with a reduction in water depth and therefore elements should be concentrated above the berm in order to model the wave Tield over it accurately. There is no unique mathematical relationship to determine the element size from the wave length but from experience investigators such as Chen and Mei (1974) have found that the ratio of the maximum element dimension to wavelength should generally be less than about 0.1 in order to obtain a reasonably accurate solution. iii) If the slope of the berm is not very mild, then partial reflection from the berm, which is not taken into account in the present formulation, may become significant Booij (1983) investigated the accuracy of the mild slope equation for a region of constant seabed slope separating regions of constant and different depths. He found that for waves propagating in a direction parallel to the depth contours, the mild slope equation produces accurate results for bottom slopes up to about 1:1. For waves propagating normal to the depth contours the mild slope equation is less accurate but can still be used for bottom slopes of up to about 1:3, where the reflection coefficient is about 4%. However, this investigation was restricted to a channel where the slope is varying only in one direction and the water depth is constant before and after the slope. The effect of slope changes in two directions and the absence of a constant depth region on one side of the sloping region, as encountered in artificial islands, is little known and the range of slopes for which results are accurate may not be as high as in the case he investigated. iv) The validation of the numerical values of wave decay factor K' and the stable wave factor T for prototype conditions is still far from being complete. Dally, Dean and Dalrymple (1984) suggest that their recommended values, obtained from 54 laboratory tests on beach slopes varying from 1/30 to 1/80, should only be taken as a guide in unknown situations rather than as authentic values. They also suggest caution when applying these data to steeper slopes. Moreover, the validity of the values which were obtained for plane 6eaches should be verified by experimental data measured for circular shoals. v) Although the selected incipient wave breaking criterion, based upon H/d = 0.78, is common in coastal engineering practice, Walker and Headland (1982) have shown that it often leads to an underestimation of the breaking wave height when compared to the breaking height and depth as determined by empirical curves presented in the Shore Protection Manual (1977). The graphs presented by Weggel (1972) show the breaking height as a function of the bottom slope, wave period and depending upon them the incipient breaking criterion in terms of H/d can vary over the range of 0.63 to 1.67. The breaking criterion is also dependent upon the type of breaker, plunging or spilling type, generated. Apart from this, for the artificial islands wave reflection from vertical walls and refraction significantly change the breaking behaviour, which will generally be much more complicated. In the light of the above, the criterion H/d = 0.78 has the merit of being simple but should only be used as a first trial and in an average sense. Moreover, for the circular island refraction can be another important factor in determining changes in wave height. The effect of refraction can be obtained by a method of ray tracing over the circular shoal. vi) The problem of wave transformation around artificial islands is a nonlinear problem if treated in a proper perspective. Small amplitude linear wave theory is usually not valid in shallow water near the range of incipient breaking. In addition, the energy dissipation factor has been considered to be a function of position only, which is not strictly true and thus can be a source of another nonlinearity. However, 55 no analytical solution exists to treat the nonlinearities to the full extent although the effect of nonlinear terms can be quite significant. vii) The difference between the measured and theoretical runup emphasises the need for an improved wave breaking model. In order to develop this, a series of elaborate experiments are needed in order to modify the wave decay and stable wave factors as well as the incipient breaking criterion. A reasonable cause of the difference can also be attributed to the increasing nonlinear behaviour of the waves in shallow water. On the experimental side, the ratio of the wave tank width and the base diameter of the island is only 1.973 and may not be sufficient for neglecting the side wall effects on the measured results. Considering all these aspects the difference between the experimental and theoretical results seems to be quite natural. viii) Although the Finite element method developed here has been applied to a simple island configuration, it can also be applied directly to cases where the slope of the berm may not be mathematically defined or for islands of other practical shapes such as hexagons, octagons etc. Thus the finite element method has the capability to handle a general class of problems where analytical methods fail and this makes the finite element technique a very powerful tool of analysis. 56 CHAPTER - VI CONCLUSION AND RECOMENDATIONS Although the calibration of the present numerical model is still not perfect and the results show discrepencies with laboratory observations, it can be safely concluded that wave breaking has a significant effect on the transformed wave field around the artificial islands and should be investigated in more detail in order to make accurate predictions of runup. The physical observations obtained so far can at best be described as indicative, and more thorough observations at a greater number of points, with systematic changes in parameters such as wave height, water depth, berm slope and wave period, are needed. Information about phase lag and instanteneous wave elevations should also be obtained and compared carefully before drawing any definite conclusions. One direct application of the present study in the design of the artificial islands is in determining the freeboard so that overtopping onto the work area can be avoided or kept to a minimum. Both experimental and theoretical studies are recomended for other practical island shapes such as octagons, hexagons, etc. There are also other design aspects as discussed below in which the present work can be helpful. i) Once the total velocity potentials are obtained at different nodes the water particle velocities at any depth and location can be found by calculating the corresponding space derivatives. Knowing the velocities, Bernoulli's equation can also be applied to yield the pressure at any location. If the pressures at different portions of 57 a vertical wall of the island are integrated over the whole area, the total force and moment acting upon the island for a particular wave condition may be determined. Calculation of forces and moments are important from the point of view of structral design and a realistic transformed wave "field helps to predict better design loads. ii) Pressures calculated over the berm are also required in the slope stability analysis of the berm, carried out by geotechnical engineers. Repeated wave motion and corresponding induced pressures increase the pore water pressure in the soil and reduce the effective shear stress. When this falls below a certain value the soil is no longer able to retain the slope and failure takes place. Thus a realistic pressure distribution is also important from the foundation engineering aspect iii) A third important design consideration pertains to the wave induced erosion around the berm. From the velocity potential the wave induced flow velocities may be obtained. Taking account of longshore currents, the sediment transport rate across any section around the island in the surf zone may then be determined by using a suitable sediment transport model such as that described in the Shore Protection Manual (1977) or such as that developed by Ackers and White (1973). After estimating the sediment transport rate, the erosion or accretion at any one section or at any location can be determined by using a one or two dimensional continuity equation. Overall it can be inferred that the present model has a wide range of possibilities in design application and can lead to better design criteria. The present work is only the first step in that direction and future work to overcome the present deficiencies will finally make it successful. 58 REFERENCES 1. Ackers, P. and White, W.R. (1973). Sediment Transport: New Approach and Analysis. Journal of Hydraulics Engineering, ASCE, Vol. 99, No. HY11, pp. 2041-2060. 2. Arthur, R. S. (1946). Refraction of Water Waves by Islands and Shoals with Circular Bottom Contours. Transactions, American Geophysical Union, Vol. 27, No II, pp. 168-177. 3. Berkhoff, J. C. W. (1972). Computation of Combined Refraction Diffraction. Proceedings, 13th International Coastal Engineering Conference, Ch. 24, pp. 471-490. 4. Bettess, P. and Zienkiewicz, O. C. (1977). Diffraction and Refraction of Surface Waves Using Finite and Infinite Elements. International Journal for Numerical Methods in Engineering, Vol. II, pp. 1271-1290. 5. Booij, N. (1981). Gravity waves on Water with Nonuniform Depth and Current Ph.D Thesis, Technical University of Delft, at Delft, The Netherlands. 6. Booij, N. (1983). A Note on the Accuracy of the Mild-Slope Equation. Coastal Engineering, Vol. 7, pp. 191-203. 7. Chen, H. S. and Mei, C. C. (1974). Oscillations and Wave Forces on a Man Made Harbour in the Open Sea. Massachusetts Inst of Technology, Report No. 190, pp. 573-595. 8. Cook, R. D. (1981). Concepts and Applications of Finite Element Analysis. 2nd Ed. John Wiley and Sons. 9. Dally, W. R., Dean, R. G. and Dalrymple, R. A (1984). Modelling Wave Transformation in the Surf Zone. Report No. CERC-84-8, US Army Corps of Engineers. 10. Dalrymple. R. A, Kirby, J. T. and Hwang, P. A (1984). Wave Diffraction due 59 to Areas of Energy Dissipation. Journal of Waterway, Port, Coastal and Ocean Engineering, ASCE, Vol. 110, No. 1, pp. 67-79. 11. Homma, S. (1950). On the Behaviour of Seismic Sea Waves Around Circular Island. The Geophysical Magazine, XXI, pp. 109-208. 12. Houston, J. R. (1981). Combined Refraction and Diffraction of Short Waves using the Finite Element Method. Applied Ocean Research, Vol. 3, No. 4, pp. 163-170. 13. Huebner, K. H. (1975). The Finite Element Method for Engineers. 1st. Ed, John Wiley and Sons. 14. Isaacson, M. (1985). Wave Transformation in the Coastal Zone. 1985 Ausgtralasian Conference on Coastal and Ocean Engineering, Vol. II, pp. 25-32. 15. Jong, J. J. A. (1976). Using Man Made Islands as Drilling Platforms in Arctic Seas. Ocean Industry, October, pp. 100-109. 16. Kirby, J. T. (1984). A Note on Surface Wave-Current Interaction Over Slowly Varying Topography. Journal of Geophysical Research, Vol. 89, No. CI, pp. 745- 747. 17. Kirby, J. T. and Dalrymple, R. A. (1983). A Parabolic Equation for the Combined Refraction-Diffraction of Stokes Waves by Mildly Varying Topography. Journal of Fluid Mechanics, Vol. 136, pp. 453-466. 18. Lautenbacher, C. C. (1970). Gravity Wave Refraction by Islands. Journal of Fluid Mechanics, Vol. 41, Part 3, pp. 655-672. 19. Liu, P. L F. (1983). Wave-Current Interactions on a Slowly Varying Topography. Journal of Geophysical Research, Vol. 88, No. C7, pp. 4421-4426. 20. Lozano, C. J. and Liu, P. L F. (1980). Refraction-Diffraction Model for Linear Surface Gravity Waves. Journal of Fluid Mechanics. Vol. 101, pp. 705-720. 21. Ocean Industry. (1985). Arctic Operations with Mobile Steel Caisson Rig. June, p. 39-41. 22. Offshore. (1985). Arctic Technology Awaiting Big Discovery. Feb, pp. 41-43. 23. Offshore Engineer. (1985). US Study Reviews Beaufort Needs. April, pp. 108-116. 60 24. Radder, A. C. (1979). On the Parabolic Equation Method for Water Wave Propagation. Journal of Fluid Mechanics. Vol. 95, pp. 159-176. 25. Sarpkaya. T. and Isaacson, M. (1981). Mechanics of Wave Forces on Offshore Structures. 1st. Ed, Van Nostrand Reinhold Company, (a) pp. 258-260. (b) pp. 153-155. 26. Segerlind, L J. (1976). Applied Finite Element Analysis. 1st Ed, John Wiley and Sons. 27. Shore Protection Manual. (1977). Coastal Engineering Research Center, US Army Corps of Engineers. 28. Tsay, T. K. and Liu, P. L F. (1982). Numerical Solution of Water Wave Refraction and Diffraction Problems in the Parabolic Approximation. Journal of Geophysical Research, Vol. 87, No. CIO. pp. 7932-7940. 29. Tsay, T. K. and Liu, P. L F. (1983). A Finite Element Model for Wave Refraction and Diffraction. Applied Ocean Research, Vol. 5, No. 1, pp. 30-37. 30. UBC SPARSPAK. (1982). Solving Sparse System of Linear Equations. Computing Centre, University of British Columbia, Vancouver. 31. Vastano, A C. and Reid, R. O. (1967). Tsunami Response for Islands: Verification of a Numerical Procedure. Journal of Marine Research, 25, pp. 129-139. 32. Walker, J. and Headland, J. (1982). Engineering Approach to Nonlinear Wave Shoaling. Proceedings, 18th. International Coastal Engineering Conference, Ch. 34, pp. 523-542. 33. Weggel, J. R. (1972). Maximum Breaker Height for Design. Proceedings, 13th International Coastal Engineering Conference, Ch 21, pp. 419-431. 34. Zienkiewicz, O. C. (1981). The Finite Element Method. 3rd Ed, McGraw-Hill Publishing Company Ltd. 61 APPENDIX - I USER'S MANUAL The finite element program ISLAND should be used in conjunction with the data generating program DATGEN. DATGEN generates the data for the analysis of a circular island with a linearly sloping bottom. The output file of DATGEN serves as the input file for ISLAND. The input data file for DATGEN should be prepared as follows. A-I.1 DESCRIPTION OF INPUT DATA FILE The input data file for DATGEN, read on unit 5, must be specified as described below. The format statements are indicated in parenthesis but a free format can also be used. The corresponding output file (which serves as the input file for ISLAND) should be assigned unit 6. HD, HS, AM, DIA (4F10.4, 1 line) HD = Deep water depth HS = Shore water depth AM = Slope of the berm. For slope = l:m, the input is m only. DIA = Top diameter of the island. NTT, SC (14. F8.4, 1 line) NTT = Plotting code for mesh. 0 = No plot, 1 = Plot 62 SC = Scale of plot. For max dimension to the output width = s:l, the input is s only. ANINC (F10.4, 1 line) ANINC = Angular increment in degrees, i.e. angle subtended by each element at the centre of the island. NR. NA (214, 1 line) NR = Number of radial increment from the vertical side of the island. NA = Number of angular increments in 180°, i.e. 180°/ANINC DR (F10.4, NR lines) DR = Values of the radial increments starting from the shore of the island Jl, J2, J3 (314, 1 line) Jl = Print option for the input data while executing ISLAND. 0 = No print, 1 = Print J2 = Breaking option in analysis. 0 = No wave breaking considered, 1 = wave breaking considered. J3 = Option for preparing input files for Graphics. 0 = No preparation, 1 = Prepare input files. T. H, G (3F12.4, 1 line) T = Wave period H = Incident wave height G = Acceleration due to gravity. 63 A-I.2 EXECUTION OF DATA GENERATING PROGRAM For execution of the data generating program, the following command has to be given after compilation. R -LOAD 5=*input file', 6='output file' If the plotting code NTT = 1, the plot of the mesh will be automatically generated. For a copy of the plot in the printronix terminal, the command is R 'PXPLOT 0=-PLOT For a copy from the QMS plotter, the command is R 'QMSPLOT 0=-PLOT A-L3 EXECUTION OF FINITE ELEMENT PROGRAM After compilation of the finite element program ISLAND, the following command should be given for execution. R -LOAD+NICL:NEWSPARSPAK 5=XDATGEN output file' 6=¥E output file' If the input option for graphics data J3 = 1, three more output files are to be specified in the execution command. These are the files which will contain island geometry data, node numbers and corresponding wave heights and the nodes and elevations to be used in the execution of the graphics program. The files should have units 7, 8 and 9 assigned to them respectively and the complete execution command will then be R -LOAD+NICL:NEWSPARSPAK 5 =*DATGEN output file' 6=TE output file' 7=*Geometry file' 8=*Wave height file' 9=xWave elevation file' 64 A-1.4 TRANSFER OF GRAPHICS FILES T O V A X / V M S Since the coloured graphics program is in Civil Engineering VAX system, the graphics input files from MTS are to be transferred to VAX before executing the graphics program. After signing on a VAX terminal the following commands should be given for this purpose. ALLOC TXAO (i.e. zero) KERMIT SET LINE TXAO CONN After these commands, dial up UBC (228-1401) on modem. Then give the following commands. G SIG "ccid" "password" R "KERMIT SEND "filename" Ctrl ]c RECEIVE "filename" CONN (Repeat from SEND for transferring other files) EXIT SIG EXIT LO 65 A-I.5 EXECUTION OF T H E GRAPHICS PROGRAM To execute the graphics program WAVEPLOT, the following commands are to be given from a VAX terminal after signing on. UNIX WAVEPLOT.SEIKO The coloured graphics program is interactive and most of the prompts are self explanatory. When a prompt asks for a displacement file either a wave height or a wave elevation file name should be specified depending upon the plot required. It should be noted that in answer to prompts, the values assigned to the directions of a normal vector (coordinate of the point to where a person is looking at) and to a up vector (coordinate of a point that specifies the orientation of a plot) should not be such that they become colinear because in such a case the program will be aborted. For example, the graphs presented in this thesis correspond to a normal vector of (0, -3, -2) and a up vector of (0, 0, 3). However, the interactive nature of the program allowes for any number of possible views and the selection of a suitable one is entirely the user's discretion. A-1.6 LIMITATION OF T H E FINITE E L E M E N T PROGRAM The maximum number of nodes used in this program so far is 1467. It has the capability to go up to 4000 nodes but problems may arise due to the unavailablity of sufficient storage space required to solve a matrix of the corresponding size. The program also issues a warning if the incident wave is too high and exceeds the wave breaking criteria in the deep water.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Transformation of waves around artificial islands
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Transformation of waves around artificial islands Talukdar, Kushal K. 1986
pdf
Page Metadata
Item Metadata
Title | Transformation of waves around artificial islands |
Creator |
Talukdar, Kushal K. |
Publisher | University of British Columbia |
Date Issued | 1986 |
Description | A numerical method based upon the finite element technique is described to determine the transformed wave field around artificial islands for regular unidirectional waves. The present work takes into account various aspects of wave transformation, including wave refraction, reflection, diffraction, shoaling and wave breaking, in a computer model which finally predicts the transformed wave heights and elevations at different locations. A modified version of the original two-dimensional mild slope equation, to include the energy dissipation from wave breaking has been used for this purpose along with criteria for the onset and cessation of wave breaking. Results are presented first for the cases where analytical solutions are available in order to check for accuracy and convergence of the present computer model. Since no previous results are available for the proper island geometries, the computer predictions are compared with measurements of the wave field around a model island in a laboratory' basin. In particular, the wave heights at specified points and the runup around the island are compared. The study indicates the need for including wave breaking in such an analysis and also suggests possible ways to calibrate certain features of the computer model for a better numerical prediction. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-07-11 |
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.0062930 |
URI | http://hdl.handle.net/2429/26331 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1986_A7 T34.pdf [ 4.44MB ]
- Metadata
- JSON: 831-1.0062930.json
- JSON-LD: 831-1.0062930-ld.json
- RDF/XML (Pretty): 831-1.0062930-rdf.xml
- RDF/JSON: 831-1.0062930-rdf.json
- Turtle: 831-1.0062930-turtle.txt
- N-Triples: 831-1.0062930-rdf-ntriples.txt
- Original Record: 831-1.0062930-source.json
- Full Text
- 831-1.0062930-fulltext.txt
- Citation
- 831-1.0062930.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062930/manifest