Citation |

- Permanent Link:
- http://ufdc.ufl.edu/UF00075474/00001
## Material Information- Title:
- Water wave interaction with porous structures of irregular cross sections
- Series Title:
- UFLCOEL-TR
- Creator:
- Gu, Zhihao, 1957- (
*Dissertant*) University of Florida -- Coastal and Oceanographic Engineering Dept Dean, Robert G. (*Thesis advisor*) Sheppard, D. Max (*Reviewer*) Kurzweg, Ulrich H. (*Reviewer*) - Place of Publication:
- Gainesville, Fla.
- Publisher:
- University of Florida
- Publication Date:
- 1990
- Copyright Date:
- 1990
- Language:
- English
- Physical Description:
- xvii, 201 leaves : ill. ; 29 cm.
## Subjects- Subjects / Keywords:
- Berms ( jstor )
Best fit ( jstor ) Boundary conditions ( jstor ) Breakwaters ( jstor ) Damping ( jstor ) Energy dissipation ( jstor ) Ocean floor ( jstor ) Reflectance ( jstor ) Velocity ( jstor ) Waves ( jstor ) Coastal and Oceanographic Engineering thesis Ph. D Dissertations, Academic -- Coastal and Oceanographic Engineering -- UF - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Abstract:
- A general unsteady porous flow model is developed based on the assumption that the porous media can be treated as a continuum. The model clearly defines the role of solid and fluid motions and henceforth their interaction. All the important resistant forces are clearly and rigorously defined. The model is applied to the gravity wave field over a porous bed of finite depth. By applying linear wave theory, an analytical solution is obtained, which is applicable to the full range of permeability. The solution yields significantly different results from those of contemporary theory. The solution requires three empirical coefficients, respectively representing linear, nonlinear, and inertial resistance. Laboratory experiments using a standing wave system over a porous seabed were conducted to determine these coefficients and to compare with analytical results. The coefficients related to linear and non.linear resistances were found to be close to those obtained by previous investigators. The virtual mass coefficient was determined to be around 0.46, close to the theoretical value of 0.5 for a sphere. The analytical solution compared well with the experiments. Based on this porous flow model and linear wave theory, two numerical models using boundary integral element method with linear elements are developed for permeable submerged breakwaters and berm breakwaters, respectively. Due to the establishment of a boundary integral expression for wave energy dissipation in a porous domain and the application of the radiation boundary condition on the lateral boundary (ies), the numerical models are highly efficient while maintaining sufficient accuracy. The numerical results show that the wave energy dissipation within a porous domain has a well defined maximum value at certain permeability for a specified wave and geometry condition. The nonlinear effects int he porous flow model are clearly manifested, as all the flow field properties are no longer linearly proportional to the incident wave heights. The numerical results agreed reasonably well with the experimental data on the seaward side. On the leeward of the breakwater, despite the appearance of higher order harmonics, the numerical model produces acceptable results of energy transmission based on energy balance.
- Thesis:
- Thesis (Ph. D.)--University of Florida, 1990.
- Bibliography:
- Includes bibliographical references (leaves 197-199).
- General Note:
- Typescript.
- General Note:
- Vita.
- Funding:
- This publication is being made available as part of the report series written by the faculty, staff, and students of the Coastal and Oceanographic Program of the Department of Civil and Coastal Engineering.
- Statement of Responsibility:
- by Zhihao Gu.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. Â§107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact Digital Services (UFDC@uflib.ufl.edu) with any additional information they can provide.
- Resource Identifier:
- 24160971 ( OCLC )
## UFDC Membership |

Full Text |

