|
UFL/COEL-89/014
WAVE-CURRENT INTERACTION OVER A
SUBMERGED BAR FIELD
By
Thomas Richard McSherry
1989
Thesis
WAVE-CURRENT INTERACTION OVER A SUBMERGED BAR FIELD
By
THOMAS RICHARD MCSHERRY
A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF MASTER OF SCIENCE
UNIVERSITY OF FLORIDA
1989
ACKNOWLEDGEMENTS
In the course of completing this thesis I have received abundant help in many forms
from a variety of people. The work required the extensive help of my Advisor, Dr. James
Kirby, and I am truly grateful for his assistance. I consider myself very fortunate to have
had the opportunity to study under Dr. Kirby. I also used the COEL extensively on two
separate occasions, and depended on the help of several people there. I thank Jim Joiner
and Sydney Schofield for their patient assistance, and the inputs from the entire group. A
great source of help at the COEL was in the way of mental release in the form of the daily
volleyball game, and although I jammed a few fingers in the less-than-benevolent net play,
I am grateful for the enthusiastic games.
I also thank Dr. Robert Dean for his help in every way, and I thank the secretaries who
have told me the logistical requirements for getting out from under the deadline hammer.
Also, thanks to the students that offered friendship and release, and provided help when I
really needed it. These people I wish the greatest of success.
TABLE OF CONTENTS
ACKNOWLEDGEMENTS ...........
LIST OF FIGURES ...............
ABSTRACT ..............
CHAPTERS
1 INTRODUCTION ..........
1.1 Problem Statement .......
1.2 Literature Review .......
2 THEORETICAL REVIEW ....
2.1 Introduction ..........
3 GOVERNING EQUATIONS .. .
3.1 Introduction ..........
3.2 Circulation Model ......
3.2.1 Governing Equations .
3.2.2 Radiation Stress Terms
3.3 Wave Model ..........
3.3.1 Governing Equation
4 FINITE DIFFERENCING ... .
4.1 Introduction ........ .
4.2 Circulation Model .......
4.2.1 Method of Solution .
4.2.2
4.2.3
Boundary Conditions .
Radiation Stress Subroutine
. . . . . . . .
4.3 Wave Model ...... ........ ........................ 35
4.3.1 Method of Solution .................... ........ ... .. 35
4.3.2 Finite Differencing .. ......... ................ .. 35
4.3.3 Boundary Conditions ......................... 37
4.3.4 Lateral Smoothing ..... ......... ........... .. 37
5 LABORATORY WORK ...... .... .... .............. 48
5.1 Introduction .. ... .... ...... ..... .. ... ... ... ...... 48
5.2 Equipment .... ...... ... .................. .. 48
5.2.1 Basin .............. ........... ........ 48
5.2.2 W avemaker .... ............ ............... 50
5.2.3 Beach and Cart System ........................ .. 50
5.2.4 Electronics .. . ........................ 52
5.2.5 Bar Field .................... .. .......... .. 54
5.3 Setup, Procedure and Data Analysis . . . . ... 55
5.3.1 Session One ....... ........... ............. .. 56
5.3.2 Session Two ... ... ... ... ... .. .. ... ... .... .. 62
6 RESULTS AND CONCLUSIONS .......................... 76
6.1 Introduction .... ........ ..... ....... ........... 76
6.2 Laboratory Work ................. .. ..... ........ .. 76
6.2.1 Reflection Results ..... ........ ............ .. 76
6.2.2 Bathymetry ...... ......... ..... ...... .... 86
6.3 Numerical Results .......... .. ................ .. 87
6.4 Conclusions and Status Report ........................ 92
BIBLIOGRAPHY ............... .................. 100
BIOGRAPHICAL SKETCH .. ............................ 102
LIST OF FIGURES
2.1 Bar field in the presence of a depth-uniform current. From Kirby (1988). 6
2.2 Contour of Tr v. IFx v. w/Wro. ....................... 9
2.3 Surface projection of transmission Tr v. frequency and current. . 10
4.1 Flow chart for model MCSERRY. . . . ..... 25
4.2 Hyperbolic tangent startup function. . . . . ... 26
4.3 Grid for circulation model .......................... 28
4.4 Differencing of radiation stress gradients in x-sweep . . .... 30
4.5 Differencing of advective acceleration in x-sweep. . . ... 32
4.6 Wave calculated from parabolic equation on domain of width W. 39
4.7 Spectrum of A(y) without smoothing. Cutoff at nA = 9.051. . 41
4.8 Spectrum of y-direction driving force before smoothing. . ... 42
4.9 Spectrum of y-direction driving force after smoothing. . ... 43
4.10 IAI without smoothing. Bars are located on the bottom of this domain. 44
4.11 IAI after smoothing. Bars are located on the bottom of this domain... 45
4.12 y-direction forcing without smoothing. . . . ... 46
4.13 y-direction forcing after smoothing. . . . ........ 47
5.1 Plan view of large basin. ......................... .. 49
5.2 Plan view of cart track system. ........................... 51
5.3 Side view of profiling cart. ......................... 52
5.4 Sample profile. ....... ... ........ . . ... 53
5.5 Sample calibration curve for wave gage. .. . . . ... 54
5.6 Bar frame design. ............................... 55
5.7 Gage locations for two-gage reflection method. . . ... 57
5.8 Reconstructed bar profile for Session One. . . . ... 60
5.9 Plan view of barfield for Session One. . . . ..... 61
5.10 Gage locations for Session One. . . . ..... 62
5.11 Three gages inline with x-direction waves . . . ..... 63
5.12 Sample envelope, Session Two. Retrieved digitally from one of the three
on-offshore transects using the moving cart gage. . . .... 68
5.13 Sample envelope height versus X. . . . ..... 69
5.14 Sample Kr with offshore distance x from Session Two. . ... 69
5.15 Plot of Kr with transect location, Session Two. . . ... 70
5.16 Plan view of barfield, showing data collection locations. . ... 71
5.17 Sample time series of magnitude of current meter data, Session Two. 73
5.18 Profile for Session Two ............................ 74
5.19 Instrument setup for Session Two. . . . ..... 75
6.1 Plot of K, with wave period, Session One data and theoretical prediction
from 1-D model Kirby (1987) ......................... 78
6.2 Plot of Kr with period, data and theoretical for the case of no discernable
rip current. ....................... ............... 81
6.3 Plot of Kr with period for the case of a discernable shore rip current.
Data and theoretical .............................. 83
6.4 Two profiles, dotted line referring to profile used in the second and third
sets of testing, where a weak current formed. . . . ... 84
6.5 Plot of Kr with period for case of discernable shore rip current, moving
toward resonance ................... ........... 85
6.6 Bathymetry 1. Shoreward direction with increasing . .... 88
6.7 |A[ for initial amplitude of .03 meters and perfectly resonant waves. .. 89
6.8 IBI for initial amplitude of .03 meters and perfectly resonant waves,
without currents.......... ........... ......... 90
6.9 Steady state mean surface, ....................... 91
6.10 No-current wave model, incident amplitude of .02 meters. K, = 0.81. 93
6.11 No-current wave model, incident amplitude of .03 meters. Kr = 0.81. 94
6.12 No-current wave model, incident amplitude of .04 meters. Kr = 0.81. 95
6.13 No-current wave model, incident amplitude of .03 meters. K, = 0.81,
and bar half-length is 3 meters ....................... 96
6.14 No-current wave model, incident amplitude of .03 meters. Kr = 0.843,
and bar half-length is 2.25 meters. . . . ..... 97
6.15 No-current wave model, incident amplitude of .03 meters. Kr = 0.88,
and bar half-length is 1.75 meters. . . . ..... 98
6.16 No-current, tapered bars with .03 meter incident wave. Kr = .774... 99
Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science
WAVE-CURRENT INTERACTION OVER A SUBMERGED BAR FIELD
By
THOMAS RICHARD MCSHERRY
December 1989
Chairman: Dr. James T. Kirby
Major Department: Coastal and Oceanographic Engineering
People are moving to the coast in great numbers that seem to increase each year. Devel-
opment along the shoreline is therefore also on the uprise, and the public often assumes that
the shore environment is stable. Recognizing the appeal of the beach, and understanding
that this environment is anything but stable, the coastal engineer must attempt to curb the
natural destructive forces from the sea that would normally hinder man's use of the beach
and it's surroundings.
One present goal in this attempt is to achieve shore stabilization with a minimum
impact on the dynamic processes in the nearshore region. One candidate receiving much
recent investigation is the shore parallel low profile bar field. The bar field acts with the
incoming wave to actually reflect energy back out to sea. This action thus protects the
beach from some of the destructive wave energy that often erodes massive amounts of dune
material in heavy weather. If the concept could be perfected, the low profile bar field could
conceivably be a huge shield against any wave possible.
The way the bar field reflects energy is by having bar spacings that are one half the
incident wave wavelength. With this satisfied, the incident energy serves to form another
wave that propagates opposite to that of the incident wave. Knowing that the bar field
would be of finite longshore length, a low spot in energy would then exist shoreward of
the bar field. What would then be expected is to have a circulation begin by this energy
differential, much in the same manner as circulation into the lee of breakwater, or along
groins. It is also expected that the current would flow outward through the middle of the bar
field. Clearly, then, the current would alter the character of the incident wave, changing it's
wavelength. Then the resonant condition between the incident wave and the immoveable
bars would be upset. If the bar field is to considered for practical implementation, the
behaviour of the rip current driven by the energy differential will need to be investigated.
The most important aspect of the wave-driven current is how it responds to slight
shifts off resonant incident wave frequency. The main hypothesis of this paper is that
the bar-wave system seeks resonance, and therefore if the wave were to be initially out of
resonance a current would be required to tune the wave. By incrementally lowering the
initial wave frequency, and observing how the system responds and whether reflection is
enhanced or reduced, the hypothesis is proven or repudiated. Another factor in this scheme
is the available energy from the incident wave. If the waves are not big enough to create
sufficient energy differentials, then the current will never reach the strength required to
tune the waves that are out of resonance. The thesis thus looks at how current strength
responds to wave heights, and how reflection behaves during the process.
This paper studies results from laboratory data taken from a finite length bar field, and
also findings from a numerical model that incorporates coupled wave equations with the
momentum equations. The numerical scheme iterates between the two parts until a final
steady state is achieved. Mean values in surface elevation, lateral depth integrated currents
and wave amplitudes are solved for through a semi-implicit algorithm.
CHAPTER 1
INTRODUCTION
1.1 Problem Statement
The need for effective shore protection seems to be increasing with each year. The
worldwide dependence upon the coastline coupled with such chronic problems as rising sea
level has placed a responsibility upon coastal engineers to investigate truly effective methods
of shore stabilization. Effective means that it must reduce erosion, hold up under severe
conditions, be economically feasible, be environmentally sound, and not impede the regular
use of the coastline. Historically, structures meant to curb erosion have had some success,
but the side effects caused by the structures were oftentimes worse than the original problem.
Today, the main thought is to develop a device that works elegantly and efficiently. One
of the promising candidates is the low profile bar field that reflects incident wave energy
by use of a phenomenon now known in water wave theory as Bragg reflection, which was
named after the analogous process in crystallography.
By placing bars along the sea bottom parallel to the shore, and making their spacing
one half of the incident wavelength, a resonant condition occurs between these bars and
the surface waves. Specifically, a reflected wave is generated through the bottom boundary
condition. This reflected wave grows in amplitude with the number of bars placed until
theoretically complete reflection can occur if there were an infinite number of bars. The
findings to date have been very promising, with upwards of 80% reflection for waves travel-
ling over four sinusoidal bars. Previous laboratory data also confirm the expectations, and
these encouraging findings are fostering continued research into the low profile bar field.
These bars would necessarily be of finite longshore length, and if part of the incoming
2
energy were reflected, a low spot in energy would exist shoreward of the bar field, with
full amounts at the sides of the bars. This differential in the mean surface elevation along
the coast could possibly generate a large enough circulation to create a fairly strong rip
through the barfield and back offshore. What this means is that the incident wavefield
will be distorted over the barfield, and the resonant condition will be upset. An important
question that has been answered by Mei (1985), is how far off resonance the waves can be
before the bars quit reflecting energy. With the introduction of a current, however, the
question is whether reflection will be enhanced or reduced.
This thesis focuses on the wave-current interaction over the bar field, and what happens
when initial wave frequency is incrementally reduced. What the thesis tries to answer is
what happens to reflection given the current, and how the current changes as wave height
changes. Also, the thesis tries to prove that the wave-bar system seeks a resonant condition
by tuning off- resonant waves closer to resonance through the rip current. These questions
are answered by means of laboratory data and a numerical model. Some theoretical review
introduces the main content, and the problem is fully developed. Then the equations
necessary to the problem are presented, followed by a numerical solution to the equations.
After this, the laboratory work is covered, and the methods of data analysis are reviewed.
Finally, the results from the laboratory work and the numerical model are given, through
which the hypothesis will be proven or repudiated. Before all of this, however, the past
research that this thesis depends upon will be reviewed.
1.2 Literature Review
This paper is a direct extension of various avenues of previous research. The radiation
stress concept developed by Longuet-Higgins and Stewart (1964), is obviously important
in any study of waves and currents. These concepts were extended to terms involving
complex surface amplitude r' by Mei (1972). The expressions in Mei's book (1983) are used
here directly,and are the driving forces in the time averaged, depth integrated momentum
3
equations that govern the circulation and mean water surface.
The other terms appearing in the momentum equations have been included in the
model, and most of those formulations are based on previous work. The bottom friction
as developed by Longuet-Higgins (1970a) is represented in the model as a function of the
linear water wave maximum orbital velocity and the mean current U. In a companion
paper, Longuet-Higgins (1970b) looked at the transition of longshore current at the breaker
line, and developed an expression for the lateral mixing. This term also is included in the
present model.
The advective acceleration terms in the momentum equations are included in the model,
and are formulated as a four grid-point average of the current values, which is done by
Winer (1988). The momentum equations are solved by an Alternating Direction Implicit
(ADI) method which was introduced by Sheng and Butler (1982). The matrix set up for
this method can be seen also in Winer (1988). Although the matrix implicitly refers to
linear terms, the model introduces limitations in time step by including the nonlinear terms
explicitly. Iteration between wave model and circulation model was used by Noda et. al.
(1974), where the importance of wave current interaction in determining the final current
field was noted.
There has been much work on Bragg reflection in the last decade, as elegant and en-
vironmentally sound methods are sought for shore stabilization. A non-resonant theory
was developed by Davies and Heathershaw (1984), and partially successful comparisons to
laboratory data were made. The theory broke down at resonant frequencies, however, and
the reflection coefficients predicted were significantly larger than the data. Mei (1985) ex-
tended the theory to include resonant frequencies, and successfully predicted the reflections
seen in previous experimentation. Kirby (1986a) allowed for the presence of a mild slope in
the theory. The theory was further extended to include currents flowing over the bar field,
although on a flat bottom by Kirby (1988). Kirby (1986b) attempted to include currents
over a mild slope without the bar field. The work on flat topography noted that the pres-
4
ence of a current enhances reflection, and presented the background that this paper uses
as a starting point. The presence of a current allows incident waves of frequencies slightly
off resonance to tune with the bar field and actually achieve greater reflections than with
perfectly resonant waves. Transmission coefficients always drop in the presence of a current.
Thus with the joint action of a depressed energy zone driving a rip current that can alter the
incident wavelength, and transmission reduced against the current which further reduces
the low energy zone, the system may tend to want to go toward perfect tuning. The scheme
depends largely, however, on the incoming available energy to drive the current. Laboratory
data, will demonstrate cases where incident wave heights were too small to sufficiently create
a rip current strong enough to keep the reflective plateau.
Various methods are used to measure reflection from the laboratory data, and the
previous work that made these methods available will be cited in Chapter Five.
CHAPTER 2
THEORETICAL REVIEW
2.1 Introduction
To begin the task of answering some of the questions raised in the last chapter, specific
wave theory should be reviewed. It has been conjectured already that reflection will increase
when a current flows over the bar field against the waves.. The thesis will seek to prove this,
and the more important concept that the entire system will try to reach a resonant condition.
The mathematical basis for these hypotheses was developed by Kirby (1988).
Referring to the following figure 2.1 for the case of a current over a sinusoidal bar field
of length L, with bar wavenumber A = 27r/Lb, Kirby showed that if one considers incident
waves at or very near resonance, the ratio of transmission with a current to transmission
without a current is
SA(L)current cosh QL (2.1
T A(L)nocurrent coshQcL
The subscripts c refer to properties in the presence of a current. The bar field has length
L, where A(L) is the value of the transmitted wave amplitude at the back edge of the bars.
The lateral bar shape is described as Db. The other terms are defined as
Q = Q (2.2)
c,
Qc = (2.3)
(C1C2)
gk2Db
o gkDb (2.4)
4wcosh2kh
Figure 2.1: Bar field in the presence of a depth-uniform current. From Kirby (1988).
( U 0)2 (2.5)
Q 2,,2 )
where C1 = Cgl + U and C2 = Cg2 U. Qc is defined later by (3.24). Note that when a
strong enough current opposes either wave, the wave can be stopped when Cg. + U = 0.
To simplify the expressions and to show what happens when the current strength changes,
consider a shallow water wave where the expressions reduce to the following;
ADb
Q--
8h
ADb 2-l
Qc 8 -, F)-
where the Froude number in the x-direction is F, = U/vg-h. It appears that Qc has a
singularity as IFxI -- 1, which corresponds to the stopping Froude number. When IFx\ < 1,
from the above simplifications it is seen that Qc > Q, and therefore TA is less than unity.
What this means is that in the presence of a current, transmission for waves which are
perfectly resonant with the bars with an opposing current is less than the transmission
X
-~ x
---- U
x=L
7
for perfectly resonant waves without a current. The current is thus perceived to enhance
reflection, and it is interesting to note that this is independent of the direction of the current
along the x-axis.
Now for the shallow water limit consider resonant frequencies, where wr is in the pres-
ence of a current, and wro is without an imposed current. These are
Wr = )A(1-F (2.6)
1ro = IVghA (2.7)
which from the ratio wr/wro = 1 Fx describes a parabola. Refer to figure 5 in Kirby
(1988) for a picture of this ratio for various values of Ah, where Ah = 0 is also shown as the
shallow water limit. Now consider any frequency w shifted off the resonant frequency wr by
an amount 0, so that w = Wr + There is a cutoff condition for this shift, above which
for a given current, the resonance is lost between the wave and the bar field. This cutoff is
described by
S2 ; = ffg- C= h. (2.8)
cutoff = 64C2
The behaviour of the wave field depends'on the relation between Q and Qcutoff. Mei (1985)
and Kirby (1988) denoted four cases where the solutions to the coupled equations described
different wave fields, and these solutions for A(x) and B(x) can be put in terms of the
transmission coefficients, Tr = IA(Lb)/A(O)I. These expressions for the four cases are
case 1: 0 > Qcutoff
PCI
T, = 1 (2.9)
[(PC1)2cos2PLb + (S')2sin2PLb]/2
case 2: t = 2,cutoff
C1
T,= (2.10)
[C2 + (fLb)211/2
case 3: Q < cutoff
T =QC (2.11)
[(QC1)2cosh2QLb + (QI)2sinh2QLb]1/2
case 4: Q = 0 (perfect tuning)
1
Tr = (2.12)
ScoshQLb
These relations may be evaluated close to shore using the approximations:
Ci = vgh(l + F)
C2 = Vh(-F)
gADb 1
8^ _( F.2
1- Fx
p2 = 2[ 02 2 f
r ~ \2gh(l F)
Q2 = -P2 when Q < ,cutoff
These expressions can be put into a program to plot the surface projection of transmission
Tr for any combination of w/wro and IFxI. Figure 2.2 shows the contour of the transmission
coefficient, and figure 2.3 is a surface projection plot of the same result.
Considering an incident frequency that is not far off wro, say 85% of it's value, then
by no-current resonant theory the transmission coefficient will be 0.8. The mean surface
elevation in the surfzone shoreward of the bar field would then be lower than surrounding
areas, and the differential could be enough to drive a current. By following this line of
constant wl/ro along the surface projection, clearly as the current increases the line heads
9
0.8
JR FROM 1.5B000E-1 TO 1.8101 CONTOUR INTERVAL OF 8.S1100E-11 PT13.31- 1.99934 LABELS SCALED BY 196B.
IF.1
Figure 2.2: Contour of T, v. IFI v. /Wro.
0.5
Figure 2.3: Surface projection of transmission Tr v. frequency and current.
11
down the bank of the large trench toward the perfect resonance ratio Wr/wro, which again is
illustrated in figure 5 of Kirby (1988). As this minimum transmission coefficient is reached,
any more current would drive the transmission coefficient up the other side of the trench.
This will decrease the energy differential and reduce the current, allowing the system to fall
into the base of the trench again.
This assumes that there is enough available energy from the wave to drive the currents
necessary to reach the trench base, but clearly shows the process that the system seeks
an equilibrium value that is as close to the perfect resonance ratio as it can get. This
is the main point of this thesis, and Chapter Six will present results which are aimed at
investigating this hypothesis.
CHAPTER 3
GOVERNING EQUATIONS
3.1 Introduction
This chapter addresses the governing equations that apply for the case of water waves
over mild topography, and the resulting mean currents and mean surface elevation. The
problem as stated must be simplified to be solved numerically, and therefore the equations
governing momentum balance in the domain between the mean values must be reduced to
two dimensions. Also, the wave equations that govern wave amplitude over the bathymetry
are presented as a parabolic equation. These implications inherently take away from the
full process, but without them the solution would be very difficult to achieve. First to be
presented in this chapter are the equations and terms that apply to the circulation. Then
the equations for the wave field will be presented.
3.2 Circulation Model
3.2.1 Governing Equations
The governing equations are presented in equations (3.1)-(3.2). They are the depth-
integrated and time averaged equations of motion, where time averaging is done over a wave
period. The derivation can be found in a number of references, most recently in Appendix
A of Winer (1988).
aU aU oU 0- 1 1
+ U- + V +g + b- s-
ft 9x 9y Dx pD pD
13
1 (dSId BS 1 drj
+ X +- - + 0 (3.1)
pD ozx dy p dy
OV OV dV Or 1 1
+ U + gV + + f y bi7s
9t 9x 9y 9y pD pD Y
+_1 + SYY) + o 0 (3.2)
pD \ Ox Oy px 0
The continuity equation completes the needed three equations for the three unknowns ?, U
and V.
+ (U D) + (V D) = (3.3)
Tt Ox By
These equations are simplified to the spatially two-dimensional case to render a realistically
solvable set of equations. The assumption does, however, take away from the complete
process which might include strong effects from vertically variable current fields. In these
equations, the terms are defined as
U = x component of mean current
V = y component of mean current
S= mean water surface elevation
p = vertically constant water density
ho = still water depth
D = total water depth = ho + 7
rl = lateral shear stress
Tbx = x component of bottom shear stress
Tby = y component of bottom shear stress
r,, = x component of surface shear stress
rsy = y component of surface shear stress
14
Sxx = x component of flux due to x-propagating wave
Sxy = y component of flux due to x-propagating wave
Syy = y component of flux due to y-propagating wave
Note that the equations will retain the nonlinear advective acceleration, gradients in the
radiation stresses, bottom friction, and lateral mixing. The numerical representation of
these terms are the same in this thesis as was presented in Winer (1988). Referring to that
dissertation, the bottom friction term is represented by equations (2.15) and (2.16), which
is a direct application of Longuet-Higgins (1970a). The lateral mixing term is in equations
(2.26) and (2.27), and also in Longuet-Higgins (1970b). The advective acceleration is rep-
resented in finite difference form in equations (4.24) and (4.25). The surface shear stress
is not included in the paper, but could be added easily enough. It was not included for
it remains out of the scope of the present problem dealing with energy differentials due to
reflection.
3.2.2 Radiation Stress Terms
The original model by Winer (1988) did not address the reflected wave that is now
present. Since a reflected wave will be expected, the radiation stress terms must be altered
to include this wave to maintain a proper balance. Also, when considering waves on a
current, the incident and reflected wavenumbers will be different. If one looks at how the
waves interact with each other, there will be terms that will oscillate spatially at a spatial
frequency related to the difference of these wavenumbers. Therefore, deriving radiation
stresses would necessarily be involved. For demonstration of what types of terms are present
in a two-wave system at zero angle of propagation, consider the situation of zero current
and thus equal wavenumbers. To begin, the linear velocity potential for a two-wave system
is given by
(xytz,) ig coshk(ho + z) A(x,y) ei(kx-wt)
2 w cosh kho
ig coshk(h + z)B(, y) ei(-k-wt)+ +
2 w cosh kh,
(3.4)
and the instantaneous surface displacement 4' is given by
'(x,y,t) = 1 [A(x,y) eke + B(x,y)e-ikx] e-'t + *
,'(xY,t) = -
(3.5)
where denotes the complex conjugate. It has been assumed that the complex amplitudes
relate to a scaling parameter c defined as E = k lA, the steepness of the wave. The amplitudes
vary with x at O(e2), and in y at O(c). In keeping with the parabolic simplification to the
wave model, the wave angles are then assumed to be at or very near zero. The wave induced
velocities are then, by (3.6-8)
u'- 0 -
0 _
.-=^ =
Oy
w' =- z
z9
g cosh k(h, + z) kA ei(kx-wt) i Axei(kx-wt)
2 w cosh kho I
+ kA*e-i(kx-wt) + i A e-i(kx-wt) kBei(-kx-wt)
iBxei(-k-t) kB*e-i(-kx-wt) + iBe-(-kx-w) (3.6)
g cosh k(ho + z) [i A, ei(kx-wt) + i A e-i(kx-wt)
- ~ e(2hkh ) i e(-x-)] (3.7)
- iBye i(-kx-wt) + i B* e-i(-kx-wt) (3.7)
i gk sinh k(ho + z) i(k-wt) +A* e-i(k-wt)
2 w cosh kho e-(-
B ei(-kx-wt) B* e-i(-kx-wt)
These wave induced velocities go in to form the radiation stress terms, along with the mean
surface elevation and wave surface displacement. The general expression for each of the
radiation stress terms are presented as a review from Mei (1983) as
(3.8)
s, = p i 2 dz + p g(- z) dz + (p /"-- d( dz -
J-ho J-ho -ho \ Jz X /
p /w'2 dz (h, + )2 + (77)2 (3.9)
J-ho 2 2
SXy = p u'v' dz (3.10)
3- ho
Syy = p v'2dz + pI g( z) dz + p d( dz -
.-ho -ho J-ho\ Jz O y
7 9 nP9
f w' dz P (ho + )2 + ) (3.11)
-ho 2 2
From these the derivation will now be completed for the two-wave system without currents.
First note that for any two complex number R and Q, that RQ* + R*Q = 2R[RQ*] and
RQ* R*Q = 2ij[RQ*]. Now consider each term in the radiation stress expressions
separately;
TERM ONE
PfhoU'2 dz
p f 2dz = 2k 2|AI2 + 2k2BI2 1 4k2R [AB*e2ikx]
+ 21A 12 + 21Bx12 + 4 R[A B e2i"kx]
4 k[AA*] 4k [ABj e2ikr] 4 k [A B*e2ikx]
+ 4 k [B B*]} p (gz)2 dz (3.12)
2 w -h cosh kho
The integral in this expression has the solution,
S(cosh k(ho dz = 2 sinh 2kh, wCg
-ho cosh kho g Vk2 sinh 2kho gk
17
TERM TWO
pf2ghog(i z)dz
p g(-- z) dz= -(h, + i7)2=p
-ho 2 2
TERM THREE
P fh f'l7uW' d dz
Pjfho z ax --5--
p 7 o7- dC dz =
J-hoJz 1 OX
Pg { 2 IA2 + 2 IB|12 + 2 R[A Al] + 4 [A B e2ikx]
+ 2 Z[B BL] + 2 R[A BL e2ikx] 8k S[A Bd e2ik. ]
- 8 k A[Ax B* e2ikx] + 2 R[Axx B* e2ikx]
-82 [A ik sinh 2k(ho + ) d( dz
S8 -ho J sinh 2kho
(3.15)
The integral is evaluated as,
f f sinh 2k(ho + () 1
h ik-h + ( dC dz = [2kho coth 2kho 1]
-ho Jz sinh2kho 4k2
(3.16)
TERM FOUR
P f-h--2 dz
(3.14)
S 2dz = ( sinhk2kh ) { IA2 + IB12 + 2R[AB*e2ikx]} (3.17)
TERM FIVE
-P(ho + )2
TERM SIX
Pg (/)2 = P.g { A2 + IB12 + 2,[AB*e2ikx]} (3.18)
Then writing the final expression for the radiation stress component Sxx,
4 (1 + sin h, k2A2 + k2 BI2 2k2 [AB*e2ik]
+ IA, 2 + iBx12 + 2R[AxxBe2ikx]
2 k S[A A*] 2 k .A[A B e2ikx] k Q[AB B* e2ikx] + 2 kQ [B B]}
+ -g(2khcoth2kh 1) { Ax12 + IB 12 + R[AA*x] + R[B B*]
+ [[A BZ e2ikx] 4 k S[A B e2ikx] 4 k b [Ax B* e2ikx]
+ R[AZ B* e2ik] 4 k2 R[A B* e2ikx] + 2 R[Ax B e2ik] }
sinh2kho) {IA2 + IB|2 + 2 R[AB*e2i]}
+ {A|12 + B|12 + 2R[AB*e2]ik} (3.19)
19
Now for the terms in the S,, expression, only the first applies.
P f hu' v' dz
and the expression for Sx follows as
sY = 4 1 sinh2k h) {[A A A] + [A, B; e2ik]
+ ~2[B Aye2ikx] + [Bx B] kQ [AA;]
k B[A Be2ikx] k [A B* e2ikx] + k3[B B*] (3.20)
For the Syy expression, which is similar to the Sxx term in the integral representations,
the final expression is
Syy = 1 + 2 k h ) {Ay2 + By2 + 2 [AB*e2ikx]
4k2 *sinh 2 k h,2
+ g(2khcoth2kh 1) {IAy\2 + IByI2 + 2l [AB*e2ik ] + R[A Ay]
+ ~[B B,] + R[AB*;e2ikx] + R[Ay B* e2ik]}
p sinh2k ho { A2 + IB12 + 2 R[AB*e2ik]}
+ {jAI2 + IB12 + 2 R[A B* e2ikx] (3.21)
The terms arising because of the two-wave system couple nicely when not considering
currents. After doing some preliminary testing of the model, it was found that the interac-
tion between the two waves in the radiation stress terms would give rise to small currents.
Since the only concern for the modelling effort was currents arising from differentials in
overall energy, these small disturbances did not matter. Thus it was decided to calculate
the wave field for a given current field and let the bar coupling extract the energy from
20
the incident wave, then set the reflected wave field to zero in the grid prior to calculating
the radiation stresses. This would not only remove the currents over the bar field due to
interplay of the partial wave, but also allow a simple form of the subroutine that included
the incident wavenumber. The equations are the same as (3.19)-(3.21), but contain only
incident A terms. The final expressions used in the model will be presented in Chapter
Four.
3.3 Wave Model
3.3.1 Governing Equation
The wave model emulates the wave equation presented in this section. It governs the
energy balance of a wave travelling over a mild bottom with fast undulations, and riding
on a current. The parent equations come from two sources, one being Kirby (1988) and
the other being a University of Florida Technical Report, Kirby (1986b). The first paper
developed coupled equations for waves on a current over a bar field on a flat bottom. The
other paper handled such a case, but on a mild slope. The forcing from Kirby (1988) is used
in the present thesis, and the main left hand side of Kirby (1986b) is also used to complete
the full equation. As a review, the coupled equations for the incident and reflected waves
from Kirby (1988) are the refraction approximation,
{TT, + (C., cosOl + Uo)Tx, + (Cg, sin0i + Vo)Ty,} = -c R (3.22)
2 {RT + (Cg cos02 + Uo)Rx, + (-C, sin2 + Vo)Ry} = fnT (3.23)
where the following terms are defined;
Qc = 'oc cos(01 + 2) -+ O1c,
(3.24)
gklk2Db
o- = (3.25)
oc 4wcoshklhcoshk2h (3.25)
fl, = (4wa coshAh)-'I {AUo(12a1 lr2)Db
[gkik2 cos(81 + 2) + 02 ( 2 + (AUo)2)]Db} (3.26)
g g
a = AhF tanhAh
2r
A = Barfield wavenumber =
Lb
Uo
Fx = Froude number =
Ii = kicosOi;
ai = intrinsic frequency = V/gki tanh kih
02 = 02
27r
w = absolute frequency =
wave period
Db = Uniform bar amplitude
The amplitudes in these equations are surface displacement amplitudes divided by intrinsic
frequency, or T = A/la and R = B/a2. From the technical report, the parabolic coupled
equations are written as,
2ia(Cgi + Uo)T + 2oi(Cgi + Uo)(kl k1)T + 2iiTVTy
+ i {[l(Cg + Uo)]. + (aiV),} T + {(CCgl V2)(T),}
=coupling with R (3.27)
2iO2(C,2 U,)R, 2U2(Cg2 Uo)(k2 k2)R 2iC2VR,
+ i {[2(Cg2- Uo)] (02V)} R {(CCg2 V2)(R)y}
= coupling with T (3.28)
By dividing by 2i from (3.27-28), and multiplying w to the coupling from Kirby (1988), a
set of equations result for the case of a bar field on a mild slope with currents. These are
then solved using a Crank-Nicolson method for each wave until convergence is achieved. A
term w that handles the wave breaking energy dissipation is also included. This is added
to the incident wave only, since the reflected wave will not be expected to break anywhere.
This term is obtained from Dally et al. (1985). The concept is that if wave height exceeds
a criterion based on the local water depth, written as
2[A|
2AI> (3.29)
D
then breaking commences, where K = 0.78. The term is then defined as
K [ 72h2
w = 2 1 A (3.30)
where 7 = 0.4 and K = 0.15 for the present paper. The full parabolic equations used in
the model are then
al(Cgi + U,)T. ioi(C,. + Uo)(ki k)T + olVT
+ {[al(C1 + Uo)]. + (aV),~ } T V)(T)y
+ -c7(Ci + Uo)T
= c R (3.31)
02(Cg2 Uo)R, + ia2 (Cg2 U,)(k2 k2)R 92VRY
+ [2(Cg2 U0)] (2V)y} R + f (CCg2 V2)(R)
= -w T (3.32)
23
This is recognized as the equations governing waves that are at perfect resonance with
the bars. If there were a slight shift Q, then the forcing would have a different form.
CHAPTER 4
FINITE DIFFERENCING
4.1 Introduction
The incorporation of the equations developed in Chapter Three involves finite differenc-
ing, decisions on initial and boundary conditions, and placement within a coherent algorithm
that performs efficiently. The latter requirement is best described by means of a flow chart.
So before moving to the details of the coding, an overall view of the model MCSHERRY
would be timely. Figure 4.1 provides the logic behind the main program.
After inputs of bathymetry and wave information, the model iterates between the wave
model WAVEMOD and circulation model CIRC. A hyperbolic tangent startup function is
multiplied to the initial offshore incident amplitude, which slowly increases it's value to the
final after a specified number of iterations. This is done so that the circulation model will
not be shocked by suddenly large forcing terms. The form of the startup function is
C(t) = 0.5[1 + tanh(t/cl c2)]
where the parameters cl and c2 are chosen to alter the offset and slope of the curve. For
the present case, cl = 3000 and c2 = 3 to render a curve as shown in figure 4.2.
The subroutine FLOOD allows the grid area to fluctuate depending in the mean surface
elevation. It is called at every iteration, and alters values at wetted or dried grids to conserve
mass flux in those grids. The subroutine CHECK simply checks for convergent values, and
sets a flag for the model to quit if the condition is satisfied at each grid point for U,V and
Figure 4.1: Flow chart for model MCSHERRY.
1.0
0.8-
0.6-
4-
0.4
0.2
0.0
0 50 100 150 200
time (sec)
Figure 4.2: Hyperbolic tangent startup function.
27
4.2 Circulation Model
4.2.1 Method of Solution
The equations presented in equations (3.1)-(3.3) are solved using an implicit Alternating
Direction method that solves for the new U value and interim 7Y value in the x-sweep, and
the new 7 and new V in the y-sweep. This method applies to the matrix setup as suggested
by Sheng and Butler (1982). The detailed treatment of this method as used here can be
seen in Chapter 4 of Winer (1988).
The two unknowns in each sweep are applied at different places in the grid, where J
is at the grid centers and the mean velocities apply to the grid edges. Note that the first
offshore value is i7, at i = 1, and the first mean current U is at i = 2. The value of U and
V at i = 1 do not exist for the grid edge values, and are not used directly in the code. As
stated before, however, the wave model uses averages of the currents and for the i = 1 value
uses U, V at i = 2. The grid domain is shown in figure 4.3.
The x-sweep occurs in subroutine XCOEF, which sets up the coefficients, then passes
them to the double-sweep solver. At the lateral boundaries, the special subroutines XCOJ1
and XCOJN are called to handle the special coefficients when j = 1 and j = N. The
y-sweep is handled by subroutine YCOEF, except at the shoreline, where the variable grid
domain determined by FLOOD is deciphered by the subroutine YCVAR.
There is no time restriction for ADI method, but some of the terms are expressed
explicitly and added to the knowns of the implicit scheme, and restrictions are introduced.
The finite difference form of each of the terms in the momentum equations will now be
presented.
The radiation stress forcing is a function of the square of the amplitudes, thus being
nonlinear. Also, since the wave amplitudes are determined at grid centers in the wave
model, the radiation stresses also apply to grid centers.
The radiation stress gradients are put into the nth time level of the scheme, which
- offshore
1 4* 4
shoreline
Figure 4.3: Grid for circulation model.
(1,1)
1* II
(M.N)
29
referring back to the last chapter for the x- and y-sweeps are
1 (OSxx Sy)
+ -Y
pD Ox Dy
1 (OSX+, OSYY
pD Ox Dy
These terms are represented in finite difference form respectively by
1 SXXij SXXilj .SXYij+1 SXYij- + SXY + SXY j_
S'-"- "- .5 sXy,+ 2A -i-Ij-1
pD AX 2AY
1 SYYij- SYYi, .SXY+j SXYi-_,j + SXYi+lj SXYl,j_
pD AY b 2AX
The symbols in the code for each variable are hopefully clear. The subroutine for the
radiation stresses is RADIAT, which is called just after the wave model and prior to the
subroutine CIRC.
Referring to figure 4.4, the terms are used to find the mean velocities, and are therefore
difference to straddle the known (i, j) velocity location. Note differencing is done over one
grid space, and is accurate to second order.
The bottom friction pFlUIU is nonlinear, but uses mean currents from unknown and
known time steps, and is thus represented implicitly. The terms are represented based on
the work of Longuet-Higgins (1970a) as
bx = pF luoblU + pF UIU
Tby = pF luoblV + pF IUIV
These are finite difference according to the sample x-component,
Figure 4.4: Differencing of radiation stress gradients in x-sweep.
= Unew 2fH +
x T sinhkh U2 + )
where the linear formulation for maximum orbital velocity has been used, T being the wave
period. Since the reflected wave is considered nonexistent when performing the radiation
stress subroutine, this also applies for this subroutine. Therefore, the linear formulation
of maximum orbital velocity for the one wave applies. The subroutine that calculates this
term is called SHEAR2, and is called just prior to CIRC.
The lateral mixing is expressed in the momentum equations as
p y y a U
P 9y Dy *Y 1y9
J-1 J J+1
1-1
I -- -- -- --
Il~l - -
31
D Orl a V\
p Ox Ox X x)
These are added to the knowns at the nth time level explicitly, and thus introduce a serious
timestep limitation. In the complete finite difference form they are expressed as
t AtEY U,j+1 2Ui, + Uij-1
Oy AY2
atr AtE +1j 2Uij + Ui-,j
It AtEY
Ox AX2
For the solution to remain in front of the fastest travelling error, the timestep must satisfy
At
< 1
AX2
These values are calculated in subroutine EXCOEF just prior to CIRC. Note that the
coefficient Ey is set to the highest value of Ex, which occurs at the breakerline.
The terms represented in the momentum equations by U and V2J in the x direction
are represented in finite difference form respectively as,
Uj Ui+l,j Ui-1,J
u'3 2AX
and
(Vi,j + Vi-,i + V-l J,+ + V1i-l,1+) U+l Ui-l,j
4 2AY
The same form applies to the y direction equations. A four point average is used for the
cross-derivative of the four surrounding velocities as shown in figure 4.5
Figure 4.5: Differencing of advective acceleration in x-sweep.
4.2.2 Boundary Conditions
The solution for the equations governing mean currents and mean surface elevation
require boundary conditions for each boundary in the grid. Winer (1988) used a fixed lid
condition at the offshore edge, with no flow into the shoreline, and free flow at the lateral
boundaries. Each condition had to be reviewed for the present case, and changes were made
where appropriate.
The original model constrained the mean elevation j to be zero at the offshore grid edge.
This was a fine assumption for the one wave system, but was found to cause instabilities
if a reflected wave occurred. The instabilities were due to a severe mean surface gradient
at the first two grid rows which gave rise to large velocities. For the standing wavefield, if
the grid edge happened to occur on the envelope antinode, the problem was accentuated.
Obviously a radiation condition was required to free that edge to allow disturbances within
J-1 j J+l
I-1
the domain to escape. The condition
+ V + 1 - = 0 (4.1)
ot Ox pD a x
was applied to the x-sweep of the model, where the interim 7 is solved for at the grid center
i = 1. After making the decision to not include the reflected wave, however, the original
fixed lid condition was again used and is what is specified now. If the model were to be
modified to include the reflected wave, an appropriate radiation condition would need to be
developed. Also, if a strong current were to be flowing across the offshore edge, an energy
set down condition might suffice.
The lateral conditions are reflective, that is the y-sweep sets V = 0 at j = 1 and j = N.
This is to remain consistent with the wave model.
There is no flow into or out of the boundary, so the x-sweep condition sets U = 0 here.
Also, there is no longshore current, so the y-sweep specifies that V = 0.
The model MCSHERRY begins by setting all unknowns to zero in subroutine INPUT.
Then, as the waveheights are ramped to their full value, the forcing slowly increases in the
circulation program. The model iterates until the difference between unknowns in successive
time steps is smaller that a prescribed tolerance.
4.2.3 Radiation Stress Subroutine
The finite differencing of the radiation stress components that were developed in Chap-
ter Three will now be presented. The calculation of these terms is performed in the sub-
routine RADIAT. This immediately follows the wave model and precedes the circulation
model. Since the wave model is only accurate to O(E2) in spatial derivatives, there are
terms based on derivatives in the radiation stresses that should be dropped as being too
small. Referring to the expressions for S.x, S.y and Syy from Chapter Three, the terms
that should be dropped are
02
8z2
Ox2 4
0 3
Ox Oy
The terms that remain are difference to remain at grid centers where the amplitudes are
calculated. As an example, the derivatives are expressed as
A Ai+l,j Ai-,l,j
2AX
A Aij+j Ai,j-1
Ay ,
2AY
A Aij+l 2Ai, + Ai,j-1
dyy -= y2 --
As was stated in Chapter Two, the reflected waves have been considered nonexistent, since
they do not physically affect the circulation landward of the bar field. To incorporate only
the differences in wave energy in the grid, and not wave-wave interaction effects, the final
expressions that exist in subroutine RADIAT therefore neglect terms consisting of reflected
wave B. The final forms used in the model are
+ pg [(A2)}] (4.2)(kh)'
SXYX 1+ (k2A 2) 2(AA)) (A A())
S4(k2)i sinh2(kho) i
S2(kho)})
Y pg ( 2(kh[ ) [( (A (l2)
4(l sinh2(kho)]
+ [(|A|2)i (4.2)
pg ( 2(kh))'
4(k2) sinh2(kho)
4 sinh2(kho)j
35
pg 2(ko)
4 sinh2(kh)
+ [(IA12i]
+ 8(k2) tanh 2(Ikh) A) [ 2+ R(A Ay,)] (4.4)
4.3 Wave Model
4.3.1 Method of Solution
The wave equations have been presented by (3.31) and (3.32). The finite differencing
follows a Crank-Nicolson method, where a tridiagonal matrix is formed for the three forward-
row amplitudes based on three known values at the present row. The solution requires only
an initial condition at the offshore edge, and lateral conditions to complete the matrix.
For the two waves, two loops make up the subroutine WAVEMOD. The first loop solves
for the incident wave A up to the shoreline. Then, the second loop starts at the shoreline
and solves for B out to the offshore edge. Due to the coupling through the bottom boundary
condition, the two sweeps must be iterated until a convergent solution is obtained.
4.3.2 Finite Differencing
Now the incident equation will be difference to serve as an example. Again, the
amplitudes T and R refer to amplitudes divided by intrinsic frequency. Using a Crank-
Nicolson method leads to the finite-difference approximation of the wave equation of
1 ([ai(Ci9 + IU)]+1 + ai(C91 + U)]} (Tj+1 J}
2AX j Tj
{[(Cgi + U)(ki k)]t+1 + [air(Cg, + U)(ki k)]}} {Tj+1 + Tj}
+ 1 [r(Cgi + U+1 r1(Ci + U)] T +1 + I'
+ A yI(' + 7(C [j+, ;_,
4AX '- I- J
36
+ 1Y {[V]i+l [rlV]-J + [rlV]}+l [ClV]_1I {Ti+1 + T}
16AY 31 )j ][+ I T;1]
2 {[(CC, V2)- vji (CCgl V2 +1 +1 [ +1
8Ay2 3 +1 j
-[(ccgl V2)+l + (CC,_ V2)+][+1 -_ T
+ [(c V2)j1 + (CCig V2)1[T1 T]
-[(CC, V2)i + (CC9 V2) -l][T T_-i]
+ 1 {[ia(Cg ) ] + U)wT] [al(C91 + U)wT]}}
fI{[wR]' + [QRi }
2 3
The reflected form is similar, with sign changes where appropriate. These are then put into
three coefficients for the i and i +1 rows, and put into a complex tridiagonal matrix double-
sweep solver. Since currents exists at grid edges, the values to place in the WAVEMOD
need to be averaged straddling the grid center. Thus Uj is really ((Uj+' + UJ), and likewise
for V. The grid domain is set so that the first calculated value of U is at i = 2, so then the
averaged value for velocities are set to U2 and Vj2 when marching from i = 1 to i = 2.
In WAVEMOD, a preamble defines the coded meanings of each variable. Then state-
ment functions are defined for the coefficients in the equations for the two waves. Before the
first loop, subroutine WAVNUM is called to determine the incident and reflected wavenum-
bers, average velocities, celebrities, group velocities, intrinsic frequencies, and forcing terms
Qoc and 1,c as defined by (3.25) and (3.26). After this is done, the iterations begin, with
a maximum of four iterations, but convergence is usually complete within 0.1% after three.
If the wave is perceived to meet the breaking criteria and start dissipating energy at some
grid point, then w is defined something other than zero, and the row is iterated again,
letting the energy out. The wave model uses depths of the uniform bottom, not including
the barfield. Therefore, shoaling over the bars is not included in the driving forces.
4.3.3 Boundary Conditions
The boundary conditions are that the lateral boundaries behave as walls. This requires
that no reflected wave disturbance reach the walls and bounce back into the domain, which
would be unrealistic.
4.3.4 Lateral Smoothing
When doing some preliminary testing on the wave model, the radiation stress gradients
that would be put into the circulation model would be plotted simply for observation.
It was found that the plots exhibited much irregularity in the driving forces, more than
was expected. In fact, the circulation model had stability problems, and the irregularities
were thought to be at fault. After some review, suggestion was made that the parabolic
approximation had altered the nature of wave properties in the transverse y-direction.
Physically, a wave can have a maximum wavelength for given frequency and depth. This
is uniform in all directions. But when assuming small incident angles, thus primarially x-
direction propogation, which is done to get from a hyperbolic equation to a feasibly solvable
parabolic equation, a modification is introduced. Consider a surface described as
7' = Aeikx (4.5)
where A is the complex wave amplitude. For a = IA|, the surface in the y- direction is
77 = aeiky (4.6)
If one considers this wave in the y-direction, the radiation stresses are proportional to
the wavenumbers, with the Sxy term being proportional to k1, and the term S,, being
proportional to k2. And then when considering the gradients of these terms as they appear
in the momentum equations, the term would be proportional to k3. Therefore, there
is a large sensitivity to the wavenumber in the transverse direction. By equating (4.5) and
(4.6), the amplitude A could be expressed as
A = aeie-ik = (y)e-ik (4.7)
Now consider a general parabolic approximation equation of form
2ikAx + Ayy = 0 (4.8)
and by making the substitution for A, the ODE is
Ay, + 2k2A = 0, (4.9)
which has the solution
A = aei\Vky (4.10)
Therefore, a limit has been placed on the maximum physically meaningful wavenumber in a
parabolic wave equation for the y-direction. Instead of the normal L, = 2, the new limit
is L, = '-- Any wavenumber in the transverse direction that is larger than this will have
severe effects on the driving forces in the transverse direction.
By showing this, consider a numerical grid with lateral width W, and some calculated
wave field signal along a transect in the y-direction, as shown in figure 4.6.
A(y) can be written as a Fourier series
00
A(y) = ~anein (4.11)
n=O
where A = J is the base spatial frequency. By taking the FFT of the signal, the spectrum
will show energy distribution over each component nA. For the parabolic model, however,
any energy that exists beyond frequency nA = VZ2k is not physically meaningful.
Now for the wavefields calculated by the numerical model WAVEMOD, a transect in y
was FFT'd, and was found to have a very small amount of energy beyond what was physical.
This is shown in figure 4.7, where the cutoff frequency component was nA = 9.051. Then a
y=O y
RA(y)
SHORELINE
Figure 4.6: Wave calculated from parabolic equation on domain of width W.
40
transect of the driving forces ( + a ) were likewise decomposed spectrally, and a large
amount of energy existed beyond what was allowed. Shown in figure 4.8 for the y-direction,
this shows how a small amount of error in amplitudes can expressly alter what the circulation
model sees. The leakage of miniscule energy into high transverse wavenumbers in the
wavefield translated into major amounts of energy beyond the cutoff frequency component
for the driving forces. The process of the radiation stress gradients being proportional to
k3 shows the drawback when making a parabolic approximation. After taking the wave
field and doing the FFT row by row, and truncating the energy beyond nA = 9.051 and
inversing, a smoothed wave field was obtained. Figure 4.9 is the spectrum of the y-direction
driving after smoothing. By figures 4.10 and 4.11, the slight improvements are seen in the
wavefield, and in figures 4.12 through 4.13 the difference is remarkable. The variance in the
energy spectrum of the y-direction driving before the smoothing was calculated as 0.00804,
while after smoothing the total variance was 0.0030. In comparison to this 60% reduction
in variance, the wave field had an unsmoothed variance of 4.065 x 10-7 and a smoothed
variance of 3.80 x 10-7, about a 6% reduction.
0.000000008
0.000000006
0.000000004
0.000000002
0.000000000
5 10 15
frequency component
Figure 4.7: Spectrum of A(y) without smoothing. Cutoff at nA = 9.051.
20
0.0015
0.0010
0.0005
0.0000
0
20 40 60
frequency component
Figure 4.8: Spectrum of y-direction driving force before smoothing.
80
0.0015
0.0010
0.0005
0.0000
0 10 20 30
frequency component
Figure 4.9: Spectrum of y-direction driving force after smoothing.
40
Figure 4.10: IAI without smoothing. Bars are located on the bottom of this domain.
Figure 4.11: JA[ after smoothing. Bars are located on the bottom of this domain.
Figure 4.12: y-direction forcing without smoothing.
Figure 4.13: y-direction forcing after smoothing.
CHAPTER 5
LABORATORY WORK
5.1 Introduction
The laboratory work was done to get a better idea of the process of Bragg reflection,
where effects not represented by the wave model were shown. The following chapter treats
the laboratory phase of this project specifically. Since the work was completed at two
separate times, the procedures and setup unique to each will be presented separately. The
results from both sessions, however, will be presented together, with comparisons to a
simple numerical model.
The work was done at the Coastal and Ocean Engineering Laboratory (COEL) of the
University of Florida. To begin the discussion, the basin and the other necessary equipment
will be described.
5.2 Equipment
5.2.1 Basin
The modelling basin is illustrated in figure 5.1. The internal walls are moveable, and
were positioned at right angles to the wave crests to minimize side wall reflection. The
basin can be filled to a depth of 60 centimeters, and allows for waves up to heights of
several centimeters. The main testing region lies inside the internal wave vanes, to minimize
intrusion of backscatter from the outer basin walls. A wavemaker sits at one end, with a
sloping beach at the other.
Figure 5.1: Plan view of large basin.
5.2.2 Wavemaker
The wavemaker has 80 paddles, each being 23.5 centimeters in width. Each paddle can
be adjusted for stroke and phase. The entire system has a frequency range of 0.5 Hz to 1.25
Hz. Since the beach was oriented at 10 degrees relative to the wavemaker, the paddles were
adjusted to generate an oblique wave at a 10 degree angle, which then arrives normally at
the beach. This is done by calculating the wavenumber for the wavefield over the desired
depth using a Newton-Raphson technique. Then the relation of paddle wavenumber kp to
the wavefield wavenumber vector kw is, by (5.1)
kp = k, sinO (5.1)
Then the paddle wavelength is Lp = -. For paddle width Wp, and paddle number np, the
phase of each paddle in degrees from 0 to 360 is, by (5.2)
(360)(np)(W,)
p = L (5.2)
5.2.3 Beach and Cart System
The beach gently slopes to a maximum height of 60 centimeters above the floor. A
structure has been placed over the beach to aid in various types of testing. It allows a cart
to ride on tracks and be placed at any (x, y) coordinate in an 8 meter by 17 meter area on
the beach. For a plan view of this track system, see figure 5.2.
The track sits several centimeters above the greatest water elevation allowable in the
basin. The cart is positioned in the y-direction by the wire connected to the handcrank at
the side of the basin. The variable speed motor fitted onto the cross-member for the cart
controls the position of the cart in the x-direction.
The cart can carry wave gages, current meters, and whatever else is necessary in an
experiment. See figure 5.3 for a side view of the cart. Since profiling is important in the
conduct of moveable bed studies, the cart is fitted with a profiling system. The profiler arm
Figure 5.2: Plan view of cart track system.
Figure 5.3: Side view of profiling cart.
extends underneath the cart, dragging over the beach surface. A potentiometer produces
a signal which is converted to readings of arm angle, while another fitted to the trailing
wheel measures in voltage how many revolutions that wheel has made. With these two
instruments, the beach (x, z) position can be traced digitally. A sample profile is depicted
in figure 5.4.
5.2.4 Electronics
The computer for the basin is a PDP-11 microcomputer. It receives all basin data
which is amplified at a box next to the basin. A remote terminal at the basin runs all of the
programs. The main program is ATOD.FOR which takes multiple channels simultaneously
at any rate for any length of time. The PDP- 11 is located in an instrument room about
200 feet from the basin.
Wave gages were the main source of data. The capacitance gages were built at the
HORIZ.
POTENTIOMETER
60 11 -TT I I I
40
N 20
0 j
0 200 400 600 800
x (cm)
Figure 5.4: Sample profile.
COEL. Figure 5.5 is a sample of the calibration points transposed on the calculated second
order curve. The response appears linear over the operational range. Calibration was done
daily, or if the temperature changed drastically. Along with the gages, a strip chart took
the signal from the cart gage in session one. It was used to trace wave envelopes from the
moving gage.
A Marsh-McBirney current meter took data for the second session. Since the calibration
was 10ft/sec = 1 Volt, and since the currents were small, an amplification of twenty
provided the necessary range to increase resolution. The data were quite erratic, and will
be addressed in a following section.
Lastly, a video camera was mounted above the beach. The camera was fitted with a re-
mote control to zoom, focus, point, and contrast. A taping device was included. Rhodamine
dye transport was recorded as a way to show current magnitude. The only drawback was
that the camera is black and white, and the dye does not show up as well as it could if the
8
6
-) 4
)-
2
0
-20
-2 I I I I I I I
1000 1500 2000 2500
digital
Figure 5.5: Sample calibration curve for wave gage.
camera filmed in color.
5.2.5 Bar Field
The bars were made of wood and metal. The same bars served both sessions, and are
shown in figure 5.6 without the metal covering. They were sized according to the local
water depths at their position of placement. The planform dimensions were 3.048 meters
long and 0.318 meters wide, and each bar was made to be 40% of the local water depth.
Once the height h was determined by inspection, the bar radius in centimeters followed
equation (5.3).
(19.05)2 h2
(19.05)2 + h = R (5.3)
2h
M
R 3.048
Figure 5.6: Bar frame design.
The full pattern was traced on paper, then transferred to the plywood frames and con-
structed. Since the main material was wood, the bars would float. They were anchored in
the first session using wires connected to the concrete floor. In the second session, the un-
dersides were filled with sand which was kept in place with a panel covering. Both methods
worked fine, but the first proved to be better in time. Another problem was wave orbital
velocities flowing under the bars, scouring large holes. Strips of sheet metal 4 inches wide
were extended off the bar edges as aprons sunk into the sand. This effectively stopped the
scouring.
5.3 Setup, Procedure and Data Analysis
This next section discusses the actual testing methods. Since different methods were
used in each session, each will be discussed individually. The first session will be addressed,
with it's purposes and procedures, followed by the second session.
NOTE
DIMENSIONS ARE
IN METERS 0.381
0 LL
h /\ \
, V \
rllll
5.3.1 Session One
The first session ran from September of 1987 to February of 1988. Since this was the
first work done in the lab on the currents attributed to Bragg reflection, the goals were
fairly loose. The purposes were specifically:
1. To measure Bragg resonant reflection.
2. To compare this reflection to theory.
3. To describe qualitatively the properties of the wave-induced current field.
4. To observe how reflection and current magnititude varied as wave period varied.
5. To observe sand movement during the experiment.
The first goal received most of the attention, with two methods being developed. The
first was to use two stationary gages, in line, and take data simultaneously. The data from
each then can be used to determine the reflection coefficient, and this method is from the
work of Goda, et. al. (1976). To review that method, for a wavefield which can be described
by a surface displacement
= a cos(kx wt + 61) + b cos(-kx wt + 62) (5.4)
the unknowns are a, b, 61, and 62. Four values for y from the data will therefore be needed.
Let two gages be located at x = 0 and x = 1, as in the figure 5.7.
If a data value is taken for 17 from each gage at times t = 0 and t = r, four values will
result. Each one can be expressed through (5.4) according to (5.5-8) as
71 = r(0,0) = a cos 61 + b cos 62 (5.5)
72 = (1, 0) = a cos(kl + 61) + b cos(-kl + 62) (5.6)
773 = 7(0, r) = a cos(61 wr) + b cos(62 wr) (5.7)
774 = 7(1, r) = a cos(61 + kl wr) + b cos(62 kl wr) (5.8)
Figure 5.7: Gage locations for two-gage reflection method.
Let the following relations hold;
c, = cos61 (5.9)
C2 = coS62 (5.10)
S1 = sin 6i (5.11)
s2 = sin b2 (5.12)
d = coskl (5.13)
e = sinkl (5.14)
f = cos wr (5.15)
g = sinwr (5.16)
Expanding out (5.5-8) and substituting (5.9)-(5.16) gives
?71 = a cl+ b 2 (5.17)
12 = (aci + b 2)d- (as b s2)e (5.18)
773 = (a c + bc2)f + (as + bs2)g (5.19)
7q4 = (aci -+bc2)f d + (bs2 asi)fe
+ (a s1 + b 82)g d + (a ci b c2)g e (5.20)
By substituting sequentially from (5.17) down to (5.20) and eliminating terms, the remaining
expressions are,
2acl = A = 171i+ 4 72 1fd (5.21)
ge
2as = B = (5.22)
g e
2bc2=C = q1 4 713 712 + d (5.23)
2 b C2 = 7 ?i 1 --- ---- (5.23)
ge
2b2 = D = 73 71f +2- d (5.24)
g e
From these expressions, the incident and reflected wave amplitudes are, by (5.25) and (5.26),
A2 + B2
a A2 B2 (5.25)
a C2 + 2
b C2 D2 (5.26)
4
The reflection coefficient is then Kr = b/a. The problem in this method is that if g or
e approach zero, the expressions for A, B, C and D blow up. Thus one must not choose
r = n( ) or 1 = n( ) where T = wave period, L = wavelength, n = integer.
To determine the accuracy of the method, a test was run where the reflection was from
the beach. To compare, a moving gage retrieved the partial standing wave envelope. It was
discovered that for various time increments r, the reflection coefficients were inconsistent.
This inconsistency was enough to scrap the two gage system, and use the envelopes as
Table 5.1: Bar heights and positions for Session One.
bar no. height (cm) x-position
1 7.62 420.
2 9.50 495.
3 11.20 570.
4 12.60 645.
the major source of reflection data. This required measuring by hand the minimum and
maximum envelope amplitudes, 7?min and r1max. Then the reflection coefficient is,
Kr = i]max r7min (5.27)
?7max + l7min
Envelope transects across the barfield were taken at four locations during the testing to get
a lateral average. The positions of these transects are presented later in this section. The
results are then presented and compared to the circulation model in the next session.
A camera was mounted above the testing site which filmed the movement of rhodamine
dye across the barfield at each wave period. The film is available for viewing through the
author.
The tests were run at several wave frequencies around the designed resonant frequency.
Each run was made from a still water condition. With reflection and current observations
made at each run, a feel of the dependence was gained. The different runs are listed in the
next section for the periods and reflection coefficients.
Observations were made in regard to sand movement during the course of the testing.
These observations are listed as part of the results in Chapter Six.
The setup consists of where the bars were oriented. By using the equation for the bars
from the local depths, the bars were built according to Table 5.1.
Plots of the testing bathymetry were done as testing commensed, and are available as
of this writing, but the actual data that was written onto magnetic tape is not immediately
I
60
40
N o20
0 2 4 6 8
x (m)
Figure 5.8: Reconstructed bar profile for Session One.
available to the author. There is no other recourse than to reconstruct what the bathymetry
was at the time from the existing plots. This is the only handicap from losing the data,
since the envelopes are printed on strip chart paper. Figure 5.8 shows the lab bathymetry
with the bars placed. The depth was 55 centimeters, and the breakerline in the testing was
at about x = 3 meters. Figure 5.9 gives an overhead view of the barfield, with four lines in
the x-direction offshore which represent the transect lines where envelope data was taken.
The profile is available as a test bathymetry in the numerical model, although the bars need
to be sinusoidal.
Once the bars were placed, the wave gages were arranged. The barfield would create
a reflected wave that would propagate and diffract once beyond the bars. Two gages were
mounted on tripods in the center of the bar field longitudinally, and directly offshore of the
last bar. As a control, one gage was put off to the side, away from the influence of the
barfield. Another gage was on the moveable cart. Figure 5.10 shows a plan of the beach
1365 1300 1175 1110
16.0 16.0 ----16.0-
24.0 "0-4.0
32.0 0 __- 32.0 ---- 32.0
n f____
Figure 5.9: Plan view of barfield for Session One.
-7
Figure 5.10: Gage locations for Session One.
region with the gage placements.
This concludes the discourse on session one. The data from the experiment will be
presented later, with the test run parameters, in Chapter Six. Now the second session will
be addressed with the procedures and setup.
5.3.2 Session Two
With the experience from the first session, the second experiment ran from September
of 1988 to December of the same year. The goals were more defined in this second session.
1. To repeat Bragg resonant reflection and a resultant current.
2. To extend the reflection to frequencies off the no-current resonant value.
3. To observe differences between moving away from resonance and shifting toward
it.
Figure 5.11: Three gages inline with x-direction waves
4. To quantify current values.
In this session, two methods were again employed. The first was from the work of
Funke and Mansard (1980). Here, three in line wave gages took simultaneous data, and
the reflection coefficient came out of the spectra. The work covered waves travelling over a
flat bottom, with no current involved. Their results therefore needed to be revised for the
addition of current and variable topography. For the setup in figure 5.11, the algorithm was
modified in the following manner;
The wave profile for any probe p = 1,2,3 can be described as a Fourier series sum of
each component k by
N / 2rkt \
qp(t) = Ap,k sin + ap,k (5.28)
k=1
64
where Ap,k is the Fourier coefficient for frequency -, T is the time length of the data run,
ap,k is a random phase shift, and N is the desired number of components for the profile.
The Fourier coefficient is gotten from the transform of the time series at each probe p, rp(t),
which for any probe would be, by (5.29)
7 {lp(t)}= Bp,k = CI,kexp i 2(Xlk +iOk
+ CR,k exp {i 27(X1 2XR1 X1P) ik + O(k)}
I Lk
+ Yp,k exp{i(pp,k)} (5.29)
These equations are actually (1) and (10) in Funke and Mansard (1980). They point out
that the transforms represent functions of the complex amplitudes for the incident and
reflected waves. The three expressions for the FFT's of the three probes then lead to a
least squares determination of the incident and reflected amplitudes, Zi and ZR.
B1,k = ZI,k + ZR,k + ZN,1,k (5.30)
B2,k = Ks,I,2 ZI,k exp[qI,i,2] + Ks,R,2 ZR,k exp[--R,1,2] + ZN,2,k (5.31)
B3,k = Ks,I,3 ZI,k exp[~I,1,3] + Ks,R,3 ZR,k exp[-'Ra,,3] + ZN,3,k (5.32)
The terms in these three equations are defined as follows;
ZI,k = CI,kexp i27r +L iOk (5.33)
SLk .2X1+2XR1)
ZR,k = CR, exp {i 27(X1 Lk X +i(O, k)} (5.34)
ZN,p,k = Yp,kexp {iPp,k} (5.35)
The shoaling coefficients in the incident and reflected directions for each probe, and group
velocity Cg are,
Ks,I,2 = (5.36)
K,3 = (5.37)
C9,3
V ^.2
sR,2 = (5.38)
V C9
Ks,,3 = Cg, (5.39)
The phase accumulations are obtained according to
yX2I
PI,1,2 = j kidx (5.40)
,1,3 = kldx (5.41)
R,1,2 = k dx (5.42)
JX2
/X1I
RR,1,3 = 3 kR dx (5.43)
JX3
(5.44)
These modifications reflect the different case of waves over a variable topography on a
current. Funke and Mansard (1980) addressed waves over a flat bottom where the incident
and reflected wavenumbers were equal, and wavelength Lk for each frequency was constant
in the domain. The reflective and incident wavenumbers k are found from the dispersion
relation that includes the current effect, using an iterative technique. Setting the ZN terms
to a small error parameter E, the least squares method is followed and the resulting equations
for unknowns ZI and ZR are
Zi {2 + 2K',I,2exp[2ii,1,2] + 2K',I,3exp[2i~P,i,3]}
66
+ ZR {2 -: /s,,,2AKs,R,2exp[i(II,1,2 R,1,2)] + 2Ks,I,3KS,,R,3exp[i(JI,1,3 T'R,1,3)]}
= 2B1 + 2Ks,I,2exp[i'@J,1,2]B2 + 2Ks,i,3exp[i',1 ,3]B45.45)
and then
ZI {2 + 2Ks,I,2Ks,R,2exp[i(PI,1,2 TR,1,2)] + 2K,I,,3Ks,,R,3exp[i('I,1,3 @R,1,3)]}
+ ZR 2 + 2I,R,2exp[-2iQR,1,2] + 2K ,R,3exp[-2iR,1,3]}
= 2B1 + 2Ks,R,2exp[-i'R,1,2]B2 + 2K',I,3 exp[-iQR,1,3]B3 (5.46)
These can be solved for the two complex amplitudes, Zi and ZR, and then Kr = ~IZR
The modifications were done numerically according to the following forms. For the
dispersion relation with opposing current,
w = -Uo k + V/gk tanh kh (5.47)
a Newton-Raphson iterative technique was used to solve for k, given depth, frequency, and
current Uo. It is assumed here that the incident wave angle 0 = 0, so that k cos 0 = k.
The simple integration scheme of the phase in space is performed by,
Wp = (kp + k~p-1)A -'-1 + Xp-1 (5.48)
2
where again p = the gage number, either 1, 2, or 3. Ax p-.p-1 is the interval size between
gages p and p 1. The intervals were on the order of 15 centimeters, so that the bottom
could be considered a plane sloping beach in these intervals, allowing the simple integration.
The group velocity used for the shoaling coefficient K, is, for each gage p,
Cgp P (1 + sinh2(h)p (5.49)
This method was tested for accuracy for reflection from the beach, being compared to
envelope data. It was found that the reflections differed between the two methods by as
67
much as 60%, and the development of the equations thus leaves something in error. Rather
than an in-depth search for an accurate form, it was again decided to use the envelope
method, as in Session One. Again, partial standing wave envelopes were retrieved from
a moving gage across the barfield, but were recorded digitally, and analyzed as per the
methods outlined below. The envelopes were taken at seven transects over the barfield. It
was found that the three center transects represented the highest consistent reflection for
a test. These locations are shown in the plan view in figure 5.16. After determining the
reflection from each transect, it was plotted with location in figure 5.15. As is apparent,
sections 850, 900, and 950 give a consistent representation of reflection. The method of
analyzation was to take each envelope, as shown in figure 5.12, and condense it into envelope
height versus x-position. A sample is shown in figure 5.13. The reflection coefficient was
calculated by taking a local maximum and immediately following minimum, using equation
(5.27). Thus a plot of Kr with x was obtained, as seen in figure 5.14. The maximum
reflection was taken from each envelope, and the three were averaged. This was done for
each wave period. The results are presented in the next chapter and show a repeatability of
the process, which is a good sign. Further similarity with theoretical output gives confidence
in the values themselves. From figure 5.14, there is a variability in measured reflection along
the bars, and the ending minimum should be disregarded since the cart may not have traced
the full wave envelope before stopping. The variability within the bar region, however, could
be attributed to focusing effects over the bar field.
This data was taken using program ATOD.FOR for a sampling frequency of 26.0485
IIz, and a run of 2048 points. The gage was the cart mounted gage, moved slowly through
the water over the barfield, insuring that the wave excursions stayed in the middle portion
of the gage wire.
The second priority was to get quantified results for current induced by Bragg reflection.
At each period, a line of current data was taken across the front edge of the barfield. The
locations were at the same locations used for the envelope measurements, except they were
I
0
CU)
*-
2
0
-2
-4
200
400 600
800
x (cm)
Figure 5.12: Sample envelope, Session Two. Retrieved digitally from one of the three
on-offshore transects using the moving cart gage.
4)
,c)
Cd
OJ
6
4
2
O
0-'-
200
Fi
0.5
0.4
0.3
0.2
0.1
0.0
200
Figure 5.14:
I I I I I I I I
I I I I I
I I I I I
400
600
I I1
x (cm)
Sample Kr with offshore distance x from Session Two.
400 600
x (cm)
gure 5.13: Sample envelope height versus X.
800
800
-
. I
I ...
' I I I
II IL~
0.5 -
0.4
0.3
0.2
0.1
6-0
0.0
600
700 800 900
1000 1100
section number
Figure 5.15: Plot of K, with transect location, Session Two.
16.0 6 16.0 16.0
2 4.8 0_ ---10
Figure 5.16: Plan view of barfield, showing data collection locations.
stationary just inshore of the first bar. Figure 5.16 is a plan view of the session two barfield
with the locations used for current meter measurements and envelope data. A side view is
shown in figure 5.18, where the water depth was 50 centimeters, and the breakerline was at
x = 320 centimeters. The current meter was placed at 750, 800, 850, 900, 950, 1000, and
1050.
Data was taken at 26.0485 Hz for 512 points. The data proved very erratic from the two-
channel Marsh-McBirney current meter. The calibration for the instrument was 10 ft/sec =
72
1 Volt, and since the expected currents would be on the order of a few centimeters per
second, an amplification was done. An amplifier was built in house to magnify the signal
by a factor of twenty. The raw results from such a time series is given in figure 5.17, which
represents the velocity magnitude. Because of the nature of the signal, any trustworthiness
of discharge calculations would be in doubt. The irregular signal could be due to interference
from the metallic bars that were about 10 centimeters away from the probe. Still water
readings taken in a container away from the bars showed a fairly stable signal. Other reasons
could be excessive low frequency water displacements, but nothing like that was observed.
Coupled with a serious time limitation as of this writing, the choice was made to disregard
laboratory quantification of the rip current. A major goal of Session Two therefore is yet
to be realized, but the reflection results weakly imply the presence of a current, and the
dye experiments from Session One that are on film definitely show the rip current. Indeed,
observations made at the time were that the currents were very mild, and did not get past
the first bar.
Unlike session one, session two was a continuous test between wave period changes. The
intention was to extend the reflective peak as long as possible, thus keeping the induced cir-
culation from previous periods going through the enhanced reflection. Once it was observed
that circulation had dropped considerably, the runs began at rest from a high nonresonant
period, and moved at slow increments toward resonance. Thus a test was made to see if
reflection commenced upon return at the same place where it slacked when moving away
from resonance. As said, the current magnitudes were so small that it was difficult to see
a definite drop in reflection or current strength. Nevertheless, the results are presented in
Chapter Six and still offer some useful lessons for future study.
The bars were the same from session one, but some of the bars were buried slightly to
accommodate the new conditions. Table 5.2 lists the heights and positions of the bars in this
15
0
r--
Q)
10
5
5 10 15
time (sec)
Figure 5.17: Sample time series of magnitude of current meter data, Session Two.
part.
For an overhead view of the barfield, refer back to figure 5.16. A side view of the
barfield is shown in figure 5.18.
The instrument setup consisted of everything being put on the cart. The three gages
were positioned inline on the cart, along with the current meter. The overhead of the
Table 5.2: Session Two bar heights and positions
Sbar no. height (cm) x-position (cm)
1 6.25 436
2 7.00 500
3 11.00 565
4 12.50 631
20
60
40
20
N 20
0
0 200 400 600 800
x (cm)
Figure 5.18: Profile for Session Two.
instrument array is shown in figure 5.19. With the points made earlier, the three-gage
system was dropped in favor of the one moving gage to retrieve envelopes across the barfield
in the x-direction.
This concludes the description of Session Two. The next chapter describes the results
obtained from the laboratory work.
Figure 5.19: Instrument setup for Session Two.
CHAPTER 6
RESULTS AND CONCLUSIONS
6.1 Introduction
The results from the work done in the past year and a half are presented in this chapter.
The results are from laboratory data, where the methods that were used to obtain this data
are outlined in Chapter Five. The other source of results is the numerical modeling effort.
This consists of the wave- current model MCSIIERRY that as of this writing handles a
limited set of conditions. The preliminary findings could be construed as inconclusive, but
definitely suggest the need for more attention. The model should be running at full speed
in the next year (1989).
6.2 Laboratory Work
The data collected consisted of reflection data from wave gages, current meter data,
bathymetric data from the profiler, and some data from the camera. The results will be
presented in tabular form and pictorally where possible.
6.2.1 Reflection Results
The Session One reflection data is presented in table 6.1. Each run represents a Static
test, i.e., startup from a still basin. The averaged data over the middle transects is also
plotted in figure 6.1, transposed upon the theoretical prediction for a one dimensional
barfield without currents. This one-dimensional model is from Kirby (1987) that excludes
currents. This has been placed over the data to simply give an idea of the effect the
currents have on the data. The incident wave had a resonant period of 0.98 seconds, and
was 7 centimeters in height. The waves arrived normally onto the beach, with a breakerline
at x = 3 meters, referring to the profile in figure 5.8.
76
I
Table 6.1: K, and comments from Session One.
period (sec) 1110 1175 1300 1365 comments on rip
0.980 0.357 0.667 0.529 0.362 weak rip current
0.910 0.347 0.189 0.235 0.091 no current
0.963 0.172 0.579 0.368 0.238 no current
1.050 0.368 0.550 0.368 0.538 rip current
1.110 0.294 0.692 0.458 0.514 significant current.
1.180 0.333 0.529 0.300 0.176 no rip current
1.110 0.333 0.474 0.647 0.428 again strong rip repeated.
1.085 0.290 0.412 0.579 0.375 slightly weaker than T=1.ll.
1.020 0.474 0.438 0.400 0.405 weak rip, dye is stagnant.
1.000 0.333 0.474 0.538 0.333 weak current
0.942 0.273 0.444 0.350 0.170 no current
The predicted reflection from the one-dimensional model underestimates the reflection
for the given bathymetry, for the envelope data describes 50% reflection. The data were not
carried to sufficiently large periods so that the reflection clearly died away. From the films
of the current, the magnitude increases in strength until a period of 1.11 seconds where it is
quite strong, after which it dies off at 1.18 seconds. This would tend to suggest that while
the waves are moving increasingly off no-current resonance, the rip current increases to
keep the waves at resonance. This supports the hypothesis that the system seeks a resonant
condition.
From Session Two, the data were more abundant, but perhaps not as telling. Much
of the setup was the same as in Session One, except that the depths over the bar field
were less. This is due possibly to the state of the beach when testing began, in that the
slope was less than in Session One. To offset this the bars were buried into the sand about
10% of their heights, which was typically about a centimeter. Also, due to the testing that
had been going on before this work, some of the sand had been removed from the main
testing region, so water depths could not be as high. In Session One, the testing depth was
55 centimeters, and for Session Two the depth was 50 centimeters. To further offset this,
0.6
0.4
0.2
0.0
0.5 1.0 1.5 2.0
period (sec)
Figure 6.1: Plot of IK with wave period, Session One data and theoretical prediction from
1-D model Kirby (1987).
79
wave heights were decreased from 7 centimeters offshore to 5.5 centimeters offshore. This
seems to be the main value that has the greatest effect on the results, for a rip current never
developed as strongly as the one appearing in Session One. The reflections were high, but
it is conjectured that there was not enough available energy to drive the circulation.
The reflection data are presented in three tables, where the first set refers to reflection
with absolutely no discernible current throughout the course of the runs. The next set was
made after digging sand out from the region in between the bars and the land. This allowed
a weak current to form in some of the periods. The third set refers to runs that started at
low frequencies and slowly approached the no-current resonant case. In all the testing in
Session Two, the runs were continuous, without stopping the wavemaker.
The averaged data from the first set with no current are tabulated in table 6.2, and
illustrated in figure 6.2, which has the theoretical result transposed. Again, the theoretical
values are from the no-current, one- dimensional model of Kirby (1987). Three or less runs
were done at each period to get an idea of the repeatability of the experiment results.
Like the first type, the second type started at resonance and moved slowly outward. A
weak rip appeared after removing some sand from the foreshore. The tabulated results are
given in table 6.3, with the plot in figure 6.3. The two profiles from set one and two are
shown in figure 6.4.
From these plots, the data follows the theoretical fairly well, but the data drops away
at later periods. The sand did not move onto the first bar at all, whereas in Session One,
the first bar had to be repeatably uncovered.
The third type started off resonance, then moved toward resonance continuously until
reflection commensed and a weak current was excited. Table 6.4 gives the results, also
plotted in figure 6.5.
The difference between this test and the other is again not clearly evident, and may be
due to the lack of dominance of the current.
The one-dimensional algorithm from Kirby (1987) is over a variable mild topography.
Table 6.2: Kr. Reflections for case of no discernable rip current, Session Two.
Period (sec)__ Kr_
1.014 0.419 0.515 0.423
1.020 0.381 0.418
1.030 0.433 0.503
1.050 0.507 0.471
1.065 0.511 0.508
1.080 0.530 0.505
1.100 0.494 0.526 0.525
1.125 0.516 0.481
1.150 0.470 0.454 0.473
1.175 0.518 0.493
1.200 0.425 0.386
1.225 0.323 0.341
1.300 0.385 0.310
1.350 0.374 0.414
1.400 0.396 0.372
1.450 0.386 0.367
1.500 0.352 0.339
0.6 I I
000
1.0 1.5 2.0
0 0
A
period (sec)
Figure 6.2: Plot of K,- with period, data and theoretical for the case of no discernable rip
current.
1.0 1.5 2.0
Table 6.3: Kr measurements from case of a discernable shore rip current. Session Two.
period (sec) i K,
1.014 0.425 0.464
1.030 0.391 0.442
1.050 0.521 0.536 0.539
1.080 0.466 0.510 0.483
1.100 0.483 0.503 0.514
1.125 0.482 0.521 0.482
1.135 0.447 0.486
1.150 0.532 0.520
1.160 0.461 0.490
1.200 0.402
1.250 0.349
1.300 0.349
1.350 0.231 0.289
1.400 0.215
1.450 0.257
1.500 0.259
1.550 0.323
1.600 0.310
1.650 0.272_
0.6
0
o o_
0.4 0
0.2
0.0
1.0
Figure 6.3: Plot of K, with period for the
theoretical.
1.5 2.0
period (sec)
case of a discernable shore rip current. Data and
1 I I 1
1 I I I I I
I I I I I
I I I I I I I I
x (m)
Figure 6.4: Two profiles, dotted line referring to profile used in the second and third sets
of testing, where a weak current formed.
Table 6.4: Kr measurements for case of discernable shore rip current, moving toward reso-
nance.
period (sec) IKr
1.650 0.216
1.600 0.330
1.550 0.323
1.500 0.316
1.450 0.316
1.400 0.356
1.350 0.249
1.300 0.302
1.200 0.436 0.460
1.175 0.464 0.459
1.160 0.518 0.512
1.150 0.529 0.520
1.135 0.487 0.556
0.6
0.5
0.4
0.3
0.2
'''''''''`'
I1 I 1 1
I I I I I
"C1~,
0.6 I
0.4-
0.2
0.0
1.0 1.5 2.0
period (sec)
Figure 6.5: Plot of K, with period for case of discernable shore rip current, moving toward
resonance.
86
It again excludes current effects. The plots do not show the kind of reflection extension
that could occur because the currents were not as strong in the second session as the first.
Thus, only so much tuning was possible. Some notes about the comparisons of the first and
second sessions are worth noting to explain some problems with the execution.
1. The waves were slightly higher in the first session, thus involving more current
driving energy.
2. The anchoring system kept the bars very still in the first session, while there was
slight movement in the second.
3. Not enough time was allowed between period shifts in the second test to allow a
proper current maintainance.
The first reason seems the main cause of the decreased current generation. As is stated in
the theory, if a wave field does not have enough available energy, the rip current will not
be strong enough to achieve tuning of the nonresonant incident waves. Thus the efficiency
of the bar field is comparable to the case with no current, and affects only closely resonant
waves. For practical considerations, the most damaging waves will probably carry enough
energy to drive a sufficient current, and the bar field will be effective when most needed.
6.2.2 Bathymetry
This last section deals with the movement of sand during the experiments. There are
profile data on backup tape from Session One that is presently not available to the author.
Observations conclude that the current was sufficient to move sand over the barfield. In
fact, it was necessary to uncover the first bar periodically, which had been buried and
which lessened the effectiveness of the barfield. In the second Session, the current was too
weak to transport sand, and the profile shown in figure 5.18 was the same throughout the
testing.
87
6.3 Numerical Results
After much tinkering, some results are available from the numerical model. As of the
writing, however, the model requires another two months of work. Specifically, lateral
mixing was left out because the time step limitations would make the model very slow, too
slow for the given time allotment. Therefore there are circulation variations in the domain
that would probably not exist if the cells were allowed to influence on another. However, it
is thought that the overall pattern is similar for both cases.
The results presented now are from the no-current wave model coupled with the cir-
culation model. The process to achieve results was to run the wave model once and get
the final radiation stresses. These terms were then multiplied by the startup function and
inserted into the circulation model. The tests were run on a common bathymetry which is
shown in figure 6.6.
Without the wave-current equations to govern the wave field, the thesis cannot answer
the question about whether the system tends toward resonance through numerical results.
However, the no-current model can answer questions about how the magnitude of the current
is affected by wave height. Therefore, three wave inputs were chosen for the bathymetry,
where the period was the no- current resonant period. As a point of interest, the wave field
and mean surface elevation that are calculated by the wave model for an input of a .03
meter wave are shown in figures 6.7-6.9.
Then the following plots show how the current pattern varies as the amplitude is in-
creased. Velocity times total depth is plotted in the figures, which describes the discharge
across a vertical line. It is seen as expected that the rip current discharge increases as
wave height increases. The circular cell over the barfield ends is thought to be attributed
to a slow leaking of energy from the incident wave over the bar field. The action stated
in Chapters Three and Four of removing the reflected wave probably account for this lost
energy. The reflection coefficients remain the same for each wave height, being calculated
simply by the maximum reflected wave amplitude in the grid divided by the initial incident
0.1
0.0
-o. --
-0.2
0 .3 i i 1 1 1 1 1 1 1 1 1 1 1 1
0 2 4 6 8 10
x (m)
Figure 6.6: Bathymetry 1. Shoreward direction with increasing x.
.os1-N
Figure 6.7: jA[ for initial amplitude of .03 meters and perfectly resonant waves.
90
.03
Figure 6.8: IBI for initial amplitude of .03 meters and perfectly resonant waves, without
currents.
|
RSS
HIDE
