Numerical modeling of current and wave interactions of an inlet-beach system

UFL/COEL-TR/073

NUMERICAL MODELING OF
CURRENT AND WAVE INTERACTIONS
OF AN INLET-BEACH SYSTEM

BY
YIXIN YAN

1987

NUMERICAL MODELING OF CURRENT AND WAVE INTERACTIONS OF AN
INLET-BEACH 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 co-adviser 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 Depth-Averaged Equation of Continuity ........ .......... 7

2.2.2 Depth-Averaged 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 NON-SLIP BOTTOM BOUNDARY CONDITION APPLIED IN WAVE EQUA-

TION ....................................... .......149

C NUMERICAL SCHEME FOR CONTINUITY-MOMENTUN EQUATION ... 152

D DOUBLE-SWEEP 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 Longuet-Higgen's Longshore Current Profiles . . . ... 41

2.6 Variations of Eddy Vicosity Coefficients. = 0. ; *** =0.001; -
=0.01; Von Khrmhn's Model with upper-limit =0.001. unit= m2/s 43

3.1 Grid Scheme ................... ................ 48

3.2 Definitions ofXj and (Y,). ....................... 49

4.1 Set-up 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 Cross-Section Views of Wave Tank . . . . .

5.3 Cross-Sections 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 x-Component 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 x-Component 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 x-Component 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 x-Component 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 x-Component 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 x-Component 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
INLET-BEACH 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 inlet-beach

system is developed in this study. The theoretical formulation is based on the time- and

depth-averaged equations of continuity and momentum, and a 3rd-order wave energy equa-

tion that retains the important nearshore effects of refraction, diffraction, current-wave

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 3rd-order 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 Crank-Nicolson finite difference scheme is adopted for the model. In each spatial

sweeping, a recently developed two time-level 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 man-made 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 inlet-beach system.

The hydrographic and shoreline changes of an inlet-beach 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 inlet-beach 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

non-uniform 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

non-linear 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

non-stationarity, 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

non-uniform 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 semi-infinite 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 beach-structure 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 current-wave

interaction problem and proposed a weekly non-linear 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

current-wave 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 wave-induced 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 Depth-Averaged 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 i--h + 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 f-h

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 Depth-Averaged Equations of Momentum

The horizontal depth-averaged momentum equations used here were originally derived

by Ebersole and Dalrymple(1979). A brief summary was given here. The Navier-Stokes

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
I-h Ot a- t J.h at

P Ou, u, 8 8s7 Oh
7) L dz = -__ uiuYdz uiu I 97 uiuj 1-h

i.ii f_7 9uiw
S-) -i-W dz = uiW 1 j -ut, I-h

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 I-h (j l-h + W I-h)

Recognizing the last two brackets of LHS are the KFSBC and BBC, respectively, we have

LHS = /uidz + -a uiujdz
9t T-h x9j f-h
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) "f-h 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 J-h

a UD)+ a (Ui'iD)+ a(UfiD)+ I a uFu.0dz +

a r
ax, J-h'

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 J-h 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 I-h
P daz J-h 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 I-h consists of two parts: hydrostatic pressure

and wave induced dynamic pressure, viz.

PI-h = -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 + - I-h
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 I-h + -e,
p xxi J-h ,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 I-h + -1.i + -T,, uiutdz (2.18)
2 p 8zi p p 9x J-h

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(pu-U.) +
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

X-Direction

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

Y-Direction

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 wave-like 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 x-direction which is true when the angle of wave propagating is small with

respect to the x-axis, i.e., it is assumed that no fast wave-like 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 y-direction. The problem is further considered to be

quasi-steady; 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 = C-2,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 C-2 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+ f-A- 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 + e2Uf-Am + eC2Vf A-Y
CJlJlr-- "x + 12z.!Zzt2 "2 + 2U Ax an1 1 ,Am22VAfiArnyl

+e2 2flVh, VhA + eC2A-Vh, 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

I|U.

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 well-known 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 + 3-lfiArA7 + 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 n-mv n-n
+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

S--f 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 = -ig--e (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 ae-i (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 nA-xx 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 + f1A1A-1 + fl '21lt A + t~l1A1t] + ff5[t71AXA
17 = A IA + AA1 1+ AlA t]+ iA

+rj1lAxl)}

Here the condition of fi = fi-1 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 set-down 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 non-zero. 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,e-i (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,

-f-Lialx = 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 + flf2-A2Ax + fL' 1-D
1 2 1 21 -D1
+ft0 + flz2r + fIfz(2AxAx + A127-1) + 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 + .a-I 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 wave-current 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(ko-fz kdz) (2.91)

where ko is the reference wave number. Equation (2.91) implies that

A" = a" ei(k'x-koz)

(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 one-dimensional 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 semi-nonlinear 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 non-dimensional 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 f-h S 2

I

36
The form was simplified by Longuet-Higgins 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 = P-f (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 Longuet-Higgins (1970) suggested a value about 0.01

for long-shore 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, semi-empirical 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 (Longuet-Higgins, 1970 a.b). The common

assumption is that turbulent mixing is proportional to the local mean velocity gradient. Fol-

lowing Longuet-Higgins' 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, Longuet-Higgins proposed a model for erw, that is proportional to the distance

offshore, z, multiplied by velocity scale

Ce = Nzx-g

(2.124)

0

SWithout Mixing
With Mixing
I-
-a

z
-JJ

Breaker Line
DISTANCE OFFSHORE

Figure 2.5: Longuet-Higgen's Longshore Current Profiles

where N is a dimensionless constant whose limits Longuet-Higgins 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 upper-limit 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 upper-limit 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)

Figul-e 2.6: Variations of Eddy Vicosity Coefficients. = 0. ; ... =0.001; =0.01; -
Von Khrman's Model with upper-limit =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 current-wave-topographic 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:

Depth-Integrated Continuity Equation

aB aUD 8VD
+ + --- 0 (2.131)
at 8x ay

Depth-Integrated Momentum Equation in X-Direction

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

Depth-Integrated Momentum Equation in Y-Direction

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 = Nzvg-D 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 grid-center variable and the grid-side vari-

able. Practically all the quantities of interest are represented by the grid-center form which

expresses the value at the center of the grid and has the following notation:

where X is the grid-center variable, the subscripts i, j represent the ith grid side in the

x-direction and the jth grid side in the y-direction, respectively, and the superscript n

represents the nth time step. The grid-side 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 grid-center 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 x-direction on the left of the vector and j denotes the jth grid point in

the y-direction on the right of the vector.

In the present study, the grid-side velocity vector was applied in the continuity and

momentum equations and a grid-center velocity vector was applied in the wave equation.

The definitions of grid-side and grid-center variables are shown in Fig.3.2

3.2 Finite Difference Scheme for Combined Continuity- Momentum Equation

The Crank-Nicolson 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

M-1

M-2

I

AX

Y
y

YI,I+1

---
Y

Figure 3.2: Definition of X.. and (Y,),

I

Yij

aQv a-9r
a = -=gD (3.3)
at ay
where Q, = UD and Qy = VD. U and V in this section denote the grid-side 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)

x-sweep
S- + z(Un+1D) + St (VD) = 0

(3.15)
(1+ Atc))Un"+ U + g6* + Atmi = 0(
y-sweep

+ -i +s 6(Vn+D) ,(VD)= 0
(3.16)
(1 + Aty)V"+1 V + A-gY6,"+ + 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 up-dated

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

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,j-1)] = 0 (3.17)

X-Direction Momentum Equation

I-(-" t- U.-) + (U,+l Ui-,,) + (Vi+1 + V,.1 + Vi-,i+ +- Vi-y)
At 2Ax 8Ay
(Uj+1 Ui,-1) = ,i -*,) Di-,) z ) ( -

~p(Di, Di-lj) 2Ay [(Sy)ij+l (Say)i,i-1 + (Syv)i-lj+l (Sy)i--,j-1]

(Di, + Di-,i)1 (ut)ij ,i ij + (FLy)id + (FLys),i (3.18)
where

I (ut)i, I, = [(um)i,j + (um)i-1,,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, + ( ,J-1

+(my)i-1,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)i-l,,j+) + 2((Ecz)i,j+l + (eez)i-1,+l)] (Vi,. V-i,,j)[((ews)i,i

+(Cew)i-_,i + (ew);,,-_ + (Ews)i-,._1) + 2((e~),i + (Ecc)i-l)]}

b) Y-Sweep
Continuity Equation

1
( )+^ --) [ v1 [ (Dj+1i + 2,,), V^(Di,, + +,,-i)]

[ayVj+i(D.j+l + Dj)* Vj(D.,1 + .,j-i).] = 0 (3.19)

Y-Direction 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,j-1) = i iD 1(D,-, ,). 1(Sy)i+,
AY2 (Dj Di-1)
-(SOy)i-lj + (S.y);+1,,-1 (S.A),i,j-1i 2 ((Sy)i,
P(Dii + Dij-l1)* 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,. Vi-l,4)[((E~w),,, + (Cw,)i-l,j

+(cw.),,j- + (eCw).-1,i-1) + 2((e~.),, + (ecs).-i)]}