UFL/COEL/TR-083
WATER WAVE INTERACTION WITH POROUS STRUCTURES OF IRREGULAR CROSS SECTIONS by Zhihao Gu Dissertation 1990 WATER WAVE INTERACTION WITH IRREGULAR CROSS POROUS STRUCTURES OF SECTIONS By Zhihao Gu A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1990 ACKNOWLEDGEMENTS It is almost impossible to fully express the author's appreciation in words to everyone who has contributed to the completion of this dissertation. First of all, the author wishes to express his deepest appreciation to his advisor and the chairman of his advisory committee, Dr. Hsiang Wang, for his advice, friendship, encouragement and financial support which enabled the author to come to the United States to pursue his higher education. The four years of graduate work with Dr. Wang have been a challenging and enjoyable experience in the author's life. The author would also like to thank Dr. Robert G. Dean, Dr. D. Max Sheppard and Dr. Ulrich H. Kurzweg for serving as the members of his doctoral advisory committee; Dr. Daniel M. Hanes for revising the dissertation and attending the final exam. Thanks are also due to all other teaching faculty members in the department during the author's graduate study: Dr. Michel K. Ochi, Dr. Ashish J. Mehta, Dr. Peter Y. Sheng and Dr. James T. Kirby (now at University of Delaware), for their teaching efforts, time and being role models to the author. The author has benefitted substantially from their invaluable experience, knowledge and continuous inspiration throughout the four years' study. The author is especially grateful for the valuable discussions with Dr. Dean and Dr. Kirby on numerous topics directly or indirectly related to this work. The continuous financial support throughout the Ph.D. program by the Coastal and Oceanography Engineering Department is also greatly appreciated. The author is deeply indebted to his former advisors and supervisors in China, Profs. Wen-Fa Lu, Zhizhuang Xieng and Weicheng Wang. Without their generous support and recommendation, it would have been impossible for the author to study abroad in the first place. The continuous support and understanding by them and all other faculty and staff members in the author's former working unit are highly appreciated. Very special thanks are given to Prof. Alf Torum, Head of Research in Norwegian Hydrotechnical Laboratory, Trondheim, Norway, and Dr. Hans Dette in LeichtweiB Institut ffir Wasserbau, Technische Universitft Braunschweig, West Germany, for providing computer, office and accommodation facilities while the author was visiting the two institutions in late 1988 with his advisor. The initiation and the provision of technical information by Prof. Torum for the topic of berm breakwater which later became a chapter in this dissertation are gratefully acknowledged. Thanks are given to Mr. Sidney Schofield, Mr. Jim Joiner and other staff members in the Coastal Laboratory for their help during the experimental phases of the study. The assistance provided by Mr. Subarna Malakar in the use of computer facilities is also appreciated. Thanks are extended to Helen Twedell for the efficiency and courtesy in running the archives and for the parties, cookies and cakes; and to Becky Hudson and all other secretarial personnel in the department for the great job they have done related to this work. The support of many friends and fellow students is warmly appreciated. Special thanks go to Mrs. Jean Wang for her moral support, constant encouragement and hospitality. Deep appreciation is given to fellow students Rajesh Srinivas and Paul Work for proofreading the manuscripts and to Byunggi Hwang for assisting with rock sorting during the experiments. The friendship shared with Yixin Yan, Steve Peene, Yuming Liu and Richard McMillen and other friends, made the author's stay in Gainesville a delightful period of time. The author's gratitude towards his family is beyond words. The author would like to thank his best friend and wife, Liqiu (also a Ph.D. student at the time), for her support, both moral and physical, patience and useful technical discussions. She was always the one to count on for assistance on weekend laboratory work, graphics work and more often on the housework beyond her share. Without her support, the road leading to this goal would have been much more tortuous. Finally, this work is dedicated to the author's wonderful mother and father. Their affection and early family education are reflected in between the lines of this work and the author's daily conduct. With all these efforts, thank god it's Phinally Done! TABLE OF CONTENTS ACKNOWLEDGEMENTS ............................ LIST OF FIGURES ................................ viii LIST OF TABLES ................................ xv ABSTRACT .................................... xvi CHAPTERS 1 INTRODUCTION ............................... 1 1.1 Problem Statement ............................ 1 1.2 Objectives and Scope ........................... 3 2 LITERATURE REVIEW ........................... 5 2.1 Porous Flow Models ........................... 5 2.2 Wave-Porous Seabed Interactions .................... 10 2.3 Modeling of Permeable Structures of Irregular Cross Sections .... 13 3 POROUSFLOW MODEL .......................... 20 3.1 The Equation of Motion ......................... 20 3.2 Force Coefficients and Simplifying Assumptions ......... 23 3.3 Relative Importance of The Resistant Forces ............. 28 4 GRAVITY WAVES OVER FINITE POROUS SEA BOTTOMS ..... 31 4.1 Boundary Value Problem ........................ 31 4.2 The Solutions of The Complex Dispersion Equation ......... 35 4.3 Results ........................ .. ......... 41 5 LABORATORY EXPERIMENT FOR POROUS SEABEDS ....... 5o 5.1 Experiment Layout and Test Conditions ................ 50 5.2 Determination of The Empirical Coefficients . 52 5.3 Relative Importance of The Resistances in The Experiment 55 5.4 Comparison of The Experimental Results and The Theoretical Values 55 6 BOUNDARY INTEGRAL ELEMENT METHOD . 64 6.1 Basic Formulation . . . 64 6.2 Local Coordinate System for A Linear Element . 68 6.3 Linear Element and Related Integrations . . 70 6.4 Boundary Conditions . . . 74 7 NUMERICAL MODEL FOR SUBMERGED POROUS BREAKWATERS 76 7.1 Governing Equations . . . 76 7.2 Boundary Conditions . . . 78 7.2.1 Boundary Conditions for The Fluid Domain . 78 7.2.2 Boundary Conditions for The Porous Medium Domain 81 7.3 BIEM Formulations . . . 82 7.3.1 Fluid Domain . . . 82 7.3.2 Porous Medium Domain . . 86 7.3.3 Matching of The Two Domains . . 87 7.4 Linearization of the Nonlinear Porous Flow Model . 88 7.5 Transmission and Reflection Coefficients . . 91 7.6 Total Wave Forces on an Impervious Structure . 92 7.7 General Description of The Computer Program . 94 7.8 Numerical Results . . . 95 8 LABORATORY EXPERIMENTS OF A POROUS SUBMERGED BREAKWATER . . . . 107 8.1 General Description of The Experiment . . 107 8.2 Wave Transmission and Reflection . . 109 8.3 Pressure Distribution and Wave Envelope Over The Breakwater 121 9 NUMERICAL MODEL FOR BERM BREAKWATERS............ 128 9.1 Mathematical Formulations............................ 128 9.1.1 CFSBC for The Free Surface Inside Porous Medium....... 130 9.1.2 BIEEM Formulations............................ 131 9.2 Linearization..................................... 133 9.3 Numerical Results .. .. .. ... ... ... ... .... ... ....134 10 SUMMARY AND CONCLUSIONS. .. .. ... ... ... .... .....146 10.1lSummnary. .. .. .. ... ... .... .. .... ... ... ... ..146 10.2 Conclusions. .. .. .. ... ... ... ... ... ... .... .....148 10.3 Recommendations for Future Studies...................... 151 APPENDICES A BOUNDARY INTEGRAL FORMULATION FOR ENERGY DISSIPATION IN POROUS MEDIA................................... 153 B EXPERIMENTAL DATA FOR POROUS SEABEDS.............. 159 BIBLIOGRAPHY........................................ 197 BIOGRAPHICAL SKETCH................................. 200 LIST OF FIGURES 3.1 Definition sketch for the porous flow model . 21 3.2 Regions with different dominant resistant forces . 30 4.1 Definition Sketch .......................... 32 4.2 Progressive wave case: h, = DS = 5.0 m and h = DW = 2.0 6.0 m. (a) Nondimensional wave number kr/(aCr2/g), (b) Nondimensional wave damping rate ki/(oa2/g) . 45 4.3 Maximum nondimensional damping rate (kg),.a/(a2/g) and its corresponding permeability parameter R as functions of nondimensional water depth h. (a2/g) . . 46 4.4 Standing wave case. (a) Nondimensional wave frequency or/(L/g) ; (b) Nondimensional wave damping rate ag/(L/g)-. . 48 4.5 Solutions based on four porous flow models. (a) Nondimensional wave frequency u,/(L/g)1, (b) Nondimensional wave damping rate o, /(L/g)1 ...... .. ................. .... 49 5.1 Experimental setup ......................... 51 5.2 Typical wave data: (a) Averaged nondimensional surface elevation (rt/H), (b) Nondimensional wave heights (H/H1) and the best fit to the exponential decay function . 54 5.3 The Measurements and the predictions vs. R. for L = 200.0 cm, ho = 20.0 cm, h = 25.0 cm. (a) Wave frequency ra,, (b) Wave damping rate ai ........................... 58 5.4 Theoretical values by the present model vs. experimental data of Table 5.5, Table 5.6 and Table 5.7. (a) Wave frequency o,, (b) Wave damping rate ai....................... 61 5.5 Theoretical values by the model of Liu and Dalrymple vs. experimental data of Table 5.5, Table 5.6 and Table 5.7. (a) Wave frequency ar, (b) Wave damping rate ai. . 62 6.1 Auxiliary coordinate system . . 68 7.1 Computational domains. .. .. .. .. ... ... ... .... ..77 7.2 Flow chart of the numerical model for porous submerged breakwaters .. .. .. .. ... ... ... ... ... .... ... ....96 7.3 Wave envelopes for (a) 'Transparent' submerged breakwater; (b) Impermeable step .. .. .. .. .. ... .... ... ... ... ..98 7.4 Porous submerged breakwater: (a) Wave form and wave envelope; (b) Envelopes of pressure and normal velocity .. .. .. ....99 7.5 Wave field over submerged breakwaters: (a) Impermeable; (b) Permeable .. .. .. .. ... ... ... ... ... ... ... ..101 7.6 Transmission and reflection coefficients vs. stone size for different wave periods. (a) Transmission coefficient; (b) Reflection coefficient. .. .. .. .. .. ... ... ... ... ... ... ....103 7.7 Transmission and reflection coefficients vs. R for different wave periods. (a) Transmission coefficient; (b) Reflection coefficient. 104 7.8 Transmission and reflection coefficients vs. R for different wave heights. (a) Transmission coefficient; (b) Reflection coefficient. .105 7.9 Wave forces and over turning moment for a impermeable submerged breakwater: (a) Wave forces; (b) Overturning moment. .106 8.1 Experiment layout .. .. .. .. ... ... ... ... ... ....108 8.2 Typical wave record. (a) Partial standing waves on the up wave side, (b) Transmitted waves on the down wave side. .. .. .. ..111 8.3 The wave spectrum of the transmitted waves .. .. .. .. ....112 8.4 The predicted Kt and K, versus the measured Kt and K,. (a) Kpvs. Kt,; (b) Kt vs. Ki.... .. .. .. .. .. .. .. ... ...117 8.5 The predicted and measured Kt and K, versus Hj/gT. (a) Kt and Ktm..; (b) Kp, and Kr.. .. .. .. .. .. .. .... ... ...118 8.6 The predicted and measured Kt and K, versus H,/gT. (a) Ks,, and Kt m; (b) K7, and K .. .. .. .. .. .. .. .... ... ...119 8.7 The predicted and measured Kt and K7 versus H,/gT. (a) K1,, and Ktm; (b) Kp, and Krm,.. .. .. .. .. ... ... .... ...120 8.8 Transmitted and reflected wave heights versus the incident wave heights. (a) Transmitted waves; (b) Reflected waves .. .. .. ...122 8.9 Transmitted and reflected wave heights versus the incident wave heights. (a) Transmitted waves; (b) Reflected waves .. .. .. ...123 8.10 The envelopes of wave and pressure distribution for T = 0.858 sec.; non-breaking wave case. (a) Wave envelope; (b) Envelope of pressure distribution. .. .. ... ... ... ... ... ....126 8.11 The envelopes of wave and pressure distribution for T = 0.858 sec.; breaking wave case. (a) Wave envelope; (b) Envelope of pressure distribution. .. .. .. .. .... ... ... ... ....127 9.1 Definition sketch for berm breakwaters. .. .. .. ... ... ..129 9.2 Flow chart of the numerical model for berm breakwaters .135 9.3 Berm breakwaters of vertical face: (a) Zero permeability; (b) Infinite permeability. .. .. .. .. .. ... ... ... .... ...136 9.4 Berm breakwaters of inclined face: (a) Zero permeability; (b) Infinite permeability. .. .. .. .. .. ... ... ... .... ...137 9.5 The Cross Section of The Berm Breakwater. .. .. ... ....138 9.6 Permeable berm breakwater of model scale with H = 5.0 cm: (a) Wave envelope; (b) Envelopes of pressure and normal velocity distribution .. .. .. .. ... ... ... .... ... ... ....140 9.7 Permeable berm breakwater of model scale with H = 10.0 cm: (a) Wave envelope; (b) Envelopes of pressure and normal velocity distribution .. .. .. .. ... ... ... ... .... ... ....141 9.8 Permeable berm breakwater of model scale with H = 20.0 cm: (a) Wave envelope; (b) Envelopes of pressure and normal velocity distribution .. .. .. .. ... ... ... .... ... ... ....142 9.9 Permeable berm breakwater of prototype scale with H = 2.0 m: (a) Wave envelope; (b) Envelopes of pressure and normal velocity distributions. .. .. .. .. ... ... .... ... ....143 9.10 Permeable berm breakwater of prototype scale with H = 4.0 m: (a) Wave envelope; (b) Envelopes of pressure and normal velocity distributions. .. .. .. .. ... ... .... ... ....144 9.11 Permeable berm breakwater of prototype scale with H = 8.0 m: (a) Wave envelope; (b) Envelopes of pressure and normal velocity distributions. .. .. .. .. ... ... .... ... ....145 A.1 Geometric relations between the vectors .. .. .. ... .....156 B.1 Caseof L =200 cm, h =DW =30cm, h, =DS =2Ocm and d4o = DD = 0.72 cm. (a) Averaged nondimensional surface elevation (?71H,), (b) Nondimensional wave heights (IT/RI) and the best fit to the exponential decay function .. .. .. .. .. ...161 B.2 Case of L = 200 cm, h = DW = 30 cm, h, = DS = 20 cm and dso = DD = 0.93 cm. (a) Averaged nondimensional surface elevation (q/H1), (b) Nondimensional wave heights (H/-I) and the best fit to the exponential decay function . 162 B.3 Case of L = 200 cm, h = DW = 30 cm, h, = DS = 20 cm and ds50 = DD = 1.20 cm. (a) Averaged nondimensional surface elevation (,l/H1), (b) Nondimensional wave heights (H/--) and the best fit to the exponential decay function . 163 B.4 Case of L = 200 cm, h = DW = 30 cm, h, = DS = 20 cm and dso = DD = 1.48 cm. (a) Averaged nondimensional surface elevation (q7/H1), (b) Nondimensional wave heights (H/7-) and the best fit to the exponential decay function . 164 B.5 Case of L = 200 cm, h = DW = 30 cm, h, = DS = 20 cm and dso = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (q7/HI), (b) Nondimensional wave heights (H/--) and the best fit to the exponential decay function . 165 B.6 Case of L = 200 cm, h = DW = 30 cm, h. = DS = 20 cm and ds0 = DD = 2.84 cm. (a) Averaged nondimensional surface elevation (t7/Hi), (b) Nondimensional wave heights (H/Rl) and the best fit to the exponential decay function . 166 B.7 Case of L = 200 cm, h = DW = 30 cm, h, = DS = 20 cm and dso = DD = 3.74 cm. (a) Averaged nondimensional surface elevation (ti/Hi), (b) Nondimensional wave heights (H/H) and the best fit to the exponential decay function . 167 B.8 Case of L = 200 cm, h = DW = 25 cm, h, = DS = 20 cm and dso = DD = 0.72 cm. (a) Averaged nondimensional surface elevation (qt/Hi), (b) Nondimensional wave heights (H/-) and the best fit to the exponential decay function . 168 B.9 Case of L = 200 cm, h= DW = 25 cm, h, = DS = 20 cm and dso = DD = 0.93 cm. (a) Averaged nondimensional surface elevation (t/Hi), (b) Nondimensional wave heights (H/7--1) and the best fit to the exponential decay function . 169 B.10 Case of L = 200 cm, h = DW = 25 cm, h, = DS = 20 cm and dso = DD = 1.20 cm. (a) Averaged nondimensional surface elevation (77/Hi), (b) Nondimensional wave heights (H/HIV) and the best fit to the exponential decay function. . 170 B.11 Case of L = 200 cm, h = DW = 25 cm, h, = DS = 20 cm and dso = DD = 1.48 cm. (a) Averaged nondimensional surface elevation (77/Hi), (b) Nondimensional wave heights (H/--) and the best fit to the exponential decay function . 171 B.12 Case of L = 200 cm, h = DW = 25 cm, h, = DS = 20 cm and dso = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (t/Hi), (b) Nondimensional wave heights (H/H-) and the best fit to the exponential decay function . 172 B.13 Case of L = 200 cm, h = DW = 25 cm, h, = DS = 20 cm and dso = DD = 2.84 cm. (a) Averaged nondimensional surface elevation (7/Hi), (b) Nondimensional wave heights (H/R-1) and the best fit to the exponential decay function . 173 B.14 Case of L = 200 cm, h = DW = 25 cm, h, = DS = 20 cm and dso = DD = 3.74 cm. (a) Averaged nondimensional surface elevation (?/H1), (b) Nondimensional wave heights (H/h) and the best fit to the exponential decay function . 174 B.15 Case of L = 200 cm, h = DW = 20 cm, h, = DS = 20 cm and dso = DD = 0.72 cm. (a) Averaged nondimensional surface elevation (7/H1), (b) Nondimensional wave heights (H/I7) and the best fit to the exponential decay function . 175 B.16 Case of L = 200 cm, h = DW = 20 cm, h, = DS = 20 cm and dso = DD = 0.93 cm. (a) Averaged nondimensional surface elevation (?/H1), (b) Nondimensional wave heights (IH/Y) and the best fit to the exponential decay function . 176 B.17 Case of L = 200 cm, h = DW = 20 cm, h, = DS = 20 cm and ds50 = DD = 1.20 cm. (a) Averaged nondimensional surface elevation (n/H1), (b) Nondimensional wave heights (H/I) and the best fit to the exponential decay function . 177 B.18 Case of L = 200 cm, h = DW = 20 cm, h, = DS = 20 cm and dso = DD = 1.48 cm. (a) Averaged nondimensional surface elevation (n/HI), (b) Nondimensional wave heights (IH/HI) and the best fit to the exponential decay function . 178 B.19 Case of L = 200 cm, h = DW = 20 cm, h, = DS = 20 cm and ds50 = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (q/Hi), (b) Nondimensional wave heights (H/-) and the best fit to the exponential decay function . 179 B.20 Case of L = 200 cm, h = DW = 20 cm, h, = DS = 20 cm and dso = DD = 2.84 cm. (a) Averaged nondimensional surface elevation (q/H), (b) Nondimensional wave heights (F/R) and the best fit to the exponential decay function . 180 B.21 Case of L = 200 cm, h = DW = 20 cm, h, = DS = 20 cm and ds50 = DD = 3.74 cm. (a) Averaged nondimensional surface elevation (n/H1), (b) Nondimensional wave heights (H/-I) and the best fit to the exponential decay function . 181 B.22 Case of L = 200 cm, h = DW = 30 cm, h, = DS = 15 cm and dso = DD = 1.48 cm. (a) Averaged nondimensional surface elevation (77/H1), (b) Nondimensional wave heights (H/H-1) and the best fit to the exponential decay function . 182 B.23 Case of L = 200 cm, h = DW = 25 cm, h, = DS = 15 cm and ds0 = DD = 1.48 cm. (a) Averaged nondimensional surface elevation (7/Hi), (b) Nondimensional wave heights (H/--- and the best fit to the exponential decay function . 183 B.24 Case of L = 200 cm, h = DW = 20 cm, h, = DS = 15 cm and dso = DD = 1.48 cm. (a) Averaged nondimensional surface elevation (i7/Hj), (b) Nondimensional wave heights (H/H--) and the best fit to the exponential decay function . 184 B.25 Case of L = 200 cm, h = DW = 30 cm, h, = DS = 10 cm and dso = DD = 1.48 cm. (a) Averaged nondimensional surface elevation (i7/Hi), (b) Nondimensional wave heights (H/-) and the best fit to the exponential decay function . 185 B.26 Case of L = 200 cm, h = DW = 25 cm, h, = DS = 10 cm and dso = DD = 1.48 cm. (a) Averaged nondimensional surface elevation (rl/Hi), (b) Nondimensional wave heights (H/7--) and the best fit to the exponential decay function . 186 B.27 Case of L = 200 cm, h = DW = 20 cm, h, = DS = 10 cm and dso = DD = 1.48 cm. (a) Averaged nondimensional surface elevation (q/Hi), (b) Nondimensional wave heights (H/H) and the best fit to the exponential decay function . 187 B.28 Case of L = 225 cm, h = DW = 30 cm, h, = DS = 20 cm and dso = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (7l/Hi), (b) Nondimensional wave heights (7/Wj) and the best fit to the exponential decay function . 188 B.29 Case of L =225 cm, h = DW = 25 cm, h, = DS = 20 cm and dso = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (7/Hi), (b) Nondimensional wave heights (7/-- and the best fit to the exponential decay function . 189 B.30 Case of L = 225 cm, h = DW = 20 cm, h, = DS = 20 cm and dso = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (77/H1), (b) Nondimensional wave heights (H/-1) and the best fit to the exponential decay function . 190 B.31 Case of L = 250 cm, h = DW = 30 cm, h, = DS = 20 cm and dso = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (7/H1), (b) Nondimensional wave heights (H/-i and the best fit to the exponential decay function . 191 B.32 Case of L = 250 cm, h = DW = 25 cm, h, = DS = 20 cm and dso = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (ql/H1), (b) Nondimensional wave heights (7/--) and the best fit to the exponential decay function . 192 B.33 Case of L = 250 cm, h = DW = 20 cm, h, = DS = 20 cm and ds50 = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (iq/H1), (b) Nondimensional wave heights (H/-)- and the best fit to the exponential decay function . 193 B.34 Case of L = 275 cm, h = DW = 30 cm, h, = DS = 20 cm and ds50 = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (q/H,), (b) Nondimensional wave heights (H/H-) and the best fit to the exponential decay function . 194 B.35 Case of L = 275 cm, h = DW = 25 cm, h, = DS = 20 cm and ds50 = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (/HI1), (b) Nondimensional wave heights (7H/) and the best fit to the exponential decay function . 195 B.36 Case of L = 275 cm, h = DW = 20 cm, h, = DS = 20 cm and ds50 = DD = 2.09 cm. (a) Averaged nondimensional surface elevation (q/H1), (b) Nondimensional wave heights (H/-1) and the best fit to the exponential decay function . 196 LIST OF TABLES 3.1 Illustration of Dominant Force Components Under Coastal Wave Conditions: [a O(1)rad/sec, V(p/y) 0(10-1)] ........... 30 4.1 Comparison of The Predictions and The Measurements ..... ..43 5.1 Material Information ...... ........................ 51 5.2 Test Cases ....... .............................. 52 5.3 Measured a, and ai for h, = 20 cm and L = 200 cm....... ...53 5.4 Comparison of The Resistances. h, = 20 cm, L = 200 cm ... 56 5.5 Comparison of Measurements and Predictions for h, = 20 cm, L = 200 cm...................................... 57 5.6 Comparison of a, and ap for do = 1.48 cm, L = 200 cm..... 59 5.7 Comparison of a,, and ap for d50 = 2.09 cm, h, 20 cm..... 59 5.8 Comparison of a, and ap for do = 0.16 am, h, = 20 cm, L = 200 cm ..................................... 60 5.9 Comparison of a, and ap for h = dso ................... 63 8.1 Test Results of Non-breaking Waves .................... 114 8.2 Comparison of Km and Kp ........ ..................... ..115 8.3 Test Results for Breaking Waves ....................... 121 8.4 Normalized Pressure Distribution ..................... 125 B.1 Test Cases ....... .............................. 159 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy WATER WAVE INTERACTION WITH POROUS STRUCTURES OF IRREGULAR CROSS SECTIONS By Zhihao Gu December. 1990 Chairman: Dr. Hsiang Wang Major Department: Coastal and Oceanographic Engineering A general unsteady porous flow model is developed based on the assumption that the porous media can be treated as a continuum. The model clearly defines the role of solid. and fluid motions and henceforth their interactions. All the important resistant forces are clearly and rigorously defined. The model is applied to the gravity wave field over a porous bed of finite depth. By applying linear wave theory, an analytical solution is obtained, which is applicable to the full range of permeability. The solution yields significantly different results from those of contemporary theory. The solution requires three empirical coefficients, respectively representing linear, nonlinear and inertial resistance. Laboratory experiments using a standing wave system over a porous seabed were conducted to determine these coefficients and to compare with analytical results. The coefficients related to linear and nonlinear resistances were found to be close to those obtained by previous investigators. The virtual mass coefficient was determined to be around 0.46, close to the theoretical value of 0.5 for a sphere. The analytical solution compared well with the experiments. Based on this porous flow model and linear wave theory, two numerical models using boundary integral element method with linear elements are developed for permeable submerged breakwaters and berm breakwaters, respectively. Due to the establishment of a boundary integral expression for wave energy dissipation in a porous domain and the application of the radiation boundary condition on the lateral boundary(ies), the numerical models are highly efficient while maintaining sufficient accuracy. The numerical results show that the wave energy dissipation within a porous domain has a well defined maximum value at certain permeability for a specified wave and geometry condition. The nonlinear effects in the porous flow model are clearly manifested, as all the flow field properties are no longer linearly proportional to the incident wave heights. The numerical results agreed reasonably well with the experimental data on the seaward side. On the leeward of the breakwater, despite the appearance of higher order harmonics, the numerical model produces acceptable results of energy transmission based on energy balance. Xvii CHAPTER 1 INTRODUCTION 1.1 Problem Statement Among the existing shore protecting breakwaters, a large number of them are rubble-mound structures made of quarry stones and/or artificial blocks. They can be treated as structures of granular materials. A breakwater is called subaerial when its crest is protruding out of the water surface and submerged when its crest is below the water level. Sometimes a breakwater is subaerial at low tide and submerged at high tide. These types of structures are usually termed as low crest breakwaters. When a beach is to be protected from wave erosion, a detached shore parallel submerged breakwater may provide an effective and economic solution. The advantages of submerged breakwaters as compared to subaerial ones are low cost, aesthetics (they do not block the view of the sea) and effectiveness in triggering the early breaking of incident waves, thus reducing the wave energy in the protected area. With the increasing interest in recreational beach protection, where complete wave blockage is not necessary, submerged breakwaters may find more and more applications. On the other hand, if wave blockage is the main objective such as for harbor and port protection, a subaerial breakwater may be more effective. The traditional design of such a structure is a trapezoidally shaped rubble mound with an inner core covered with one or more thin layers of large blocks to form the armor layer(s) to protect the core. Owing to the demand for deeper water applications, the armor sizes have become larger and larger. This would greatly increase not only the cost but also the structural vulnerability. A relatively new type of structures, called 2 berm breakwater, is attracting more and more attention. For a berm breakwater, the armor layer(s) is replaced by thicker layer(s) of blocks of much smaller sizes. The seaward face of the structure has a berm section instead of the traditional uniform slope. The berm section is intended to enhance the structural stability and to trigger wave breaking. Wave attenuation behind a porous breakwater is affected by three mechanisms: wave reflection on the seaward face of the structure, waves breaking over the structure and flow percolation inside the porous structure. The former one is conservative and the latter two axe dissipative. The main focus of this study is on the third mechanism, that is the dissipation due to flow percolation. For such a purpose, a numerical solution is sought since an analytical solution for such irregularly shaped structures is nearly impossible. The main advantage of numerical methods is the flexibility of handling complex geometries and boundary conditions. The disadvantage of any numerical approach is the lack of generality for the solution, since it is usually implicit in terms of the variables so that the influences of the parameters can only be examined case by case. In order to examine the validity of the model and the nature of the dissipative force, the porous flow model is applied to an infinitely long flat seabed of finite thickness subject to wave action. This condition can also be viewed as a submerged breakwater with infinitely long crest. Because of the geometrical simplicity, an analytical solution is attainable. This case is used to compare with the existing solutions proposed by other investigators, to examine the nature of wave attenuation as a function of flow and material properties and, more importantly, to guide the design of an experiment so that the empirical coefficients in the porous model can be determined. Successful determination of these coefficients is crucial to the validity of the model. The experiment which is to be carried out subsequently must demonstrate that the results are stable (or the experiments are repeatable) and that 3 they cover, at least, an adequate range of intended application. A numerical solution with computer code is then to be developed for arbitrary geometry. Boundary Integral Element Method (BIEM) is proven to be very efficient for boundary value problems with complicated domain geometries. With this method, the solution is expressed in boundary integrals and no interior points have to be involved in the solution procedure. The aim of this study is to develop an efficient computer code based on such a method. Finally, the validity of the solution is to be assessed by a set of experiments. 1.2 Objectives and Scope Specifically, the objectives of this study are listed as follows: 1. Develop a percolation model suitable for unsteady and turbulent porous flows, 2. Verify the model through an analytical solution and laboratory experiments for the case of a flat porous seabed subject to linear gravity waves, 3. Based upon the porous flow model, develop a numerical solution for submerged and berm breakwaters of arbitrary cross sections using the boundary integral element method. To achieve these goals, the research is carried out in the following steps: 1. Derive the porous flow model, 2. Examine and interpret the relative importance of the various terms in the model to narrow the scope of the study, 3. Obtain the analytical solution of wave attenuation over a flat porous seabed, compare the solution with the existing ones and examine the wave energy dissipation process, 4. Conduct a laboratory experiment for the flat porous seabed case to determine the empirical coefficients in the porous flow model and to verify the analytical solution, 5. Develop a numerical model using the BIEM method for porous submerged breakwaters, 4 6. Verify the numerical model with laboratory experiments of a submerged rubblemound breakwater, 7. Modify the numerical model to represent porous berm breakwaters. CHAPTER 2 LITERATURE REVIEW This chapter is a brief literature review of the past efforts made in studying the flow in porous media with and without water waves. The review is separated into three groups: one on porous flow model, one on analytical solutions for wave and porous-seabed interactions and the other one on the modeling of interactions between waves and porous structures of irregular cross sections. Since the porous media encountered in coastal engineering are largely of the granular type, made of sand, gravel, quarry stones or artificial blocks, the deformation of the solid skeleton of the media is usually negligibly small as compared to that of the pore fluid. The literature review deals with granular media only; other types of porous media, such as poroelastic or poroplastic media and so on, are not included. 2.1 Porous Flow Models A large number of porous flow models have been proposed over the past few decades as an effort to quantitatively model the flows in porous media. Some of the models were based purely on experiments and some of the others on certain theoretical considerations. The derivations of the these models are based primarily on the following three approaches: 1) Simple element approach which views a porous medium as the assembly of simple elementary elements such as a bundle of tubes and so on, 2) Microscopic approach which recognizes the midcrostructure of porous media in derivations. The final formulations of the models usually have to be obtained by taking the spatial average of the microscopic equations, 3) Phenomenological approach which homogenizes the porous media as a continuum and the flow within 6 the medium is assumed to be continuous. The third approach is the most popular for porous flows in granular materials since the microscopic structure is not well understood and the corresponding constitutive equation not well established. The first porous flow model-Darcy's law-was proposed by the French engineer Darcy based on his experiments more than a century ago. In his model, the characteristic properties of porous media are lumped into one parameter-permeability coefficient-and the model has the form - Vp= (2.1) where Vp is the gradient of the pore pressure, it and q- are, respectively, the dynamic viscosity and the specific discharge velocity of the pore fluid and K. is called the intrinsic permeability coefficient, which reflects the collective characteristics of porous media, such as the porosity, the roughness of the pore walls, the tortuosity, the connectivity etc. This model is very simple in form and reasonably accurate for steady porous flows within media of low permeability, generally of the order of Kp = 10-9 10-12 m2, where the flows are normally laminar dominant. However, for porous media with relatively higher permeability, the porous flow is no longer laminar dominant, turbulence plays a larger and larger role with increasing pore size. In such cases, Darcy's model tends to over estimate the discharge velocity for a given pressure gradient. Dupuit and Forchheimer (cited in Madsen, 1974) modified Darcy's law by adding a term quadratic in velocity to take into account the turbulent effects: - Vp = (a + b I 4I)q" (2.2) where a = A/KP, and b is coefficient for turbulent resistance defined by Ward (1964) as b = pC! (2.3) 7 where p is the density of the fluid and C1 is a nondimensional coefficient. The two coefficients a and b were further expressed in terms of particle diameter by Engelund (1953) and Bear et al. (1968) in the forms (1 n)v (2.4) a = ao n2d2 1-n r b = bo 1n (2.5) where ao and b0 are nondimensional coefficients and n is the volumetric porosity of porous media. The above two models, Eq.(2.1) and (2.2) are both for steady porous flows, although some applications have been made for certain unsteady flows, such as wave induced porous flows in seabeds (Putnam, 1949; Liu, 1973; etc.). These applications were primarily on low permeability media like fine sand. In general, for unsteady flows, the inertial force arises due to the acceleration of the pore fluid and such a force can be quite large and can even be dominant when the permeability is high. To model the unsteady porous flow, a number of models have been suggested. The first model found in the literature containing the inertial force term was the one by Reid and Kajiura (1957). It is the direct extension of Darcy's model: - Vp = P8t (2.6) KP nat The inertial term in this model is in fact the local acceleration term in Navier-Stokes equation with the fluid velocity expressed as a specific discharge velocity q- It is clear that this term contains only the inertia of the fluid but not the inertia induced by the fluid-solid interaction, or the added mass effect. The extension of Dupuit-Forchheimer model leads to Polubarinova-Kochina's model (Scheidegger, 1960; McCorquodale, 1972). It has the form - Vp = a+b2 + (2.7) 8 where a, b and c are empirical coefficients. Murray (1965), when studying the viscous damping of gravity waves over a permeable seabed, proposed - Vp = !6 +pn- (2.8) where k is the permeability of porous media (defined differently from Kp) and 6 is the spatially averaged microscopic velocity of the pore fluid, related to 4* by 4 = n6Z. Comparing to Eq.(2.7), the coefficient c is equal to p in Murray's model. McCorquodale in 1972 further modified the expression of this coefficient as P (2.9) n and his model (McCorquodale, 1972) becomes - Vp= (a +bT) pa j (2.10) n at At about the same time as McCorquodale, Sollitt and Cross (1972) also extended Eq.(2.2) to include the effects of unsteadiness. In their model, the inertia resistance consists of two components, one is the inertia of the pore fluid, and the other one is the inertia induced by the virtual mass effect due to the fluid-solid interactions; the model reads + z~pC+ 1-n Ca a Vp = ~ (n + n'P) (2.11) where C is again the spatially averaged microscopic velocity vector and Ca is the virtual mass coefficient. The value of the virtual mass coefficient C. in Sollitt and Cross's model was assumed to be zero in their actual computation because it was unknown at the time. It has been so assumed in almost all the later applications of this model (Sulisz, 1985, etc.). When comparing to the experimental data, Sollitt and Cross found that the correlation improved by taking nonzero values for Ca and asserted that "one 9 cannot predict the magnitude of this coefficient a priori because the virtual mass of densely packed fractured stone is not known" (Sollitt and Cross, 1972, p1842). They also pointed out that "Evaluation of Ca, however, may serve as a calibrating link between theory and experiment in future studies" (same page). Hannoura and McCorquodale (1978) made an attempt to determine the value of C0 by laboratory experiment. They measured the instantaneous velocity and the pressure gradient in the experiment, and calibrated the data with their semitheoretical porous flow model: - Vp = (a + b ap(l+ Co)Lq (2.12) The results so obtained for C. scattered in a range of -7.5 +5.0. Dagan in 1979, based on the microscopic approach, arrived at a generalized Darcy's law for nonuniform but steady porous flows (Dagan, 1979): -Vp = V2q (2.13) where -y is a constant coefficient depending only on the geometry of the media and it was defined as dZL (2.14) In applying this model to the wave-porous-seabed interaction, Liu and Dalrymple (1984) added an acceleration term to the Dagan's model and it becomes - Vp- = j -ppq+nP-q V2- (2.15) K,, n at K, It was found, by performing a dimensional analysis (Liu, 1984), that the last term in Dagan's model is a second order quantity in comparison to the other terms. The inertial term in this model is the same as that of Reid and Kajiura and that of Sollitt and Cross where C. is set to be zero. Based on the phenomenological approach, Barends (1986) added another porous flow model to the list: - V = C aq gq (2.16) = at K(q) 10 and K(q) was defined by Hannoura and Barends (1981, see Barends, 1986) as /pd, 2.7 K(q) = g (2.17) with f Ca (1 n) Pins where C is the drag coefficient related to the porous Reynolds number (R =1 q do/v), a and /3 are two constants and d, is the relevant grain size. With so many models, it is difficult to decide which one is most appropriate for the porous flow problem in this study without further analysis. 2.2 Wave-Porous Seabed Interactions The computation of gravity wave attenuation over a rigid porous bottom has been performed by a number of investigators. The classic approaches by Putnam (1949) and Reid and Kajiura (1957) was to assume that the Laplace equation is satisfied in the overlying fluid medium and that the bottom layer can be treated as a continuum following Darcy's law of permeability. In Reid and Kajiura's paper, although the inertial term was included in the porous flow model, they later neglected it and concluded that Darcy's law is adequate for sandy seabeds. Under the assumption of low permeability, Reid and Kajiura found, for infinitely thick seabeds, that a2 gk,. tanh krh (2.18) ki = 2(agp/v)k,. (2.19) 2kh + sinh 2krh for progressive waves, where a is the wave frequency, k, and k are, respectively, the wave number and the damping rate, h is the water depth and v is the kinematic viscosity of the fluid. Hunt (1959) examined the damping of gravity waves propagating over a permeable surface using Reid and Kajiura's porous flow model and retained the inertial term to the end. The thickness of the permeable beds was infinite and water above 11 the permeable surface was assumed viscous as opposed to the inviscid approach by the previous investigators. The tangential velocities at the surface of the permeable bed were set to be zero and the continuity of the normal velocity and the pressure were observed on the same boundary. The resulting dispersion equation was quite lengthy. Under the consideration of sand beds of low permeability, the results for kr and k, are k = ko + 2k (2.20) r 2o 2koh + sinh 2koh 2k0 aK k- = 2ko (- + ko ) (2.21) S 2koh+sinh2koh (2.21)2 with a' = gko tanh koh (2.22) Murray (1965) investigated the same problem by the use of stream function and the porous flow model given in Eq.(2.7). Instead of letting the tangential velocities at the interface be zero as done by Hunt, Murray assumed that the tangential velocities were finite and continuous. The first order solution of the dispersion equation for the spatial damping was k + 2koN S2koh + sinh 2koh (2.23) where k = ko + k, with ko being the initial wave number for the impervious bottom as defined in Eq.(2.22) and N= naKp For fine sand, N -+ 0 and k, = v~ a(2.24) 2kh + sinh 2krh and the temporal damping for known k in such case is a = i/uk (2.25) vAsinh 2kh 12 Liu (1973) added an interface laminar boundary layer in his solution for the same problem. In his paper, the water in the fluid domain was regarded as inviscid and the porous flow was assumed to follow Darcy's law. The imaginary part of k, after the boundary layer correction, was given by 2k,+ L) , 2kh + sinh 2kh (aKF2/v + ,7 a 2 (2.26) with k, remaining the same as that in Eq.(2.18). It is noted that this is the same expression as that given by Hunt (1959). The correction term of Eq.(2.26) to Eq. (2.19) is found to be insignificant unless Kp is extremely small. Dean and Dalrymple (1984) later obtained the expressions for shallow water waves over sand beds: k,- =a[1 1- 4 h/ 211/2 (2.27) k, (1 hh/g) K = ( --) (2.28) All the solutions above are for infinitely thick seabeds. Liu (1977) solved the problem for a stratified porous bed where each layer has a finite thickness and different permeability. It was found that the wave damping rate was insignificant to the permeability stratification while the wave induced pressure and its gradients are affected significantly by stratification. The porous flow model employed in the analysis also obeys Darcy's law only. Liu and Dalrymple (1984) further modified the solution, for a bed of finite thickness, by replacing Darcy's law with Dagan's unsteady porous flow model, which partially included the effects of unsteady flows. The dispersion relationship from the homogeneous problem-without laminar boundary layer corrections-was found as R gk gk (R + i)(1 tanhkh) +Rtanhkh, tanhkh(1 -k-ctnhkh) = 0 (2.29) n a2 72 or equivalently Cr2 gk tanh kh=- tanh kh, (gk r2 tanh kh) (2.30) Rn where h, is the seabed thickness and R is the permeability parameter defined as R = crK' It has been shown that the corrections by the laminar boundary layers were not significant since "the damping is largely due to the energy losses in the porous medium rather than the boundary layer losses" (Liu and Dalrymple, 1984, p47). Therefore, it is reasonable to believe that if an appropriate porous flow model is adopted, the solution for a homogeneous problem, without the correction of laminar boundary layers, will be sufficient for engineering purposes. Besides the theoretical studies, the interaction of waves and porous seabeds was also investigated experimentally by Savage (1953). The experiment was carried out with progressive waves over sand beds of median diameters, d, = 0.382 cm. and d, = 0.194 cm, respectively. The permeability coefficient of the sand bed was measured as 44.9 x 10-10m', and the dimension of the seabed was 0.3 meters thick and 18.3 meters long. The wave heights at the beginning and the end of the sand bed were recorded and adjusted to eliminate the effects of the side friction. Most of the theories listed above were claimed to agree well with the data in this experiment. 2.3 Modeling of Permeable Structures of Irregular Cross Sections The ability to predict the performance of rubble-mound breakwaters under the attacks of ocean waves is critically important in designing such structures. Significant amounts of effort have been devoted to this subject. In the early stage of development, research was basically confined to laboratory experiments and only the empirical formulae extracted from the experiments were available for designs. In recent years, with the continuing efforts by researchers, more and more theoretical models for different types of breakwaters are becoming available. For subaerial 14 breakwaters, in addition to a large number of laboratory experiments, a multitude of theoretical models have been developed. The analysis of such structures has been approached both analytically and numerically. Analytically, Sollitt and Cross (1972) investigated the problem of wave transaission and reflection by a vertical face (rectangular cross section) porous breakwater. The porous flow model used in the formulation was Sollitt and Cross's model although the virtual mass coefficient C. was assumed to be zero in the computation. In their solution, the whole computation domain was divided into three regions, two fluid regions and one porous region. The velocity potential functions in the three regions were assumed to be the summations of the fundamental mode and the evanescent modes according to the wave maker theory. The nonlinear term in the flow model was linearized according to Lorentz's hypothesis and the normal velocities and the pressure were matched at both vertical faces of the breakwater. The transmission and the reflection coefficients for long waves, simplified from the complete solution of complex matrix form, were given by Kb = i (2.31) 1+ --n -(S if + n2) K, = (2.32) S-if +n2-i2nab where S is the coefficient for the inertial resistance of porous media, i is the imaginary unit, f is the linearized coefficient for velocity related resistances, b is the width (or thickness) of the breakwater and all the other symbols have the same meanings as before. Madsen (1974) later re-examined the same problem beginning directly with the linear long wave theory. He adopted the Dupuit- Forchheimer model and introduced the inertial resistance into the momentum equation. The resulting expression of Kt and Kr are almost identical to Eq.(2.31) and (2.32). By carrying out the volume 15 integration over the breakwater for the energy dissipation, Madsen obtained an explicit expression for the linearized friction factor f, that is n KpIa. Kpla. 16b 1 f = n[-(1 pa) + (, + -pa-)2 + 16b ai (2.33) K,,1 2a 2a 37rh where a and b are the coefficients in Dupuit-Forchheimer's porous flow model. A similar solution for this problem was also given by Scarlatos and Singh (1987). Madsen et al. (1976, 1978) further extended the long wave solution to a trapezoidal porous breakwater with, again, Dupuit-Forchheimen's model. The solution was, however, much more complicated than that for the crib type. Since the shore protection breakwaters are usually irregular in shape and very often with several layers of stones of different sizes, analytical solutions become impractical. For such complicated geometries, numerical approaches have to be adopted. The most commonly used numerical schemes are finite element and Boundary Integral Element Methods (BIEM). The first finite element model was developed by McCorquodale (McCorquodale, 1972) using the McCorquodale porous flow model, for computing the wave energy dissipation in rockfalls. In his model, the entire cross section of the structure was divided into small triangle-elements and the variation of the physical quantities was interpolated by a time dependent element function: = (1 + 2X + 3y)t + 34 + SX + 6y (2.34) Since the numerical computation was only carried out in the porous domain, the interaction between the porous domain and the fluid domain was not modeled, although the free surface within the porous region was well predicted. Due to the large amount of work for data preparation in using a finite element model, the Boundary Integral Element Method (BIEM) became popular since the mid-1970s. With BIEM, the discretization is only on the boundaries as opposed to the entire domain with the finite element method. Ijima et al. (1976) applied 16 this method to porous breakwaters with constant elements, which assumes that the physical quantities remain constant over each individual element. The whole computation domain was divided into three regions, two fluid regions and one porous region. In the porous domain, Darcy's law was used and in the two fluid regions, two artificially defined vertical boundaries were placed at the offshore and inshore ends. The patterns of the vertical distribution of the velocity potential at these two boundaries was assumed undistorted by the presence of a structure. The reflected and the transmitted potential functions outside the computational domain were given by ,.(X, z) = [ekz'+Ae'z-] cosh k(z +h) (2.35) -i())Iv r cosh kh Ot(x,z) = BOe-iA(,+,,)coshk'(z + h) (2.36) cosh k'h' where I and 1 are the distances from the left and the right vertical boundaries to the origin, k and k' are the wave numbers in the reflection and the transmission regions, respectively. The reflection and transmission coefficients are then K, =1 Ao 1 (2.37) Kt =l0Bo 1 (2.38) The agreement of the transmission coefficient for impermeable floating structures with the experiment data was fairly good, but the comparison of Kt and K, for a sloped-face permeable breakwater was not as satisfactory. Finnigan and Yamamoto (1979) further added the modulated standing wave modes to and Ot, while keeping the constant element and Darcy's law unchanged in their model. The two potential functions were, respectively, _,k.z_,.cosh k(z + h) +cosk,(z + h) ,(XZ) = [e,+c oecosh kh + Amek+(z) cos k,+h) m(1 cos2kh (2.39) B~-'''+ l~hI C cosh kz + h') + 0 M -a''cskC~o 'h t(x,z) = B0eik(z+1)I + B e- cos k'(z + h) (2.40) cosh k'h' m1cos kmlh' with gk tanh kh = -gkm tan kmh = 2 (2.41) = -gk' tan knh = a2 (2.42) where Ai and Bi are unknown complex coefficients. In their model, the series were truncated at a certain point to make the number of unknown coefficients equal to the number of the elements. Owing to the introduction of the series, the procedure became much more complicated and the computation much more time consuming than that of Ijima et al. The advantage of a constant element is its simplicity and efficiency. But it has an inherent problem-awkward behavior around "corner points". It has been found that significant errors arise around such points even for the simplest case-sinusoidal wave propagating in a constant water depth over a impermeable bottom. Sulisz (1985) employed linear elements in his BIEM model for subaerial porous breakwaters. Eqs.(2.39) and (2.40) were adopted for the lateral boundary conditions. As a result, the formulation became awfully complicated. In modelling the porous media, Sollitt and Cross's model was adopted. However, the inertia term due to virtual mass effect was dropped as the corresponding coefficient C. was not known. For submerged breakwaters, the modelling is still largely based on laboratory experiments and empirical formulae. For impermeable submerged breakwaters, quite a few empirical formulae for transmission coefficient have been established from experimental data by different authors (see Ba.ba, 1986). For example, the equation given by Goda (1969) for the transmission coefficient due to overtopping is K= = 0.5[1 sin F (2.43) 2a H, 18 with a = 2.2 and /3 = 0.0 0.8 for vertical face breakwaters, where F is the depth of submergence of the breakwater crest and Hi is the incident wave height. Another equation for impermeable submerged breakwaters was given by Seelig (1980): = C(1 F (1- 2C)F (2.44) where R is the wave run up given by Franzius (1965, cited in Baba, 1986) R = HIC(O.123-HI)(C2vTH7U+c3) (2.45) with C, = 1.997, C2 = 0.498 and C3 = -0.185. d is the total water depth and C = 0.51- 0,11B (2.46) h where B is the crest width and h is the height of the breakwater. Several other empirical equations and methods are also available for designs. The latest theoretical approach for impermeable submerged breakwaters was given by Kobayashi et al. (1989) by using finite- amplitude shallow-water equations. The numerical results were found to agree well with the experimental data by Seelig (1980) for such structures. However, for permeable structures, no parallel empirical formula could be found in the literature besides the nomograph by Averin and Sidorchuk (1967, cited in Baba, 1986). The research for such structures is still in the stage of physical experiments. Dick and Brebner (1968) carried out an experiment on solid and permeable breakwaters of vertical faces. In the experiment, it was discovered that, over a certain wave length range, the permeable breakwater was much better than the solid one of the same dimension in terms of wave damping, and that the permeable breakwater had a well defined minimum value for the coefficient of transmission. It was also found that a substantial portion, 30% 60%, of the wave energy of the transmitted waves was transferred to higher frequency components. The equivalent 19 wave height for the transmitted waves in his paper was defined as Heq = 2 vf2 jT2 (t t(2.47) Due to the smallness of the submergence of the breakwater crest, it is felt that the majority of the waves in the experiment were breaking waves. Dattatri et al. (1978) tested a number of submerged breakwaters of different types, permeable, impermeable, rectangular and trapezoidal. One of the conclusions was that the important parameters affecting the performance of a submerged breakwater are the crest width and the depth of submergence. The transmitted waves were found "to be irregular though periodic". They also reported that porosity and wave steepness did not have significant influence on the transmission coefficient, which contradictory with the observations by Dick et al. (1968). Seelig (1980) conducted a large number of tests on the cross sections of 17 different breakwaters for both regular and irregular waves. Most of the breakwaters tested were rubble-mound porous structures with multi-layer designs. Beyond the experiments for impermeable breakwaters by which Eq. (2.44) was obtained, he also tested these permeable structures as submerged breakwaters. Since "no generalized model was available" at the time for such breakwaters, the experimental data with regular waves were compared to Eq.(2.44). It was found that the formula is quite conservative in estimating the transmission coefficient for permeable submerged breakwaters. The phenomenon of wave energy shifting to higher order harmonics in the transmitted waves was again reported without further explanation. Baba (1986) conducted an experiment on concrete submerged breakwaters and compared the data with the four widely used computational methods for wave transmission coefficient for impermeable submerged breakwaters. He concluded that the formula given by Goda (1969) was the most suitable one in the case of a shore protecting submerged breakwater. CHAPTER 3 POROUS FLOW MODEL Examining the existing porous flow models listed in the literature review, it is noted that most of them contain one, two or all of the following components: the linear velocity term, the quadratic velocity term and the term proportional to the acceleration of the fluid. For the first component all the models give the same definition while for the rest two terms, different models have different forms. To be able to model the porous media with confidence, it is necessary to formulate the porous flow from fundamental principles. 3.1 The Equation of Motion The equation of motion for a solid body placed in an unsteady flow field can be expressed in the following general terms: FP +FD+FI+Fb-FF =ma (3.1) where m = mass of the solid; a = acceleration of the body; Fp = pressure force; FD = velocity related force also known as the drag force; F = acceleration-related force also known as inertial force; Fb = body force; and Ff = frictional force due to surrounding solids. We now apply this equation to a control volume of a mixture of fluid and solids, as shown in Fig. 3.1, and restrict our discussion to the twodimensional case. The x-axis coincides with the horizontal direction and the z-axis points vertically upward. The spatially averaged equation of motion for the solid skeleton in the x-direction becomes a- dx(1 nAZ)dz + FD. + P-' + Fb.z FfZ = p.(1 n) ii,, dxdz (3.2) where Lp = pressure gradient in the x-direction; nA = averaged area porosity aw p w + -dz p + -d 8z Oz U P Bu u-dx p + "axdx Op p+- dz Figure 3.1: Definition sketch for the porous flow model defined as the ratio of Aid to Atotal with Avo being the area of the voids and Atomt being the total control area; n = volumetric porosity defined as the ratio of Void to Vtota~; Tb, is the body force of the solid; p, = density of solid; u, and ti, = velocity and acceleration of the solid. The subscript z denotes the x- direction and the overbar denotes spatial average. Since, from this point on, we will be dealing exclusively with spatially averaged values, the overbar will be dropped. The force balance on the pore fluid can be established in a similar manner: aP - dzx na. dz F, FI, + Ffy, = pn dz dz (y, i,) ax (l.-f. with p being the density of the fluid and Fb&1 being the body force of the fluid in z direction; u and tif equal, respectively, the actual spatially averaged velocity and acceleration of the pore fluid, defined as u = ua dv (3.4) od = V old (3.3) 22 where ua = actual pore velocity. If the fluid and solid mixture is now being treated as a porous medium, and thus a continuum, we require that 1. All the field variables, such as that defined by Eq.(3.4), be independent of the volume of integration. 2. 6 < L, where 5 and L are respectively, the length scales of the pore and the system. For the wave attenuation problem, Eq.(3.3) is of special interest. We introduce here yet another field variable q representing the apparent spatially averaged fluid velocity, which is related to Uf by q = 1 uadv "f uadv = n.u (3.5) here V is the total volume Vtota. This velocity also known as the discharge velocity, is related to the actual discharge over a unit surface area, or, Q/A. These apparent properties are of final interest in engineering. Since we are now dealing with a continuum, we have Dqq + V =t = n" -(3.6) ijt at q*V at here it has been assumed that the convective acceleration is negligible. The various force terms are established as follows: (A). Drag Force The drag force consists of two terms: that due to laminar skin friction and that due to turbulent form drag; the former is proportional to the velocity and the latter is proportional to the velocity squared. Their functional form, when expressed in terms of the final field variables, may be written as FD2 = [Ax(qz nu..) + B. q nu, I (q- nu.)Idx nA dz (3.7) The coefficients A, and B. are to be determined later. (B). Inertia Force The inertial force can be treated as an added mass effect and expressed in the following form: Fj = C(6 nt,,)dx nAz dz (3.8) The coefficient C, will be determined later. (C). Body Force In the pore fluid, there is no horizontal body force, and the vertical body force is balanced by the static pressure gradient and, thus, vanishes. In the solid skeleton, the vertical body force is the net weight. If the solid skeleton is subjected to an unsteady vertical motion, this term should be included; otherwise, it can be ignored. Substitution of Eqs.(3.5) through (3.8) into Eq.(3.3) results in the following differential equation: aP - -- Ax(q< nu,2) + Bz I q nu, I (q. nu,) ax + ( +Cz)(. -nik..) (3.9) This is the basic equation of pore fluid motion. In a general case, it is coupled with Eq.(3.2) and must be solved simultaneously. Only a special case with no movement of the solid skeleton will be analyzed here. For such a case, Eq.(3.9) reduces to ap = Aqz + Bx, I q I q. + (' + C. )4. (3.10) axn The force balance in the z-direction leads to aP az -. A~q + B. I qI q. + (' + C,)5k (3.11) 3.2 Force Coefficients and Simplifying Assumptions The evaluation of force coefficients involves a great deal of empiricism. However, successful engineering application relies heavily upon successful estimation of these coefficients. The part of the flow resistance which is linear in q clearly will lead to 24 a Darcy-type resistance law. The coefficient A is, at least, related to four foundational properties, one of the fluid-the dynamic viscosity-and three of the porous structure- -the porosity, the tortuosity (one that defines effective flow length) and the connectivity (one that defines the manner and number of pore connections). Obviously, the more those factors can be specified explicitly, the more accurately the value of A can be determined. However, it is also generally true that the more factors explicitly introduced the more restrictive the range of application. If none of these four fundamental properties is expressed explicitly, A is simply the inverse of Darcy's hydraulic permeability coefficient. If the dynamic viscosity, tL, is factored out, we have the empirical law: A=KA (3.12) Kp where Kp is known as the intrinsic permeability. A number of investigators including Engelund (1953, cited in Madsen, 1974) and Bear et al. (1968, cited in Madsen, 1974) attempted to relate A to porosity as well and came up with the relationship of the following type: (1i-n)3 it A = a0 n 2 (3.13) where d, is a characteristic particle size of the pore material. The coefficient a0 obviously still contains the other properties such as tortuosity and connectivity. Since both of these factors have directional preference, a0 should also be a directional property and, in general, is a second-order tensor. For isotropic material ao becomes a scalar. Engelund (1953, see Madsen, 1974) recommended ao = 780 to 1500 or more (3.14) with the values increasing with increasing irregularity of the solid particle. Attempts were also made to sort out the effects of tortuosity and connectivity (Fatt, 1956a,b). The conditions invariably become more restrictive and one is required to specify soil characteristics beyond the normal engineering properties. 25 The characteristics of the coefficient B can be determined by examining the total form drag resistance acting on a unit volume of granular material of characteristic size d,. For this case, we have FDN. = 6, P-,CDAP I U1 I ,U, (3.15) 2 where AP is the projected area of individual particles; CD is the form drag coefficient of the individual particles; b is a correction coefficient for CD accounting for the influence of surrounding particles; and m is the total number of particles per unit volume. Since the total projected area should be proportional to (1 n)dx dz/d., we have, after substituting q for uf: FDN. = CD nI q dxdz (3.16) 2 nd. the constant of proportionality is absorbed in 6. Comparing Eq.(3.16) with the second part of Eq.(3.7), we obtain 2 DnAn2d, If we let: b0 = -5- (3.18) Eq.(3.17) can be expressed as B.n (3.19) Pboflf2d, This equation is very similar to that recommended by Engelund (1953), with the exception that nA was replaced by n, the volumetric porosity, in his formula. Again, because b0 has directional preference, it should also be a tensor of second order. Engelund recommended b0 = 1.8 to 3.6 or more (3.20) for granular material of sand-sized particles. 26 The coefficient C. is the overall added mass coefficient for the porous medium. Its characteristics can be determined in much the same manner as those of B,, or m Fl. = (3.21) where C. is the added mass coefficient of each individual particle and V is the volume of each particle. One may then readily obtain, by comparing Eq.(3.21) with Eq.(3.8), C. = pC.(1-n) (3.22) nAz n Since C. is related to the volume of the solid skeleton, one does not expect significant directional preference, unless the geometry of the element deviates significantly from a sphere shape. A new coefficient, Cm, can be defined such that C. = p + n = p[n + Co(1 n) (3.23) nAz nA, which has the same meaning as the mass coefficient in hydrodynamics. Clearly the force coefficients in other directions can be defined in a similar fashion. Equations.(3.10) and (3.11) can now be expressed in the following general form: A4+B q1 + Cm =-Vp (3.24) Here, the carat indicates a second order tensor. This equation can be generalized to include the motion of the solid simply by replacing q with q Uo. We now proceed to assume that the porous material is isotropic and that all the coefficients reduce to scalar quantities and nAz = Ax =- n Recall that the convective acceleration in q has been assumed, in Eq.(3.6), to be negligible, i.e. q= T 27 Substituting these conditions into Eq.(3.24), we obtain (A + B I q 1)4 + C,, -- -Vp (3.25) at This equation is similar to the well-known equation of motion in a permeable medium such as given by Reid and Kajiura (1957) and others, with the exception that the added mass effect is now formally introduced. To seek a solution to this equation, the common approach is to ignore the inertial term and linearize or ignore the nonlinear term. Under such conditions, Eq.(3.25) can be reduced to the Laplace equation by virtue of mass conservation. To retain the inertia term, the most convenient approach is to assume the motion to be oscillatory, as is the case under wave excitation. Now we let q= qtx, z)e cit (3.26) where Fis a spatial variable only, and a is a generalized wave frequency, which could be a complex number. Substituting Eq.(3.26) into Eq.(3.25) leads to - Vp = (A- iaCm, + B I q 1)" (3.27) A new set of nondimensional parameters are defined as follows: Permeability parameter, R, R- pgp (3.28) with K = nd; (3.29) A a o (1 -n) Inertia parameter, /) =n + C.(1 -n) (3.30) n2 Volumetric averaged drag coefficient, Cd Cd = B = b01 n (3.31) p n~d, 28 Substituting these parameters into the above equation reduces it to the following familiar form: - Vp = po( 1 ig + Iq I) (3.32) R GC This is essentially the same equation used by Sollit and Cross (1972) and others for porous breakwater analysis. 3.3 Relative Importance of The Resistant Forces Before solving Eq.(3.32) with an overlying wave motion, it is useful to examine the importance of the resistance forces and to determine the nature of the fluid motion in the porous medium. By taking ratios of the three respective resistance forces, three nondimensional parameters can be established: Inertial Resistance f1 n + C(1 (3.33) Laminar Resistance f = ao(1 n)3 Turbulent Resistance f, b0 Laminar Resistance fJ = aon(1 n)21 (3.34) Turbulent Resistance f,, b0(1 n) Rf (3.35) Inertial Resistance f = Can(1 n + n R, C. where R, and R are two forms of Reynolds number defined as R_ q Ido (3.36) and II d! (3.37) where v is the kinematic viscosity of the fluid. Obviously both Reynolds numbers signify. the relative importance of the inertial force to the viscous force. The origins of the inertial forces are, however, different. In Rf, the inertia is of a convective nature and the resistance arises due to change of velocity in space (fore and aft the body) whereas in Ri, the inertia is of a local 29 nature and the resistance arises due to the rate of change of velocity at a specific location. The ratio of Ri and R! is the Strouhal number, which clearly identifies the different origins of the two inertial forces. In the range of common engineering applications, the magnitudes of various coefficients can be estimated as follows: n s 0.3-0.6 ao ~0(1) bo 0(1) Therefore, Eqs.(3.33) through (3.35) reduce to t I ~ 10-2R1 (3.38) I-_ I ~ 10-2Rf (3.39) f I _R (3.40) When Eqs.(3.38) to (3.40) are plotted on a R1 and Ri plane as given in Fig. 3.2, we identify seven regions where the three resistance forces are of varying degrees of importance. There are three regions where only one resistance force dominates, that is, the dominant force is at least one order of magnitude larger than the other two. There is one region where the three forces are of equal importance. Then, there are three intermediate regions where two out of three forces could be important. Since both Reynolds numbers are flow-related parameters, accurate position of a situation within the graph cannot be determined a priori. However, a general guideline can be provided with the aid of the graph or with Eqs.(3.38) through (3.40). We give here an example of practical interest. In coastal waters, we commonly encounter wind waves with frequency in the order of 0(1)(rad/sec) and the nondimensional bottom pressure gradient defined as V(p/y) in the order of O(10-1). 106 IN Kegion 10 > lo f 10o- > lOf, 102 A, f, f-, A 1 L Region 10-1- f I Region ft > 10f, 0A > lof, > 10f 10-6R, 10-6 10-4 10-2 1 10 10 106 Figure 3.2: Regions with different dominant resistant forces Under such conditions, the importance of the three resistance components for bottom materials of various sizes can be assessed. Table 3.1 illustrates the results. Table 3.1: Illustration of Dominant Force Components Under Coastal Wave Conditions: [o O(1)rad/sec, V(p/q) 0(10-1) Description Size Disch. Vel. Rf R Dominant S Range (m/sec) Force Coarse sand <2 mm. < O(10-3) < O(1) < 0(1) laminar or finer Pebble, laminar or small 1 cm. O(10-2) O(102) O(102) turbulence gravel inertia Large gravel 10 cm. O(10-1) O(104) O(104) turbulence crusted stone inertia Boulder 0.3 1.0 m O(100) O(106) O(0) turbulence crusted stone inertia Artificial turbulence blocks, > 1.0 m > O(100) > O(106) > O(106) inertia large rocks III CHAPTER 4 GRAVITY WAVES OVER FINITE POROUS SEA BOTTOMS We consider here the case of a small amplitude wave in a fluid of mean depth h above a porous medium of finite thickness h.. The bottom beneath the porous medium is impervious and rigid. The subscript s will be used here to denote variables in the porous bed. The basic approach is to establish governing equations for different zones separately and to then obtain compatible solutions by applying proper matching boundary conditions. In the pure fluid zone, it is common to assume the motion essentially irrotational except near the interface. Most of the investigators neglected the influence of the boundary layer with the exception of Liu (1973) and Liu and Dalrymple (1984) who included in their solutions two laminar boundary layers at the mud- line. Since 'the damping is largely due to the energy losses in the porous medium rather than the boundary layer losses (Liu and Dalrymple, 1984) ', the boundary layer effect will be ignored in this study to simplify the mathematics. 4.1 Boundary Value Problem The governing equation for the velocity potential function in the fluid domain is v2D= 0 -h< z - VP. = pufoq h+ h.) <<-h (4.2) where fo is the linearized resistance coefficient. {z T 7 AV x 1k1 I I V {\ [x, h Figure 4.1: Definition Sketch If the porous medium is rigid, thus, incompressible, the continuity equation for the discharge velocity, V i = 0, leads to VP=0 ( A + h < h 4.3) The boundary conditions of entire system are r(x, t) = = a ei(k-o0 at z = 0 (4.4) g 8t . g. 0-=0 at z = 0 (4.5) 0*o :; ;*o'o'ono 04 o oooo0o0oo0 0O 00 Figure 4.1: Definition Sketch If the porous medium is rigid, thus, incompressible, the continuity equation for the discharge velocity, V. = 0, leads to VIP, =0 (h + h.) <-h (4.3) The boundary conditions of entire system are 1 aD ei(kz._at) ,7(x, t) = gat= a eat z=0 (4.4) a,-.1) z=O at z=0 (4.5) alD p- = P, at z=-h (4.6) at 8_ 1 8P, az f z at z=-h (4.7) 8z pufo 8z 33 BP, -7 = 0 at z = -(h + h,) (4.8) 8z The potential and pore pressure functions are assumed to have the forms #(z,z,t) = [Acoshk(h + z) + Bsinhk(h + z)ei(kz-at) (4.9) P,(x,z,t) = Dcoshk(h+ h, + z)es(kz-at) (4.10) where A, B and D are unknown complex constants. Equation (4.8) is already satisfied by P, in Eq.(4.10). Introducing Eqs.(4.9) and (4.10) into Eqs.(4.6) and (4.7), A and B can be expressed in terms of D, A = i- Dcoshkh, (4.11) Pa D B sinh kh, (4.12) pefo 4 is then given by D(Xz, t) = D [icosh coshcoshk(h + z) + 1 sinhkh sinh k(h+ z)]eCi(kz-at) (4.13) Pa fo D can be obtained by the free surface boundary condition given by Eq.(4.4): Spga (4.14) cosh kh cosh kh.(1 tanh kh tanh kh.) fo Finally, Eq.(4.5) along with Eq.(4.13) gives the complex wave dispersion equation: a2 gk tanh kh = tanh kh,(gk oa2 tanh kh) (4.15) fo here either a or k, or both, could be complex, as well as the coefficient fo. In the above equation, when setting Ca = Cf = 0, it becomes the dispersion equation obtained by Liu and Dalrymple (1984). If C1 # 0, as will be the case in this study, the coefficient fo can not be a known value a priori. It becomes another unknown besides a or k, and therefore the procedure of solution will be different from that employed by Liu and Dalrymple. Before solving Eq.(4.15), a few limiting cases are examined here. 34 For the case of infinite seabed thickness and low permeability (laminar skin friction dominates, fo -- ), Eq.(4.15) reads a2 gk tanh kh = -iR(gk a2 tanlh kh) (4.16) which is the same as obtained by Reid and Kajiura (1957). If R 0 as with an impervious bed, Eq. (4.16) reduces to the ordinary dispersion relationship for a finite water depth h. On the other hand, if the permeability approaches infinity, we have R --+oo, n-- 1.0, #--- 1.0, fo--.-i The dispersion relationship expressed in Eq.(4.15) can be shown as 02 = gk tanh k(h + h,) (4.17) Physically, this is the case of water waves over a finite depth h + ha, or, the solid resistance in layer h, vanishes. Another limiting case is where the water depth approaches zero and waves now propagate completely inside the porous medium. The dispersion relationship from Eq.(4.15) becomes a i tanh kh, (4.18) It can be easily shown that the same dispersion relationship can be obtained by directly solving the problem of linear gravity waves in a porous medium alone. Again, for the case of laminar resistance only, Eq.(4.18) becomes a2 = -iRgk tanh kh (4.19) which states the Darcy-type resistance law. On the other hand, when R -- co, we obtain from Eq.(4.18) a2 = gk tanh kh (4.20) or, as expected, the dispersion relationship in a pure fluid medium. 35 When the virtual mass and the drag coefficients are set to be zero, Eq.(4.15) gives the dispersion equation for the homogeneous solution by Liu and Dalrymple (1984) which is u2 gk tanh kh= 1 tanh kh,(gk a2 tanh kh) (4.21) R 4.2 The Solutions of The Complex Dispersion Equation As mentioned in the previous section, Eq.(4.15) is an equation involving both of the complex unknown variables of either a or k, and the complex coefficient fo. The other necessary equation can be obtained from the linearization process (This step is not necessary if C1 = 0, as in Liu and Dalrymple's solution). The common method of evaluating the linearized resistance coefficient, fo, is to apply the principle of equivalent work. This principle states that the energy dissipation within a volume of porous medium during a time period should be the same when evaluated from the true system or from its equivalent linearized system (ED)I = (ED)a (4.22) where ED is the energy dissipation in a controlled volume during one wave period and the subscripts I and n1 refer to linearized and nonlinear systems, respectively. It can be readily shown (Appendix A) that such energy dissipation (considered as a positive value) can be expressed in the form of a boundary integral ED = d d (4.23) where ED is a complex energy dissipation function and is the complex energy flux normal to S, with the real parts of them being the corresponding physical quantities. Here S is the closed boundary of the computation domain. The complex energy flux function 7,, can be expressed as (Appendix A) 1 t e 2i t . w(t + UnP*) (4.24) 2 36 where u, is the normal velocity at the boundary S and p" is the conjugate of pore pressure p, both of them are complex quantities; a, is the wave frequency, areal value, and the subscript r is used to distinguish the complex a. Physically, there are two classes of problems: standing waves of a specified wave number, and progressive waves of a specified wave frequency. In the former case, k is real, a is complex, and the solution has the following form r7(xt) = ae'i' cos kxe- t o, < 0 (4.25) where ar is the imaginary part of the complex a and a, is the wave frequency. In the latter case, a is real, k is complex and 7 becomes ti(x,t) = ae-kze(ktz-at) ki > 0 (4.26) where k, and ki are, respectively, the real and imaginary parts of k with k, being the wave number. For standing waves, the pore pressure function is P, = PC-'" = Dcoshk(h + h, + z) coskxeite-a't (4.27) Since P, is periodic in x, the boundary curve S for the contour integral in Eq.(4.23) can be chosen as x = 0, z = -(h + h,), x = L and z = -h with L being the wave length. As x = 0 and x = L are the antinodes of the standing wave, there is no normal velocity, i.e. u,, = 0, on these two vertical boundaries and also at the impermeable boundary z = -(h + h.). Then f rt+r lL( p-it ED jt 0 + wo p) dx dt (4.28) with kD Wo= U, I,=h= k- sinh kh, cos kx eait (4.29) paf and f = fo for linearized system (4.30) 1 --i )3+- C IwoI R a - fA + f2 wO j for nonlinear system (4.31) with f = -- i = Cd (4.32) R a In Eq.(4.31), the magnitude of w0 is approximately taken as wo | (fo) sinhkh, (4.33) af with (fo) = (4.34) cosh kh cosh kh,(1 tanh kh tanh kh,) where the decaying wave amplitude ae By introducing Eq. (4.29) into Eq.(4.28) and carrying out the contour integration regarding f as a constant, the energy dissipation within that one-wave-long portion of a seabed is obtained as kLTsih oD'_(ff ) | [D(f) I ED = kLT sinh2kh,( T(t) + T2(t)] (4.35) 8p h f f where Ti (t) and T (t) are the nondimensional functions of t generated by the integration with respect to time. Substituting the linearized and the nonlinear resistance coefficients given in Eqs.(4.30) and (4.31) for f in Eq.(4.35) respectively, equating the energy dissipations 38 by the two systems according to Eq.(4.22) and assuming that a and k are the same for both systems, we have D(fT) T (t) + 72I D (t) = 10 fo D'(fA + f2 Iwol)1) 12 D(T = m) fi f2 Iwo I + D(f + f2 Iwo 1) 2(t) (4.36) 2fl +AIwo I fl+f2 IWoI Since this equation has to be satisfied at any time instant, we obtain the following equations: D 2(fo) D 2(f + f2 IO 1) (437) - (4.37) fo A + f2 IwoI I D(fo) j I D(f + f2 Iwo 1)12 (4.38) O fl + f2 IWOI Therefore, fo = fh + f2 IwO 1 (4.39) Substituting the expressions for fl, f2 and wo in the above equation, it becomes fo = 1 i + Cd I kD(fo) sinh kh, I (4.40) R C Ia Ufo It is clear, from the definitions of the parameters, that all the terms in Eq.(4.40) are complex quantities in the standing wave case. For progressive waves, the pore pressure function P, is given by Eq.(4.10) and the contour S for the integration is the same as that for standing waves. The energy dissipation is t+T E = ds dt t+ = f (- 7- +1 F dz + f n d)dt (4.41) t =0 ==L =-h By substituting the expression of ', derived from Eq.(4.10) into the above equation, the summation of the first two integrals is found to be O(I ki 1) while the last 39 one is 0(1). If we assume that I ki j< 1, which is generally true in coastal waters (see the results for the progressive wave case), then the energy dissipation becomes ED = t 7,% dx dt (4.42) Eq.(4.42) will clearly lead to the same relationship as that given by Eq.(4.40) except that the averaged wave amplitude d in D, in this case, is the spatially averaged value of ae- kiz over [0, L]. In the progressive wave case, we assume that the wave period is given and that there is no time dependent damping, i.e. a is real. Therefore, the first and the last term on the R.H.S. of Eq.(4.40) are all real numbers. It is then obvious that (fo)s--/3 (4.43) and 1 Cd I kD(k, fo) I (4.44) R alfol where the subscripts i and r are for imaginary and real parts, respectively. To actually compute fo requires specification of two fundamental quantities of n and d, (see Eq.(4.40) and the definitions for R, Cd and 6). It is often more convenient to specify the permeability parameter R = aK,/v (R' = a, Kp/v for standing waves), as opposed to specifying the actual granular size, d,. Under these circumstances, we could replace Cd with the following expression (Ward, 1964): C, Cf f (4.45) Cd =- a with - bo(1 n) (4.46) 6 O =nao(1 n)3 Now C1 is a nondimensional turbulent coefficient. To be consistent with the suggested values of ao and b0 given in Eqs.(3.14) and (3.20), C1 should be in the range of 0.3 1.1 for a porosity of n = 0.4. Equation (4.40) now becomes: 1 i/+ C1 I kD(k, fo) sinh kh (47) f = a I afo1 40 This is the final form of the equation used for computing fo for both cases. The dispersion equation, Eq.(4.15), which is coupled with Eq.(4.47), is solved iteratively. The procedures to obtain the solution for the cases of standing waves and progressive waves are different. a) For standing waves a = oa, + icai is complex, k is real and known. Equation (4.15) can be rewritten as gk(tanhkh tanhkh,) a2' = o = Q(fo) + iQ,(fo) (4.48) 1 tanh kh, tanh kh fo with Qr(fo) and Qi(fo) being the real and imaginary parts of a2, respectively. The iteration procedures are summarized as follows: ,(n+1) 1 [ VQ(f'n')+Q;)(fO(-) + Qr(f(n))] (4.49) - [Q(f))+Q2(fn)) Q,(fo) ] (4.50) 1 C k I r( n), (4.51)f V") i i(n)f() (4.51) f 1) 1 if (4.52) 0 R where the superscript n denotes the level of iteration, and the criterion of convergence is f a) fo(-1 and o") '(n-) j< E (4.53) I f(n) a_, n (n) with e being a pre-specified arbitrarily small number. It is set to be 1.0% in this model. b) For progressive waves k = k, + iki is complex and a is real and known. In this case, the dispersion equation can be written as i F(k, fo) = a2 gk tanh kh + -tanh kh,(gk a tanh kh) = 0 (4.54) and the iteration was carried out as follows: and the iteration was carried out as follows: (4.55) r(k"+), f")) = 0 41 F,(k(n+1),fon)) = 0 (4.56) f(n) _1 j,3+ 1 C1 I kD(k(n),fo"))l (4.57) R VraV' lafoCn)I fo i) (4.58) R where F, and F are the real and imaginary parts of F and n indicates the iteration level. The criterion of convergence for such case is fn) f O I-) and I k(-) ) ( (4.59) fofl) k(n) with e being a pre-specified arbitrarily small number. It is again set to be 1.0% in this model. 4.3 Results The predicted damping rates kip from the solution of Eqs.(4.15) and (4.47) are first compared to the laboratory data ki, by Savage (1953) for progressive waves propagating over a sandy seabed. The experiment was conducted in a wave tank of 29.3 meters (96 ft) long, 0.46 meters (1.5 ft) wide and 0.61 (2 ft) meters deep. The porous seabed was composed of 0.3 meters thick of sand with the medium diameter of 3.82 mm. The water depths are h = 0.229 m. and h = 0.152 m. The data for the water depth of h = 0.102 m was ignored for the reasons given by Liu and Dalrymple (1984). The wave conditions and the comparison of the damping rates are listed in Table 4.1. The parameters used in the solution of Eqs.(4.15) and (4.47) are: n = 0.3, ao = 570, b0 = 2.0 and C. = 0.46 and the average wave height H in the table was calculated according to Ho HLg) In Ho (1 Ho HLg where Ho and HLg are the wave heights measured at two points of Lg apart with Lg = 18.3 m (60 ft). The values computed by Liu and Dalrymple (1984) are also 42 listed. The relative error in the table is defined as A =1kim ,IX100% (4.60) kim with ki being either ki or kiLD. In this case, the errors are of similar magnitude between the present model and the model by Liu and Dalrymple (1984). This is because the experimental values fall in the region where the inertial resistance due to the virtual mass effect and the turbulent resistance, neglected in Liu and Dalrymple's model, are unimportant. In Fig. 4.2 through Fig. 4.5, the solutions of the complex dispersion equation are illustrated graphically. The case of progressive waves are demonstrated first. The specified conditions for the progressive waves are: H = 1 m, T = 4 seconds, h, = 5 m and n = 0.4. Figure. 4.2 plots the values of k, (wave number) and kic (damping rate) against R (nondimensional permeability parameter) for three different water depths, h = 2, 4, and 6 meters. The equivalent particle sizes for the range of R values are also shown in the figure; they cover a range from 2.3 mmn to 2.3 mx. The thick dlash line is the solution of the dispersion equation given by Liu and Dalrymple (1984, Eq.(4.5)) for the case of h -4 m. The correction term of the laminar boundary layer was not included in this curve since it is negligible in this case. From these results, a number of observations can be made: 1. As expected, the wave number decreases (or wave length increases) monotonically with increasing R, from one limiting value k,,, corresponding to the case of an impervious bottom at depth h, to the other limiting value k,2, corresponding to the case of a water depth equal to h + h, when the lower layer becomes completely porous. 2. Compared with the linear resistance, the nonlinear and the inertial resistances dominate for the complete range of R values displayed. In the region where only linear resistance dominates (Rf < 1, R,~ < 1), the bottom effects are relatively small and often negligible (outside the R range shown). Table 4.1: Comparison of The Predictions and The Measurements Run No. H (cm) T (sec.) ki,(m-') k1p(m-1) A,% kLDF m) ALD% h = 0.229 m 1 6.74 1.27 0.0379 0.0313 17.4 0.0338 10.8 2 4.45 0.0318 0.0337 6.0 0.0338 6.3 3 5.26 0.0303 0.0328 8.3 0.0338 11.6 4 1.93 0.0317 0.0369 16.4 0.0338 6.3 12 6.65 0.0334 0.0313 6.3 0.0338 1.2 13 4.36 0.0283 0.0338 19.4 0.0338 19.4 14 1.90 0.0374 0.0370 1.1 0.0338 9.6 5 6.25 1.00 0.0411 0.0350 14.8 0.0390 5.1 6 4.66 0.0357 0.0372 4.2 0.0390 9.2 7 2.17 0.0397 0.0411 3.5 0.0390 1.8 8 2.08 0.0393 0.0413 5.9 0.0390 0.8 15 5.48 0.0383 0.0360 6.0 0.0390 1.8 16 4.52 0.0397 0.0374 5.8 0.0390 1.8 17 1.93 0.0373 0.0415 11.3 0.0390 4.5 9 5.29 0.80 0.0296 0.0314 6.1 0.0343 15.8 10 3.55 0.0297 0.0334 12.5 0.0343 15.5 11 2.28 0.0335 0.0351 4.8 0.0343 2.4 18 7.05 0.0264 0.0296 12.1 0.0343 29.9 19 4.08 0.0320 0.0328 2.5 0.0343 7.2 20 2.14 0.0305 0.0353 15.7 0.0343 12.5 h = 0.152 m 21 4.00 1.27 0.0552 0.0574 3.9 0.0600 8.7 22 1.83 0.0379 0.0640 68.9 0.0600 58.9 23 3.83 0.0555 0.0579 4.3 0.0600 8.1 24 1.36 1.00 0.0449 0.0770 71.5 0.0731 62.8 25 3.40 0.80 0.0658 0.0700 6.4 0.0767 16.6 ki, the experimental damping rate given by Savage (1953); kip the theoretical values by the present model; k LD the theoretical value by Liu and Dalrymple (1984) 44 3. The wave attenuation, and hence the wave energy dissipation, shows a peak. This peak occurs when the magnitude of the dissipative force (velocity related) equals to that of the inertial force (acceleration related). Depending upon the water depth, this attenuation could be quite pronounced under optimum R values. For instance, for the case h = 4 m, Fig. 4.2 shows the wave height is reduced to about 74% of its original value (or about 46% energy dissipation) over approximately 2 wave lengths. 4. The locations of the peak damping are quite different from that of Liu and Dalrymple (illustrated here for the case of h= 4 m); they occur at a higher permeability. The magnitude of the peak damping is generally smaller than that of Liu and Dalrymple's solution. The values of k,. from the two solutions are also different. In Fig. 4.3 the maximum damping rate and the corresponding permeability parameter R are plotted as functions of nondimensional water depth. The curves from Liu and Dalrymple (1984) are also plotted for comparison. The trends of (kj),,'s are similar but the corresponding R behaves quite differently from the two solutions. The case of a standing wave system is also illustrated here in Fig. 4.4. The behavior is very similar to that of the progressive waves. The grain size in the figure is calculated according to a,. for the case of h = 4.0 m. Finally, the solutions of the dispersion equation for standing waves based on Darcy's model (DARCY: = 0, C1 = 0), Dagan's model (DAGAN: P3 = .1, C1 = 0), Dupuit-Forchheimer's model (D F: = 0 and C1 follows Eq.(4.46) and Sollitt-Cross's model (S C: 0 follows Eq.(3.30) with C1 following Eq.(4.46) and n = 0.4, ao = 570, b0 = 3.0 and C. = 0.46) for standing waves are compared in Fig. 4.5. Under the same wave conditions, the differences among the various solutions are seen to be very pronounced. for the permeability range displayed here. When the permeability parameter becomes less h h ............. h h PARTICLE SIZE (CM) 2.32 23.20 = 4.0 m Liu and Dalrymple H = 1.0 m = 2.0 m T = 4.0sec = 4.0 m h, = 5.0 m = 6.0 m I i I I V m -1.0 0.0 1.0 2.0 PERMEABILITY PARAMETER LOG(R) 2.32 23.20 -2.0 -1.0 0.0 1.0 2.0 PERMEABILITY PARAMETER LOG(R) 232.00 3.0 Figure 4.2: Progressive wave case. (a) Nondimensional wave number k,/(a2/g), (b) Nondimensional wave damping rate ki/(a'/g). 0.23 4 I.8 - 1.4 - 232.00 -2.0 0.23 0.10 0.08 0.06 0.04 0.02 0.00 .... .. . . . . . . .. .. .. .. .. .. .. .. .. I m 46 4.0 (L-D): Liu and Dalrymple (1984) 3.0 2.0. R .0 o ~ (k,) ma:(L-D) 0.0 " sR (L-D) -1.0 -2.0 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 Normalized water depth log(h a'2/g) Figure 4.3: Maximum nondimensional damping rate (ks).,=/(a2/g) and its corresponding permeability parameter R as functions of nondimensional water depth h ('/g). 47 than 10'2, all solutions converge to a single curving following Darcy's law. For very highly permeable seabeds, the damping rate based on D F model is very high and tends to increase with the permeability as oppose to approaching to zero according to common sense. Also in the high permeability region, the wave frequency based on Darcy's model is much larger than the correct value (determined by Eq.(4.17)). The reason for such large error is that the force balance of the pore fluid in highly permeable media is now mainly between pressure gradient and the inertia rather than the velocity related frictions. PARTICLE SIZE (CM) 2.28 22.69 ------ A = 4.0 m Liu and Dalrymple h = 2.0 m ............. h = 4.0 m ------ h = 6.0 m 225.44 H = 1.0m L = 20.0 m h, = 5.0 m 1.60 1.50 1 .40 I .30 1 20 1.10 -2.0 -1.0 0.0 PERMEABILITY i i 1.0 PARAMETER (a) 2.0 LOG (RI 0.23 2.28 22.69 225.44 0.01 0.0 -2.0 -1.0 0.0 1.0 2.0 3.0 PERMERBILITT PARAMETER LOG(RI Figure 4.4: Standing wave case. (a) Nondimensional wave frequency ar,/(L/g)1; (b) Nondimensional wave damping rate aj/(L/g)i. 0.23 - .---... -...--------..... .. ... 1.00 0.90 0.06 0.05 0.04 0.03 0.02 . . . . PARTICLE SIZE (CM) 2.28 22.69 S C ............. ODAGAN DARCY 0 F H = 1.0 m L = 20.0 m h, = 5.0m ---------------------- / / h=4m /' ..-/: ... ... ... ... ... ... ................................ " ....----------.. -1.0 0.0 PERMEABILITY 2.28 ~I ~ I I 1.0 2.0 PARAMETER LOG(R) 22.69 -2.0 -1.0 0.0 1.0 2.0 PERMEABILITT PARAMETER LOG(R) Figure 4.5: Solutions based on four porous flow models. (a) Nondimensional wave frequency o,./(L/g)i, (b) Nondimensional wave damping rate ai/(L/g)i. 0.23 .i. 225.44 1.4 1.3 1.2 -2.0 0.23 225.44 0.12 0. 10 0.08 0.06 0.04 0.02 0.00 S 5 C ,............. ORGAN /' \ .O... ARCY . \ .... DARCY /\ ./ .. - I t -. . CHAPTER 5 LABORATORY EXPERIMENT FOR POROUS SEABEDS The experiment was carried out in a wave flume in the Laboratory of Coastal and Oceanographic Engineering Department of University of Florida. The flume is about 15.5 meters long, 0.6 meters wide, 0.9 meters high and equipped with a mechanically driven piston-type wave maker. All the tests were conducted with standing waves. 5.1 Experiment Layout and Test Conditions Figure 5.1 shows the experiment arrangement. A porous gravel seabed was constructed at the end of the wave flume in the opposite side of the wave maker. A sliding gate was positioned at one wavelength from this end of the tank to trap the standing wave after a sinusoidal wave system was established by the wave maker. The decay of the freely oscillating standing wave was then measured by a capacitance wave gage mounted at the center of the compartment. The damping rate of the porous seabed for each particular test condition was determined by applying a least squares fit to the data according to the following equation: 1 Hi (o,)i = mCHln ) (5.1) Ti Hj-i where j refers to j-th wave. The corresponding wave frequency was obtained by averaging individual waves. The contribution to the wave damping due to the side walls and the bottom was subtracted from the data according to the following equation 1- ln[1I + e2aieT e2a;.T (5.2) 2T wave maker hs o **S L/2 L/2 Figure 5.1: Experimental setup Table 5.1: Material Information d5o(cn) f0.72 0.93 1.20 1.48 2.09 2.84 3.74 porosity 0.349 0.349 0.351 0.359 0.369 0.376 0.382 where u, is the actual seabed damping rate, a; is the gross damping rate of the seabed-wall system, ai,,, is the damping rate by the side walls and the bottom, and T is the wave period. The bed material used in the experiment was river gravel of seven sizes, ranging from d50 = 0.72 cm to 3.47 cm, with all sizes having a fairly round shape and smooth surface. The porosity of the material increases slightly with the diameter, as given in Table 5.1. With the grain sizes being determined by the material selection, the adjustable independent parameters left in the dispersion equation are the water depth h, the seabed thickness h,, and the wave length L. A total of 36 different cases was tested Table 5.2: Test Cases d5o (cm) h, (cm) h (cm) L (cm) 0.72 20 20,25,30 200 0.93 20 20,25,30 200 1.20 20 20,25,30 200 1.48 10,15,20 20,25,30 200 2.09 20 20,25,30 200,225,250,275 2.84 20 20,25,30 200 3.74 20 20,25,30 200 and these test conditions are summarized in Table 5.2. The damping effect due to side walls alone was found to be very small, and an average value of uri,, = -0.01 sec-I was used for all the corrections. 5.2 Determination of The Empirical Coefficients The experimental results for the conditions of h, = 20 cm and L = 200 cm are given in Table 5.3 (complete results are given in Appendix B). Each data point represents an averaged value of 10 to 20 tests. Figure 5.2(a) shows a typical example of measured waves and the exponential fitting to the wave heights. The solid line in Fig. 5.2(a) is the ensemble average whereas the two dash lines are the envelopes of plus and minus one standard deviation from the mean value at each time instant. From the wave form, it is noted that the waves in the experiment were more or less nonlinear waves as oppose to what was assumed. The effect of the nonlinearity to the calculation of the measured damping rate can be minimized by using the ratio of wave heights instead of wave amplitudes. In Fig. 5.2(b), the dashed curve is the exponential decay function with the a, obtained from the data ensemble by the least square analysis. Based on the measured a, and a, listed in Table 5.3, the empirical coefficients ao, bo and Ca were then determined by multi-variate linear regression analysis such Table 5.3: Measured a, and ai for h, = 20 cm and L = 200 cm h (cm) dso (cm) R; (cm) Y" (cm) o,(s-') a',g(s-1) o,(s-1) 0.72 11.08 6.49 4.9054 -0.0686 -0.0570 0.93 10.93 6.35 4.9402 -0.0678 -0.0562 1.20 10.47 5.52 4.9370 -0.0814 -0.0694 30.0 1.48 9.03 4.71 4.9620 -0.0878 -0.0757 2.09 9.75 4.85 4.9673 -0.0948 -0.0824 2.84 9.80 4.95 5.0111 -0.0931 -0.0808 3.74 9.80 5.71 5.0341 -0.0658 -0.0543 0.72 8.70 5.15 4.6218 -0.0799 -0.0678 0.93 8.51 4.75 4.6856 -0.0905 -0.0781 1.20 8.00 4.46 4.6793 -0.0984 -0.0858 25.0 1.48 7.28 3.91 4.7352 -0.1162 -0.1030 2.09 6.46 3.57 4.7422 -0.1105 -0.0975 2.84 7.37 4.34 4.7893 -0.0940 -0.0816 3.74 7.37 4.26 4.7865 -0.0813 -0.0693 0.72 6.70 4.19 4.2716 -0.1184 -0.1047 0.93 6.28 3.55 4.3426 -0.1334 -0.1192 1.20 6.07 3.60 4.3565 -0.1320 -0.1179 20.0 1.48 4.49 2.57 4.4027 -0.1484 -0.1337 2.09 5.19 3.12 4.4325 -0.1546 -0.1396 2.84 5.41 3.33 4.4728 -0.1456 -0.1311 3.74 5.56 3.58 4.4646 -0.1313 -0.1173 * H1 is the average value of the first wave heights of the data group and H is the average height of the complete train. 54 0.60 z H = 5.29CM o - 0.40 L = 200. CM > .T = 1.27 5 , 0.20 U U- 0.00. o -0.20 = -0.40 O z -0.60 0.0 4.0 8.0 12.0 16.0 20.0 TIME (SEC) 1.20 ODW= 30.0 CM 1.00 S.. = 15.0CM 00= 1.48CM LI-i = 0.80 . 2: 0.60 .. ... C LIN4 ; 0.40 ...-... .. S0.20 . 0.00 0.0 4.0 8.0 12.0 16.0 20.0 TIME (SEC) Figure 5.2: Typical wave data: (a) Averaged nondimensional surface elevation (7l/Hi), (b) Nondimensional wave heights (1/H) I and the best fit to the exponential decay function. that the error function defined as = 1 r p)2 + (Uim iP)2] (5.3) _1=1 arm im was minimized, where the subscript m represents the measured values and p denotes the predicted values; M is the number of data points. The best fit was found to exist when: ao = 570 bo = 3.0 C, = 0.46 The added mass coefficient obtained here, C, = 0.46, is close to the theoretical value of 0.5 for a smooth sphere. The values of a0 and b0, to an extent, re-confirm those given by Engelund (1953) and the others. 5.3 Relative Importance of The Resistances in The Experiment Introducing a0, bo, C., and the averaged wave height F into Eq.(4.29) for I qI and using Eqs.(3.33) through (3.35), the relative importance of the various resistances of the tested cases can be established precisely. The results are given in Table 5.4. As we can see from these ratios, for the first four grain sizes, all three resistances are about equally important. Whereas for the last two larger diameter materials, the turbulent and the inertial resistances are evidently dominant over the linear resistance, with the inertial force approximately double that of the nonlinear resistance. 5.4 Comparison of The Experimental Results and The Theoretical Values Table 5.5 shows the comparisons between experimental results and the theoretical values of ai and a,. The relative error, defined as A% =1 a I x100% (5.4) am Table 5.4: Comparison of The Resistances. h, = 20 cm, L = 200 cm. d5o (cm) q I (cm/s) -c I R IR, R If,/f, If.if,'f, /fi h = 30.0 cm 0.72 0.57 4.788 37 225 0.93 1.32 1.42 0.93 0.69 4.800 58 379 1.56 2.06 1.32 1.20 0.75 4.819 81 630 2.63 2.88 1.10 1.48 0.78 4.842 104 960 4.18 3.71 0.89 2.09 0.96 4.869 181 1933 8.90 6.48 0.73 2.84 1.09 4.889 281 3577 17.12 10.10 0.59 3.74 1.33 4.900 453 6240 30.90 16.34 0.53 h = 25.0 cm 0.72 0.53 4.526 34 213 0.88 1.21 1.38 0.93 0.62 4.545 52 358 1.48 1.85 1.25 1.20 0.70 4.567 76 597 2.49 2.71 1.09 1.48 0.75 4.596 101 915 3.99 3.60 0.90 2.09 0.84 4.639 159 1833 8.44 5.70 0.67 2.84 1.10 4.657 284 3407 16.31 10.21 0.63 3.74 1.18 4.679 400 5959 29.51 14.43 0.49 h = 20.0 cm 0.72 0.50 4.180 32 196 0.81 1.14 1.41 0.93 0.55 4.206 46 332 1.37 1.64 1.20 1.20 0.66 4.232 72 553 2.31 2.56 1.11 1.48 0.61 4.286 81 849 3.70 2.89 0.78 2.09 0.85 4.322 160 1707 7.86 5.73 0.73 2.84 1.00 4.356 258 3187 15.26 9.27 0.61 3.74 1.15 4.381 393 5579 27.63 14.18 0.51 Table 5.5: Comparison of Measurements and Predictions for h, = 20 cm, L = 200 cm. d,,) (cm)] 1? Urm( .. o'.(-)A% (S )1 aip(S- ) AjZ _______ h = 30.0 cm __ _ 0.72 0.1791 4.9054 4.7880 2.39 -0.0570 -0.0541 5.07 0.93 0.3021 4.9402 4.8001 2.84 -0.0562 -0.0628 11.60 1.20 0.5110 4.9370 4.8187 2.40 -0.0694 -0.0706 1.66 1.48 0.8448 4.9620 4.8419 2.42 -0.0757 -.0.0751 0.76 2.09 1.8755 4.9673 4.8685 1.99 -0.0824 -0.0734 11.00 2.84 3.7429 5.0111 4.8891 2.43 -0.0808 -0.0661 18.26 3.74 6.9543 5.0341 4.9001 12.66 -0.0543 -0.0613 12.85 h = 25.0 cm___0.72 0.1687 4.6218 4.5264 2.06 -0.0678 -0.0709 4.49 0.93 0.2866 4.6856 4.5445 3.01 -0.0781 -0.0840 7.49 1.20 0.4843 4.6793 4.5666 2.41 -0.0858 -0.0935 8.97 1.48 0.8095 4.7352 4.5963 2.93 -0.1030 -0.1000 2.94 2.09 1.7819 4.7422 4.6389 2.18 -0.0975 -0.0956 1.94 2.84 3.5772 4.7893 4.6566 2.77 -0.0816 -0.0912 11.76 3.74 6.6123 4.7865 4.6795 12.24 -0.0693 -0.0781 12.67 h = 20.0 cm______0.72 0.1559 4.2716 4.1799 2.15 -0.1047 -0.0907 13.34 0.93 0.2656 4.3426 4.2061 3.14 -0.1192 -0.1109 6.96 1.20 0.4509 4.3565 4.2318 2.86 -0.1179 -0.1228 4.16 1.48 0.7496 4.4027 4.2857 2.66 -0.1337 -0.1319 1.32 2.09 1.6656 4.4325 4.3217 2.50 -0.1396 -0.1308 6.33 2.84 3.3408 4.4728 4.3560 2.61 -0.1311 -0.1200 8.42 3.74 6.1676 14.4646 14.3805 11.88 -0.1173 -0.1076 18.24 is uniformly very small for a, (wave frequency). For a, (wave damping), the range of error is larger, with the maximum being 18.26% of the measured value. The agreement between the predictions and the measurements shown in the above table is demonstrated by Fig. 5.3 where the values of a7, and ai in Table 7 for the case of h = 25 cm are plotted against the permeability parameter. 10.0 -1.0 0.120 0. 100 0.080 0.060 0.040 0.020 0.000 -1.0 -0.5 0.0 0.5 PERMEABILITY PARAMETER LOG(R) -0.5 0.0 0.5 PERMEABILITY PARAMETER LOG() Figure 5.3: The Measurements and the predictions vs. R. for L = 200.0 cm, h, = 20.0 cm, h = 25.0 cm. (a) Wave frequency at, (b) Wave damping rate oj. WAVE LENGTH = 200 CM BED THICKNESS = 20 CH WATER DEPTH = 25 CH a0 & 0 + + * PREDICTED x MEASURED -I Table 5.6: Comparison of a, and ap for do = 1.48 cm, L = 200 cm. h (cm) IH (cm) I rm.,(s-) I o7.(s-1) IA,.% Iaim(s- ) Icip(S- )A% h. = 15.0 cm 30.0 5.29 4.9268 4.8291 1.98 -0.0525 -0.0600 14.13 25.0 4.00 4.6715 4.5836 1.88 -0.0697 -0.0794 13.98 20.0 3.12 4.3173 4.2565 1.41 -0.0928 -0.1054 13.56 h. = 10.0 cm 30.0 5.97 4.9163 4.8145 2.07 -0.0368 -0.0416 13.13 25.0 4.61 4.6310 4.5629 1.47 -0.0492 -0.0551 11.92 20.0 3.48 4.2688 4.2309 0.89 -0.0724 -0.0729 0.70 Table 5.7: Comparison of 0, and oap for d5o = 2.09 cm, h, = 20 cm. h (cm) H (cm) I ) crpC5T I A% Uim(s ) I i(s- ) IZ% L = 225.0 cm 30.0 4.42 4.5114 4.4401 1.58 -0.0831 -0.0766 7.73 25.0 3.91 4.2913 4.1991 2.15 -0.0977 -0.1006 2.92 20.0 2.96 3.9699 3.8995 1.77 -0.1270 -0.1286 1.21 L = 250.0 cm 30.0 4.40 4.1524 4.0725 1.92 -0.0752 -0.0790 5.02 25.0 3.79 3.9119 3.8358 1.94 -0.1010 -0.1002 0.81 20.0 3.32 3.5862 3.5361 1.40 -0.1384 -0.1282 7.34 L = 275.0 cm 30.0 4.18 3.8371 3.7592 2.03 -0.0783 -0.0788 0.68 25.0 3.32 3.6285 3.5339 2.61 -0.0993 -0.0962 3.15 20.0 2.90 3.3298 3.2497 2.41 -0.1264 -0.1211 4.16 We now proceed to compare the theoretical computations with those data which were not included in the calibration. These results are given in Table 5.6 and Table 5.7. The agreement is just as good, if not better. In Fig. 5.4, the theoretical values of a,'s and a,'s by the present model are plotted, respectively, against the measurements for all the experimental data (36 cases). The overall least mean square error defined in Eq.(5.3) for the data of all 36 cases is about 0.008. While the solution of the dispersion equation (Eq.(4.21) based on the linear porous flow model by Liu and Dalrymple (Liu and Dalrymple, 1984) Table 5.8: Comparison of a,, and op for ds0 = 0.16 rm, h, = 20 cm, L = 200 cm. h (cm) R a,.m(s ) arp(-S) A,% aim(s5 ) 30.0 5.4 x 10. 4.8908 4.7638 2.60 -0.0070 -4.5 x 10.7 25.0 5.1 x 10-7 4.6312 4.4957 2.93 -0.0090 -5.6 x 10-7 20.0 4.8 x 10-7 4.2824 4.1428 3.26 -0.0122 -6.8 x 10-7 is plotted in Fig. 5.5. In obtaining the solution, the only empirical coefficient ao appeared in their dispersion equation was calibrated to the data in Table 5.3 with the same routine as that for the solution of Eq.(4.15), and the nonlinear regression yielded a0 = 2200. The overall least mean square error for the complete data set (36 cases) is as high as 0.327. It can be seen from the figure that considerable discrepancies occur for the damping rate, and the maximum relative error is about 87% of the experimental values. The prediction of wave frequency by their solution is about the same as that by this study. A test over a sand bed (d50 = 0.16 mM) was also conducted. Active sand movement was observed along with the rapidly formed ripples. As a result, theoretical predictions are several orders of magnitude smaller than those actually measured (Table 5.8). The basic assumption of no particle movement was violated. However, the predictions of wave frequencies were as good as the coarse-grained cases. The same test, with dye injection into the sand bed, disclosed that the flow inside the sand bed was very slow, and indicated that the assumption of impermeability invoked in the linear wave theory is justified in such situations. When the thickness of a gravel bed is reduced to one grain diameter, the model also failed to yield good results as shown by the values in Table 5.9. In this case, the effect of the boundary layer at the interface, which was not included in the model, obviously cannot be ignored. The nature of a turbulent boundary layer over a thicker seabed was qualitatively investigated for the case of h,= 20 cm, ds0 = 1.48 cm, h =20 cm and L =200 cm, 2.0 3.0 4.0 5.0 EXPERIMENTAL FREDUENCT 0.25 0.20 0.15 0. 10 ., 0.05 " 0.00 p 0.00 0.05 0.10 0.15 EXPERIMENTAL DAMPING 0.20 0.25 RATE Figure 5.4: Theoretical values by the present model vs. experimental data of Table 5.5, Table 5.6 and Table 5.7. (a) Wave frequency at, (b) Wave damping rate a,. I 5.0 4.0 3.0 -1 1 t I 1 6 i 2.0 3.0 4.0 5.0 EXPERIMENTAL FREQUENCY 0.25 0.20 0.15 0.10 0.05 0.00 0.00 0.05 0.10 0.15 0.20 0.25 EXPERIMENTAL DAMPING RATE Figure 5.5: Theoretical values by the model of Liu and Dalrymple vs. experimental data of Table 5.5, Table 5.6 and Table 5.7. (a) Wave frequency or, (b) Wave damping rate a0. ~1' A. * .' * A .** . .4 Table 5.9: Comparison of crn. and ap for h, = d5o d50 (cm) H (cm) crn(S -) Ir.5 -) I r% UI(,(s-) rp(s -) TA,% h = 30.0 cm 0.72 6.05 4.9157 4.7666 3.03 -0.0024 -0.0031 29.66 1.20 5.88 4.9374 4.7723 3.34 -0.0040 -0.0041 3.04 2.09 5.58 4.9559 4.7827 3.50 -0.0075 -0.0042 43.69 3.74 5.21 4.9808 4.8002 3.63 -0.0097 -0.0049 48.99 h = 25.0 cm 0.72 4.64 4.6495 4.4991 3.23 -0.0031 -0.0041 31.58 1.20 4.53 4.6675 4.5068 3.44 -0.0055 -0.0056 0.96 2.09 4.28 4.6822 4.5205 3.45 -0.0100 -0.0056 43.88 3.74 4.02 4.7318 4.5438 3.97 -0.0130 -0.0064 50.59 h = 20.0 cm 0.72 3.24 4.2936 4.1470 3.41 -0.0057 -0.0053 7.19 1.20 3.23 4.3134 4.1571 3.62 -0.0094 -0.0075 20.32 2.09 3.07 4.3339 4.1755 3.66 -0.0147 -0.0074 49.62 3.74 2.88 4.3950 4.2061 4.30 -0.0170 -0.0081 52.29 by dye studies. Turbulent diffusion was spotted over virtually the entire porous domain and up to approximately 1.0 cm above the interface. The effects of turbulent boundary layers to the wave damping are important in wave interaction with porous seabeds and more theoretical and quantitative experimental research is needed. CHAPTER 6 BOUNDARY INTEGRAL ELEMENT METHOD The Boundary Integral Element Method (BIEM) is an efficient numerical method for conservative systems such as those governed by Laplace equation. It is efficient because the computation is carried out only on the boundaries rather than over the entire domain as it would be in finite difference and finite element methods. Considerable amount of computation time and storage space can be saved without losing accuracy. In addition to the efficiency, the data preparation becomes much simpler as compared to the other two methods. For non-Laplace equations, the application of this method may become difficult and sometimes impossible. By the efforts of many scientists, the scheme has been extended to more and more problems governed by non-Laplace equations (Brebbia, 1987). In this chapter, we restrict ourselves to the two dimensional Laplace equation only. 6.1 Basic Formulation Let U and V be any two continuous two dimensional functions, twice differentiable in domain D which has a close boundary C. According to the divergence theorem (Loss, 1950; Franklin, 1944, cited in Ligget and Liu, 1983) for continuity in a volume fD(V ifC6 d (6.1) where V is the vector operator, V is any differentiable vector and n- is the unit outward vector normal to C. If we define that 6 = UVV and then V = VVU in Eq.(6.1), two integral equations can be generated in terms of U and V, f(VU VV + UV2V)dA = UVV i.s ( C (6.2) 65 (VV VU + VV2U)dA= VVU- iids (6.3) Subtracting (6.3) from (6.2), we have (UV'V VV'2U)dA = (UVV VVU) d (6.4) Introducing the expressions VV au V n = VU n = n (6.5) and assuming that both U and V satisfy Lapalce equation, i.e. V2U V2V = 0 (6.6) Equation (6.4) becomes r av au (U V an )ds = 0 (6.7) To apply this equation, we choose U as the velocity potential 4 and V as a free space Green's function, G. Both of them satisfy Lapalce equation. Equation (6.7) can then be rewritten, in terms of 1 and G, as 8G(PQ) 8(Q)ld. C[.D(Q) aG(PQ) G(P,Q) ]ds = 0 (6.8) where P is a point in the domain DnC and Q is a point on C. One of the free space Green's functions for two dimensional problems is G(P,Q) = In r (6.9) where r is the distance between point P and point Q and it can be expressed by r = V(zp ZQ)2 + (zp zq)2 on the z z plane. Substituting Eq.(6.9) into Eq.(6.8), it becomes #(Q) Br 8-(Q) [ r In r (Q) ]ds = 0 (6.10) r 5n an 66 Despite the fact that the Green's function has a singularity at r = 0, the contour integral in Eq.(6.10) exists and can be worked out by removing the singular point from the domain. In general, Eq.(6.10) becomes C In r( ]d (6.11) cz(P) TfQ)ar8 after removing the singularity at point P. Here a is 27r if P is an interior point and is equal to the inner angle between the two boundary segments joining at P if it is on the boundary. Generally, point P can be anywhere in the domain DnC, while Q is always on C due to the contour integration. In BIEM, since only the boundary values of D and a are solved, P has to be kept on the boundary C in an the process of computation. When 4 and an on C are solved, the interior values can be derived with Eq.(6.11) by placing P at the point of interest inside D. The closed boundary C in Eq.(6.11) is the same as that in Eq.(6.10) except that it does not pass point P as it does in Eq.(6.10). The singular point of the Green's function r = 0 has been excluded by a circular arc of infinitesimal radius from the domain DnC. Thus there is no singularity on the new contour C. Eq.(6.11) is the basic equation for BIEM formulations. Discretizing the boundary C into N segments, and breaking the contour integration into N parts accordingly, Eq.(6.11) reads N 4D ar, ao N =-I ri5)ds=Ei (6.12) 3=1 arn an1 The curved segment C is usually replaced by a straight line to simplify the integration. This simplification generally does not introduce significant error provided that the segment is small enough. Before carrying out the integration, the type of element has to be determined. The formulations are different for different types of element. The commonly used elements in two dimensional BIEM problems are: constant element, linear element, 67 quadratic element, special element and so on so forth. The classification of elements is based on the type of the function used to interpolate the values over an element. For example, on a constant element, the physical quantities and their normal derivatives are assumed to have no change over the element; while on a linear element, the quantities and their normal derivatives are interpolated by a linear function between the values at the ends of the element, i.e. = N1( )D-F g(C )4y+1 (6.13) o-- (C) = Ng (e)( a ) i+ N2( alp )j+l (6.14) Cj C ej+i where N, and N2 are linear functions with the feature of N, (C.) = 1 N1(ej+1) = 0 g2( ) = 0 N2(e,+1) = 1 In the same manner, higher order elements can be defined by replacing N1 and N2 with higher order functions. When the variation of the quantities is known, the variation function can then be chosen as the interpolation function to obtain a more accurate element. Such elements are usually classified as special elements. Generally, it is difficult to judge which element is superior over the others without analyzing the particular problem. As a rule of thumb, the higher the order of an element, the more accurate the approximation would be for the same element size. But for a higher order element, as a trade-off for precision, the formulation would be much more complicated and tedious, and the computation might be much more time consuming. The constant element works fairly well for problems with smooth boundaries and continuously changing boundary conditions, but could generate considerable errors at 'corner points' where the boundary conditions are not continuous. For boundary with such 'corners', the linear element is usually a better choice. / / pi Figure 6.1: Auxiliary coordinate system 6.2 Local Coordinate System for A Linear Element To formulate a problem with linear elements, an auxiliary local coordinate system, as shown in Fig. 6.1, has to be established to facilitate line integrations over the element. The two axis C j7 are perpendicular to each other with one lying on the element and pointing in the direction of integration, from Pj to Pj+1, and the other one being in the same direction as the outward normal vector ii. Based on this auxiliary coordinate system, the integration over each segment can be completed analytically and be expressed in terms of 4., 4j+j and t7, the distances from P, P+i and P to the local origin defined in Fig. 6.1. The values of 4, e,+i and j7i are determined by the global coordinate information of three points. Since S+ 171 = ? (+ AX) + 171 = rt," +i where 4,,+4 = vi Xi+z) + (i zi+) then . -"<+ (6.15) 2Ay + = fi + A i (6.16) The value of mi is therefore r/,I= ,/ -C (6.17) However, the numerical test showed that Eq.(6.17) could lead to significant error of j iri due to runoff errors. A more accurate expression for I 7i j can be obtained by computing the distance from Pi to the straight line PyP+ directly from the global coordinates, i.e. i ne I + z zi (6.18) I r/, I-- +A, where A = zj+1 zy zy+1 zy The sign of ri depends on the relative position of Pi to the boundary line element PjPi+1I If P is in the same side of the element with ii, 7i is positive, otherwise, it is negative. Mathematically, it can be expressed by the following two equations for vertical and nonvertical elements: sign (77i) = (sy Xi) Az +y4 = xy (6.19) (xj x,) Az sign (7i) = AI (Azio X+,) Az 1 x (6.2019) Ax Azio sign (ri) = AxzI x1+1 x (6.20) where Ax = xY+ Xi Az = zy+1 zy Az Azio = (X- ) +z -z Therefore (6.21) ri = sign (7i) |i, I 70 6.3 Linear Element and Related Integrations By the definition stated in the last section, on a linear element, a quantity and its normal derivative are interpolated by a linear function between the values at the two nodes, P and P+1. For the velocity potential function 4 and its normal BA derivative we have, in terms of 9 coordinates, an ( +(i+- i) + ~ ii+) fi< C i+1 (6.22) &+1 as as as as =[(--)j+l (a-)i]C + [+ e(a+1- ) 5-n M ej(6.23) ej+1 ej+ It is obvious that a= ( D) ( e. ) = Il ) ; = 1 CC+1) =C )+ Substituting Eqs.(6.22) and (6.23) into Eq.(6.12), and noting that r, = 7, + 2 (6.24) Bri ar= r (6.25) an= a = ri the integration over segment j (between Pj and Pi+1) can be carried out analytically in the local coordinate system (Ligget and Liu, 1983): le (r In r, )a d [1 ap 2 ( = Ill + H*2 >+1 Kix( )t K ( (6.26) where the superscripts and + refer to the positions immediately before and after the nodal point denoted by the corresponding subscript. And (6.27) H = -I1- + e1+2 71 I = I 22. Ki4 = -, + (j+1llj Ki, = Ii2j -( l 2 _ 1 fi+1 r = fi+ -fifi ri 8n i- i,7- 2+1 2As 17+?+ = 1~j fJ +, ti 1 ari d = fi+ ri 8n 2Ai fiJ 1 [t an- (+1) Yi Vi (6.32) - tan-'( )] 'i S 1 f Inri+d2 fi+1 fi fi - {(,? + C+1)[1n(7 + C+1) 1] + + -(,7, + C)[11n(, + I) 1]} 1 f +In ri d( -~ 1nn + e;-C 'i7~+~ 2~+ j 1 2 s 2 = ~ {.i(77? + ex)-6n(to + ) 2(e6 1- e) 2A( +2rjj[tan-l(e7+l) -7t7-(,)]} th 77- with (6.28) (6.29) (6.30) -i 2A, I CJ+ i e ;dC r. (6.31) d xi r? (6.33) (6.34) 72 for all ij and i#j+l and Ai = fi+x fi If i = j, the above integrals have to be re-evaluated to avoid the singularity of the Green's function, 1 I j = lim 1 cI = 1 lin ;j+1 i e0 :i+: 'o 1 = lim ij+1 ieo J i+i 1 r d Ji+1 1 ar 0 + dg = 0 ei+, ri a fI +1 1 B i+, In r d fi+ A(2InA- 1) 1 l i + d = lirm+;- In r de + 1 ~ -0 44+. = lnA- 1 Similarly, when i = j + 1, - fiEfz-+ a ri 1lim la d( = 0 fi+ ~ iC-0 t i r8n S +1 i e-o 1 ri an 1 lim I nried( j +1 i e- = (1-2lInAj) 4 (6.35) (6.36) (6.37) I1+11 1+11J (6.38) (6.39) (6.40) (6.41) -; 1i ur ++1 -i = InA -1 lnri d( (6.42) Summarizing all the If's according to Eq.(6.12) for node points i = 1 -- N, and assuming that (6.43) (6.44) 80 at all the nodes, a system of linear equations with unknowns of P and a can be anobtained: obtained: N = L(8n)i j= 1 P4H1= ,N + THil bi, 1 C A?j = Hi... + Hj - 1 ifi=j 4.- = 0 if i # j Li, = K, N + K,1 LiE = K_ i1 + Ke ai j = 2, 3, ..., N Equation (6.45) can also be expressed by a matrix equation where (6.45) with (6.46) (6.47) (6.48) = (-) 74 where 4' and q,, are the vectors of order N x 1 containing the unknown T and L respectively. The coefficients in R and L are determined solely by the boundary geometries. 6.4 Boundary Conditions In Eq.(6.45), there are N linear equations and 2N unknowns. The additional N equations necessary for obtaining a unique solution are usually introduced by the application of the boundary conditions. There are essentially three types of boundary conditions for boundary value problems. I i = T (6.49) II = -- (6.50) an an III = k5 (6.51) n where D, .'n' f and k are known functions on the boundary. In water wave problems, the first relation usually means that the pressure distribution is known, and the second describes a given flux distribution on the boundary. They are also referred as 'essential' and 'natural' conditions (Brebbia,1989), respectively. The third one some times takes place when none of the quantities in the first two equations are specified but only the relationship between them can be found. A typical example of this is the free surface boundary condition in linear wave theory. For nonlinear wave problems, the right hand side of the above equations may include other known functions. As might be noticed in Eq.(6.45) or Eq.(6.48) there are two unknowns for each nodal point. Since we have obtained one equation for every node on the boundary, only one of the three relations needs to be specified at each node, if the assumptions in Eqs.(6.43) and (6.44) are strictly satisfied. However, for some nodal points where the boundary changes directions, ('90)- may not be equal to (P)+. For such cases, more than one boundary conditions may be necessary for one node. The treatment 75 of such points will be discussed in the next chapter. By introducing the boundary conditions at all N boundary nodal points, another system of N equations with the same unknown in Eq.(6.48) is established. The number of unknowns is now equal to the number of the equations. It is usually convenient to eliminate N unknowns with the 2N equations, and the resulted matrix equation can be expressed as AX=b (6.52) in which A is a known matrix of order N x N, X is the unknown vector containing 4 or &D or D for some part of the boundary and a for the rest of it, depending on which one is not specified by the boundary condition; and b is a known vector resulted from boundary conditions. CHAPTER 7 NUMERICAL MODEL FOR SUBMERGED POROUS BREAKWATERS In this chapter, a numerical model of BIEM for submerged porous breakwaters is developed by using the unsteady porous flow model given in Chapter 3. The basic function of the numerical model is to compute the wave flow field and related quantities such as the wave form, the dynamic pressure and the normal velocity along all the boundaries and so on. With these quantities, the wave transmission and reflection coefficients and wave forces can then be obtained. 7.1 Governing Equations The computation domain of the problem, as shown in Fig. 7.1, consists of two sub-domains, the fluid domain and the domain of the porous medium. In the fluid domain, the water is considered inviscid and incompressible. The flow induced by gravity waves is assumed irrotational. Thus, the governing equation in this domain, for the velocity function 41, is the Laplace equation, vl.(D 0(7.1) with fluid velocities being defined as U = (7.2) ax W, = -aD (7.3) While in the porous medium domain, the viscosity of the fluid cannot be ignored since the flow is largely within the low Reynolds number region. The flow is 0 0 0 0 0 0 0 0 0 O Figure 7.1: Computational domains described by the linearized porous flow model, pafof= -VP (7.4) where P is the pore pressure, fo is the linearized resistance coefficient and q is the discharge velocity vector in the porous medium, q = iu + kw with iand j being the unit vector in z and z directions. Taking curl of Eq.(7.4), pafoV x f= -Vx VP = 0 (7.5) leads to Vx= 0 (7.6) This means that the homogenized porous flow is irrotational, and hence a velocity potential function 0, in the porous domain exists such that T= V , (7.7) According to Eq.(7.4) and Eq.(7.7) P -po 0 (7.8) here p, a and fo are all treated as constants. Substituting Eq.(7.4) into the continuity condition for the porous flow in terms of the discharge velocity, V. 0 (7.9) the governing equation for the porous medium domain becomes the Laplace equation in terms of the pore pressure, V2P=0 (7.10) In this study, the waves are assumed to follow the linear wave theory = Oeit (7.11) P = Pe iat (7.12) It is to be noted that the time function is now expressed as eiat, as opposed to eit used in the seabed problem presented in previous chapters. 7.2 Boundary Conditions The whole computation domain has four types of boundaries, free surface, impermeable bottom, permeable interface of different sub-domains and the artificial lateral boundaries. These boundary conditions are discussed here. 7.2.1 Boundary Conditions for The Fluid Domain 1. The free surface The combined linear free surface boundary condition in the linear wave theory is a = a (7.13) wz 9 where z is the vertical coordinate and g is the gravity acceleration. 2. Impervious bottom The boundary condition on an impermeable bottom is simply nonflux condition: a_ = 0 (7.14) n where n is the direction normal to the boundary which can be of any shape-sand bar, or soil trench and so on, as long as it can be considered impermeable. 3. The permeable interface The permeable interface is the common boundary either between the fluid domain and the porous domain, or between two different porous domains. For the sake of simplicity, we consider only one porous sub-domain in this chapter. The same formulation can be easily extended to multi-porous-domain configuration. The boundary conditions on such a boundary are the continuity of pressure and mass flux. They are p = P (7.15) at or equivalently iap'- p (7.16) with the linear wave assumption, and .a = 1 ap (7.17) an, pafo an~i where n1 is the outward normal of the fluid domain and nn is the outward normal of the porous domain. However, these two notations are only used when confusion could occur, otherwise, n without subscript always denotes the outward normal vector of the domain in discussion. 4. The lateral boundaries There are two vertical lateral boundaries, one is on the offshore side and one is on the lee side of the porous domain, as shown in Fig. 7.1. In this study, the radiation boundary condition for these two boundaries is adopted. Such a boundary 80 condition assumes that the waves at these two boundaries are purely progressive, and that the decaying standing waves generated by the object inside the domain are negligible at these boundaries. Comparing with the matching boundary condition with the wave maker theory employed by Sulisz (1985) and some other authors, the radiation boundary condition offers significant simplicity in programming and provides sufficient accuracy with much less CPU time. In applying this boundary condition, the two lateral boundaries have to be placed far enough from the structure. Numerical tests show that a distance of about two wave lengths from the toe of the structure is more than adequate for this purpose. At these two lateral boundaries, the potential functions for the transmitted and reflected waves are assumed to have the forms OrC(,z) = e'k(z+'l)R(z) (7.18) 't(x,z) = e-k'('-")T(z) (7.19) where T(z) and R(z) are two unknown functions of z; k and k' are the wave numbers at the two boundaries respectively, and the subscripts t and r here refer to the transmitted and reflected waves, respectively. On the lateral boundary of up-wave side, the potential function for the incident wave is known: 01(x, z) = e-'Ck(+')A(z) (7.20) with A(z) = Hg cosh k(z + h) (7.21) 2a cosh kh Therefore, = ki+ 01 at x-- (7.22) 0= 4 at x=l' (7.23) where I and are, respectively, the distances from the origin of the coordinates to the offshore and onshore lateral boundaries and 0 is the unknown potential function. 81 On the lateral boundary of transmission (onshore side), x = 1', (g -ik'', = -ik'4 (7.24) an ax ax in which k' is determined by gk'tanh k'h' = a2 (7.25) where h' is the water depth at x = 1'. While on the lateral boundary of reflection (offshore side), x = -1, ,. = 4 01 (7.26) ao a 8 a4,. ik, iko, an -x ax ax - ik4' ik(O 0') = 2ik, ikO or a= 2iko' ikO (7.27) an in which k is the incident wave number determined by gk tanh kh = a2 (7.28) where h is the water depth at x = -1. 7.2.2 Boundary Conditions for The Porous Medium Domain In this domain, there are two types of boundaries, one is the common boundary with the fluid domain and the other one is the impervious bottom. On the interface, the boundary condition is the same as the one specified by Eq.(7.16) and (7.17). On the impermeable bottom, the no-flux condition, in terms of the pore pressure, is .p= 0 (7.29) n= Again, this boundary can have any shape, it can be the surface of the domain for core materials if it is considered as impermeable, etc. 82 7.3 BIEM Formulations The formulations in the two domains are slightly different because of the different boundary conditions. In the fluid domain, due to the radiation boundary condition on the lateral boundary of offshore side, the terms containing the incident wave potential are introduced, which will form the RHS vector of the matrix equation. On the free surface, the normal derivative are expressed in terms of owing to the CFSBC. In the porous sub-domain, the matrix equation, after applying the no-flux condition, has to be manipulated to match with the fluid domain. 7.3.1 Fluid Domain According to the formulae given in chapter 6, Eq.(6.12) can be expanded, in terms of node values, assuming that 4- = + for all the nodes, as: jOf = (H,1 + HN)41 + (Ha + H2)02 + (H?2 + Hg3)s... + (HIN,1 + H N)ON, + ... + (~U.2 1 + H.)+N2N . + H~3..1+ 4 3)ON3 + + (H 1 +X HNb + (H, .-~ + H ..)ON.. + ... + (H + H1)qN - KS,+1- Ks,2- Kg+f- K320 K,35+ -... - KFN-11, K ,1 N +, -...- K~.2--4_f K + % ,N3-_10;- K,, + .. K,?,,,, n- K ++,-. - F K 2 - K~NN..-1;N. K, N,.4N,,. ...- KN-1N K.,N;l (7.30) i = 1, 2 ......, N where ai is the inner angle of node i and the subscripts 1, NI, N2, N3, Nes and Ne,, refer to the nodes of 'corner points' as shown in Fig. 7.1. The superscripts and + refer to the positions immediately before and after the nodal point in the direction of contour integration, and the subscript n refers to the normal derivative. Here 0is not necessarily equal to 0+ for all the nodes. |

Full Text |

xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID ES9FOVIO7_7Q93EH INGEST_TIME 2017-07-14T23:05:50Z PACKAGE UF00075474_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES xml record header identifier oai:www.uflib.ufl.edu.ufdc:UF0007547400001datestamp 2009-02-20setSpec [UFDC_OAI_SET]metadata oai_dc:dc xmlns:oai_dc http:www.openarchives.orgOAI2.0oai_dc xmlns:dc http:purl.orgdcelements1.1 xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.openarchives.orgOAI2.0oai_dc.xsd dc:title Water wave interaction with porous structures of irregular cross sectionsUFLCOEL-TR dc:creator Gu, Zhihao, 1957-University of Florida -- Coastal and Oceanographic Engineering Deptdc:subject Coastal and Oceanographic Engineering thesis Ph. DDissertations, Academic -- Coastal and Oceanographic Engineering -- UFdc:description b Statement of Responsibility by Zhihao GuThesis Thesis (Ph. D.)--University of Florida, 1990.Bibliography Includes bibliographical references (leaves 197-199).Typescript.Vita.Funding This publication is being made available as part of the report series written by the faculty, staff, and students of the Coastal and Oceanographic Program of the Department of Civil and Coastal Engineering.dc:date 1990dc:type Bookdc:format xvii, 201 leaves : ill. ; 29 cm.dc:identifier http://www.uflib.ufl.edu/ufdc/?b=UF00075474&v=0000124160971 (oclc)dc:source University of Floridadc:language Englishdc:rights All rights reserved, Board of Trustees of the University of Florida xml version 1.0 encoding UTF-8 REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EOO4BHO6A_JSVJUF INGEST_TIME 2017-11-28T19:59:15Z PACKAGE UF00075474_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES |