UFL/COEL/TR083
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. WenFa 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 WavePorous 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/H1) and the
best fit to the exponential decay function . . .. 54
5.3 The Measurements and the predictions vs. R. for L = 200.0 cm,
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.; nonbreaking 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/Hi) 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/Hi) 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/H1) 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/7H1) 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(101) . ... 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 Nonbreaking 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
rubblemound 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 porousseabed 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 modelDarcy's lawwas 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 parameterpermeability
coefficientand 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' ~
1012 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)
1n
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 NavierStokes
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 fluidsolid interaction, or the added mass effect.
The extension of DupuitForchheimer model leads to PolubarinovaKochina'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 fluidsolid 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 waveporousseabed interaction, Liu and Dalrym
ple (1984) added an acceleration term to the Dagan's model and it becomes
Vp = pq+pq 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 WavePorous 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 problemwithout laminar boundary layer correctionswas 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 1010m2, 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 rubblemound 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 reexamined 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 DupuitForchheimer'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, DupuitForchheimen'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 triangleelements 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
mid1970s. 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(')+ Aoek()] (z +h) (2.35)
cosh kh
et(x,z) = Boei*(+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 slopedface 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'eik'(z+1 ) csh k'( + h') + Bm ek'') '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 problemawkward behavior around "corner points". It has been found
that significant errors arise around such points even for the simplest casesinusoidal
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 shallowwater 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 rubblemound porous structures with multilayer 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 = accelerationrelated
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 xaxis coincides with the horizontal direction and the zaxis
points vertically upward. The spatially averaged equation of motion for the solid
skeleton in the xdirection 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 xdirection; 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 zdirection 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 Darcytype resistance law. The coefficient A is, at least, related to four founda
tional properties, one of the fluidthe dynamic viscosityand 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 secondorder 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 1n
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 sandsized 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 wellknown 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= qx, z)eit (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 1n
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 ~ 102Rf (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 flowrelated 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(101).
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,
106
106 104 102 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(101)1
Description Size Disch. Vel. Rf Rl Dominant
Range (m/sec) Force
Coarse sand < 2 mm. < 0(103) < 0(1) < 0(1) laminar
or finer
Pebble, laminar
or small 1 cm. 0(102) 0(102) 0(102) turbulence
gravel inertia
Large gravel 10 cm. 0(101) 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)=gat=a at z=0 (4.4)
,+.1) az=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(kt) (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 Darcytype 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 e2it (4.24
= (unp e2i+ 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) = aekizei(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 eio" = Dcosh k(h + h, + z) cos kxe'ieit'' (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 pe2' + 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 onewavelong 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 aek'" 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 prespecified 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 prespecified 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(m1) 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), DupuitForchheimer's model (D F:
p = 0 and Cf follows Eq.(4.46) and SollittCross'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
(LD): Liu and Dalrymple (1984)
3.0
 2.0.
R
N.
S.0
h
0.0 0
"., R (LD)
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 102, 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 pistontype 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) = In( ) (5.1)
Tj Hji
where j refers to jth 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
seabedwall 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 sec1 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 multivariate 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(s1)
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
t4
= 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, reconfirm
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(s1) 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(s1) 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(s1) aop(s) A,% im() p(s1)
30.0 5.4 x 107 4.8908 4.7638 2.60 0.0070 4.5 x 10
25.0 5.1 x 107 4.6312 4.4957 2.93 0.0090 5.6 x 107
20.0 4.8 x 107 4.2824 4.1428 3.26 0.0122 6.8 x 107
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 coarsegrained 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(s1) 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 nonLaplace 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 nonLaplace 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 tradeoff 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
(+) = ( )+ ( fiCi+) ,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) tanl()]}
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 reevaluated 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 eo
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 C0tI ri 8n
1 ,i+l* 1 ar,
= li~+ od = 0;
j T (e#40 ri an
= lim Inri d4
Tj+1 i T ei
= (12 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
subdomains, 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 ei'"
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 subdomains 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 shapesand
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 subdomain in this chapter. The
same formulation can be easily extended to multiporousdomain 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) = ei'(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 upwave
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 noflux 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 subdomain, the matrix equation, after applying the
noflux 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 + ...
+ (HiNu1 + 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.NlNiNYnN ~ N
KI,1 _1 KON+y+, ... K,,.,_,N K..
F Kl 2 n (7.30)
K~NN..14;N,. Ki,N,.4+N,,. ... K.,N1, 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.
