• TABLE OF CONTENTS
HIDE
 Front Cover
 Title Page
 Acknowledgement
 Table of Contents
 List of Figures
 List of Tables
 Abstract
 Introduction
 Literature review
 Porous flow model
 Gravity waves over finite porous...
 Laboratory experiment for porous...
 Boundary integral element...
 Numerical model for submerged porous...
 Laboratory experiments of a porous...
 Numerical model for berm break...
 Summary and conclusions
 Appendix A. Boundary integral...
 Appendix B. Experimental data for...
 Bibliography
 Biographical sketch






Group Title: Technical report – University of Florida. Coastal and Oceanographic Engineering Program ; 83
Title: Water wave interaction with porous structures of irregular cross sections
CITATION PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
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
Physical Description: xvii, 201 leaves : ill. ; 29 cm.
Language: English
Creator: Gu, Zhihao, 1957-
University of Florida -- Coastal and Oceanographic Engineering Dept
Publication Date: 1990
 Subjects
Subject: Coastal and Oceanographic Engineering thesis Ph. D
Dissertations, Academic -- Coastal and Oceanographic Engineering -- UF
Genre: bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )
 Notes
Thesis: Thesis (Ph. D.)--University of Florida, 1990.
Bibliography: Includes bibliographical references (leaves 197-199).
Statement of Responsibility: by Zhihao Gu.
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.
 Record Information
Bibliographic ID: UF00075474
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved, Board of Trustees of the University of Florida
Resource Identifier: oclc - 24160971

Table of Contents
    Front Cover
        Front Cover
    Title Page
        Title Page
    Acknowledgement
        Acknowledgement 1
        Acknowledgement 2
        Acknowledgement 3
    Table of Contents
        Table of Contents 1
        Table of Contents 2
        Table of Contents 3
    List of Figures
        List of Figures 1
        List of Figures 2
        List of Figures 3
        List of Figures 4
        List of Figures 5
        List of Figures 6
        List of Figures 7
    List of Tables
        List of Tables
    Abstract
        Abstract 1
        Abstract 2
    Introduction
        Page 1
        Page 2
        Page 3
        Page 4
    Literature review
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
    Porous flow model
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
    Gravity waves over finite porous sea bottoms
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
    Laboratory experiment for porous seabeds
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
    Boundary integral element method
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
    Numerical model for submerged porous breakwaters
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
    Laboratory experiments of a porous submerged breakwater
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
    Numerical model for berm breakwaters
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
    Summary and conclusions
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
    Appendix A. Boundary integral formulation for energy dissipation in porous media
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
    Appendix B. Experimental data for porous seabeds
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
        Page 173
        Page 174
        Page 175
        Page 176
        Page 177
        Page 178
        Page 179
        Page 180
        Page 181
        Page 182
        Page 183
        Page 184
        Page 185
        Page 186
        Page 187
        Page 188
        Page 189
        Page 190
        Page 191
        Page 192
        Page 193
        Page 194
        Page 195
        Page 196
    Bibliography
        Page 197
        Page 198
        Page 199
    Biographical sketch
        Page 201
        Page 203
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 T0rum, Head of Research in Nor-

wegian Hydrotechnical Laboratory, Trondheim, Norway, and Dr. Hans Dette in

LeichtweiB Institut ffir Wasserbau, Technische Universitit Braunschweig, West Ger-

many, 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 .................... ........ ii

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 POROUS FLOW 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 ....... 50

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 BIEM Formulations ........ ................. 131

9.2 Linearization ............................... 133

9.3 Numerical Results .......... ..... ............ 134

10 SUMMARY AND CONCLUSIONS ..................... 146

10.1 Summary .. ............... ............ .. 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/(a2/g), (b) Nondi-
mensional wave damping rate k/(o2/g) . . . 45

4.3 Maximum nondimensional damping rate (ki),na/(a2/g) and its
corresponding permeability parameter R as functions of nondi-
mensional water depth h. (a2/g) . . . .. 46

4.4 Standing wave case. (a) Nondimensional wave frequency o,/(L/g) ;
(b) Nondimensional wave damping rate ao/(L/g) . . 48

4.5 Solutions based on four porous flow models. (a) Nondimensional
wave frequency r,/(L/g), (b) Nondimensional wave damping
rate a,/(L/g). ... .. ... . .. .... 49

5.1 Experimental setup ......................... 51

5.2 Typical wave data: (a) Averaged nondimensional surface eleva-
tion (rt/H1), (b) Nondimensional wave heights (I/H-1) and the
best fit to the exponential decay function . . .. 54

5.3 The Measurements and the predictions vs. R. for L = 200.0 cm,
h, = 20.0cm, h = 25.0cm. (a) Wave frequency a,, (b) Wave
damping rate a. ........................... 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 a,,
(b) Wave damping rate ai....................... 61

5.5 Theoretical values by the model of Liu and Dalrymple vs. ex-
perimental data of Table 5.5, Table 5.6 and Table 5.7. (a) Wave
frequency ra, (b) Wave damping rate . . . ... 62

6.1 Auxiliary coordinate system . . . ...... 68









7.1 Computational domains ....................... 77

7.2 Flow chart of the numerical model for porous submerged break-
w aters . . . . . . . . 96

7.3 Wave envelopes for (a) 'Transparent' submerged breakwater; (b)
Impermeable step ......................... 98

7.4 Porous submerged breakwater: (a) Wave form and wave enve-
lope; (b) Envelopes of pressure and normal velocity. . 99

7.5 Wave field over submerged breakwaters: (a) Impermeable; (b)
Perm eable. ............................. 101

7.6 Transmission and reflection coefficients vs. stone size for differ-
ent 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 sub-
merged 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)
Kp vs. Kt; (b) Kp vs. Km ... .................. 117

8.5 The predicted and measured Kt and K, versus Hi/gT2. (a) Ktp
and Ktm; (b) Kp and Krm. ...................... 118