4AxA1 l
(FLy)iJ = 4AxAy{(U +: U+ -)[((Cwu)i+l + (uW)i+l,j-1 + (Ey)i,j
+( ) ) + 2(()+,, + ( ,c)i+ (- 1)[((

+( )i,-1 + ( ).i-1,j + (cw),-1,f-1) + 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 up-dated 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 x-sweep and
y-sweep. 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

X-Sweep

f -j + Pi ,i UiRj Pi-j1,iU = (3.21)
(3.21)
Aij U," + RI V+1,j RI jij = Bij
Y-Sweep

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 = F2i-1 + 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 = idry-1, is the last grid at

56
which water depth has non-zero 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 x-direction.

The lateral boundary conditions basically have three kinds for the concerned region:

the symmetric free boundary conditions; the no-flow 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 = F4i-I + E4- l,1
where

Zij QijF4j
F3i =
1 + Qi,,E4j

E3j = Qij-1
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 parabolic-hyperbolic

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 Crank-Nicolson 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 --gk-2
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 grid-center

values in this calculation.

The Crank-Nicolson 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 -j-I)] [(C)ji+1,,(A.i+j+1 2Ai+1,, + Ai+1,-1) + (Cl)j,j

1 ((C2),+, (C21,s) (A1,
(Ai~+1 2Aj + Aj-I)] + ((2), (A+i+ + Ai,)
2Ax ((C2)i,+1 + (C2)idj)
1 ((C3)+1-s),+ (C(3),j-1)
+ 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 = R2j-1 + R1, Ai4j-1 (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)
1-Q
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(Ai-la+) = 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 (on-offshore) with Ax = 0.25m and n = 4

in the y-direction 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: Set-up on a Plane Beach. Eq.(2.138); Linear Wave; A Data
o (
o-

Figure 4.1: Set-up 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 out-flow from

downstream would not be balanced by the in-flow 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 open-boundary-condition 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 x-direction with Ax = 5.0 m and n=27 in y-direction 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 4-c,)/ 4-4 -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 y-axis.

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 non-dimensional 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 fore-edge 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 mid-section 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 10-4 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. Ismail-Awadalla

(1981) conducted laboratory experiments in wave-current interaction on a flat bottom with

a nozzle injecting dye against the normal incident waves. The results demonstrated the

physical phenomenon of current-wave interaction but were not sufficiently quantitative to

verify numerical models.

In the present experiments, we focus on a beach-inlet system, where the phenomena of

wave transformation over a gentle slope, coupled with the effects of a non-uniform 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 x-direction and

40 in y-direction. 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

Cross-Section A-A

Screen

Trench
4-*
Diffuser

F- IIter
^-

Cross-Section B-B

Figure 5.2: Cross-Section Views of Wave Tank

T
0.8m
I -

Flow

~ T

Inlet
17

m
IDry
IWET

Figure 5.3: Cross-Sections of Measurement

U

II III l l llH\ xkl l~

- I I I : I IFl I l l l l l l l
F-l 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 cross-section 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 Marsh-Mcburney type electro-magnetic 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 6-channel 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 time-integrated

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

MILLISECOND	CLASS.METHOD	MESSAGE
0	sobekcm_page_globals.constructor
0	sobekcm_page_globals.constructor	Application State validated or built
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.constructor	Navigation Object created from URI query string
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.display_item	Retrieving item or group information
0	sobekcm_page_globals.get_entire_collection_hierarchy	Retrieving hierarchy information
0	sobekcm_assistant.get_entire_collection_hierarchy
0	cached_data_manager.retrieve_item_aggregation
0	cached_data_manager.retrieve_item_aggregation	Found item aggregation on local cache
0	item_aggregation_builder.get_item_aggregation	Found 'all' item aggregation in cache
0	system.web.ui.page.page_load (ufdc.page_load)
0	sobekcm_page_globals.constructor.on_page_load
0	html_echo_mainwriter.add_style_references	Adding style references to HTML
0	html_echo_mainwriter.add_text_to_page	Reading the text from the file and echoing back to the output stream
37	html_echo_mainwriter.add_text_to_page	Finished reading and writing the file