UFL/COELTR/073
NUMERICAL MODELING OF
CURRENT AND WAVE INTERACTIONS
OF AN INLETBEACH SYSTEM
BY
YIXIN YAN
1987
NUMERICAL MODELING OF CURRENT AND WAVE INTERACTIONS OF AN
INLETBEACH SYSTEM
By
YIXIN YAN
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
1987
ACKNOWLEDGEMENTS
The author would like to express his deepest appreciation to his adviser, Dr. Hsiang
Wang, Professor of Coastal and Oceanographic Engineering, for his continuing advice and
encouragement throughout this study, also to coadviser Dr. James T.Kirby, Assistant
Professor of Coastal and Oceanographic Engineering, for his guidance and support. Appre
ciation is extended to Dr. Robert G. Dean, Graduate Research Professor of Coastal and
Oceanographic Engineering, Dr. D. Max Sheppard, Professor of Coastal and Oceanographic
Engineering, Dr. Ulrich H. Kurzweg, Professor of Engineering Science and Dr. James E.
Keesling, Professor of Mathematics, for their valuable suggestions in this research.
The author also wants to thank Dr. Y. Peter Sheng and Dr. Joseph L. Hammack for
their Suggestions and help.
The author is indebted to all staff of the Coastal Engineering Laboratory at University
of Florida, specially to Marc Perlin, Sidney Schofield, Vernon Sparkman, Charles Broward
and Jimmy E. Joiner, for their cooperation and assistance with the experiments.
Thanks also go to Subarna B. Malakar for the helping with programming and Ms.
Lillean A. Pieter for the drafting of figures.
Finally, the author would like to thank his wife Ran Liu for her support, encouragement
and patience.
TABLE OF CONTENTS
ACKNOWLEDGEMENTS ................................ ii
LIST OF FIGURES .................................... vi
LIST OF TABLES .................................... x
ABSTRACT ..................... ................... xi
CHAPTERS
1 INTRODUCTION .................... ............... 1
1.1 Statement of Problem .. .... ......... .... ......... 1
1.2 Past Study .................... ......... ....... 2
1.3 Scope of Study ... ......... ..... ............... 3
2 THEORETICAL EQUATIONS OF THE MODEL ................ .... 5
2.1 Introduction ..... ......... .... ..... .......... .... 5
2.2 Governing Equations ....... ....... ...... ... ........ 5
2.2.1 DepthAveraged Equation of Continuity ........ .......... 7
2.2.2 DepthAveraged Equations of Momentum ............... 10
2.2.3 Wave Energy Conservation Equation ................. 14
2.2.4 Dispersion Relation ...... ....... ........ ...... .. 30
2.2.5 Radiation Stresses, Lateral Mixing Stresses and Bottom Stresses 35
2.2.6 Wave Breaking Criterion ................. ...... 42
2.3 Summary ..................................... 44
3 NUMERICAL SCHEME .. ....... ................... 46
3.1 Grid System and Definition of Variables . . . . ... 46
3.2 Finite Difference Scheme for Combined Continuity Momentum Equation .47
3.2.1 Linearized Implicit Model ........................
3.2.2 Numerical Scheme for the Full Equations . ... . ...
3.2.3 Method of Calculation and Boundary Conditions ........
3.2.4 Stability Criteria ....................... .......
3.3 Finite Difference Scheme for the Wave Equation . . . .
3.4 M ethod of Solution ................................
4 PERFORMANCE OF THE MODEL ........................
4.1 Wave Setup on Plane Beach ...........................
4.2 Open Lateral Boundary Conditions and Longshore Current Distributions on
a Plane Beach ... ... ... .. ........ .. .. ...
4.3 Nearshore Circulation ..............................
4.4 Combined Refraction, Diffraction and Current Interaction Case . .
4.4.1 The Circular Shoal on a Flat Bottom . . . . .
4.4.2 Circular Symmetric Shoal on Plane Beach . . . .
5 EXPERIMENTAL COMPARISONS ........................
5.1 Purpose of Experiments .............................
5.2 Facilities and Instruments ............................
5.2.1 The Laboratory Facilities . . . . . .
5.2.2 Instrumentation .............................
5.3 Experimental Arrangement and Measurements . . . .
5.4 Experimental Results and Numerical Comparisons . . . .
5.4.1 Wave Shoaling on Plane Beach . . . . . .
5.4.2 Flow Field due to Inlet on Sloping Beach with no Waves . .
5.4.3 Inlet on Sloping Beach with Waves. . . . ...
6 CONCLUSIONS AND RECOMMENDATIONS . . . .. .
APPENDICES
A WAVE EQUATION DERIVED BY USING PERTURBATION METHOD
S63
64
66
66
70
81
81
82
82
86
89
90
91
94
95
134
. 137
B NONSLIP BOTTOM BOUNDARY CONDITION APPLIED IN WAVE EQUA
TION ....................................... .......149
C NUMERICAL SCHEME FOR CONTINUITYMOMENTUN EQUATION ... 152
D DOUBLESWEEP METHOD IN SOLVING NUMERICAL MODELS ...... 155
BIBLIOGRAPHY ........... .... ...................... 159
BIOGRAPHICAL SKETCH ........................ ...... 164
LIST OF FIGURES
2.1 Outline of Studied Area............................ 6
2.2 Variations of Dispersion Relation with U/co = 0.1 and c, ...... Linear;
. Hedges; Whitham; Composite . . . ... 33
2.3 Variations of Dispersion Relation with U/co = 0.1 and *.....* Linear;
. Hedges; Whitham; Composite . . . ... 34
2.4 Velocity Profiles for different friction models (2.107); (2.117) 39
2.5 LonguetHiggen's Longshore Current Profiles . . . ... 41
2.6 Variations of Eddy Vicosity Coefficients. = 0. ; *** =0.001; 
=0.01; Von Khrmhn's Model with upperlimit =0.001. unit= m2/s 43
3.1 Grid Scheme ................... ................ 48
3.2 Definitions ofXj and (Y,). ....................... 49
4.1 Setup on a Plane Beach. Eq.(2.138); Linear Wave; A Data 62
4.2 Longshore Current Generated by Obliquely Incident Waves.  no mix
ing; breaking mixing coeff.=0.01m2/s; breaking mixing coeff.
=0.01m2/6 and eddy vis. coeff.=0.005m/s . . ..... 65
4.3 Velocity Field of Periodic Bottom Topography . . .... 67
4.4 Wave Amplitude Variations along Centerline Aoko = 0.16; 
Aoko = 0.32 .................................. 69
4.5 Wave Amplitude Variations across Centerline Aokco = 0.16; 
Aoko = 0.32 ............. ... ........ ............ 71
4.6 Wave Amplitude Variations across Centerline Aoko = 0.16; 
Aok0 = 0.32 .................... ...... ....... ... .. 72
4.7 Bottom Topography : a Circular Shoal with Parabolic Configuration on
the Plane Beach ............................... 73
4.8 Description of Water Depth . . . . . 74
4.9 Velocity Vector Field for Circular Shoal Located on the Plane Beach (no
inlet) . . . . . . . . . .
4.10 Longshore Current Distributions at Lowest Boundary. No Inlet; 
with Inlet .. . ............... .........
4.11 Velocity Vector field for Circular Shoal Located on the Plane Beach with
Current from Inlet ................... ........
4.12 Contour Plot of Wave Amplitudes on the Circular Shoal Located on the
Plane Beach (no inlet) ............................
4.13 Contour Plot of Wave Amplitudes on the Circular Shoal Located on the
Plane Beach with Current from Inlet . . . ..... ...
5.1 Plan View of Wave Tank ..........................
5.2 CrossSection Views of Wave Tank . . . . .
5.3 CrossSections of Measurement . . . ..... ...
5.4 Velocity Distribution at I=Idry . . . . . .
5.5 Instrument Calibrations ...........................
5.6 Instrument Connections ...........................
5.7 Velocity Variations along the Beach, Y7(1,j) = 0; o (1,j) =
a2/(2 sinhkh) .. ... .. .. ... .. .. .. ... ... ..
5.8 Wave Height Variations along the Beach, i((1,j) = 0; o r7(1,j) =
a2/(2sinhkh) .. .. ... ... ... .. ... ... .. .. ..
5.9 Comparisons between m = 43 with Az = 0.12 m, ; and m = 83 with
Ax = 0.06 m, at I1 = 15,30 for m = 43 and 12 = 29,59 for m = 83
5.10 Wave Height Variations along the Beach (no current), Energy Eq.
(2.138); Linear Shoaling; x Measured Data . . .
5.11 Velocity profiles for Qo = 0.0056m3/sec. (no waves) . . .
5.12 Velocity profiles for Qo = 0.0076m3/sec. (no waves) . . .
5.13 Variation of xComponent Velocity. Case a: Qo = 0.0056m3/s; T = 0.87
sec.; Ho = 0.97 cm. Numerical Model; o Measured Data . .
5.14 Comparisons of Calculating Velocities between with and without Waves.
Case a: Qo = 0.0056m3/s, T = 0.87 sec.; Ho = 0.97 cm. with Waves;
no W aves ...... .............. ............
5.15 Variation of xComponent Velocity. Case b: Qo = 0.0056m3/s; T = 1.16
sec.; Ho = 0.88 cm. Numerical Model; o Measured Data . .
79
80
83
84
85
87
88
89
.92
.92
93
5.16 Comparisons of Calculating Velocities between with and without Waves.
Case b: Qo = 0.0056m3/s, T = 1.16 sec.; Ho = 0.88 cm. with Waves;
no W aves .. .......................... 107
5.17 Variation of xComponent Velocity. Case c: Qo = 0.0056mS/s; T = 1.36
sec.; Ho = 0.97 cm. Numerical Model; o Measured Data. ... 108
5.18 Comparisons of Calculating Velocities between with and without Waves.
Case c: Qo = 0.0056m5/s, T = 1.36 sec.; Ho = 0.97 cm. with Waves;
no W aves ................... ............. 109
5.19 Variation of xComponent Velocity. Case d: Qo = 0.0076mS/s; T = 0.87
sec.; Ho = 1.23 cm. Numerical Model; o Measured Data . ... 110
5.20 Comparisons of Calculating Velocities between with and without Waves.
Case d: Qo = 0.0076m3/s, T = 0.87 sec.; Ho = 1.23 cm. with Waves;
no W aves ......... ........................ 111
5.21 Variation of xComponent Velocity. Case e: Qo = 0.0076m3/s; T = 1.16
sec.; Ho = 0.95 cm. Numerical Model; o Measured Data . ... 112
5.22 Comparisons of Calculating Velocities between with and without waves.
Case e: Qo = 0.0076m3/s, T = 1.16 sec.; Ho = 0.95 cm. with Waves;
no W aves ................... ............. 113
5.23 Variation of xComponent Velocity. Case f: Qo = 0.0076m3/s; T = 1,36
sec.; Ho = 0.97 cm. Numerical Model; o Measured Data ...... 114
5.24 Comparisons of Calculating Velocities between with and without Waves.
Case f: Qo = 0.0076m3/s, T = 1.36 sec.; Ho = 0.97 cm. with Waves;
no W aves ................... ..... ........ 115
5.25 Velocity Field and Contour Lines in Flat Bottom. f = 0 . .. 116
5.26 Velocity Field and Contour Lines in Flat Bottom. f = 0.01 . 117
5.27 Contour Plot of Velocity. Case i: Qo = 0.0056mS/s (no waves) . 118
5.28 Contour Plot of Velocity. Case a: Qo = 0.0056m3/s; T = 0.87; Ho =
0.97 cm . . . . . . . . . 119
5.29 Contour Plot of Velocity. Case b: Qo = 0.0056m3/s; T = 1.16; Ho =
0.88 cm . . . . . . . . . 120
5.30 Contour Plot of Velocity. Case c: Qo = 0.0056m3/s; T = 1.36; Ho =
0.97 cm ....................... ...... ...... 121
5.31 Contour Plot of Velocity. Case ii: Qo = 0.0076mS/s (no waves) . 122
5.32 Contour Plot of Velocity. Case d: Qo = 0.0076m3/s; T = 0.87; Ho =
1.23 cm ... ... .... ... ... .. ... .. ... . 123
5.33 Contour Plot of Velocity. Case e: Qo = 0.0076m3/s; T = 1.16; Ho =
0.95 cm .. . . .. . . . . .. 124
5.34 Contour Plot of Velocity. Case f: Qo = 0.0076m3/s; T = 1.36; Ho =
0.97 cm . . . . . . . .. 125
5.35 Bottom Friction Coefficients f as Function of um/Ui . . : 126
5.36 Wave Amplitude Variations. Case a: Qo = 0.0056mS/s; T = 0.87 sec.;
Ho = 0.97 cm. Numerical Model; A Data . . . ... 127
5.37 Wave Amplitude Variations. Case b: Qo = 0.0056m3/s; T = 1.16 sec.;
Ho = 0.88 cm. Numerical Model; A Data . . . ... 128
5.38 Wave Amplitude Variations. Case c: Qo = 0.0056mS/s; T = 1.36 sec.;
Ho = 0.97 cm. Numerical Model; A Data . . ..... 129
5.39 Wave Amplitude Variations. Case d: Qo = 0.0076ms/s; T = 0.87 sec.;
Ho = 1.23 cm. Numerical Model; A Data . . . ... 130
5.40 Wave Amplitude Variations. Case e: Qo = 0.0076mS/s; T = 1.16 sec.;
Ho = 0.95 cm. Numerical Model; A Data . . .... 131
5.41 Wave Amplitude Variations. Case f: Qo = 0.0076m3/s; T = 1.36 sec.;
Ho = 0.97 cm. Numerical Model; A Data . . . ... 132
5.42 Wave Amplitude Variations along Centerline . . . ... 133
LIST OF TABLES
Test Condition .................. ....
The values of ec and f in the Tests ...........
Bottom Friction Coefficients Used in the Model ..
Incident Wave Height ...................
. .. . 90
. .. .. 96
. . . 98
.. 99
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
NUMERICAL MODELING OF CURRENT AND WAVE INTERACTIONS OF AN
INLETBEACH SYSTEM
By
YIXIN YAN
August 1987
Chairman: Dr. Hsiang Wang
Major Department: Engineering Sciences
A numerical model for the prediction of nearshore hydrodynamics of an inletbeach
system is developed in this study. The theoretical formulation is based on the time and
depthaveraged equations of continuity and momentum, and a 3rdorder wave energy equa
tion that retains the important nearshore effects of refraction, diffraction, currentwave
interaction and energy dissipation. Variational principle as well as a parallel method by us
ing perturbation technique is applied to arrive at the same wave equation. The application
of 3rdorder wave equation is restricted to small incident wave angle and slowly varying
topography.
A numerical scheme is developed to solve the above governing equations. The double
sweep CrankNicolson finite difference scheme is adopted for the model. In each spatial
sweeping, a recently developed two timelevel implicit method by Sheng is applied. The
numerical model is tested against existing numerical models of simpler cases and is found
to yield reasonable comparisons for all cases.
Laboratory experiments are conducted to verify the numerical model and it is found
that the model performs well in predicting current field but grossly underpredicts the wave
amplitude in shallow water, particular, in the center region of the inlet jet.
CHAPTER 1
INTRODUCTION
1.1 Statement of Problem
Tidal and river inlets are prominent coastal features. They play a vital role in the
evolution process of coastlines. Along the east coast and gulf coast of the United States,
inlets are numerous because of the extensive barrier island system. In these regions, the
rules of inlets are particularly significant owing to the heavy development pressure and
the competing interests of various human activities. To maintain navigation and water
commerce, for instance, often requires inlets to be stabilized with manmade structures.
The disruption of littoral processes owing to the presence of inlet structures, on the other
hand, poses a major threat to down drift beaches, thus diminishing the value of recreational
beaches and endangering waterfront structures. Many coastal communities spent millions
of dollars installing coastal protective structures or simply transporting sand from one place
to another. However, these measures do not always render satisfactory results owing partly
to our poor understanding on complicated interactions of an inletbeach system.
The hydrographic and shoreline changes of an inletbeach system are the consequence
of sediment movement which, in turn, is caused by the hydrodynamic forces exerted on the
system. Therefore, to better predict the behavior of inletbeach interaction one must first
be able to predict the hydrodynamics.
Physically, the problem in hand is a complicated one. Both wave and current field are
modified by the nearshore topography and the inlet geometry. In addition, the spatially
nonuniform current will interact with the unsteady flow field induced by waves.
Mathematically, there is no unified wave theory in the first place capable of describing
waves traveling through large areas of varied depth. Furthermore, the current field that we
~
2
are dealing with has a different scale of variation from that of the local wave field. The
nonlinear interaction between the two systems is expected to be important as well as the
nonlinear behavior of the wave field.
1.2 Past Study
In the absence of an inlet, the water motion in the nearshore zone is primarily induced
by waves from offshore as modified by the bottom topography and/or shore structure.
The influence of bottom topography on wave direction, height and length has long been
recognized and as early as in the 1940s analytical methods were already being developed
to solved the so called wave refraction problem (Johnson et al. 1948; Arther et al. 1952)
by tracing wave rays based on the principle of geometrical optics. Since then, significant
progress has been made in treating the shallow water wave transformation phenomenon.
Arthur (1950) included the effect of current into the ray methl,. Dobson (1967) developed
one of the first computer codes to trace wave ray over an irregular bottom. Skovgaard et
al. (1975) extended the ray method to include the effects of wave energy dissipation due
to bottom friction. Noda et al. (1974) developed a finite difference scheme based on the
conservation of wave number and wave energy, that compute spatial distributions of wave
number and wave energy in shallow water. His method was extended by Wang and Yang
(1981) and Chen and Wang (1983) to wave spectral transformation and the inclusion of
nonstationarity, thus enabling the computation of spatial and temporal distributions of
random waves. Birkemier and Dalrymple (1976) and Ebersole and Dalrymple (1979), on
the other hand, expanded Noda's method by including the depth averaged continuity and
momentum equations to computer wave induced nearshore circulations. In the presence of
an inlet, the nearshore wave field is further modified due to the interactions with a spatially
nonuniform current and the diffraction effect of inlet jetty structures.
Traditionally, the problems of wave refraction and diffraction were treated separately.
Penny and Price (1944,1952) were among the first to apply Sommerfeld's (1896) solution
to wave diffraction around a sharp obstacle such as a semiinfinite thin breakwater. Other
3
solutions of similar nature were obtained for different geometries (Wiegel,1962; Memos,1976;
Monfetusco,1968; Goda et al.,1971).
Only recently, Berkhoff (1972) developed what is now referred to as the mild slope
equation that governs the combined refraction diffraction phenomenon. Various solutions
became available on wave field modified by different beachstructure configurations such as
waves in the vicinity of a thin groin (Liu and Lozano,1979; Liu,1982), edge waves by barred
beach topography (Kirby et al.,1981), and waves over circular shoal (Radder, 1979), etc.
Booij (1981) further introduced the ambient current element into the basic wave transfor
mation equation and derived the governing equation on variational principle. The resulting
equation is a hyperbolic type and is, in general, difficult to solve. Booij's equation reduces
to a parabolic type similar to Berkhoff's if the mild slope assumption is invoked. A number
of simple examples were given by him. Kirby (1983) systematically studied the currentwave
interaction problem and proposed a weekly nonlinear model appropriate to thr Stokes wave
expansion.
The present study extends the Booij's and the Kirby's model to the problem of inter
action between water waves and inlet flows.
1.3 Scope of Study
The primary purpose of this study is to develop a numerical model for the prediction of
wave characteristics and current field of an inlet beach system. Specifically, the mathemati
cal formulation of the problem was performed first. This task entailed the establishment of
the governing equations sufficient to solve the variables of interest. Out of necessity, simpli
fied assumption are made such that the problem is amenable to the present computational
capability.
First of all, the fluid is treated as incompressible. Second, the surface shear stress as
well as part of the turbulent shear stress (normal component) is also neglected. Third,
also one of the most major ones is that the irrotationality is assumed in the wave equation
allowing the existence of a wave velocity potential and a current velocity potential, although
4
such an assumption is in contradiction to both physical reality and mathematical results
of the combined velocity field. The equations governing the mean flow including both
continuity and momentum equations are depth averaged. The equation governing the wave
motion with due considerations of refraction, diffraction and current interaction is .quasi
three dimensional owing to the velocity potential assumption. Two parallel methods were
used to derived the wave energy equation: one based on variational principle and the other
on perturbation method. By introducing proper scaling factors, the wave energy equation
is further simplified by containing terms of only the same order of magnitude. Other
relationships including dispersion relation, radiation stresses, bottom shears and lateral
mixing compatible to the system under consideration are also formulated.
A numerical scheme is developed to solve the governing equations by double sweep
marching method. The scheme is implicit which allows for a large time step, especially the
two time level fully alternating direction implicit scheme utilized by Sheng (1981).
Attempts to verify the numerical model are made by comparing the present model
with known analytical and numerical solutions of simple geometries and with laboratory
experiments carried out in the Coastal and Oceanographic Engineering Laboratory. Due to
wave basin limitations, the experiments were carried out only for plane beach with an inlet
normal to the contours.
CHAPTER 2
THEORETICAL EQUATIONS OF THE MODEL
2.1 Introduction
In order to solve the flow field in a specified region, the governing equations with proper
initial and boundary conditions must be provided. Furthermore whether a unique solution
exists in the domain of interest or whether mathematical singularities will appear that will
cause the equations failing to give sensible solutions should also be examined.
In this chapter, the theoretical development of the governing equations for both the
mean flow field and the wave field is given in detail.
There are five basic equations employed in the computation including four flow equa
tions: one continuity equation, two horizontal momentum equations and one wave equation,
and one equation relating wave number and wave period.
2.2 Governing Equations
The basic equations, as mentioned before, are based upon the principle of conservation
of mass, conservation of momentum, and conservation of wave energy and upon the gravita
tional wave dispersion characteristics. The wave equation under consideration includes the
effects of wave refraction, diffraction and shoaling during transformation from deep water
to shallow water regions. These equations, together with initial and boundary conditions,
give a mathematical representation of the physical problem of interest.
Since the problem of interest entails waves propagating towards shore over a complex
bathymetry, wave refraction, diffraction and breaking are expected to occur as well as
currentwave interaction due to inlet discharge. Figure 2.1 depicts the physical system being
studied. The problem has two apparent time scales: a short time scale associated with the
Inlet
** / / i
i11 / *.
a 
Sa. ..
* Shoal
>. \a 
 ( 
x
y
Figure 2.1: Outline of Studied Area
7
oscillation of the incident waves, and a larger time scale over which the characteristics of
the incident wave and mean flow, such as wave height, incident wave angle and current are
perceptibly modified. The larger time scale may also include the effect of changes in wind,
tides and topography. Since our attention here is towards mean quantities,tnany of.them
may be averaged over a wave period to remove the direct effect of oscillation. In addition,
since we are mainly interested in the net transport quantities the knowledge in the detailed
structure of the velocity distribution over depth is not needed and the equations may be
averaged over depth, reducing the entire problem to two dimensional in the horizontal plane.
In the derivation of wave energy conservation equation, due to the extreme complexity
in mathematics, some restrictions have to be imposed on the orders of magnitude of wave
height, bottom slope and current strength.
In the model, the quantities to be determined are as follows: the horizontal components
of the mean current including waveinduced current and divergent current from the inlet;
the local wave amplitude under the influence of refraction, diffraction and breaking ; the
local wave angle and wave number and the mean water surface level.
2.2.1 DepthAveraged Equation of Continuity
A Cartesian coordinate system is defined with x in the onshore direction, normal to the
coastal line, y in the longshore direction and z vertically upward, as shown in Fig.2.1.
The summation convention adopted here states that in Cartesian coordinates whenever
the same letter subscript occurs twice in a term,that subscript is to be given all possible
values and the results in horizontal space add together.
The continuity equation reads
Op 8pui 8pw
p +p + = 0 (2.1)
9t 8xi az
where i = 1,2 represents x,y direction, respectively, p is the water density and ui = {u, v}
represents the {x, y} components of velocity,respectively, w is the velocity component in the
z direction.
8
In application, water density which only varies slightly is considered as constant. The
continuity equation reduces to
atii 8w
S+ = 0 : (2.2)
axi 8z
Integrating Eq.(2.2) over depth, we obtain
S+ )dz = 0 (2.3)
where rq is free surface elevation, and h is water depth. Using Leibniz's rule which states
that
a a(') a() 8lf a c
f (x, z)dz = dz + f(z, ) f (x )
T8 CW) (z) Z Tz ax
Eq.(2.3) can be written as
a f? ar} ah
_x f uidz ui 1, u ih + x I+ w 1I= o (2.4)
e U zi +zi
in which the subscripts after vertical stroke are at which the value being calculated. The bot
tom boundary condition (BBC) and the kinematic free surface boundary condition (KFSBC)
are given as
ah
i + w = 0 (BBC) at z = h(x, y) (2.5)
axi
8an a_
= w ui (KFSBC) at z= r7(x,y,t) (2.6)
Substituting (2.5) and (2.6) into (2.4) results in
a? 9 /a r
l + a udz = 0 (2.7)
Tt hzi h
We now introduce
ui = Ui + u
(2.8)
9
where Ui are time and depth independent uniform current velocities, j7 is the mean surface
elevation, u'1 and 1r' are the fluctuating components. u', can be separated further into two
parts:
u'1 = ut +u (2.9)
where u" and uf are the slowly varying wave motion and the fast varying turbulent fluc
tuations, respectively. The properties of these quantities, after taking average over wave
period, are expressed as
vi = 7i, 7 =
U' =, = o
where the upper bar represents the time average over one wave period
1 rT
f= fdt
T Jo
Therefore, (2.7) becomes, after taking time average
a+ Udz + u'dz 0 (2.10)
at 9zi h Bri fh
Defining
f u'dz .h udz
(h + ) (h+I)
the final form of the mean continuity equation is arrived at
S+ (UjD) = 0 (2.11)
Tt 8zi
where D = h+1 is total depth and Ui = U+ fi is total mean current velocity averaged
over one wave period with U, = {U,V}. Equation (2.11) represents conservation of mass
per unit surface area under the assumption of water density being constant in space and
time, and the bottom unchanged with time.
10
2.2.2 DepthAveraged Equations of Momentum
The horizontal depthaveraged momentum equations used here were originally derived
by Ebersole and Dalrymple(1979). A brief summary was given here. The NavierStokes
equation in the horizontal (x,y) directions takes the following form:
9ui aui 9ui 1 ap 1 ari+ 1 a8r
+ U + W + + = ,+ (2.12)
at a x z p axzi p ax p az
and in vertical (z) direction it is
Ow aw aw 1 ap 18rj, 1 Orz
a+t U+w + = + + + +(2.13)
9t 9x, Qz p Qz p 9x, p Qz
i = 1,2, j= 1,2
where p is the total pressure including dynamic and static pressures and riq is the shear
stress, where the first subscript is the surface that the stress acts on and the second subscript
is the direction of shear stress. Introducing the continuity equation (2.2) into (2.12) and
(2.13) we have
aui Oui+ u Ouiw 1 Op 1 Ori+ 1 arOT
+ + = + + (2.14)
7t 8j ax z p xi p Ox, p 8z
and
Ow Owu, aww 1 p 1 9rz 1 arz
7 + + = + + (2.15)
9t 9x, Qz pz p 9xj p Qz
Integrating over depth and utilizing the Leibeniz's rule, the following term by term results
are obtained:
17 a u j a 7 a 7
ii) d= uidz ui 1
Ih Ot a t J.h at
P Ou, u, 8 8s7 Oh
7) L dz = __ uiuYdz uiu I 97 uiuj 1h
i.ii f_7 9uiw
S) iW dz = uiW 1 j ut, Ih
Thus the left hand side of Eq.(2.14) becomes
OLHS = a a a ,,)
LHS = T uidz uiujdz ", (+uj 117 W1w17)
at _t, h X h ax
ah
Ui Ih (j lh + W Ih)
Recognizing the last two brackets of LHS are the KFSBC and BBC, respectively, we have
LHS = /uidz + a uiujdz
9t Th x9j fh
Again introducing the definition of
Ui = Ui + u; + u!
into the above equation and time averaging over a wave period yields term by term
i) idz = (,t + up + u)dz = (U D)
where Ui and D are total current velocity and total water depth respectively as defined
earlier, and
ii) I U dz = ( + u + )(i + "u + u)dt
Xii) "fh fui+
Noting that there is no interaction between wave orbital velocity and turbulent velocity,
Phillips (1977), the above term can be expanded as follows
a a a a uudz+
S Ui.Udz + + Uiu)dz + I uuYdz +
azj h
a dz
xi Jh
a UD)+ a (Ui'iD)+ a(UfiD)+ I a uFu.0dz +
a r
ax, Jh'
Introducing Ui = Ui + fi into the above expression, the LHS of Eq.(2.14) becomes
a(UiD) + a (UUD) ai ( D) + a J7 uudz +
I
a j uudz (2.16)
The right hand side of (2.14) is, after integrating over depth and averaging over time,
RHS = dz+ Aidz+ +r idz
p h xzi p h zx p Jh z
here water density is considered as constant. After neglecting the second term which is the
lateral shear stress due to viscosity, we obtain
R1 1 8 "r 1 8h 1 1
RHS= pdz + p I, + p h + 7,i I, Ti Ih
P daz Jh p zi p i p p
The datum of pressure is chosen as the atmosphere pressure, so that at free surface, p I,
is zero. The mean pressure at the bottom p Ih consists of two parts: hydrostatic pressure
and wave induced dynamic pressure, viz.
PIh = pg(h+) + Pd 
= pgD + d h
where pd is the dynamic pressure. Multiplying the above equation by h we have
1 ah 1 a y 1 ah
p I (gD9 ) gD +  Ih
P axi 2 8zi azi P azi
The RHS can now be written as
1 a /" 1 a 2y 1 ah 1
~~ pdz + (gDZ) gD + pd Ih + e,
p xxi Jh ,2 zx p di p
1
b (2.17)
P
where
Yei = rsi is time averaged surface shear stress
and
TU = ri \h is time averaged bottom shear stress.
13
Finally, combining (2.16) with (2.17), the time averaged and depth integrated horizontal
momentum equation is given by
a a a a r a I
(U.D) + a(U UiD) a ( ,D) +8i+dz+ _.ujtdz
18 f 18 Ia 1 h 1 1
pdz + (gDO) gD + pd h + Ti 6i
p 8zi h 2i aXi p aXi p p
or rearranged as
a(UD) + a(UUiD) = 9Dz _I i d ld+pf u dz
1 1 Ah 1 1 a t t (
P9D2i pf, pD] + Pd Ih + 1.i + T,, uiutdz (2.18)
2 p 8zi p p 9x Jh
Collectively, the second to the fifth terms on the RHS are known as the radiation stress or
the excess momentum flux due to the presence of waves, i.e.,
Sij = f p6idz + p u*uudz pgD 2i, piiujD (2.19)
h h 2
where
i1 i=j
Considering that the product of bottom slope and the mean dynamic pressure at the bottom
is negligible, we obtain
8 8 Bj 1 a8S D 8
(UiD) + (UiUjD) = gD  a(puU.) +
t 9xj 8zi p xj p xj
1 1
ri rbi (2.20)
P P
after assuming uiu is independent of depth. The term pt4u is Reynolds stresses and
their normal component (when i = j) can usually be neglected.
Writing in component form yields
XDirection
a a a I_ 1_as." 1 as_
(UD) + (UD) + (UVD) = gD +
ot xZ 9y x8 p x ':p ay;
D 97t 1 1
 + + TbZ (2.21)
P ay p p
YDirection
a a a 1 as,, 1 asyy
(VD) + (UVD)+ (D)=  gD' a
at zx Ty ay p ox p ay
D aT 1 1
a + Tz S by (2.22)
where re = pufu is the lateral friction due to Reynolds stress. In the model, Tri will not
be considered. The possible formations of T, and fT will discuss later in this chapter. They
are flow related variables.
The average continuity and momentum equations can be linearized. Noda et al.(1974)
gave a simple form of those equations after neglecting the nonlinear and time dependent
terms in (2.11) and (2.20). They are
(UiD) = 0 (2.23)
axi
1 8s,, 1
0 = gD p as (2.24)
azi p xj p
The surface shear stresses and lateral shear stresses are also neglected.
2.2.3 Wave Energy Conservation Equation
The equation considering both refraction and diffraction with current has been devel
oped by Boiij (1981) and corrections were made later by Kirby (1984). Both used varia
tional principle. In this study, a method that modified Kirby's derivation is developed. In
15
parallel, the perturbation technique is also applied. Same results are obtained and are given
in Appendix 1. Here only the derivation from variational principle is presented.
The variational principle, first applied in the irrotational flow motion by Luke (1967),
states that the variation of a certain conservative quantity, often an integral, vanishes if the
function that describes the evolution of the physical process is subjected to small variations.
It is represented as
6[ft Ldxdt] = 0 (2.25)
where x = xi + yj is in horizontal space. For irrotational flow, the energy is conserved and
L can be chosen as the Bernoulli sum:
L Jl[t + (V4)2 +gz]dz (2.26)
h 2
where
V is the three dimensional gradient operator,
a a a1
+ a3y + zk
and k are unit vectors in x, y and z directions, respectively.
The potential function consists of two parts
4 = w. + t' (2.27)
where 4, and 0c are wave potential and current potential functions,respectively. A weak
point about the expression of 0c is that the mean current field may be rotational, since
we do not know a priori whether the mean flow field derived from momentum equations
satisfies the condition of irrotationality everywhere, i.e.,
curl U = 0
However, the irrotational assumption will be used here.
16
For a slowly varying bottom topography 4 can be expressed as (Kirby, 1983):
O(x,y,z,t) = 640(jXz, py) + e,1(x, 1/2ypypzXzt)+
C'2w2 (z, 1/2y, py, px, z, t) (2.28)
where p is a scale parameter with an order of magnitude O(Vhh/kh) representing the spatial
change over the space of wave length in which Vh is horizontal spatial gradient and k = 27r/L
is wave number where L is wave length. 6 is the strength parameter of current defined as
Vhc/V 'g. c is the Stokes wave steepness parameter e = 0(ak), a is the wave amplitude.
The spatial scale parameter p is introduced to differentiate two different spatial scales
in the model: one is related to the fast wavelike functions, and the other one is connected
to the slow spatial modulations, such as bottom topography, the changes of current speed,
wave length etc. Therefore, we have
x = x+Iix=X+X2 (2.29)
y = pl/29 9+ =Y +Y2 (2.30)
and their derivatives are replaced by
a a a
= L+ p (2.31)
Tz aX +' X2
a /2 a a
P + p (2.32)
ay IY a Y2
The scales as selected imply that spatial derivative with respect to phase function occurs
only in the xdirection which is true when the angle of wave propagating is small with
respect to the xaxis, i.e., it is assumed that no fast wavelike variations occur in the y
direction. It is consistent with the forward scattering approach that the complex amplitude
A is allowed to vary as 0(p1/2) in the ydirection. The problem is further considered to be
quasisteady; time derivatives only appear in wave term.
There are three magnitude parameters 6, p and e in the formulation. Different ratios
among those parameters will yield different forms of the energy equation. In the present
17
study, we adopt Kirby's (1983) model given the relationship as follows
S1, =1 (2.33)
which represent slow varying topography, strong mean current, and small incident wave
angles.
In equation (2.28) 4,i and ,w2 are first and second order potential function, respec
tively.
1wl = A'(x, p'/2y, x, pry, t)fm(zx, py, z) (2.34)
0w2 = A'(x, pu'y, /px, py, t)fm(px, py, z) (2.35)
The Aj can be split into two parts: a modulated amplitude function and a phase function:
Am(Z, p/2y, pAy, t) = A'(p/2y, px, py)eim + c.c (2.36)
where A'i is the complex amplitude, i is the phase function and c.c is the complex conju
gate. The phase function is defined as
S= kdx t (2.37)
where k, is the wave number in main (or x) direction and w = 2r/T is the angular wave
frequency in which T is the wave period. The real wave direction which is different from
the main direction is absorbed in the amplitude function. The subscripts and superscripts
in the equations represent the order and harmonics, respectively. However the absolute
superscript values should be no bigger than the subscript values. Substituting (2.34) and
(2.35) into (2.28), yields
4 = C2,c((z,C2y) + Am(z,ey,C2x,C2y,t)fim((e2x,y,) +
e2Am(x, Cy, c2x, e2y, t)fl( 2Z, ey, z) (2.38)
where 6 and p have been replaced by C2 and c2, respectively, by virtue of Eq.(2.33).
Substituting (2.38) into (2.26), we obtain with (2.37), term by term
A = cAf +m ,C2 A f(m
9ti t T C "2t/Z
(2.39)
18
VO! = Vhe + ehfA'I+ fA i+ Cfl AVhA + c3A'Vhfl
+eCffA 2 Afmym IJ+t + ffAA4 + CA4V hf f
+E2fAX' + 2efAj + C f2VhA' + eA'Vh
+(cArf, + c'Affc)k (2.40)
where
Vh is the horizontal gradient operator,
Vh + 3.
Therefore, by definition we have
VhO = U +V (2.41)
where velocities include wave induced mean current.
The first two terms in Eq.(2.40) represents linear wave superimposed upon uniform
current over flat bottom. In Eq.(2.26), the (Vq)2 term when computed to fourth order
becomes
IxAh A _4 In An x I2 2X
(V))2 = (VhRc)2 + e flff!AXAX f+ 1f1finA1yAy + ff/2"/2"'A2Xn
+f fm A An + n cnffAA + e2UfAm + eC2Vf AY
CJlJlr "x + 12z.!Zzt2 "2 + 2U Ax an1 1 ,Am22VAfiArnyl
+e2 2flVh, VhA + eC2AVh, Vhfim + e22Uf2Ax + c32VfmAm,
+C42f2mVhc VhA + Ce2AmVhO Vhf2 + C 2fmfnAm A
4 mm in m n I, rS n m rn C n rAn m fn o
+e 2fff x2AAiX + c32fffnA2 AixAx + e 2Am Azffz (2.42)
The subscripts X, Y and z are partial derivatives. The superscripts m and n represent
interaction between two harmonics.
Free surface envelope can be considered as two parts: one is related to mean current
as t7o and another one is driven by waves as tr'. Later we will see tro has the same order as
IU.
1 (, y, t)= 7o(x, y) + 1'(x, y, t)
(2.43)
or in the stretched coordinates
rl(x, y, ez x,2y,t) = 90(622x,,2y) + '(x, fY, c2x, y, t) (2.44)
where 7' has the order 0(c) at most and it can be expressed as
17 = 07' + 2 72 + 0(e3) (2.45)
where ti and 172 are first and second order surface elevations, respectively. By virtue of
(2.45), (2.43) can be rewritten as
r7 = 7vo + Ci + c217 (2.46)
Expanding (2.26) in Taylor series about z = iro, results in
L = [t + (V)2]dz + g7 gh2
h 2 2
= lot + 1(VO)2]dz + 1'[ t + I(VO)2].=o, + 1 a[ + [ (v)]=o
1 13 a 1 1 2 1
6 zlot + (V)2],= + O2 g+ +*** (2.47)
Introducing a shifting coordinate defined as
z' = z to, h' = h + ro
and substituting into (2.47), we get, after dropping the primes on h' and z' for convenience,
L = _[t + 1 (V)']dz + 17'[O + (VO)2]=o + 10 7 2 a[tt + I (V)2],=o
1 ,s 1 1 12
+i z7 3[t + (V)2)']o + g7 gho (2.48)
where we denote ho as still water depth. The expansion is carried out to fourth term only
as the higher order terms will not affect the third order wave equation.
20
Certain integration in (2.48) that involve fir of Eqs.(2.34) and (2.35) are defined as
follows for future use:
ff"dz = Q9~
(2.49)
f" ff"dz = P."
/ J] *13
ft f.dz = Ri"
(2.50)
(2.51)
where i, j, m and n are free indexes and subscript z represents partial derivatives. From
Stokes wave theory, the wellknown expressions for the f, and fin are
fl cosh k(h + z)
1 1 cosh kh
2 _2 = cosh 2k(h + z)
cosh 2kh
The typical results in (2.49)(2.51) that will be used later are given as follows:
cosh kh
1 k sinh3 h
1 1 + 2 cosh2 kh
S3k cosh kh sinh3 kh
p22 1 h sinh 4kh
1
R12 = 4k
R2 coshsinh
3 cosh kh sinh Ich
(2.52)
(2.53)
= k2 sinh 4kh
S2sinh kh k
where c is the wave celerity, c = L;
co is the group velocity, cg = c n with
n = (1 + h .
a = /gk tanh khi is the local angular wave frequency, it will be discussed in section (2.2.4).
Substituting (2.39) and (2.42) into (2.48),with the definition (2.49)(2.51), we have
L = eQA P +xA + PAA+ A iY Ai +" R" AiA
+eQ`UAx + c2 QmVAy, + EQDVh0 c VFhAm + 'Pmn"Am An
+ PmAxA x + c3RmA + (12Ao + c(c2 + + c32ir, n) + e rL Ai
+c of Uf,"A + cA v(Vf"A A+ + cf f A A~ + e f4/rvhfAj V' A
+ cm f2nAL + e 2mfnA + e(A + e 3 kfl m A
S2 it 2 1 i ix 4 1 m
13m fi"l /I AIA tA? & + r lel"~g+e '2 f lfj'lz"x Ae 2 V
+ // i"A + C r f A + Jm mVf/inAy + 4_k mUf2nAnA n
+c i f AU f(A + _+ + CVfJAzz1
2 2
+e k7UffA n + ,E4imInzzA, + ?S3mU IzAaX (2.54)
where the condition of Vhfin = 0 has been involved. All terms are evaluated at z = 0. For
simplicity, we only keep terms containing il and A1, since as will be shown L varies with
respect to I and A4 only. In (2.54), the superscripts are free indexes (nonrepeating) which
may not assume the same value for different terms, but they give same value in one term.
The terms of 2i are defined as
S= n mntlf + m+n= k (2.55)
S= Skmn7rit l7 k+m+n=j
(2.56)
where
when m = n or k = m = n, 6mn and Skmn equal to unit, otherwise they assume the
values of 2 and 3, respectively. Those give
o = o ^ +22r%71
After introducing the following
L = Lo + L+ E + e + e4 + (2.57)
Equation (2.54) can be separated by grouping the terms of the same order:
0 (e)
1 11
La = 2p1"AmxA + 3lfiArA7 + Q3VA fi1 + 2 r/nA
After introducing the followin(2.59)
L = L0 + cL1 + 22L2 + cOLs + E4L4 +... (2.57)
Equation (2.54) can be separated by grouping the terms of the same order:
o(C)
0(e2)
L2 1 "P IAxA x 2 +1 Q11 VA1 y + 2 + 9 1 Ait
,Uf Ax (X (2.59)
O(C3)
L3 = QVhc VhA+ + Pim"A xAnX + R"nAA + g + g7 + mf2nA
+t2mfnA + r fjnm AfxAIx + fiJfifjAjAn + rlVfnAy,
+?1UfA + Uf!A + m + UfA (2.60)
+rtllvf A t + rUllvfnAur + 2 "f2inAlt 2nUflAx (2.60)
1 =rmn Y' AD + pmnAmv "AiXAnx mtfnVi," VhA .rVl.nf2nA 2y
L4 Am= !PAIIA + y + AIA"AA + fifl VAhAc + t?7'Vf2 A1y
mmk n I+lfm ft A n n nmv nn
+tl f7lZ2nAl iA2 x fim I2n AXx + z",A" 1 + 2 n I Airx "
I1 1 1
+1 /22 1 AxAx + t + t' 2 VI + 1 / ,A
+ 1 n + 2vVfnAn 1 +n + 17k U n
+ifmnAmjA, A I??V/rA~y, + ^2AUf^n + U f n, A X
1 1
+ m "Vn.LA n+ smUfnf, AUx (2.61)
6 i
There is no Lo term, as it contains neither r7l nor A, term. For each order, varying
with respect to i71 and A,, we will obtain a series of equations. The variational principle
gives
8L a aLi aL
B a(i) Vh[a( = De (2.62)
7B 7t aB, a(VhB)
where B represents either 171 or A1, D, is the dissipation term that is caused by bottom fric
tion and turbulence that eventually are converted into internal energy, which is irreversible.
(One possible formation of De is given in Appendix A by using perturbation method.)
Equation (2.62) can be written in detail as
a a a aa a a a a
S( ) ( ) ( )]Li = D, (2.63)
aB aX2 aBx, X Y1 Yi By Y2 B, t aB ~
The following discussions are based on the orders of L.
1) 0(C)
Varying (2.58) with respect to r17 yields
170= (Vh0) (2.64)
2g
which states that the surface elevation by mean current is simply balanced by kinetic pres
sure. The same result is obtained if one uses Bernoulli equation at order 0(1).
2) O(e2)
Equation (2.59) varied with respect to r7', results in
n = fi(An + UAx)
9
Sf DAi (2.65)
g Dt
Sis the total derivative related to the x direction
Da a
Dt Qt ax
For n = 0 or the zeroth harmonic, A' is independent of X and t, which results in
rl =0 for n=0 (2.66)
From Stokes linear waves, we know that
Al = ige (2.67)
201
A I = ig ei (2.68)
201
where a is a complex amplitude and al is defined as
a1 = w Ukl (2.69)
which is the apparent local angular frequency that represents the waves riding on unidirec
tional flow. According to linear wave theory, the dispersion relationship is
a = w Ukl Vk2 (2.70)
where kc = k cos 0 and k2 = k sin 0 are the wave numbers in x and y direction, respectively
and 0 is the incident wave angle as shown in Fig.2.1. From (2.69) and (2.70), we find
a1 = a + Vk2 = a + 0(c)
since for small incident wave angle, sin 0 0 = 0(c). Hence for n = 1 and n = 1,fi1 =
f I = 1 at z = 0, we have
S 1 DA y a
1 1= ae (2.71)
71 g Dt 2
1 IDA'1 aei (2.72)
g Dt 2
both are the first harmonic solutions and the same as given by linear wave theory, rl
and rij1 have the same amplitudes, but traveling in opposite directions. As mentioned
before, the complex amplitude "a" absorbs the difference between real wave direction and
its main direction, since phase function f is only defined in x the direction. Next, varying
L2 with respect to Aj, (j = m or n), we have
Rm"A1 fm ir Pr nAxx Uf rl~ix = 0 (2.73)
a) let m = 0 (the zeroth harmonic), it gives
RnOA = 0
11 l=0
Therefore,
A= 0 or fo_= 0 from (2.51) (2.74)
Later we can prove that if fiz = 0, then fo = 0.
b) Considering the first harmonic, we obtain
[cc9(kc k2)] = 0 (2.75)
9
It is important to note that (2.75) is of the order of 0(E3), since k2 k2 = k = k2 sin2 0 =
0(c2). This term later will be combined with other third order terms to form the wave
equation which is of the order of 0(c3)
3) O(0c)
Varying L3 with respect to ri and Al (j = n, m or k) respectively, we obtain
1 1
g9'n + fAf + 2/mnA"A A + I/f AFA" + VfA"AY + n Af
+UfenAn + t7UfizAlx = 0 (2.76)
R"Ag + rff/",A Vh (orVh c) PnAxxx f1nt7 ff,(Ax )x
VfW Ufn2 f flr,/1 1X = 0 (2.77)
26
we emphasize again that the superscripts only take the same values in one term; there is no
relationship among different terms. There are only two unknowns, r2 and A', in the above
two equations and so that they can be solved.
a) Collecting all the terms with zeroth harmonics in (2.76), we get
2 1 f. l Ai + f1A1A1 + fl '21lt A + t~l1A1t] + ff5[t71AXA
17 = A IA + AA1 1+ AlA t]+ iA
+rj1lAxl)}
Here the condition of fi = fi1 has been applied.
The values on the right hand side of above equation are known, so that
0gk( 9k' tanh2 kh k tanh kh
2 = gk tanh2 ktan ] (2.78)
4a2 4o2 2
Since we only consider up to 0(e), kcand al can be replaced by k and a, and vice versa.
Equation (2.78) becomes
o sinh 2k + 0(C4) (2.79)
2sinh 2kh
This is the setdown of mean free surface with reference to the shifted coordinates. From
(2.77), we obtain
A0 = V h" (QVhc) (2.80)
where m can assume any arbitary value. Hence, A can not be determined as f remains
unknown. However we will show later that A0f2 is nonzero. From (2.80), the previous
statement in case 2 is confirmed, that is, let m = 0, we must have fo = 0 for f/, = 0.
Otherwise, A will approach infinity.
b) The first order harmonics in (2.76) yield the following solution
I = 1[VA if2aiA'] (2.81)
9
From (2.77), together with (2.81), we have
A2 0 (2.82)
9
27
For m = 1, if we assume that f2 has the same form as fl, the terms in the bracket are
identically equal zero. But, if m is any value less than or equal to 1, i.e., m < 1 1, we may
select A' = 0 such that
1 = VA Iy, = iVayei (2.83)
al
which states that the surface elevation of the first harmonic of the second order of the wave
is balanced by the variation of wave amplitude modulation in longshore direction.
c) For the first negative harmonic, (2.76) and (2.77) give
A2' = 0 (2.84)
r71 = iVay,ei (2.85)
O1
d) Considering the second harmonic, after lower harmonic terms have been dropped, we
obtain from (2.76) and (2.77), respectively
9rl + C2 2 + 2 2 12 2 2 1  1 1 1
gr72 + f2A + + fl 'A1 + f2UA2x + f.rt7\A, + flUtA\x = 0
2 2
and
R12A 2 12 14_I 12 DA2 2 _1 X U2t
R1 A + fl' A7 PA12Axx x (Aix')x UV2 f,
fLialx = 0
Solving two equations for two unknowns AJ and ,2 yield
A = i. e' + O(Iei) (2.86)
16
= I 2k cosh kh
7 = cosh (2 + cosh 2kh)e"2O + 0(s) (2.87)
The8 sinhs kh
They are the same as the second order Stokes waves.
4) 0 (4)
Varying L4 with respect to ,i, gives
fnVhc VAn + fm2.A A" + f JfAn AA + 7f2.An + t A
Srm ~n mn n r+n nn 1 r
+7 f1/zAxAAx + f11fizztAA" + U fzAx + 7Uf"Anx + 2iA
1 An ,nx n 1n I I I 2X
+i2/ffIAnX + 2^ V f2AY1 2+ I zVfAAy, + Vf2"Ay = 0
For the first harmonic, we have
D A2
VhAc'VhAA 1 1 + fil A2 A + 0 A1fifA x+22zl1
V.c AVAfl + f1 A1'AA + flhf^2A1A2 + flf2A2Ax + fL' 1D
1 2 1 21 D1
+ft0 + flz2r + fIfz(2AxAx + A1271) + flfflZ
If l l1 1 ADA 1 2 ADAr1
(2AAr "1 + A1 11) + AfZ2 + f1lz2" = 0
where
DA' g
Dt 2 '
S_ _9_ae ;
Dt 2
DA2 3 2 (2eiz
Dt 8
The only unknown in the equation is f2zA, so that
1 2 2F 12
fIAo = [Vh~c Vh() + g Vh Va + .aI a 2] (2.88)
0 01 f0oo 1
in which
F = gk2 (3 2 tanh2 kh + 4 sinh2 kh sinh2 kh cosh2 kh)
16sinh' kh
Next, varying L4 with'respect to Ai, also considering first harmonic, we get
fl/22zl 2 + flJ' zf27A2 + tf22ri1) + +f \l (2A1
29
+2Ar1) + f rl('r1 2 +1 0 )1 i + 1 1 llAy
(P'A x)x2x P11(A) x Vh (1Vh,,) f1(1Ax)x (1 A x
97 1X 1 122 
rlAf)x 1 Aj + ALx) VAxy, = 0
where again fo = f1 = 1 is applied. Each term in the above equation has been calculated
before. After some tedious algebraic manipulation, we obtain
2(cc,kc + aiU)ax, + 2alVay2 i(cc V2)ayy, + a[(ccgc + U )al
+()y,]a + ikcck'l a 'a = 0 (2.89)
al
where k' = kc cosh 4kh + 8 2 tanh2 kh
where k' = k.
C9 8sinh4kh
kI and 1l in (2.89) can be replaced by k and a without introducing errors larger than 0(E4).
Equation (2.89) contains information of wave refraction and diffraction, as well as wave
amplitude dispersion and wavecurrent interaction. Finally, as we have noted before the
first harmonic for L1 varying with respect to A1 as expressed by (2.75) is third order and,
therefore, should be combined with (2.89) to reach the wave euqation. Hence, we obtain in
the (x, y) coordinates
,,ccgkl + U (w V
2(cc,ki + oaU)a, + 2oaVay i(cc, V')ay, + a cck U ) + ()v
01 U1
+icc (k' k2)]a + ikcc,k'l a Ja = 0 (2.90)
This is the final form of the wave equation without dissipation. If dissipation is to be
considered, then a dissipation term D, could be included in the left hand side. One possible
formation of De is given in Appendix 1. Following Kirby and Dalrymple (1983), a phase
shifted wave amplitude is introduced to simplify the wave number computation
a = Aei(kofz kdz) (2.91)
where ko is the reference wave number. Equation (2.91) implies that
A" = a" ei(k'xkoz)
(2.92)
30
where a" is the real wave amplitude. The final wave equation, after dropping the double
primes, takes the following form:
2(cc,k; + oU)A. + 2aVA, i(cc, V')Ay + 2[(c ).
V
+(V) + icco(k2 k2) + 2i(cck + o1U)(ko ki)]A + i(o, a2)A
(i 1)gk c2(I
csh kh )A + 2(ccgki + oaU)WA = 0 (2.93)
cosh2 kh 023
where the cubic amplitude term in Eq.(2.90) has been replaced by i(o0 o0)A (Kirby and
Dalrymple,(1986)), in which
a = gk(1 + fle2D) tanh(kh + fze)
it will be discussed in next section. The last term on the left hand side is the decay of
energy due to wave breaking and the details of which will be discussed in section 2.2.6, and
W is defined as
W = 0 before wave breaking
W = 2(1 2 I) after wave breaking, 7 = 0.39 and q = 0.4
The dissipation term De is given in appendix A. Equation (2.93) is the final form of the
wave equation. It is similar to the wave equation given by Kirby (1984) except the ad
ditional imaginary term in the second last term. (see Liu 1986). For plane beach and in
the absence of current, linearized (2.90) reduces to the conventional onedimensional wave
energy conservation equation. It reads
a = ao(coocy)(bo/b)
where quantities with subscripts 0 refer to known condition and b is the width between
adjacent wave orthogonals.
2.2.4 Dispersion Relation
When waves travel over a varying topography or, are embraced by current, the wave
number vector will change its value. The four equations established in the previous section
31
enable us to solve the four unknowns: ij, U, V and H, provided the wave number vector is
known. Since wave number is a horizontal vector, in general two equations are required to
solve its magnitude and direction. The wave number irrotationality and the wave disper
sion relationship derived from dynamic free surface condition are conmmonli used. In the
presence of a current field the mathematics could become complicated if the nonlinear be
havior is fully accounted for. In the present study, a seminonlinear dispersion relationship
is considered and is incorporated directly in the wave equation; thus the linear dispersion
relationship is used in the model. Also, by assuming small incident wave angle, the di
rectional information on wave number is incorporated into the complex wave amplitude as
presented in the previous section, thus, avoiding the inclusion of the wave number irrota
tional equation explicitly. The linear dispersion relation with current is given as (Dean and
Dalrymple, 1984)
w = a + Uki + Vk2 (2.94)
where
a = V/g ktanh kh
kl = k cos 0
k2 = k sin 0
In shallow water, tanh kh  kh approximately, Eq.(2.94) becomes
w = kvgh + Uki + Vk2 (2.95)
Several modifications have been proposed. Whitham(1967) gave a Stokes nonlinear disper
sion relation for intermediate depth region in the absence of current:
oo = or 1 + C2Dk (2.96)
where
c = ak
cosh 4kh + 8 2 tanh2 kh
Dk =snh
8 sinh4 kh
32
As waves enter shallow water, (2.96) becomes invalid as Dk approaches (kh)4
In the shallow water region. Hedges (1976) proposed a different formula
o'o = vgk tanh k(h + a) (2.97)
for the case of no current in which a is wave amplitude. For the same wave frequency, it
yields longer wave length than that given by the linear dispersion relation.
More recently, by combining both cases discussed above, Kirby and Dalrymple (1986)
constructed a composite model which satisfies both intermediate water depth, Eq.(2.96),
and shallow water depth, Eq.(2.97),
a"o = [gk(1 + f2Dk) tanh(kh + f2z)]1/2 (2.98)
where
fi = tanh5 kh
f2 = (kh/ sinh kh)'
In deep water region, (2.98) approaches (2.96) as fz  0; conversely, in shallow water
region, it becomes (2.97) as fi + 0. Now assuming there is no nonlinear current effect on
wave dispersion, the presence of current causes a Doppler shift and Eq.(2.98) becomes
w = t'o + Ukl + Vk2 (2.99)
Figures 2.2 and 2.3 shows the results of w for a constant nondimensional velocity U/co =
0.1 and varying of E from 0.01 to 0.4, where cc = vgk is the deep water wave celerity.
The dispersion equations given by (2.94), (2.96) and (2.97) are also plotted in Figs.2.2 and
2.3 for comparison.
The emprical nonlinear dispersion (2.98) can be used to modify the wave equation as we
did in Eq.(2.93), in which the Eq. (2.98) is incorporated directly in the governing equation
in the form of a term which directly produces the required amplitude dispersion (Kirby and
Dalrymple, 1986). Consequently, the linear dispersion relation (2.94) will be used in this
model.
S=0.01
0.00
%.0oo
0.50
0.50
1.00
1.00
Figure 2.2: Variations of Dispersion Relation
Hedges; Witham; Composite
1.50
1.50
2.00
2.00
2.50
2.50
with U/co = 0.1 and c, ..... Linear; .
E=0.05
\ ______
I.
.
I 
o
0
r~^I
0
0
34
o
e=0.2
0N
o
.00 0.50 1.00 1.50 2.00 2.50
Kh
0
o
e=0.3
.00 0.50 1.00i 1.50 2.00 2.50
Kh
E=0.4
o .

0
00 0.50 1.00 1.50 2.00 2.50
Kh
Figure 2.3: Variations of Dispersion Relation with U/co = 0.1 and c, ...... Linear; 
Hedges; Whitham; Composite
35
Once the dispersion equation (2.99) is brought into the model, the cubic amplitude term
in Eq.(2.90) should be dropped, since the amplitude dispersion has already been included
in the nonlinear dispersion equation (2.98).
A wave number vector which is irrotational can be defined as:
k = Vw (2.100)
where k is a wave number vector, k = {ki,k2} and w corresponds to the scaler phase
function. Using the mathematical property that the curl of a gradient is identically zero, it
is shown that
V x V = 0 (2.101)
or
Vx k = 0 (2.102)
In general, Eq.(2.102) together with Eq.(2.99) is necessary to solve the wave number
field of a domain.
However, in the present study the wave angle information has already been included in
the complex wave amplitude given by Eq.(2.92). Therefore equation (2.102) does not need
to be included explicitly. However, the assumption that wave approaching angle is small
must be observed.
2.2.5 Radiation Stresses, Lateral Mixing Stresses and Bottom Stresses
In the wave momentum equations (2.21) and (2.22), there are some undetermined terms:
radiation stresses, lateral mixing stresses and bottom shear stresses. They are not indepen
dent variables and are calculated by the following equations.
Radiation Stresses
The radiation stresses are given in (2.19):
Si = pSidz + ulu dz IpgD 2j puiitD
h fh S 2
I
36
The form was simplified by LonguetHiggins and Stewart (1962) to second order for a linear
progressive wave system. For normal incident waves, it can be approximately by
[S]= S (2.103)
[S.O S 1 [ E(2n 0.5) 0
So So o 0 E(n 0.5)
where E is local wave energy E = pgl a 12, a is complex wave amplitude and n is the ratio
of wave group velocity to the wave celerity.
In general, when waves enter beach with an angle, Sqi can be obtained by
[ Si ] = A ][ S [ A ] (2.104)
where [ A ] is directional matrix,
[A]= icos0i sin]0 (2.105)
[ sin cos 0
hence
S. SY nE[n(cos2 + 1) ] Ensin 20
[Sj]= (2.106)
Sy Syy En sin 20 E[n(sin2 0 + 1) ]
where 0 is local wave angle.
Bottom Shear Stresses
The bottom shear stresses are customarily expressed as
yM = pf uti, ut (2.107)
where ut is total velocity composing of both wave orbital and mean current velocities, and
uti is its component form. The quantity f is the friction factor. The magnitude of total
velocity I ut  is equal to Vu2" + v2. The u and v velocity components can be expressed as
S= U + u = U + u cos 0 (2.108)
v = V + uv = U + u' sin 8
(2.109)
37
where U and V are mean current components as defined before. u' and u' are wave orbital
velocities in the x and y direction, respectively and u" = [(ut)2 + (uw)2]1/2.
The total velocity ut I can be expressed as
Iut = [U2 + V2 + (u,)2 + 2Uu" cos 0 + 2Vu" sin 0]1/2 (2.110)
The local wave orbital velocity u' is purely oscillating in time
uw = ucos at (2.111)
where um is the maximum wave orbital velocity at the bottom
um aI(2.112)
sinh kh
The time averaged bottom shear stresses in the z and y directions then become
b = Pj (U +umcosO cost) j at I dt (2.113)
f 0T
Tby = P (V +umsin cos at) ut dt (2.114)
Alternative, for computational convenience, it can be rewritten as
= p (U+umcos cos t) It I d(at) (2.115)
b = Pf (U + u sin cos et)I ut d(at) (2.116)
These integrals can be calculated by applying Simpson's rule. The serious problem by
applied the above bottom friction model is the computational time. A simplified model is
introduced with
7 = pf I u't I U (2.117)
where
I u't I= V/U + V2+ U
38
since there is no interaction term between wave orbital velocity and mean current velocity.
In the component form the above equation gives
b = pf lu't U (2.118)
Yb = pfl 't V (2.119)
The difference between the two models (2.107) and (2.117) is minor for cases U > um as
shown by the example in Fig.2.4, in which the same bottom friction coefficient is used (the
testing conditions refer to Fig.5.3 in chapter 5), but the saving in computational time is
more than 60%. In fact the small difference can be offset by increasing the bottom friction
coefficient in (2.117) to account for the absence of u, term. Equation (2.117) will be used
in our model.
Another very important factor is selection of the friction coefficient f. Based on field ob
servations and laboratory experiments LonguetHiggins (1970) suggested a value about 0.01
for longshore current calculations. The bottom friction coefficient for waves alone presents
a more difficult situation. For smooth turbulent flow over fixed bed, Kamphuis(1975) gave
f = (2.120)
where Re is Reynolds number, R, = la and v is the kinematic viscosity.
For rough turbulent flow, friction coefficient no longer depends upon Reynolds number.
Investigations were conducted by Jonsson(1966)and Grant & Madsen(1979). The most
commonly used relation is
f exp[5.213(')0'14 5.977] y < 0.63
/ = (2.121)
0.15 r > 0.63
where r is the roughness of the bed and a' =1 a 1. Equation (2.121) is an approximate
explicit solution given by Swart (1974) to an implicit, semiempirical formula by Jonsson
and Carlsen (1976).
SIVJ
00!0
(33S/W31)
00 '0
n
0
0
0
C)
00'
00'02
(3IS/N3) n
00 *0
(3]S/W3)
o 0
O ,a
* N
oo
o V
0
o I
a,
o tU
0 ;fj
o
o 5
 0&
=rV
0
00'OZ
n
00 *C
O o'C
r
I
40
The situation becomes further complicated once combined wave with current, the state
of the art of modeling such flows is rather rudimentary. The total velocity in (2.107) becomes
much larger than the wave orbital velocity in the present case of very strong current. The
actual values of f used in the model will be adjusted as a function of (a, T) based on
experimental results to be discussed in chapter 4.
In application, the bottom friction coefficient may be treated differently in x and y
directions, in the present study, the bottom friction coefficient in y direction is taken as one
half of the value in x (or main) direction.
Lateral Mixing
Lateral mixing due to turbulent diffusion is a complicated flow phenomenon. Rigorous
mathematical treatment is beyond the state of art. Figure 2.5 shows the effects of lat
eral mixing on longshore current distribution (LonguetHiggins, 1970 a.b). The common
assumption is that turbulent mixing is proportional to the local mean velocity gradient. Fol
lowing LonguetHiggins' formulation on longshore current mixing, the lateral shear stress
is written as
aU av
e = P au + a (2.122)
where e, and cy are the horizontal eddy viscosities in the z and y direction respectively. We
define
ei = Ei + ei i = y (2.123)
where e,, and cci are the eddy viscosities of wave and current, respectively which are
assumed to be independent. The inclusion of Eci is based on the fact that even though
there is no wave, the flow from the inlet still generates turbulence. As the simplest possible
approach, LonguetHiggins proposed a model for erw, that is proportional to the distance
offshore, z, multiplied by velocity scale
Ce = Nzxg
(2.124)
0
SWithout Mixing
With Mixing
I
a
z
JJ
Breaker Line
DISTANCE OFFSHORE
Figure 2.5: LonguetHiggen's Longshore Current Profiles
where N is a dimensionless constant whose limits LonguetHiggins gave as
0 < N < 0.016 (2.125)
Following Ebersole and Dalrymple (1979), N is chosen as 0.01 and ew is allowed to vary
linearly with z to the breaking line. From this point seaward the coefficient remains at
this constant value. There is no elaborate consideration of ey. Ebersole and Dalrymple
assumed a constant value and here it is chosen as half the value of c,, throughout the field.
The eddy viscosity c, theoretically, will change from point to point. One of the formulas
was due to Von K&rman's similarity hypothesis:
21dV
Ce = mi j (2.126)
dU
e = m dy.  (2.127)
ay
where
dV d2V dU d U
ml = k 2 and m2 = k 
dz dzZ dy dy2
and k = 0.4 is Von Khrman's constant. In the numerical calculation, these coefficients go
to infinit near a point of inflection. An upperlimit value should be imposed. Alternatively,
both eddy viscosities can be treated as constants. Figure 2.6 give the effects of eddy viscosi
ties on the velocity profiles (the testing conditions refer to Fig.5.3 in chapter 5). It shows
that they have only minor difference between constant 0.001 m2/sec and Von Khrmhn's
mixing model with upperlimit equal to 0.001 m2/sec. Therefore in this investigation, in
stead of Eqs.(2.126) and (2.127) the constant eddy viscosity is selected.
2.2.6 Wave Breaking Criterion
When waves enter shallow water, they will eventually become unstable and break,
dissipating the energy in the form of turbulence and work against bottom friction. For a
very mildly sloping beach, spilling breaker occurs and the simplest criterion predicted by
solitary wave theory is (McCowan (1894))
H
y = ()b = 0.78 (2.128)
where D is total water depth and subscript b denotes the value at breaking.
Dally et al. (1984) proposed that the decay of energy flux with distance in the surf
zone is proportional to the excess of energy flux above a stable value, given the relation
=(Ec) = [Ec, (Ec,).] (2.129)
where q is a constant related to the rate of energy decay. The quantity (Ec,), is the stable
energy flux for a broken wave in water of total depth D. The last term in the wave equation
(2.93) which appears only after wave breaking
2(cc,ki + ~rU)WA
I ..
43
0
'O.O0 0.80 1.60 2j.40 31.20 4.00
Y (M)
o
0 1=16
C
'0.00 0.80 1.60 2.140 3.20 4.00
Y (M)
Figule 2.6: Variations of Eddy Vicosity Coefficients. = 0. ; ... =0.001; =0.01; 
Von Khrman's Model with upperlimit =0.001. unit= m2/e
I
can now be completely defined as
W = 41 A1 (2.130)
by virtue of Eqs. (2.128) and (2.129). The value of q is chosen as 0.4 following Kirby and
Dalrymple (1986).
2.3 Summary
In this chapter, the analytical model of a currentwavetopographic interaction system
is formulated under the assumption of strong current and slowly varying bottom. Basic
equations together with essential terms are developed. They are summarized as follows:
DepthIntegrated Continuity Equation
aB aUD 8VD
+ +  0 (2.131)
at 8x ay
DepthIntegrated Momentum Equation in XDirection
U + au V au ay 1 aszs 1 as la2 ) 1
 + U + V = g Da + L (2.132)
at axz y zx pD z pD ay p ay pD
DepthIntegrated Momentum Equation in YDirection
av av av aBT 1 aas., 1 as 1 a8a 1
+ U + V g vu + T 7y (2.133)
9t x ay ay pD ax pD 9y p a8 pD (2
in which we have substituted (2.131) into (2.21) and (2.22), where
'S, Sf E[n(cos2 0 + 1) ] Ensin 20
[Sii= (2.134)
S,, Su Ensin20 E[n(sin2 0 + 1) ]
9U 8V
e = P'eya + Pav (2.135)
ay az
4i = eWi + Cci; ci is a constant and ez is given by
wx = NzvgD with N = 0.01;
1
twy = 2wz
T = pfI u' I Ui (2.136)
u' 1= i U2 + V2+ + U. .(2.137)
Um=
Ssinh kh
Wave Equation
2(cc,ki + olU)A, + 2a1VA, i(cc, V2)Ayy +O [(c1 kl U1)z
al
V
+(V) + icc,(ki k2) + 2i(cckl + alU)(ko kl)]A + i(o a2)A
01
(i 1)gk2 i.
(S h)2 V()A + 2(ccgkl + o'U)WA = 0 (2.138)
cosh2 kh 2w
where
oO = [gk(l + fle2Dk) tanh(kh + fze)11/2
fi = tanh5 kh
f2 = (kh/ sinhkh)4
Scosh 4kh + 8 2 tanh kh
8 sinh' kh
S W = 0 before wave breaking
W = (1 I) after wave breaking, y = 0.39 and q = 0.4 (2.139)
When waves enter the surf zone, the last term in (2.138) will be brought into calculation
till waves reach stable condition.
The dispersion relation used in this model is
w = a + Uki + Vk2 (2.140)
and the wave angle is included in the complex wave amplitude A.
CHAPTER 3
NUMERICAL SCHEME
In this chapter, the numerical schemes for solving the continuity, momentum and wave
equations are presented, along with boundary and initial conditions.
The numerical scheme requires (1) the definition of a grid system with a system which
provides a systematic method for identifying variables of interest and (2) the conversion of
the governing equations into finite difference forms.
3.1 Grid System and Definition of Variables
A rectangular grid mesh was established over the area of interest as shown in Fig.3.1,
where x and y denote the onshore and longshore direction, respectively. Two kinds of
variables are employed in the computation: the gridcenter variable and the gridside vari
able. Practically all the quantities of interest are represented by the gridcenter form which
expresses the value at the center of the grid and has the following notation:
where X is the gridcenter variable, the subscripts i, j represent the ith grid side in the
xdirection and the jth grid side in the ydirection, respectively, and the superscript n
represents the nth time step. The gridside variable expresses the value of the quantity
along any grid side. The only variable that appeared in this form in the present study is the
velocity vector; it is denoted by a subscript s to differentiate it from the gridcenter value
and has the following form:
47
Y, is the vector quantity along a grid side: the subscripts are such that i denotes the ith
grid point in the xdirection on the left of the vector and j denotes the jth grid point in
the ydirection on the right of the vector.
In the present study, the gridside velocity vector was applied in the continuity and
momentum equations and a gridcenter velocity vector was applied in the wave equation.
The definitions of gridside and gridcenter variables are shown in Fig.3.2
3.2 Finite Difference Scheme for Combined Continuity Momentum Equation
The CrankNicolson finite difference scheme in double sweep form is applied here to
solve the continuity and momentum equations and later also to solve the wave equation.
3.2.1 Linearized Implicit Model
A variety of numerical schemes are presently available to solve the equations for nearshore
hydrodynamics similar to those developed in Chapter 2. Ebersole and Dalrymple (1979)
presented an explicit model for nearshore circulation which requires small time steps due
to the stability criterion, thus relatively long computational time. Vemulakonda (1982)
presented an implicit model allowing for much larger time steps. Both models adopt a
three time level schemes. Butler (1981) and Sheng(1981)proposed a two time level implicit
model. This method not only saves computational time but also shortens the development
of the difference scheme for the governing equations. To demonstrate the structure of their
model, the simplified linearized equations are discussed here. The depth averaged linearized
continuity and momentum equations are
aO ae9Q aoQ
at + Z + a = (3.1)
9t ( 2 a y
r= gDL (3.2)
at ax
X
I
I
I
I
I
I
I
I
!
I
Y
I
I
I
I
I
I
I
i
I
I
I
I
I
I
I
i
!
!
I
I
I
I
I
n+1 n
2 1
Figure 3.1: Grid Scheme
M
M1
M2
I
AX
Y
y
YI,I+1

Y
Figure 3.2: Definition of X.. and (Y,),
I
Yij
aQv a9r
a = =gD (3.3)
at ay
where Q, = UD and Qy = VD. U and V in this section denote the gridside values. In
matrix form, linearized equations can be rewritten as
Wt + AW, + BW, = 0 (3.4)
where
\ 0 1 0 0 0 1\
W= Qz A= gD 0 0 B= 0 0 0
Qv 0 0 gD 0 0
Sheng (1981) implemented a two time level fully implicit scheme in computing the external
model in his three dimensional mode.
+11 1
,,j + I r(Q)n+l ti
at a + [ i+ (Q.)! ] + [(Q) (Qy)1] = (3.5)
(Q,)! ,, (Q))" + D [(T,,tl
At AX 0
At AQ + d ] 0 (3.6)
\"t (Qy)?, gD _n+l n+li. (3.7)
At + Z~yh ,
or in the matrix form
t (Wn+ Wn) + A 6 + 6Y)W"n+ = 0 (3.8)
where &S and 8u are central spatial difference operators which provide the differentiated
quantities at the same points as the time difference quantities. Equation(3.8) involves only
the (n+1) time level in the last two terms, i.e., the whole equation is staggered in time.
Rearranging (3.8) by time levels, we obtain
(1 + 2fI + 2/y)Wn+1 = W (3.9)
where
At At
fzt AS,; fl MB6
2Ax 2Ay
51
By adding the quantity of 4zflBy(Wn+l W"), which has a second order truncation error
as
asw
(At)2AB
8t8z8y
we have
(1 + 2P,)(1 + 2py)W"+' = (1 + 4PzP)W" (3.10)
This factorized finite difference equation is still a second order approximation to the differ
ential equation (3.4). The advantage of using (3.10) lies in the fact that the solution of the
factorized form can be split into two separate one dimensional operations. To achieve this,
we introduce an intermediate value, W*, such that
(1 + 2).)W* = (1 2py)W (3.11)
(1 + 2P8)W"+1 = W* + 26yW (3.12)
From now on, quantities without the time level denoted are the values at time level n.
Equation (3.10) is recoverable from (3.11) and (3.12) by eliminating W*.
A double sweep method is used to solve (3.11) and (3.12). At the first sweep, the values
at time level n are known we solve for W*. The full solution W"+1 then is obtained based
upon known W* values. After some manipulations, the first sweep in x direction gives
y + ,9z + Y6vQY = 0
(3.13)
Q+ Q. + =At
Q X Q.+ gD6S. = 0
and the second sweep in the y direction yields
r+ r + ~ TY6Q = 0
(3.14)
Qn+ Qv + tgD6,S"+ = 0
On the first sweep, Un+1 and i* are found, and on next sweep the values of Vn+1 and ri"+1
are obtained. Therefore, after one loop, the quantities at time level n + 1 are known. The
52
advantages of this scheme are (a) the finite difference equations are brief, at second sweep
(3.14) there is no Q,, involved; (b) a significant savings in computation time is realized with
large time steps and less computations.
3.2.2 Numerical Scheme for the Full Equations
Following the above model, the full equations (2.131), (2.132) and (2.133) in chapter 2
can be written in numerical form as (see Appendix C)
xsweep
S + z(Un+1D) + St (VD) = 0
(3.15)
(1+ Atc))Un"+ U + g6* + Atmi = 0(
ysweep
+ i +s 6(Vn+D) ,(VD)= 0
(3.16)
(1 + Aty)V"+1 V + AgY6,"+ + Atm2 = 0
where c.U and CyV are bottom shear stresses in x and y direction, respectively. By (2.136)
in chapter 2, we have
f f
cZ = u  and Cy = uI ,y
where I ut [I is the total velocity defined in (2.137) and here is specified in each direction.
mi and m2 are the remaining terms in Eqs. (2.132) and (2.133).
Again, on the x sweep the solutions for T* and U"+1 are obtained; these updated
values are brought into the y sweep to obtain i~i+ and V"+'. Writing the above two set
of equations in detail, we have
a) XSweep
Continuity equation
1 1
Aty(, ,i) + 2Ax [U, *^3.(Di+l,y + D.,) Ui'(Di, + D,)]
+ 2A [Vi,+l (Di,j+ + Di,j) V,,j(Di + Di,j1)] = 0 (3.17)
XDirection Momentum Equation
I(" t U.) + (U,+l Ui,,) + (Vi+1 + V,.1 + Vi,i+ + Viy)
At 2Ax 8Ay
(Uj+1 Ui,1) = ,i *,) Di,) z ) ( 
~p(Di, Dilj) 2Ay [(Sy)ij+l (Say)i,i1 + (Syv)ilj+l (Sy)i,j1]
(Di, + Di,i)1 (ut)ij ,i ij + (FLy)id + (FLys),i (3.18)
where
I (ut)i, I, = [(um)i,j + (um)i1,,l + 2(ut)i,j
um is wave orbital velocity as defined in (1.112) in chapter 2.
(,U),,i = {u(i + [0.25(Vi + j +1 i,, +V; ,+1)2)1/
(FLyy)i,i =
(FLuz)jj =
1
4(Ay)2 {(Uj+l Ui,[((EvY)i,, + (e),,)+1 + (')l, + (E'y)1ij+1
+2((Ec),,j + (e)i,,+1)] (U1,, U,.)[((,Ew)i, + ( ,J1
+(my)i1,i + (cy),i,,1) + 2((ey),, + (Ic),_)]}
1
^4az*y{O(V+1 Vi1.,i+)1)[((ew,)i.,+ 1+ (ew1),i + (ew)i,i,j
+(e'w)il,,j+) + 2((Ecz)i,j+l + (eez)i1,+l)] (Vi,. Vi,,j)[((ews)i,i
+(Cew)i_,i + (ew);,,_ + (Ews)i,._1) + 2((e~),i + (Ecc)il)]}
b) YSweep
Continuity Equation
1
( )+^ ) [ v1 [ (Dj+1i + 2,,), V^(Di,, + +,,i)]
[ayVj+i(D.j+l + Dj)* Vj(D.,1 + .,ji).] = 0 (3.19)
YDirection Momentum Equation
( .".* v)Un+l+ U"+ + q_1t + Un+) ( ,
(v.1n+l vj 1 J +,,,i
,tj,j 8+ (V+,)+ ( +ij V_,j,) +
zi 8A 2Ay
(V,,+1 V,j1) = i iD 1(D,, ,). 1(Sy)i+,
AY2 (Dj Di1)
(SOy)ilj + (S.y);+1,,1 (S.A),i,j1i 2 ((Sy)i,
P(Dii + Dijl1)* a
(Syv)i,1) 2(Di .+ ,)J ). (u)ii Vi +FL ),+ (FLzJJ '(3.20)
2(D,, + D ij 1 1.7
where the bottom shear stress coefficient has taken as half the value of that in x the direction
as stated in (2.2.5), and
I (ut)ij I, = [(um)i,i + (um),,i,] + 2(u,)i,j
(tz),j = (v~; + [0.25(Ul1 + U,,+,1 + U,, + U",+ _1)2}1'/2
1
(FLz)idj = 4(Az)2 {(Vi+l, Vij[((ew,)il, + ( ( (Cw)i+l, 1 + (ECwz)i,1)
+2((eoz),+i,i + (E)j)] (Vi,. Vil,4)[((E~w),,, + (Cw,)il,j
+(cw.),,j + (eCw).1,i1) + 2((e~.),, + (ecs).i)]}
4AxA1 l
(FLy)iJ = 4AxAy{(U +: U+ )[((Cwu)i+l + (uW)i+l,j1 + (Ey)i,j
+( ) ) + 2(()+,, + ( ,c)i+ ( 1)[((
+( )i,1 + ( ).i1,j + (cw),1,f1) + 2((eey)i, + (Y).,,_1)]}
In (3.17) to (3.20), the quantities with no time level are at time level n; and the subscript
(.) outside the bracket in(3.19) and (3.20) denotes updated values, i.e., D. = h + *. Also
we have used the latest x direction velocity, Un+1, in the y sweep to accelerate the conver
gence. In application, the (ec,) and (eey) are considered as constants as discussed in chapter 2.
3.2.3 Method of Calculation and Boundary Conditions
A double sweep method has been developed in each direction,i.e., in both xsweep and
ysweep. The advantage of this sweep method lies in the fact that it gives an explicit formula
at each sweep to calculate the results.
55
Collecting the terms from (3.17) to (3.20) by time level, so that the unknown quantities
are to the left and the known quantities are to the right, we have
XSweep
f j + Pi ,i UiRj Pij1,iU = (3.21)
(3.21)
Aij U," + RI V+1,j RI jij = Bij
YSweep
f+1 + QiV^ Q"n 1 = Zy
+ (3.22)
Rij "nl + R2 ri1 R2ifj1 = S,
where the quatities without superscripts are at time level n. The definitions of undefined
quantities are given in Appendix D, which also shows the details of the following derivation.
The x sweep equations, after some manipulations, can be reorganized as
( ,i = Fli + EliU+l
(3.23)
Un+i = F2i1 + E2_li_j
where
SCi Pji,F2,
1 + Pjy E2i
Eli = Pi
1 + Pi,,E2,
F Bi+lj RIFli+i
F2 =
Ai+lj + R1Eli+I
R1
E2i =
Ai+lj + R1Eli+i
The boundary conditions state that at the shoreline (i=idry), the velocities should be equal
to given velocities in the inlet section and they disappear elsewhere, and the mean water
displacement V, assumes to be identically equal zero. Therefore, the first sweep will solve
coefficients Fl, F2, El and E2 starting from i=idry till offshore line (i=l). Once those
coefficients are determined, the second sweep will directly yield results of T* and Un+'
starting from i=l (offshore boundary) till i=iwet, where iwet = idry1, is the last grid at
56
which water depth has nonzero positive value. The offshore boundary condition is given at
deep water, i=l, by specifying the mean water displacement to be zero. This is approximate
assumption to maintain a steady state mean flow in the region of concern with no increase
or decrease of water mass. A more rigorous condition is to allow for water level variation
at the offshore boundary such that the change in hydrostatic force is balanced by the net
momentum flux and bottom shears in the xdirection.
The lateral boundary conditions basically have three kinds for the concerned region:
the symmetric free boundary conditions; the noflow boundary conditions and the non
symmetric free boundary conditions. They will be specified for specific problems later.
The y sweep is similar to the x sweep and (3.22) can be rewritten as
ntl = F3j + E3 V.n+l
(3.24)
VAn1 = F4iI + E4 l,1
where
Zij QijF4j
F3i =
1 + Qi,,E4j
E3j = Qij1
1 + Qi,E4i
S i+l R2F3j+l
F4j Ri,i+ + R2E3i+l
R2
E4 = R
Ri'j+1 + R2E3j+i
The first sweep will determine F3, E3, F4 and E4, then the following sweep will give the
values of "n+1 and V"+1. The shoreline boundary condition for V is set equal to zero at
i=idry; and offshore boundary condition at i=l supposes that it satisfies the mass con
servation, i.e., (VD)a = 0. The lateral boundary conditions will be discussed for different
situations later.
3.2.4 Stability Criteria
The stability criteria of these schemes can not be established theoretically. In general,
the stability criterion for explicit scheme can be guided by the following inequality:
At < ((3.25)
1 uI 1 + g@
where I ut I is the magnitude of total velocity. For application, it becomes
At < Min.[ X (3.26)
for I ut [< /g
Actually, the critical time step, At, in our study by using a two time level implicit model
is at least one order larger than that given by (3.26).
3.3 Finite Difference Scheme for the Wave Equation
The final form of wave equation (2.138) in chapter 2 is mixed parabolichyperbolic
partial differential equation. The advantage of parabolic type partial differential equation
is that the marching method can be used in step by step calculation. Booij (1981) proposed
a parabolic model following on Radder (1979). Then the wave field can be split into a
transmitted and a reflected field. In his parabolic approach, the reflected field is completely
neglected. Kirby (1983) gave another approximation, after assuming that the complex
amplitude can absorb the difference between real and image wave directions for small angle
incident waves, and also assuming that wave amplitudes are changing more rapidly in y
direction than in x direction for slowly varying topography. Both of those assumptions
have been introduced in our derivation of the wave equation.
For the parabolic partial differential equation, the CrankNicolson finite difference
scheme is applied. For simplicity, (2.138) in chapter 2 can be rewritten as
,1 +C iCA 1 aCs
A + CoA C + C A+ A + (C4 + iCs)A = 0 (3.27)
CA 2z C2 ay
where
2aV
Co =
CO= 2o1V
P
cc V2
P
P
1
V
C3 = V
O*1
C4 = R+W
C5 = (cc,(ki k2)] + (ko ki) R+ (o 2)
in which
P = 2(ccgkli+oU)
R gk2
Pcosh2 kh V2
As the waves enter the surf zone, the second term in C4 will be added into calculation. The
quantity of j A 12 in a2 is the value evaluated at last iteration, unlike Kirby who specified
this value at the same iteration. The values of velocity components U and V are gridcenter
values in this calculation.
The CrankNicolson finite difference scheme constructs Eq.(3.27) as
Ai+ A, 1
A,+1i Ai + 4 [(Co)i+ij(Ai+i,i+i Ai+1,,1) + (Co)i,
(Aij+I A jI)] [(C)ji+1,,(A.i+j+1 2Ai+1,, + Ai+1,1) + (Cl)j,j
1 ((C2),+, (C21,s) (A1,
(Ai~+1 2Aj + AjI)] + ((2), (A+i+ + Ai,)
2Ax ((C2)i,+1 + (C2)idj)
1 ((C3)+1s),+ (C(3),j1)
+ 77V[ Ai+,+lj + Ai']
4Ay (Cz2)i+,, (C2)i,j
+ [((C4)i+lj + i(Cs)i+1j)A.+I,j + ((C4),,j + i(C5),i,)A1.j] = 0 (3.28)
A similar double sweep method is used to solve (3.28), resulting in
A.j = R2j1 + R1, Ai4j1 (3.29)
The definitions of R1, R2 and other details are also given in Appendix C.
59
The boundary conditions are stated as follows:
a) at offshore, i=1, the wave amplitudes are specified in complex forms;
b) at shoreline, i=idry, the wave amplitudes are supposed to die off completely;
c) the lateral boundary conditions can be considered as
A = iksinOA (3.30)
which assumes that the waves travel mainly in the x direction on a slowly varying topography
and Snell's law is valid for first order correction. The finite difference form of (3.30) gives
1+9
Ai.,+1 = Ai, Q (3.31)
1Q
where
Q = (iksin 0A )i,j+1
The local wave propagation angle 0, as discussed in chapter 2, can be found through
the complex wave amplitude. By Eq.(2.92) of chapter 2, we have
Ai+lai
In(A'al 'i i(k ko)Ax (3.32)
In(Aila+) = i2k2Ay (3.33)
where the double primes in Eq.(2.92) have been dropped. The local wave angles are in
cluded both in kl and k2, so that either of above equations can be used to obtain 0.
3.4 Method of Solution
With the specified initial and boundary conditions, the finite difference equations (3.21),
(3.22) and (3.28) can be solved in a loop of iteration.
The initial condition is assumed to be the state of rest. The velocity field, both U and V
components as well as the mean free surface displacement, if, are initially set equal to zero. A
bottom topography is input into the model. The initial wave amplitude field is calculated
60
by solving the wave energy equation with velocities equal to zero, which corresponds to
Snell's law for plane beach before waves break.
To avoid "shock loading" the model, both wave amplitude and velocities of inlet dis
charge are built up from zero to their termination value over a number of iterations. We
choose 100 iterations to set up those values by using the hyperbolic tangent function
F = Fo tanh(IT/25) (3.34)
where F is either wave amplitude or velocity of inlet discharge; Fo is their termination
values; IT is the number of iteration. The approach declares that the termination value
will be reached as total iteration is up to 100, while tanh 4 approaches the value of unity.
The program will run till the variables U, V, 7 and A reach their steady state. The
convergency is also checked by the balance of total discharge in and out of the controlled
boundary.
There are two coefficients to be determined in a running case: one is the bottom
friction coefficient and the other is the eddy viscosity coefficient. For fixed inlet discharge,
the eddy viscosity is assumed to be constant, but the bottom friction coefficient may vary
for different combinations of wave and current. Further explanations on how to select these
two coefficients are given in the next two chapters.
I 
CHAPTER 4
PERFORMANCE OF THE MODEL
To evaluate the performance of the numerical model, a number of tests were run, mainly
utilizing existing numerical models of simple cases as bench marks. Little attempt was made
to compare directly with existing experimental data for two reasons:
(1) Experimental data that can be used to test the model to its full term, i.e., current
wave interaction due to inlet on sloping beach, does not exist.
(2) The bench mark numerical models selected here usually already did their compar
isons with laboratory data, whenever available. Since the purpose of the present model is
not to seek improvement of any specific numerical model but rather to extend the capa
bilities to new applications, it is deemed sufficient for verification purpose to compare with
numerical results only.
4.1 Wave Setup on Plane Beach
The model was run for a case of normal incident wave on a plane smooth beach. The
conditions were as follows: T=1.14 sec, deep water wave height Ho = 6.45 cm and beach
slope = 0.082. This was the same case used by Ebserole and Dalrymple (1979) to compare
with laboratory data by Bowen, Inman and Simmons (1968). Since it is a two dimensional
case, we have selected m = 52 in the x direction (onoffshore) with Ax = 0.25m and n = 4
in the ydirection longshoree) with Ay = 0.4m. The water depth at the deep water end
in the numerical model was about 1 m and the mean water level displacement due to set
down can be neglected. The results are shown in Fig.4.1. The wave height was computed
on the basis of linear wave shoaling and the wave equation as expressed by (2.138) was
used for wave height computation. The data from Bowen, Inman and Simmons (1968) were
also plotted for comparisons. Apparently, these two approaches are similar except around
61
0
Q
0
V)
0 A
O
IJ
o
0A
o
0.00 80.00 160.00 240.00 320.00 400.00 L80.00
X (CM)
Figure 4.1: Setup on a Plane Beach. Eq.(2.138); Linear Wave; A Data
o (
o
Figure 4.1: Setup on a Plane Beach. Eq.(2.138); Linear Wave; A Data
breaking line.
4.2 Open Lateral Boundary Conditions and Longshore Current Distributions on a Plane Beach
Two typical open boundary conditions can be specified for the present numerical model.
One of them is the periodic lateral boundary conditions, which requires that
Q(, 1) = Q(i, N + 1)
Q(i, 2) = Q(i,N + 2)
Q(i, 3) =Q(i,N + 3)
and so forth, where Q is any quantity of interest. This boundary condition can obviously
be applied to topographies with periodicity such as the example given by Ebsersole and
Dalrymple (1979), but can also be relaxed to apply to an arbitrary topography that has no
periodicity as long as the boundaries are chosen far enough away from the area of interest
and there is no source or sink in the area.
In the case where a source or sink exists, such as the present study with an inlet,
the periodic lateral boundary condition can no longer be applied as the outflow from
downstream would not be balanced by the inflow from the upstream. An approximate free
boundary condition in this case is to let
8Q
0
ay
at far enough distance on both upstream and downstream. The above boundary condition
implies no longshore variations far away from the area of interest.
A typical example for openboundarycondition computation is the longshore current
generated by obliquely incident waves. The problem was studied analytically by Longuet
Higgins (1970). The case illustrated here is for the conditions: beach slope = 0.075; incident
wave angle = 50; bottom friction coefficient = 0.01; wave height = 0.02 m and wave period
= 1 second. The results are shown in Fig.4.2.
U
64
If there is no lateral mixing effects at all, the velocity decreases rapidly from maximum
to zero at breakline. The apparent existence of velocities outside the breaker zone is due
to the fact. that the advective acceleration terms in the difference y momentum equation
cause velocities outside the surf zone to be contaminated by velocities inside the breaking
zone. The strong mixing effects are demonstrated in the other two cases: one only has the
lateral mixing due to breaking with a mixing coefficient = 0.01 m2/sec and the other in
cludes both mixing from breaking waves with coefficient= 0.01 m2/sec and flow turbulence
with an eddy viscosity coefficient = 0.005 m2/sec.
4.3 Nearshore Circulation
In this example, currents induced by waves over irregular bottom topography are com
puted. The bottom topography is slightly modified from that used by Noda et al. (1974)
and Ebersole and Dalrymple (1979) and is generated by the following equation:
h = a x[ 1 + Ae(z/20) sin10 (y V tanf) ] (4.1)
where
s: beach slope = 0.025
x: distance from shoreline
A: periodic beach length = 104 m
A: amplitude of bottom variation = 10 m
8: = 30*
in which V is used instead of z on the right hand side as in the original equation by Noda.
This modification is necessary to reduce the variations on the lateral boundary conditions.
We choose m=40 in xdirection with Ax = 5.0 m and n=27 in ydirection with Ay = 4.0
m. The other parameters are given as
T = 4.0sec.; ao = 0.75m; 0 = 50; f = 0.012
breaking mixing coefficient = 0.01m2/sec.
O
0
0
O
*
. %
o V.
10 <
2.20 0.00 2.20 4.40 6.60 8 80
V (CM/SEC)
Figure 4.2: Lonshore Current Generated by Obliquely Incident Waves. no mixing; 
breaking mixing coeff.=0.01m2/s; breaking mixing coeff. =0.01m2/s and eddy vis.
coeff.=0.005m2/8
66
The resulting velocity field is shown in Fig.4.3. Although this case offers no direct com
parison with the result by Ebersole and Dalrymple, the flow patterns are similar and the
current field is physically plausible.
4.4 Combined Refraction, Diffraction and Current Interaction Case
This example deals with the wave refraction problem as illustrated in the previous
three cases, but also the combined effects of refraction and diffraction under the influence
of topography and the interaction with current. A case of a submerged circular shoal with
a parabolic configuration resting on a beach is studied. Similar topography is commonly
seen around inlet where an ebb shoal is located outside an inlet. This idealized topography
was first investigated by Radder (1979) and then by Kirby and Dalrymple (1983) on a flat
bottom with normal incident waves. Therefore, the numerical model was first run at the
same condition used by Radder. Then, the case of obliquely incident waves was considered
with the shoal on a plane beach. Finally, the model was treated with the presence of an
inlet to the above configuration.
4.4.1 The Circular Shoal on a Flat Bottom
Radder (1979) suggested two configurations of shoals in his study of parabolic wave
equation. Configuration I, of a circular shoal on a flat bottom is represented by the depth
profile
{ h= h + eor2 r
S= ho r>R (4.2)
where
r2 = (Z x.)2 + (y _y)2
eo = (ho hm)/R2
Xm and ym are the coordinates of the center of the shoal, and hm is the water depth at that
point; ho is the water depth at the flat bottom area. The portion of r < R prescribes a
4 4 c I . 4  . 4C 4 C . 4. 4 ' N N
4 4 4 .  .~ 4 44 4 ,, ,C ' N4 N N 
4 4 4~  . ~ *. .
4 4 4 L 4 4 4 C 4 *4 * C4 C ) N~ N N N . *. 
 4 ~~ ~ ~ 4 4 4. . c~
.c ..4 4c,)/ 44 v.e~N .
4.c .  N 4,,4' 4 ~. .~\ C
S C4~~,//C ~ ~~4 . C .* .
c .5.5 '5. ) N N 4... .
 4 5. 5..
5..

/ / 4 
4/'
I
 4 
.*/S .*
.** * *
.*****es
.,,, 
0.596E+00
s5o MPXIHUH VECTOR
Figure 4.3: Velocity Field of Periodic Bottom Topography
68
parabolic curve. The numerical values selected are
hm/R = 0.0625; ho/R = 0.1875; Lo/R = 0.5
where Lo is the wave length in the flat bottom area.
For such a configuration, the wave rays of a linear refraction approximation are focus
into cusped caustics. Peregrine (1983) and Kirby and Dalrymple (1983) discussed this
case and have shown that two jumps in wave condition fan out in approximately the same
manner as the caustics emanating from the region of the cusp. The region between the
jump conditions, directly in the lee of the submerged shoal, would then be dominated by
waves of approximately uniform amplitude propagating in the incident wave direction.
The grid spacings used by Radder are described as
Ay/Az = 0.5; and Lo/Ax = 1,2,4 and 8.
and it is noted here that the larger the value of Lo/Az, the better the resolution. For the
shoal studied here, Lo = 0.25m is used, the corresponding Lo/Ay = 2.5. The center of the
shoal, with the radius R = 10Az(5Ay) is located along the line of symmetry with respect
to the yaxis.
Without consideration of current, the wave energy equation becomes
2ccgkiA, iccgAyV + [(kiccg). + iccg(k2 k2) + i2(kiccg)(ko ki) ]A
+i(hl (i)A i1gk )A + (2ccgkl)WA = 0 (4.3)
cosh2 kh V 2w
Following Kirby and Dalrymple (1983), the nondimensional incident wave amplitudes are
chosen as koAo = 0.16 and 0.32 in which ko = 21r/Lo. The variations of wave amplitudes
along the center line (y/R = 0) in the x direction are shown in Fig.4.4 in which x/R = 0
starts from the foreedge of the shoal with the center of shoal located at z/R = 1.0. For both
cases, in the neighborhood of the center of the shoal the wave amplitudes drop significantly
below the incident wave amplitude. The peak value of koAo = 0.32 is less than half the value
of koAo = 0.16; it indicates the effects of nonlinearity with the increase of wave steepness.
Figure 4.4: Wave Amplitude Variations along Centerline Aoko = 0.16; Aoko = 0.32
70
In other words, the wave energy is much easier to be trapped and harder to decay for the
shorter waves at the line of symmetry. Also, there is a slight phase shift in the position of
maximum wave amplitudes. Figures 4.5 and 4.6 gives the wave amplitude variations in the
longshore direction at the location of z/R = 3.0, 4.0, 5.0 and 6.0, respectively. y/R = 0 is
the symmetric center line of the shoal. These plots show that as waves travel downstream
from the shoal the influence zone becomes wider. The spreading of the influence zone
for both cases are almost the same. Since wave amplitudes for koAo = 0.16 have a large
amplification at y/R = 0, which can also be shown in Fig.4.4, they drop rapidly to below
the incident wave amplitude to keep the wave energy balance.
These results are similar to those described by Kirby and Dalrymple (1983).
4.4.2 Circular Symmetric Shoal on Plane Beach
In the case of the circular shoal located at the center of a uniform sloping beach as shown
in Fig.4.7, the water depth, referring to Fig.4.8, is prescribed by the following formulas:
h = z tan a r > Rcosa
h = xtan a h' r _< Rcosa (4.4)
where
r' = (X Xn)2 + (y y,)2
h' = h"/ cos a,
h B + VB2 4C
h' =
2
R2 2r
B = and
hn tan2 a sin a
r2 R
C=
sin2 a tan2 a
The values used in the examples are: R = 0.7521m; hn = 0.085m; Ax = Ay = 0.15m and
the center of shoal is located at grid point i=17 and j=13.
The model is first run for a case with no inlet and then for a case with the presence of
an inlet. In both cases, the input waves form an angle of 50 with the normal of the shoreline.
,,
x/R=3.0
O"s
gU,
c.
C3
In
in V
o *_
o o 1'. 30 260 3. 90 5.20 6.50
Y/R
in
x/R=4.0
o
cr
"n I/
in
.00oo 1.30 2.60 3.90 5.20 6.50
T/R
Figure 4.5: Wave Amplitude Variations across Centerline Aoko = 0.16; Aoko = 0.32
in
.r
x/R=5.0
C
A,
o
b0oo 1.30 2.60 3.90 520 6.50
Y/R
in
x/R=6.0
Ca
.
in
.oo i'.ao 2'.60 3.90 5.20 6.50
Y/R
Figure 4.6: Wave Amplitude Variations across Centerline Aoko = 0.16;  Aoko = 0.32
I
Figure 4.7: Bottom Topography : a Circular Shoal with Parabolic Configuration on the
Plane Beach
h
Figure 4.8: Description of Water Depth
The other input conditions are
T = 1.0sec.; ao = 0.02m; a = 0.075; f = 0.01
breaking mixing coefficient = 0.01m2/sec.
For the case of no inlet, the eddy viscosity coefficient is equal to zero and for the case with
inlet, a value of 0.005m2/sec is used. The inlet discharge is equal to 0.007m3/sec.
The velocity field with no inlet is shown in Fig.4.9. Waves break at the top of the
shoal. The current induced by the waves is characterized by a strong onshore flow over the
shoal consistent with the wave incident angle; at the lee side of the shoal, the flow splits
into two directions, one upstream and one downstream. The upstream flow develops into
two vortexes rotating in opposite directions; the downstream flow forms a longshore flow
system with two portions. This longshore current pattern at the downstream boundary is
sketched in Fig.4.10.
The velocity field with the inlet is shown in Fig.4.11. The effects of inlet flow are seen
75
to be quite pronounced. First of all, the two upstream vortices are strengthened and pushed
farther offshore. Secondly, the downstream longshore current system is also strengthened
with a flow reversal close to the shoreline, which is caused by the lower mean water surface
elevation at the midsection of the shoreline in the shadow of the shoal and also caused by
the lower pressure due to the higher velocities in the inlet, and the influence is extended
further offshore which is due to the fact that the momentum of the inlet discharge pushes
the flow into the lower pressure region. The mean water surface elevation is also lowered in
this region. The longshore current distribution at the downstream boundary for the case
with inlet is also shown in Fig.4.10. The location of null flow on the lee side of the shoal is
now being pushed seaward onto the shoal due to the inlet flow momentum.
Figures 4.12 and 4.13 show the contour lines of wave amplitudes. The unit in these
figures is 104 meter. For the case of no inlet, the wave refraction effect causes the waves
to focus on the shoal and this leads to wave breaking and a low region just leeward of the
shoal crest. The diffraction effect, on the other hand causes waves to focus on the leeward
of the shoal and creates a local high region there.
The effects of current are most pronounced in the center region between the inlet and the
shoal. The local high region behind the shoal as mentioned above is enhanced considerably.
The waves at the location of inlet entrance are also steepened and the heights increased.
/ / .Il *5
 / /1 / 
'I I
 '. / a *
J J I
*~ ~ a' S 
0.22/0
S0. 222+00
MAXIMUM VECTOR
5F
Figure 4.9: Velocity Vector Field for Circular Shoal Located on the Plane Beach (no inlet)
__
20.00 '10.00oo o.00o b.oo 2b.00 3.00
V (CM/SEC)
Figure 4.10: Longshore Current Distributions at Lowest Boundary. No Inlet; with
Inlet
78
I I 1 
0 216E+00
"*s *** 111 
.  I i / y /,. ,
... .. / ,, / / / ... 4.. .. ., ,, 
.' '' '
** # i /
Current* *
/ 1\ 1 1 l
** \ \ i \ \
S i l
.
O.216E+00
SMRXIMUM VECTOR
Figure 4.11: Velocity Vector field for Circular Shoal Located on the Plane Beach with
Current from Inlet
Figure 4.12: Contour Plot of Wave Amplitudes on the Circular Shoal Located on the Plane
Beach (no inlet)
Figure 4.13: Contour Plot of Wave Amplitudes on the Circular Shoal Located on the Plane
Beach with Current from Inlet
CHAPTER 5
EXPERIMENTAL COMPARISONS
In order to test the validity of the numerical model, laboratory experiments were con
ducted in a square basin in the laboratory of the Coastal and Oceanographic Department,
University of Florida. The experimental set up, procedures and results together with nu
merical simulations are reported in this chapter.
5.1 Purpose of Experiments
Although there are many theoretical studied in recent years on the subject of wave
current interaction, laboratory experiments and field data are scarce. IsmailAwadalla
(1981) conducted laboratory experiments in wavecurrent interaction on a flat bottom with
a nozzle injecting dye against the normal incident waves. The results demonstrated the
physical phenomenon of currentwave interaction but were not sufficiently quantitative to
verify numerical models.
In the present experiments, we focus on a beachinlet system, where the phenomena of
wave transformation over a gentle slope, coupled with the effects of a nonuniform current
from an inlet can be studied. There are four basic quantities to be investigated: mean
velocities: U in x direction and V in y direction; the mean water surface displacements
i and the wave amplitudes A. Since Wr is a small quantity and requires a long time to
establish a steady value, the effects of wave reflection in the wave basin preclude meaningful
measurements. Only wave amplitude and current field were investigated.
82
5.2 Facilities and Instruments
5.2.1 The Laboratory Facilities
The experiments were conducted in a square 7.22 m by 7.22 m wave basin with a depth
of 0.8 m as shown in Fig.5.1 and 5.2. A fixed plane beach was constructed of cement with
1:13.3 slope (or s = 0.075). The toe of the beach was 0.85 m from the wave generator. The
inlet was located at the center of the beach end. It was 1.44 m wide and the bottom was
horizontal. The total length from the leading edge of the inlet to the diffuser was 1.2 m
as shown in Fig.5.2. For all cases tested the water depth was kept 0.369 m measured from
the bottom of the basin. Thus, the water depth was 0.045 m in the inlet. A grid system
was established as shown in Fig.5.3 to facilitate numerical computation and laboratory
measurement. The grid size is 0.12 m x 0.18 m. There are 43 grids in the xdirection and
40 in ydirection. The inlet was located at the center of basin from j = 17 to j = 24. The
leading edge of the inlet was at i = 37. The toe of the sloping beach was located at i = 1.
Waves were generated by a flap type wave generator. The bottom of the flap was
elevated 0.2 m above the bottom to permit the excess water from the inlet to return freely
under it. The wave generator was driven by a D.C. motor. Wave period and amplitude
were controlled by adjusting the frequency of the motor and the eccentricity of the linkage.
Uniformly long crested waves can be obtained for small amplitude waves. Cross waves
would appear as the wave amplitudes became larger. Standing waves also existed under
certain circumstance due to the short distance from the toe of the beach to the wave flap.
To prevent the leakage of wave energy from behind the wavemaker, horse hair mats were
placed between the wave paddle and the wall, so that the wave motion behind the paddle
was greatly reduced. On the other hand, due to space limitation, wave energy absorber
could not be placed at the beach end. Therefore, reflection of waves was serious at times.
The current system consisted of a pipe system, a flow filter, a diffuser and a reservoir.
As shown in Fig.5.1, a 2.5 horsepower centrifugal pump located between the flow filter and
the diffuser was used to transport the water. The flow filter situated behind the wave flap
I
83
C4
,. Ia
;
0.a
wC*
,,SI'L
CCI
E E:'
a.. 0 q
EE
U.
= I V
U. *
Figure 5.1: Plan View of Wave Tank
Flap Type
Wavemaker
7.22m
CrossSection AA
Screen
Trench
4*
Diffuser
F IIter
^
CrossSection BB
Figure 5.2: CrossSection Views of Wave Tank
T
0.8m
I 
Flow
~ T
Inlet
17
m
IDry
IWET
Figure 5.3: CrossSections of Measurement
U
II III l l llH\ xkl l~
 I I I : I IFl I l l l l l l l
Fl Iit111 il
III ~ ~ ~ ~ 1 12ill! lll
l l l IlTl l l l l
T
29 27 24
14 12
86
has ten control valves, allowing the water to flow naturally to the recirculation pipe without
disturbing natural flow condition. A diffuser with a length 1.4 m was installed in a trench
at the end of the inlet. It was a section of pipe with a series of different sized holes to
distribute the flow uniformly. The gap around the diffuser was filled in with pebbles and a
filter screen made by two layer a net filter with horse hair in between was set up in front
of the trench to promote uniform flow across the inlet. However, the distance between the
diffuser and the intersection of the shoreline and the inlet entrance was only 0.6 m. It was
very difficult to attain a uniform flow condition for such a short distance. After many fine
tunings, an acceptable uniform flow condition was finally reached. Fig.5.4 shows the velocity
distribution in the inlet at the crosssection that intersects the shoreline. The adopted final
velocity distribution was 7/25/86 profile.
5.2.2 Instrumentation
The main properties measured included wave characteristics and current. Four resis
tance type wave gages were used to measured the spatial distributions of wave heights and
periods. Current was measured by a MarshMcburney type electromagnetic current meter.
In addition to wave and current, the discharge from the inlet was monitored by a commer
cial flow meter manufactured by Sitnet Scientific. The wave gage was calibrated statically
by varying the gage submergence stepwise. It was found that the gage could not maintain
linearity for the full gage length. Therefore, gage calibration was finally done for each lo
cation prior to each test. The current meter was supplied with calibration but was further
calibrated, verified and adjusted by towing test in a long flume in the same laboratory.
The flow meter was also calibrated by measuring the discharge filling in a reservoir. Fig.5.5
shows the calibration curves of the three instruments; the results were consistent and stable.
A 6channel strip chart recorder was used to display the data. The wave gage signal was
first amplified before feeding into recorder. The current meter signal was timeintegrated
to eliminate short term fluctuations. The layout of the instrument arrangement is shown in
Fig.5.6
x 7/9/86
0 7/21/86
0 7/25/86
x
LO
o o
0o o x
0 0 x 0 x 9 0 0
0 0 o
x
0
0
O
0.2 H
Xo(m)
Figure 5.4: Velocity Distribution at I=Idry
0.8 
0.6 
O
0
0.4
I
0.0
0
88
FLOW METER CALIBRATION
2 4 6 8 10
FLOW METER RECORDING (ft/sec)
> 10
00
j
> 10
0
oo
j
.UJ 4
c 10
IJ
UJ
a: o
CURRENT METER
6 8 10
RECORDING (cm/sec)
WAVE GAGE CALIBRATION
Figure 5.5: Instrument Calibrations
CURRENT METER CALIBRATION
I I I I I
LEGEND
x X Direction
O X Direction
 O Y Direction
A Y Direction
I I I I