8.6 The predicted and measured Kt and K, versus Hi/gT2. (a) Kp,
and Kt,; (b) Kp and K,. .......... . .... 119

8.7 The predicted and measured Kt and K, versus Hi/gT2. (a) Kp
and Kt,; (b) K, and K,,. ..................... 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 Case of L = 200 cm, h = DW = 30 cm, h, = DS = 20 cm
and dso = DD = 0.72 cm. (a) Averaged nondimensional surface
elevation (rl/H1), (b) Nondimensional wave heights (H/H-i) 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 (r/H1), (b) Nondimensional wave heights (H/Hi) 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 dso = DD = 1.20 cm. (a) Averaged nondimensional surface
elevation (rl/H1), (b) Nondimensional wave heights (H/H1) 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 (l7/H1), (b) Nondimensional wave heights (H/ H) 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 (rl/Hi), (b) Nondimensional wave heights (H/-7i) 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 (rt/Hi), (b) Nondimensional wave heights (H/Hl) 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 (t7/Hi), (b) Nondimensional wave heights (H/Il) 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 (rt/Hi), (b) Nondimensional wave heights (H/-IH) 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 (rt/Hi), (b) Nondimensional wave heights (H/-i-) 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 (q7/Hi), (b) Nondimensional wave heights (H1/H1I) 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 (rl/Hi), (b) Nondimensional wave heights (H/-l) 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 (tr/Hi), (b) Nondimensional wave heights (H/IH1) 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 (r?/Hi), (b) Nondimensional wave heights (H/Ti1) 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 (a/lh) 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 (r/H1), (b) Nondimensional wave heights (H/H-i) 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 (r/H1), (b) Nondimensional wave heights (IH/Hi) 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 dso = DD = 1.20 cm. (a) Averaged nondimensional surface
elevation (n/H1), (b) Nondimensional wave heights (H/Hi) 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 (r/Hi), (b) Nondimensional wave heights (H/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 dso = DD = 2.09 cm. (a) Averaged nondimensional surface
elevation (q/Hi), (b) Nondimensional wave heights (H/-1-) 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/Hj), (b) Nondimensional wave heights (H/RH) 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 dso = DD = 3.74 cm. (a) Averaged nondimensional surface
elevation (7/H1), (b) Nondimensional wave heights (H/II) 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 d0o = DD = 1.48 cm. (a) Averaged nondimensional surface
elevation (?7/Hi), (b) Nondimensional wave heights (H/-1) 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 (il/H1), (b) Nondimensional wave heights (H/IH1) 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 (ir/Hi), (b) Nondimensional wave heights (H/H1) 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 d5o = DD = 1.48 cm. (a) Averaged nondimensional surface
elevation (fl/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 (rl/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 (l/H1i), (b) Nondimensional wave heights (7H/'i) 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 (i7/Hi), (b) Nondimensional wave heights (H/H-) 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 (r]/H1), (b) Nondimensional wave heights (H/H1) 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 d5o = DD = 2.09 cm. (a) Averaged nondimensional surface
elevation (rl/H1), (b) Nondimensional wave heights (H/H-) 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 (q7/H1), (b) Nondimensional wave heights (H/H) 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 dso = DD = 2.09 cm. (a) Averaged nondimensional surface
elevation (ir/H1), (b) Nondimensional wave heights (H/-i-) 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 dso = DD = 2.09 cm. (a) Averaged nondimensional surface
elevation (qr/Hi), (b) Nondimensional wave heights (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 dso = DD = 2.09 cm. (a) Averaged nondimensional surface
elevation (rl/Hi), (b) Nondimensional wave heights (I1H/71) 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 d5o = DD = 2.09 cm. (a) Averaged nondimensional surface
elevation (rl/H1), (b) Nondimensional wave heights (I/7H-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/q) 0(10-1) . ... 30

4.1 Comparison of The Predictions and The Measurements . 43

5.1 M 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,r and ap for dso = 1.48 cm, L = 200 cm. .... 59

5.7 Comparison of a, and ap for dso = 2.09 cm, h, = 20 cm. .... 59

5.8 Comparison of a, and a, for dso = 0.16 mm, 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 K, ..................... 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 im-

portant 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 con-

temporary 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 lin-

ear 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 struc-

ture and flow percolation inside the porous structure. The former one is conserva-

tive and the latter two are 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 disadvan-

tage 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 attenua-

tion 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 effi-

cient 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 ele-

ment 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 dissipa-

tion 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 break-

waters,







4

6. Verify the numerical model with laboratory experiments of a submerged rubble-

mound 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 deforma-

tion 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 microstructure 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 char-

acteristic properties of porous media are lumped into one parameter-permeability

coefficient-and the model has the form

Vp = q (2.1)
K,

where Vp is the gradient of the pore pressure, ip and q are, respectively, the dy-

namic 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-' ~

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 dom-

inant, 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 pres-

sure 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 f )q" (2.2)


where a = AI/K,, and b is coefficient for turbulent resistance defined by Ward (1964)

as
b = pC! (2.3)
57K







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)Sv
a = ao nd (2.4)


1-n
b = bo (2.5)
n2d,

where ao and bo 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 = P+ (2.6)
K, 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 = aq+ bq2 + c (2.7)
at







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 = u+f pn (2.8)
k at

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 q by
q'= ni. 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 + + b) 4+ (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

+ pC+ 1- n ai
Vp = (n- + n U) + p(1 + Ci ) (2.11)
Kp TKp n at

where u is again the spatially averaged microscopic velocity vector and Ca is the
virtual mass coefficient.

The value of the virtual mass coefficient Ca 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 Ca by laboratory experiment. They measured the instantaneous velocity and

the pressure gradient in the experiment, and calibrated the data with their semi-

theoretical porous flow model:

Vp = (a + b q)q+ p(l + Co) (2.12)

The results so obtained for Ca 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 = q- V2q (2.13)
Kp Kp
where is a constant coefficient depending only on the geometry of the media and
it was defined as
dL (2.14)
80
In applying this model to the wave-porous-seabed interaction, Liu and Dalrym-
ple (1984) added an acceleration term to the Dagan's model and it becomes

Vp = -pq+p-q V' (2.15)
K,, nat 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 Ca is set to be zero.

Based on the phenomenological approach, Barends (1986) added another porous
flow model to the list:

V = ( + g (2.16)
n at K(q)







10
and K(q) was defined by Hannoura and Barends (1981, see Barends, 1986) as

/ pd (2.17)
K(q) = g (2.17)

with
Ca(1 n)
~ n5
where C is the drag coefficient related to the porous Reynolds number (R =q q

d,/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 pa-
per, 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 kh (2.18)
2(aK,/v)k,
ki = 2(apv)k,. (2.19)
2krh + sinh 2kh
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 perme-

able 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

kr = k + 2k (2.20)
k, = ko 2koh + sinh 2koh


2ko aK, K2
ki = 2k o (- + ko (2.21)
2koh + sinh 2koh v 2(

with
a2 = 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

2k a + 2koN
2koh + sinh 2koh

where k = ko + ki with ko being the initial wave number for the impervious bottom

as defined in Eq.(2.22) and
N= naKp
N =

For fine sand, N -+ 0 and
i/-v
k, = __k_ V a (2.24)
2kh + sinh 2krh
and the temporal damping for known k in such case is

a = -i-= kh- (2.25)
v2 = vsinh 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, )
ki- =2k (aKp/v + k, ) (2.26)
2k,h + sinh 2k,h FT2a

with k, remaining the same as that in Eq.(2.18). It is noted that this is the same ex-

pression 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, [1 1- h/g g 2 (2.27)
f T 4v h
(1 u2h/g) aK,
ki = ( ) (2.28)
2h V

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
( + i)(1 g tanhkh) + R tanhkh, tanhkh(1- s- ctnhkh) = 0 (2.29)
n a2 a2









or equivalently

2 gk tanh kh = i .tanh kh,(gk a2 tanh kh) (2.30)
R n
where h, is the seabed thickness and R is the permeability parameter defined as
uK
R= K

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-10m2, 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. Signif-

icant 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 theoreti-

cal 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 transmis-

sion 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 al-

though 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

Kt = b (2.31)
1+ (S if + n2)

S if n2
K, = (2.32)
S- if + n2 i2n-
gh
where S is the coefficient for the inertial resistance of porous media, i is the imag-

inary 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 Ka Kla 16b 1
f = [-(1- ) + (1 + )2 + bai (2.33)
K,1 2a 2 37r h

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 rockfills. 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:

C' = (01 + 2XP + 03y)t + 34 + x+ + 06y (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

,(, z) = [e(-')+ Aoe-k(-)] (z +h) (2.35)
cosh kh


et(x,z) = Boe-i*(+I,)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 (2.37)

Kt =| Bo 1 (2.38)

The agreement of the transmission coefficient for impermeable floating struc-

tures with the experiment data was fairly good, but the comparison of Kt and Kr
for a sloped-face permeable breakwater was not as satisfactory.
Finnigan and Yamamoto (1979) further added the modulated standing wave

modes to 4, and Ot, while keeping the constant element and Darcy's law unchanged

in their model. The two potential functions were, respectively,

r_(, z) = [(-) co (z + h) + [ A+.(.-,)cos km,( + h)
cosh kh + A cos kmh
(2.39)









-(x, z) = B'e-ik'(z+1 ) csh k'( + h') + Bm e-k'-') 'cosk k h (2.40)
cosh k'h' m1 cos kmh'
with
gk tanh kh = -gkm tan kmh = a2 (2.41)

gk' tanh k'' = -gk' tan kmh' = 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 condi-
tions. 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 ex-

periments and empirical formulae. For impermeable submerged breakwaters, quite

a few empirical formulae for transmission coefficient have been established from ex-

perimental data by different authors (see Baba, 1986). For example, the equation

given by Goda (1969) for the transmission coefficient due to overtopping is

Kt = 0.5[1 sin ( F) (2.43)
2 a 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):

F F
K,= C(1 ) -(1- 2C)F (2.44)

where R is the wave run up given by Franzius (1965, cited in Baba, 1986)

R = HIC,(0.123 I)(c2v d+c3) (2.45)

with Cx = 1.997, C2 = 0.498 and C3 = -0.185. d is the total water depth and

0.11B
C = 0.51- 1B (2.46)

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 exper-

iments.

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 = 2fT () dt (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 break-
water 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:

F,+ FD + FI + Fi F, = ma (3.1)

where m = mass of the solid; a = acceleration of the body; F, = pressure force;

FD = velocity related force also known as the drag force; FI = 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 two-

dimensional 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

dz(1 nAz)dz + FD., + FPI + Fbaz, fz = P.(1 n) t~, dxdz (3.2)
where = pressure gradient in the direction; A = averaged area porosity
where ap- = pressure gradient in the x-direction; nA = averaged area porosity
axz










aw ap
S+ dz p + -d
Oz Oz


U

p


Tu
u+ -dx

Op
p+ dz


Figure 3.1: Definition sketch for the porous flow model

defined as the ratio of Avoid to Atotal with Avoid being the area of the voids and

Atota being the total control area; n = volumetric porosity defined as the ratio of

Vvoid to Vtotai; Fb, 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
a- -d nA dz FD, Fi, + Fsy, = pn dz dz (ty, t,,)
ax


with p being the density of the fluid and Fb/e 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 I = uo Udv (3.4)
"void = o old


(3.3)







22

where u, = 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 6 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 = ua dv = o uadv = n uf (3.5)
V fv V .,oid

here V is the total volume Vtota.

This velocity also known as the discharge velocity, is related to the actual dis-
charge 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

Dq 8q aq
+ = Vq n" *tf (3.6)
Dt at 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

FDz = [Ax(qz nu,) + B. q nu, [ (q= nu.)]dx 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:
FI, = Cz( ni,,)dx nA 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
=- A(qf nu,,) + Bze q nu, I (q, nu,,)
ax
+ (+ Cz)(4 -n i.) (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
BP p
S= Aqz + B qz I q9 I + ( + C.)4. (3.10)
ax n
The force balance in the z-direction leads to

aP p
Arq + B, I q I q + (- + C,), (3.11)
8z n

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 founda-

tional 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, p, is factored

out, we have the empirical law:
A = (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:
(1 n)3 it
A = ao (3.13)

where d, is a characteristic particle size of the pore material. The coefficient ao
obviously still contains the other properties such as tortuosity and connectivity.

Since both of these factors have directional preference, ao 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 connectiv-

ity (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. = ], -,CDAp I uf I Us, (3.15)
2

where Ap is the projected area of individual particles; CD is the form drag coefficient

of the individual particles; 6 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:

FDNZ = 6CD I q dzdz (3.16)
2 n2d,

the constant of proportionality is absorbed in 6. Comparing Eq.(3.16) with the

second part of Eq.(3.7), we obtain

p 1-n
B, = CD --- (3.17)
2 nAzn2d,

If we let:
-CD
bo = 2 (3.18)

Eq.(3.17) can be expressed as

B = pbo (3.19)
P Azb 2d,

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 bo 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
FI, = pC.Vf (3.21)

where Ca is the added mass coefficient of each individual particle and Vp is the

volume of each particle. One may then readily obtain, by comparing Eq.(3.21) with

Eq.(3.8),

C, = p -- (3.22)
nAz n
Since Cz is related to the volume of the solid skeleton, one does not expect

significant directional preference, unless the geometry of the element deviates sig-

nificantly from a sphere shape. A new coefficient, Cm, can be defined such that

Sp n +C(1 n) 2
Cm, = p + Cz = p[ + C -- (3.23)
nAz nAz '

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+ | q 1 q+ 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 u,.

We now proceed to assume that the porous material is isotropic and that all

the coefficients reduce to scalar quantities and


nAz = Az = n

Recall that the convective acceleration in q has been assumed, in Eq.(3.6), to be

negligible, i.e.

Sq=
q= -
Bt







27
Substituting these conditions into Eq.(3.24), we obtain

(A + B q I)q+ Cm = -Vp (3.25)

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= q-x, z)e-it (3.26)

where q'is 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 iaC, + B I q 1) (3.27)

A new set of nondimensional parameters are defined as follows:

Permeability parameter, R,

R = -g, (3.28)

with
pr n2d2
Kp = (3.29)
A a- o (1- n)3
Inertia parameter, /
n + C(1 ) (3.30)
n2 (3.30)
Volumetric averaged drag coefficient, Cd

B 1-n
Cd = B = bo1n (3.31)
p nsd,







28

Substituting these parameters into the above equation reduces it to the following

familiar form:
1 Cd
Vp = pa( i3 + I q I) (3.32)
R C

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 fi n + C(1 n)33)
Laminar Resistance fG ao(1 n)3


Turbulent Resistance f, bo
Laminar Resistance f I aon(1 n)2 (3.34)


Turbulent Resistance f,, bo(1 n) R(33)
Inertial Resistance f Can(1 n + n) R
Ca

where R, and R/ are two forms of Reynolds number defined as

R q Id (3.36)


and
Ri I 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 R/, 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 R, 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

Co ~ 0(1)

ao ~ 0(103)

bo ~ 0(1)

Therefore, Eqs.(3.33) through (3.35) reduce to

t 10-'R, (3.38)

-j I ~ 10-2Rf (3.39)

LfI R- (3.40)
fi RI
When Eqs.(3.38) to (3.40) are plotted on a Rf 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 com-

monly encounter wind waves with frequency in the order of O(1)(rad/sec) and the
nondimensional bottom pressure gradient defined as V(p/y) in the order of O(10-1).











106 N Kegion

fn > 0lfi t, f.
10o I > 10f


102 ,/.n f, A


1
L Region
102 / I Region

fA > lofn A f > loft
104 f, > 10f,


10-6
10-6 10-4 10-2 1 10 104 106

Figure 3.2: Regions with different dominant resistant forces

Under such conditions, the importance of the three resistance components for bot-

tom 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)1


Description Size Disch. Vel. Rf Rl Dominant
Range (m/sec) Force
Coarse sand < 2 mm. < 0(10-3) < 0(1) < 0(1) laminar
or finer
Pebble, laminar
or small 1 cm. 0(10-2) 0(102) 0(102) turbulence
gravel inertia
Large gravel 10 cm. 0(10-1) 0(104) 0(104) turbulence
crusted stone inertia
Boulder 0.3 1.0 m 0(100) 0( 0(10 O(105) turbulence
crusted stone inertia
Artificial turbulence
blocks, > 1.0 m > 0(100) > 0(106) > 0(106) inertia
large rocks_____














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 vari-
ables 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 Dal-

rymple, 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
V2 = h < z < 0 (4.1)

In the porous medium domain, the linearization of the nonlinear pore pressure

equation, Eq.(3.32), yields

VP, = pafoq (h + h,) < z < -h (4.2)

where fo is the linearized resistance coefficient.













{z

V A / xe,








the discharge velocity, V = 0, leads to


VP =0 0 + 4.3
Ia h


TheP b dr cosditoseaned ti'e4 1 are N o
I0 0 0 0 0 0 0 0 0 0 o o 0 0 0 0
















g at
+C -000 0 at z = 0h (4.)
'o''o'o'o' o 4 o o 0 0 Oo 0 0 :




Figure 4.1: Definition Sketch

If the porous medium is rigid, thus, incompressible, the continuity equation for
the discharge velocity, V q = 0, leads to

VP, = 0 (h + h,) = -h (4.3)

The boundary conditions of entire system are
1 a ei(k._t)
,T7x,t)=ga-t=a at z=0 (4.4)


,+.1) a-z=O at z=0 (4.5)


a9
p--=P, at z =-h (4.6)


az 1 aP.
-az p=- a zat z= -h (4.7)
aZ pu7o az







33
0P,
ap =0 at z = -(h + h,) (4.8)
8z

The potential and pore pressure functions are assumed to have the forms

S(z,z,t) = [Acoshk(h + z) +Bsinhk(h + z)lei(k-'t) (4.9)

P,(x,z,t) = Dcoshk(h+ h, + z)e'(k-") (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- cosh kh, (4.11)
pa
D
B = sinh kh, (4.12)
pefo

is then given by

(X, z, t) = [icosh kh, coshk(h + z) + 1 sinhkh, sinh k(h+ z)]e(k-t) (4.13)
pa fo

D can be obtained by the free surface boundary condition given by Eq.(4.4):

pga
D (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 equa-
tion:
a2 gk tanh kh = --tanh kh,(gk a2 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 CI # 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 laminarr skin

friction dominates, fo -* ), Eq.(4.15) reads

a2 gk tanh kh = -iR(gk oa tanh 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 + h,, 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
a2 = -i tanh kh, (4.18)
fo
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

o2 = -iRgk tanh kh, (4.19)


which states the Darcy-type resistance law. On the other hand, when R -* oo, 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

a2 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)ni (4.22)

where ED is the energy dissipation in a controlled volume during one wave period

and the subscripts I and nl 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= = ds d (4.23)

where ED is a complex energy dissipation function and .7 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 e-2it (4.24
= (unp e-2i+ 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

r (x, t) = a e'' cos kze-'"'' o, < 0 (4.25)

where ai 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 r becomes

ti(x,t) = ae-kizei(k'-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, = p e-io" = Dcosh k(h + h, + z) cos kxe'ie-it'' (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,), z = L and z = -h with L being

the wave length.

As x = 0 and z = 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
1 t+T f L
ED = j (wo pe-2' + wo p') dx dt (4.28)


with
kD
Wo = u, I=-h= -- sinh kh, cos kz e"'t (4.29)
paf









and

f = fo for linearized system (4.30)
1 Cd
f = --i+ IwoI
R a
= fi + f2 IWo for nonlinear system (4.31)

with

f = -- iP a = C (4.32)
R a
In Eq.(4.31), the magnitude of wo is approximately taken as

o = Dfo) sinhkh, (4.33)
afo

with

D(fo) = (4.34)
cosh kh cosh kh,(1 tanh kh tanh kh,)

where the decaying wave amplitude aec'i has been approximately replaced by u,
which is the averaged wave amplitude in the period of [0,T].
With such approximations, f for the nonlinear system is no longer a function
of x and t. The approximations made above are not expected to cause significant
errors for the final results.
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

kLT oD'(f) t + D(f) I
ED = sinh2kh, T(t) + D( 2(t)] (4.35)
8p h f f

where Ti(t) and T (t) are the nondimensional functions of t generated by the inte-
gration 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 T(f) T(t) + I (t) =
f1 fo

D( + I+f2 I wo 1)(T 12 T
fi+ f2 wo) D(f + f2 wo) T2(t) (4.36)
l+/2(f Iwo I fl+f2 Iwo I
Since this equation has to be satisfied at any time instant, we obtain the follow-
ing equations:
D2(fo) D (fi + f2 ) (4 37)
(4.37)
fo f + f2 IwoI

I D(fo) I2 I D(fi+ f2 Iwol)12 (4.38)
fo f + f2 I I
Therefore,
fo = f + f2 Iwo (4.39)

Substituting the expressions for fl, f2 and wo in the above equation, it becomes
1 + Cd I kD(f) sinh kh (4.40)
R a |fol
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
ED = f ds dt
Jt s


= -Jt (- dz+ Fn dz + Lf 7nd)dt (4.41)
St ==0 =L =-h
By substituting the expression of 7, derived from Eq.(4.10) into the above equa-
tion, 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 k, j< 1, which is generally true in coastal waters

(see the results for the progressive wave case), then the energy dissipation becomes
= t+T ft
ED = t 7 dx dt (4.42)
t O
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-k'" 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)i = -/ (4.43)

and
1 Cd I kD(k, fo) I(4.44)
(fo)r = + (4.44)
R ao | fo
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 P). It is often more
convenient to specify the permeability parameter R = aK,/v (Rr = 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):

Cd = a- = (4.45)

with
C bo(1 n) (4.46)
bao(l n)

Now Cf is a nondimensional turbulent coefficient. To be consistent with the
suggested values of ao and bo given in Eqs.(3.14) and (3.20), Cf should be in the
range of 0.3 ~ 1.1 for a porosity of n = 0.4. Equation (4.40) now becomes:
1 + C I kD(k, fo) sinh kh I
fo = + f(4.47)
R |Ra afo







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 = a, + iai is complex, k is real and known. Equation
(4.15) can be rewritten as


gk(tanh kh tanh kh,)
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 [ Q(f(S)) + Q?(fO()) + Qr(f(n"))] (4.49)

n+1 -1 [Q(f) + Q(fn)) Qr(f") ]- (4.50)

1 C1 k Ir(o ), (4f51)
fo = i +((n) (4.51)

f ) -1 i (4.52)
R

where the superscript n denotes the level of iteration, and the criterion of conver-
gence is
fn) O-) and o'(n) I<_ (4.53)
I fo(n) I_ and a ((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

F(k, fo) = a2 gk tanh kh + tanh kh,(gk a2 tanh kh) = 0 (4.54)
fo
and the iteration was carried out as follows:


(4.55)


Fr(k"+), f")) = 0







41
F,(k(n+1),f1")) = 0 (4.56)

f n) 1 1 Cf I kD(k"(),fo) (457)
R VrV l \afo") I

fol = 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
fo") fc"-1) k ,(") f("-)
fI n)0 ) f ( and I k(-) -| (4.59)

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, bo = 2.0 and Ca = 0.46 and the average wave height H in
the table was calculated according to
Ho HLg
H = (1 i
HoIn 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
k Iko,
A% = i m | x100% (4.60)
kim

with ki being either kip 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 kI (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 mm to 2.3 m.

The thick dash 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) monoton-

ically with increasing R, from one limiting value kr, corresponding to the case of

an impervious bottom at depth h, to the other limiting value kr2, 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.) kim(m-1) k,(m-) ,% kiLDm) ALDo%
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


kim the experimental damping rate given by Savage (1953);
kip the theoretical values by the present model;
kiLD 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: 3 = 0, Cj

= 0), Dagan's model (DAGAN: 3 = ., Cj = 0), Dupuit-Forchheimer's model (D F:
p = 0 and Cf follows Eq.(4.46) and Sollitt-Cross's model (S C: f follows Eq.(3.30)

with Cf following Eq.(4.46) and n = 0.4, ao = 570, bo = 3.0 and Co = 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.0m Liu and Dalrymple H = 1.0m
= 2.0 m T = 4.0sec
= 4.0m h, = 5.0m
= 6.0 m


SI I
-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,/(o2/g), (b)
Nondimensional wave damping rate ki/(a'/g).


0.23


4I'


1.8 -


1.4 -


232.00


-2.0
-2.0


0.23


0.10



0.08



0.06



0.01



0.02



0.00


.. .. .. .. .. .. .. .. .. .. ... ...
. . .
. . - - - - - -


I .










46






4.0


(L-D): Liu and Dalrymple (1984)

3.0




- 2.0.
R


















N.
S.0
h



0.0 0
"., R (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'/g)


Figure 4.3: Maximum nondimensional damping rate (ki).,=/(a2/g) and its corre-
sponding permeability parameter R as functions of nondimensional water depth
h (a'Cg).







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


-2.0 -1.0 0.0 1.0 2.0
PERMERBILITY PARAMETER LOG(R)
(a)


2.28


22.69


-1.0 0.0 1.0 2.0 3.0
PERMEABILITT PARAMETER LOG(R)


(b)

Figure 4.4: Standing wave case. (a) Nondimensional wave frequency ar/(L/g)2; (b)
Nondimensional wave damping rate a/i(L/g)2.


0.23


225.44


h = 4.0m Liu and Dalrymple H = 1.0m
h = 2.0m L = 20.0m
........... h = 4.0 m h, = 5.0m
----- = 6.0 m


--. ----- -----..-----....................


1.60

1.50

1 .40

1 .30

1.20

1. 10


1.00

0.90







0.06


0.05


0.04


0.03


0.02


0.23


225.44


0.01

0.00 -

-2.0


.


.














PARTICLE SIZE (CM)
2.28 22.69


S5 C
............. DAGAN
DARCT
.. 0 F


H = 1.0 m

L = 20.0m

h, = 5.0m
/-" ---------------------


/
/ h=4m
/
....----------..


-1.0 0.0

PERMEABILITY


2.28


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), (b) Nondimensional wave damping rate a,/(L/g).


0.23
1 5


.i


225.44


1.4 -



I 2
1.3 -



1.2 -


-2.0


0.23


225.44


0.12


0. 10


0.08


0.06


0.04


0.02


0.00


S5 C
............. DORGAN


/ \


/


S \ h=4m


I t


-





. J













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 H-
(oi) = I-n( ) (5.1)
Tj 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

ai = In[ + e2aT e2a"'T (5.2)
2T-





















wave
h maker




L/2 L/2


Figure 5.1: Experimental setup

Table 5.1: Material Information

dso(cm) 0.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 ac is the actual seabed damping rate, ai, 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 d5o = 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

dso (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 ai, = -0.01 sec-1 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 devi-

ation 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 ai obtained from the data ensemble by the least square analysis.

Based on the measured or 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) ,(-1) ') (s-) a(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
0.40 iL = 200. CM
> I T = 1.27 S5
0.20


UL 0.00


o -0.20
t-4
= -0.40
C
z
-0.60
0.0 4.0 8.0 12.0 16.0 20.0
TIME (SEC)




1.20
DO= 30.0CM
1.00 S.. 0= 15.0CM
'I 00= 1.48CM
LI-
= 0.80 .


: 0.60 '.. .
C3
N4
S0.40 ... .. ...


S0.20 .


0.00
0.0 4.0 6.0 12.0 16.0 20.0
TIME (SEC)



Figure 5.2: Typical wave data: (a) Averaged nondimensional surface elevation
(r7/Hi), (b) Nondimensional wave heights (7'/H-) and the best fit to the expo-
nential decay function.









that the error function defined as

= 1 p)2 + (- iP)2] (5.3)
M EI(
1=1 arm aim

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 bo, to an extent, re-confirm

those given by Engelund (1953) and the others.

5.3 Relative Importance of The Resistances in The Experiment

Introducing ao, bo, Ca, and the averaged wave height H into Eq.(4.29) for I qI

and using Eqs.(3.33) through (3.35), the relative importance of the various resis-

tances 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 theoret-

ical values of ai and a,. The relative error, defined as

A% =1 m PI x100% (5.4)
am



















Table 5.4: Comparison of The Resistances. h, = 20 cm, L = 200 cm.

dso (cm) |q (cm/s) \aI Ry R f,/f, If/fi fn,/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,o (cm) R rm(1) (1) A% sm() iap(s1) AjL
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 2.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 2.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 4.4646 4.3805 1.88 -0.1173 -0.1076 8.24



is uniformly very small for a, (wave frequency). For ai (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 ar 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(R)



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 ar, (b) Wave damping rate oj.


WAVE LENGTH = 200 CM
BED THICKNESS = 20 CM
WATER DEPTH = 25 CM


S + 0
o



o PREDICTED

x MEASURED


- I












Table 5.6: Comparison of a, and a, for dso = 1.48 cm, L = 200 cm.

h (cm) H (cm) crm(.s-) (r,.s- ) Ar% ai(s-1) c,)(s-) A1%
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 a, and ap for dso = 2.09 cm, h, = 20 cm.

h (cm) H (cm) arm(s- ) p(- ) Ar% cim(s) ap(s-1) i%
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 a7'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 ,p for dso = 0.16 mm, h, = 20 cm, L = 200 cm.

h (cm) R ar(s-1) aop(s-) A,% im(-) p(s1)
30.0 5.4 x 10-7 4.8908 4.7638 2.60 -0.0070 -4.5 x 10-
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 move-
ment 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. How-

ever, 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, dso = 1.48 cm, h =20 cm and L =200 cm,




































2.0 3.0 4.0 5.0

EXPERIMENTAL FREQUENCY


0.25


0.20


0. 15

*Ag
."
? **
0. 10 ..


0.05 "


0.00 *

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 ao, (b) Wave damping rate a,.
























5.0




4.0




3.0


I- I I I 1 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


6.0


0.25


EXPERIMENTAL CAMPING 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 o,, (b) Wave damping
rate a,.


.i;Z
'j
.+
...
Q
.;
.i
.~


.*.i


.4 .. .


,










Table 5.9: Comparison of cra and ap for h, = d5o


d50 (cm) H (cm) Grmr(-1) rp (-1) A,% im(c(s-) rp(s-1) A,%
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. Consid-
erable 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 differen-
tiable 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

f D(V f u)dA= 6. ds (6.1)

where V is the vector operator, i is any differentiable vector and n' is the unit
outward vector normal to C. If we define that v = UVV and then v7 = VVU in
Eq.(6.1), two integral equations can be generated in terms of U and V,


J(VU. VV + UV2V)dA = fUVV ids (6.2)







65

(VV VU + VV2U)dA= VVU* ids (6.3)
Subtracting (6.3) from (6.2), we have

D(UVV VV2U)dA = (UVV VV) i ds (6.4)

Introducing the expressions
VV au
VV n= = VU -n= n (6.5)
on an

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 )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

G(P a)(Q)
C[(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) = nr (6.9)

where r is the distance between point P and point Q and it can be expressed by

r = V(zp Xq)2 + (zp zq)2

on the x z plane.
Substituting Eq.(6.9) into Eq.(6.8), it becomes

#(Q) Br 8(Q)-
[ In r n( 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

(P) = In r ]ds (6.11)
cz i(P) T(Q a
r an an

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
a
boundary values of D and are solved, P has to be kept on the boundary C in
an
the process of computation. When 4 and 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 inte-
gration into N parts accordingly, Eq.(6.11) reads
N 4D ar, a1 N
a = Inrn)ds = Eli (6.12)
=1 a, n an=1
The curved segment Cj is usually replaced by a straight line to simplify the inte-

gration. 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 deriva-

tives 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.

= -N(e)D- +Ng(C)4y+1 (6.13)
-- (C) = N~ ()( )i + N ()() )j+ (6.14)

Cj C ej+i

where N1 and N2 are linear functions with the feature of

N (C) = 1 N1(e+1) =

Ng2() = 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 with-

out 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 gen-

erate 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.


















/
/
-'ii
/
/


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 sys-

tem, as shown in Fig. 6.1, has to be established to facilitate line integration over

the element. The two axis C r7 are perpendicular to each other with one lying on

the element and pointing in the direction of integration, from Pj to Pjy+, and the
other one being in the same direction as the outward normal vector in. Based on this

auxiliary coordinate system, the integration over each segment can be completed

analytically and be expressed in terms of 4., 4i+I and t7i, the distances from Pj,

P,+i and Pi to the local origin defined in Fig. 6.1. The values of j, (e+i and i7, are
determined by the global coordinate information of three points.

Since

+ 17+ 2 = rI-

(i+ Ai)X +171 = ~tij+i

where

r,+, = vi x+) + (z. z+)

^ l= \{i- Za+l)2 + (z.- *'~)








then
'- Yj (6.15)
2Aj

+l = ,i + A (6.16)

The value of ri is therefore

mI= ,/'- (6.17)

However, the numerical test showed that Eq.(6.17) could lead to significant error

of I ri due to runoff errors. A more accurate expression for I ri j can be obtained by
computing the distance from Pi to the straight line PyP,+i directly from the global
coordinates, i.e.
I A(x- x) + zi (6.18)

where
A = jz+1 z)
zy+l zy
The sign of ri depends on the relative position of Pi to the boundary line element

PjPj+1. If Pi is in the same side of the element with ii, r7i is positive, otherwise, it
is negative. Mathematically, it can be expressed by the following two equations for
vertical and nonvertical elements:

sign (7) = (S zi) Az =x (6.19)
(xj x,) Az

sign (ri7) = A Azi X,) A x+ (6.20)
Ax Azio

where

Ax = xY+ Xi

Az = zy+ z
Az
Azio = (X- r) + z z

Therefore


(6.21)


rni = sign (r1i) I 7, 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,+l. For the velocity potential function 4 and its normal
derivative we have, in terms of q r] coordinates,
an
(+) = ( )+ ( f-iCi+) ,i i+1 (6.22)
Cj+1-


a = ([)i] + [C.a+l)+ e
n a a a (6.23)
j e <++1

It is obvious that
a(D) = ; (-)(a) = )(-



) = +1; ) )= (an)

Substituting Eqs.(6.22) and (6.23) into Eq.(6.12), and noting that

r, = 7, + (6.24)
ri a= r (6.25)
an a 7i r,
the integration over segment j (between P, and Pj+1) can be carried out analytically
in the local coordinate system (Ligget and Liu, 1983):
f fi+ 9ri a j +r
l (r In r, )dan

= Il 4+ H +1,-K i( ), +- K (i+ (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 -i = -ia1 + j+iy12j
cj 'ij ti






71
KI. .= I 1 2lI

Kl = -,j + 1j+4

K, = I2j- lSi2


1 1f i+1 =r
i+1 fi fi ri Kn


- ij7 +1I
2AE 171?+


= ~1
j1


fJ+


1ar dEr = r fi+l
ri 8n 2Ai


1 [t an-'( 1
A ji


(6.32)


- tan-'( )]
ri


= 1 fi+nrid.
fi+1 j 1i

= -{(, + J+1)[ln(o77 + -+1)- 1]


-(7,? + )[1ln(,? + ?) 1]}


1=In ri d(
*+ n( + -+


1 2
= {+I i(2 + 77 ) ln(7? + 2(J) 2+ )


+2rjj[tan-'l(e'+l) -tan--l()]}
th 77-


with


(6.28)

(6.29)

(6.30)


rii
2A,


I CJ+ i


;dC
r.


(6.31)


d xi
rT


(6.33)


(6.34)







72

for all i j and i j + 1 and

Ai = Ei+1 ti

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 -
Ij 1= lim
j+1 1i -
:i+: C- 'o


1
-= lim
ji+1 i e-o


Si+i1 =
Ji+1 1 ar 0

ei+, ri an

i+, In ri d
fi+
*(,+<


= (2nA 1)
4


= 1+ im t In r, de
$y~i-,<-

= lnA- 1


Similarly, when i = j + 1,


1 /f++i- 1 ari
= lim 5l dC = 0
S+1i C-0tI ri 8n

1 ,i+l-* 1 ar,
-= li~+ -od = 0;
j T (e--#40 ri an


= lim Inri d4
Tj+1 -i T ei


= (1-2 In A,)
4


(6.35)



(6.36)


(6.37)


I1+1.1



j13+11J

Ij111


(6.38)


(6.39)



(6.40)


(6.41)













= nA 1


n ri d


(6.42)


Summarizing all the Ij's according to Eq.(6.12) for node points i = 1 N, and

assuming that


(6.43)

(6.44)


a0
at all the nodes, a system of linear equations with unknowns of P and n can be

obtained:


N
Si= E Qin
j=1


P4,1 = H~,_N + Hli 6,1aj
?j., = Hijl + H ,i
^**2 HLi b ijci


1 ifi=j
.-- 0 if i # j


Sij = K, N + Kj,1

LiJ = Ki_ i- + Kaa

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)


(an







74

where 4 and 4, 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 J = (6.49)

II = (6.50)
an an
III a- k= (6.51)
on

where D, f and k are known functions on the boundary.
In
In water wave problems, the first relation usually means that the pressure distri-

bution is known, and the second describes a given flux distribution on the boundary.
They are also referred as 'essential' and 'natural' conditions (Brebbia,1989), respec-
tively. 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 the-

ory. 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, ('-)~ may not be equal to ( )+. 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 or D for some part of the boundary and for the rest of it, depending
an an
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 4, is the Laplace equation,

V2 = 0 (7.1)


with fluid velocities being defined as

S= (7.2)
ax


w = (7.3)
az

While in the porous medium domain, the viscosity of the fluid cannot be ig-

nored since the flow is largely within the low Reynolds number region. The flow is


































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, f = iu + kw with iand j being the
unit vector in x and z directions.

Taking curl of Eq.(7.4),

pafoV x 4= -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 $, in the porous domain exists such that
i I(
Io
!0






























potential function , in the porous domain exists such that


T= V (,


(7.7)









According to Eq.(7.4) and Eq.(7.7)
P
=-po- (7.8)
pefo

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 = (7.10)

In this study, the waves are assumed to follow the linear wave theory

S= Oeiot (7.11)

P = peiot (7.12)

It is to be noted that the time function is now expressed as eiat, as opposed to e-i'"

used in the seabed problem presented in previous chapters.

7.2 Boundary Conditions

The whole computation domain has four types of boundaries, free surface, im-

permeable 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)
8z 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:

S= 0 (7.14)
on

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 do-

main 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
S = P (7.15)

or equivalently
iap' = p (7.16)

with the linear wave assumption, and

a 1 p (7.17)
QnI pafo anri

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 struc-
ture. 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

r (xz) = e'k(z+')R(z) (7.18)

t(x, z) = e-i'(L-')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:

,;(z, z) = e-'k(+')A(z) (7.20)

with
Hg cosh k(z + h) (7.21)
2a cosh kh
Therefore,

S= + 4' at x= -1 (7.22)

S= 4t at x= (7.23)

where I and 1' 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 = I',
(g 8 -ik', = -ik'4 (7.24)
an aZ az

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,

r, = 4 4I (7.26)

a a a, a,._ ik, ik,
an az ax ax


= ik ik(4O 4i) = 2ik, ik4'

or
a 2ikd; 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 pres-

sure, is
ap 0 (7.29)
an
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 dif-
ferent 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 equa-
tion. On the free surface, the normal derivative are expressed in terms of 4 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- = 0+ for all the nodes, as:

aoji = (H,1 + HN)41 + (H i, + H, 2)02 + (H?2 + H3s)s...

+ (HN,-1 + HN4) + ... + (Hi=_It + H~fN)4 + ...

+ (HiNu-1 + HN.3)3 + + + (H,-~ + H.1)Nj
+ (H2..- + H~..)N + ... + (H -i + Ht`N
,+1 K,- K3 3- -

Ki11 KN1 + ...- K,~N24,lN~ -K -
S, Kfl4nNs i,Ns OnNs Ii.N-lNiNYnN ~ N-
KI,-1 _1 KON+y+, ... K,,.,_-,N K-..
F- Kl 2 n- (7.30)
K~NN..-14;N,.- Ki,N,.4+N,,. ... K.,N-1, KN; (7.30)

i = 1, 2, ......, N

where a, is the inner angle of node i and the subscripts 1, Ni, N2, N3, Ne and N,,
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 0-
is not necessarily equal to 0+ for all the nodes.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs