FINITE AMPLITUDE HARBOR OSCILLATIONS: 
THEORY AND EXPERIMENT 
by 
STEVEN R. ROGERS 
S.B., Massachusetts Institute of Technology 
( 1972) 
Submitted in partial fulfillment 
of the requirements for the degree of 
Doctor of Science 
at the 
Massachusetts Institute of Technology 
(January, 1977) 
Signature redacted 
Signature of Author ...... - ...•...... ,\~r~ ...... . 
I'"\,., Departme~t of Ph}(S)i cs, January, 1977 
Signature redacted 
Certified by . /' • , ~,,.., •. {J . . . . . .. -.-.- A 'fh~sis. S~p~r~i ;o; 
Signature redacted 
Accepted by • • • • • • • • ...-r -·/ • J ., -.-- ' ... ' .... - • • .- • . • . . . . • 
Chairman, Department of Physics 
FINITE AMPLITUDE HARBOR OSCILLATIONS:
THEORY AND EXPERIMENT
by
STEVEN R. ROGERS
Submitted to the Department of Physics
on January 13, 1977 in partial fulfillment of the requirements
for the Degree of Doctor of Science
ABSTRACT
A nonlinear theory of finite amplitude, time-periodic oscillations
in a narrow harbor is developed, starting from the Boussinesq equations
for shallow water waves. The narrowness of the harbor permits a linear
treatment of radiation damping and a one-dimensional treatment of the
nonlinear harbor response. The wave field in the harbor is expanded in
harmonics, whose spatial dependence is governed by a set of coupled,
nonlinear ordinary differential equations, subject to two-point boundary
conditions. Solutions are obtained by analytical perturbation tech-
niques and by an iterative numerical procedure. Results indicate sig-
nificant nonlinear effects at large amplitudes.
Harbor resonance experiments were carried out using a model
harbor placed in a large, shallow wave basin. The effects of boundary
layer dissipation, flow separation, and spurious reflections in the
wave basin are analyzed. The experimental observations are found to
agree reasonably well with the proposed nonlinear theory.
-2-
Finally, the nonlinear theory is applied to large-scale harbors
and some theoretical predictions are presented.
Thesis Supervisor: Chiang C. Mei, Professor of Civil Engineering
-3-
ACKNOWLEDGEMENTS
The author wishes to express his appreciation to the Fannie
and John Hertz Foundation for a generous fellowship which enabled
him to devote full time to his studies, and to the Fluid Dynamics
Branch of the Office of Naval Research (Contract N062-228) for pro-
viding funds for laboratory supplies and computer time.
Many members of the Ralph M. Parsons Laboratory for Water
Resources and Hydrodynamics gave of their time and knowledge to
assist in this endeavor. Among them are Mr. Edward F. McCaffrey
(Research Engineer) and Mr. Stanley M. White, who designed the data
acquisition system used in the experiments, and Mr. Roy G. Milley,
who helped with the construction of the experimental apparatus.
Mrs. Elizabeth A. Quivey did a superb job typing the manuscript.
I would especially like to thank Mr. Dick K.-P. Yue and Mr.
Demosthenes C. Angelides for many valuable discussions, and for their
enthusiasm and fine sense of humor.
Professor Chiang C. Mei, who supervised this work, has been
a constant source of inspiration and guidance. I am deeply grateful
for his time, his encouragement, and, most of all, his friendship.
To my parents, who have taught me the importance of education,
I am eternally grateful.
Finally, I should like to dedicate this work to my wife,
Rahel, whose patience, understanding, and devotion have really made
it all possible.
-4-
TABLE OF CONTENTS
Page
TITLE PAGE 1
ABSTRACT 2
ACKNOWLEDGEMENTS 4
TABLE OF CONTENTS 5
LIST OF FIGURES 7
LIST OF TABLES 12
CAST OF CHARACTERS 13
1. INTRODUCTION 15
1.1 Description of the Problem and Review of Literature 15
1.2 Elements of Water Wave Theory 17
1.3 Linearized Theory of the Long, Narrow Bay 26
1.4 Scope of this Investigation 30
2. WAVE SYSTEM OF THE OCEAN 32
2.1 Motivation for a Linear Treatment 32
2.2 The Radiated Wave 33
2.3 The Impedance Boundary Condition 36
2.4 On the Validity of the Linear Approximation 40
3. RESONANT OSCILLATIONS OF A LONG, NARROW BAY 42
3.1 Reduction to One Space Dimension 42
3.2 Equations Governing Time-Periodic Solutions 44
3.3 Regular Perturbation Analysis for the Lowest Resonant
Modes or for Very Small Amplitude Waves 48
3.4 An Energy Theorem and Relation to the Korteweg-de
Vries Equation 61
-5-
Page
3.5 Numerical Solution by Iteration 68
4. EXPERIMENTAL STUDY OF THE LONG, NARROW BAY 75
4.1 Equipment and Procedures 75
4.1.1 Wave Basin 75
4.1.2 Wave Gauges 77
4.1.3 On-line Digital Computer 77
4.1.4 Model Harbor 79
4.2 Multiple Reflections from Wavemaker 80
4.3 Real Fluid Effects 89
4.3.1 Viscous Dissipation in Wall Boundary Layers 89
4.3.2 Separation Loss at the Entrance 93
5. EXPERIMENTAL AND THEORETICAL RESULTS 101
5.1 Comparison of Theory with Experiment 101
5.2 Some Theoretical Predictions for a Large Scale Harbor 131
6. CONCLUDING REMARKS 162
6.1 Validity of the Various Approximations 162
6.2 On the Importance of Nonlinear and Frictional Effects 163
6.3 Suggestions for Future Research 164
LIST OF REFERENCES 167
APPENDIX A: Numerical Iteration Procedure for Solving a
Nonlinear, Two-Point Boundary Value Problem 171
APPENDIX B: Fourier Analysis of Experimental Data 182
BIOGRAPHY 184
-6-
LIST OF FIGURES
Figure Page
1.1 Geometry of long, narrow bay. 27
3.1 Regular perturbation theory: Frequency response
of large-scale harbor (h = 20 m., 2a = 100 m.,
L = 1000 m., A = .03).
(a) First harmonic amplitude at the back wall. 52
(b) Second harmonic amplitude at the back wall. 53
(c) Zeroth harmonic amplitude at the back wall. 54
3.2 Regular perturbation theory: First resonant mode
of large scale harbor (h = 20 m., 2a = 100 m.,
L = 1000 m., A1 = .03, = 1.41).
(a) Spatial variation of the harmonics. 55
(b) Surface profiles at times t = mn/4,
m = 0,l,2...7. 56
3.3 Regular perturbation theory: Second resonant mode
of large-scale harbor (h = 20 m., 2a = 100 m.,
L = 1000 m., A1 = .03, z = 4.36).
(a) Spatial variation of the harmonics. 57
(b) Surface profiles at times t = mn/4,
M = 0,l,2...7. 58
3.4 Regular perturbation theory: Third resonant mode
of large-scale harbor (h = 20 m., 2a = 100 m.,
L = 1000 m., A1 = .03, z = 7.36).
(a) Spatial variation of the harmonics. 59
(b) Surface profiles at times t = mr/4,
m = 0,1,2...7. 60
4.1 Experimental set-up. 76
4.2 Typical calibration curves of a wave gauge, at
two attenuator settings. 78
-7-
FigurePage
4.3 Theoretical study of multiple reflections from
wavemaker: Placement of images. 83
4.4 Experiment (-) vs. theory (-) for the wave ampli-
tude in the ocean (h = .5 ft., 2a = .33 ft.,
L = 1.211 ft., w = 4.067 sec-1, Lm = 31.17 ft.).
Top: harbor closed. Middle: harbor open.
Bottom: the radiated wave.
(a) A1 = .015 86
(b) A1 = .027 87
(c) A1 = .040 88
4.5 Entrance loss: Normalizd first harmonic ampli-
tude, |T/T0 1, vs. normalized friction factor, a. 100
5.1 Experiment vs. inviscid nonlinear theory (? = 0).
(A,+,x): measured, first, second, and third har-
monic amplitudes. (-): nonlinear theory. C----):
linear theory for first harmonic amplitude.
(a) L = 1.211 ft., A = .015 107
(b) L = 1.211 ft., A = .027 108
(c) L = 1.211 ft., A = .040 109
(d) L = 4.173 ft., A = .015 110
(e) L = 4.173 ft., A = .027 111
(f) L = 4.173 ft., A = .040 112
(g) L = 7.136 ft., A1 = .015 113
(h) L = 7.136 ft., A = .027 114
(i) L = 7.136 ft., A = .040 115
5.2 Experiment vs. inviscid nonlinear theory Cf = 0)
with exact dispersion relation. (A,+,x): mea-
sured first, second, and third harmonic amplitudes.
(-): nonlinear theory. (----): linear theory
for first harmonic amplitude.
(a) L = 4.173 ft., A = .015 116
-8-
Figure Page
5.2 (b) L = 4.173 ft., A1 = .027 117
Cc) L = 4.173 ft., A1 = .040 118
5.3 Experiment vs. nonlinear theory with separation
loss (f = .35). (A,+,x): measured first,
second, and third harmonic amplitudes. (-):
nonlinear theory. (----): linear theory for
first harmonic amplitude.
(a) L = 1.211 ft., A1 = .015 119
(b) L = 1.211 ft., A1 = .027 120
Cc) L = 1.211 ft., A1 = .040 121
(d) L = 4.173 ft., A1 = .015 122
(e) L = 4.173 ft., A = .027 123
(f) L = 4.173 ft., A = .040 124
(g) L = 7.136 ft., A = .015 125
(h) L = 7.136 ft., A = .027 126
(i) L = 7.136 ft., A = .040 127
5.4 Experiment vs. nonlinear theory with separation
loss (f = .35) and exact dispersion relation.
(A,+,x): measured first, second and third har-
monic amplitudes. (-): nonlinear theory.
(----): linear theory for first harmonic ampli-
tude.
(a) L = 4.173 ft., A = .015 128
(b) L = 4.173 ft., A = .027 129
Cc) L = 4.173 ft., A = .040 130
5.5 Inviscid nonlinear theory (f = 0): Frequency
response of large-scale harbor (h = 20 m.,
2a = 100 m., L = 1000 m., A1 = .03).
(a) Number of harmonics, N*, required in the
numerical solutions. 137
-9-
Figure Page
5.5 (b) First harmonic amplitude at the back wall
according to nonlinear theory (A) and
regular perturbation theory (-). 138
(c) Second harmonic amplitude at the back wall
according to nonlinear thoery (+) and
regular perturbation theory (-). 139
(d) Third harmonic amplitude at the back wall
according to nonlinear theory (x). 140
(e) Zeroth harmonic at the back wall according
to nonlinear theory (A) and regular pertur-
bation theory (-). 141
5.6 Inviscid nonlinear theory C? = 0): First reso-
nant mode of large scale harbor (h = 20 m.,
2a = 100 m., L = 1000 m., A = .03, k = 1.41).
(a) Spatial variation of the lower harmonics. 142
(b) Spatial variation of the higher harmonics. 143
(c) Surface profiles at times t = m7/4,
m = 0,1,2...7. 144
5.7 Inviscid nonlinear theory (f = 0): Second
resonant mode of large-scale harbor (h 20 m.,
2a = 100 m., L = 1000 m., A1 = .03, z = 4.36).
(a) Spatial variation of the lower harmonics. 145
(b) Spatial variation of the higher harmonics. 146
(c) Surface profiles at times t = m/4,
m = 0,1,2...7. 147
5.8 Inviscid nonlinear theory C? = 0): Third reso-
nant mode of large-scale harbor (h = 20 m., 2a =
100 m., L = 1000 m., A1 = .03, z =7.36).
(a) Spatial variation of the lower harmonics. 148
(b) Spatial variation of the higher harmonics. 149
(c) Surface profiles at times t = mn/4,
m = 0,1,2...7. 150
- 10 -
Figure Page
5.9 First harmonic amplitude at the back wall of the
large-scale harbor vs. incident wave amplitude.
(----): linear theory. (A): inviscid nonlinear
theory (? = 0). (x): nonlinear theory with separa-
tion loss (f = .35). 151
5.10 Nonlinear theory with separation loss (? = .35):
First resonant mode of large-scale harbor (h = 20 m.,
2a = 100 m., L = 1000 m., A1 = .03, z = 1.41).
(a) Spatial variation of the lower harmonics. 152
(b) Spatial variation of the higher harmonics. 153
(c) Surface profiles at times t = mnr/4,
m = 0,l,2...7. 154
5.11 Nonlinear theory with separation loss C? = .35):
Second resonant mode of large-scale harbor (h = 20 m.,
2a = 100 m., L = 1000 m., A1 = .03, = 4.36).
(a) Spatial variation of the lower harmonics. 155
(b) Spatial variation of the higher harmonics. 156
(c) Surface profiles at times t = mnTr/4, m = 0,1,2...7. 157
5.12 Nonlinear theory with separation loss (f?= .35):
Third resonant mode of large-scale harbor (h = 20 m.,
2a = 100 m., L = 1000 m., A = .03, = 7.36).
(a) Spatial variation of the lower harmonics. 158
(b) Spatial variation of the higher harmonics. 159
(c) Surface profiles at times t = mir/4,
m = 0,,2...7. 160
5.13 Response of large-scale harbor to nonmonochromatic
incident wave (h = 20 m., 2a = 100 m., L = 1000 m.,
A2 = .050, A3 = .015, = 1.41). 161
- 11 -
LIST OF TABLES
Table Page
2.1 Orders of Magnitude of Incident, Reflected, and
Radiated Waves in the Ocean 40
5.1 Composition of the Incident Waves Used in the
Experiments 102
5.2 Exact and Approximate (Boussinesq) Wavenumbers
in the Experiments 103
5.3 Resonant Modes of the Large Scale Harbor According
to Linearized Theory 131
- 12 -
CAST OF CHARACTERS
Symbol
a Half-width of narrow bay
A Twice the incident wave amplitude
en(x) Basis function in the numerical scheme
E Energy density averaged over time
fJ'f ,2
Friction parameters
g Acceleration of gravity
G Green's function
h Depth of water
H Interval size in finite difference mesh
HO Zeroth order Hankel function of first kind0
k Wavenumber
kc Minimum wavenumber for cross modes
9,L Dimensionless and dimensional harbor length
MLM Dimensionless and dimensional length of model ocean
( )n the n-th harmonic of ( )
n S A vector normal to surface S
N* Number of harmonics present in a numerical solution
NPT Number of points in the finite difference mesh
p Pressure
P R Radiated power
Pv Power dissipated is viscous boundary layer
Q Source strength for radiated wave
- 13 -
Symbol
r Radial coordinate
Rc Critical Reynolds number
Rns' 9ns Coupling constants in the nonlinear equations
t Time coordinate
T,T0  Coefficient of first harmonic, in linearized theory
u Velocity vector
U Horizontal velocity in harbor entrance
x Horizontal coordinate
x* Near field x-coordinate
x A slow scale in x, e.g. x = EX
y Horizontal coordinate
y* Near field y-coordinate
z Vertical coordinate
Z Radiation impedance
cx Normalized separation loss parameter
6 Dimensionless half-width of narrow bay
6 v Boundary layer thickness
V Two-dimensional gradient operator =
s Nonlinearity parameter
TI Surface elevation
A A relaxation parameter; also wavelength
P2 Dispersion parameter
v Viscosity
p Fluid density
Velocity potential; also a phase angle
W Circular frequency -14 -
CHAPTER 1
INTRODUCTION
1.1 Description of the Problem and Review of Literature
A harbor is a partially enclosed body of water which is joined to
the sea by a relatively narrow entrance. The entrance, which is quite
often a man-made system of breakwaters, serves to shield the harbor from
the most common disturbances, namely wind-generated waves of the open
sea. Because it limits the amount of energy exchange between the harbor
and the sea, a narrow entrance occasionally has the damaging effect of
trapping long waves, which may be originated by tsunamis, surf beats or
other sources, and of creating within the harbor large oscillations which
may persist for days. This is the phenomenon of harbor resonance or
"seiching." The deleterious effects of these oscillations are well-
known--ships are torn loose from their moorings causing damage to dock
facilities and neighboring vessels, and strong currents in the harbor
entrance prevent the safe passage of incoming and outgoing ships, cf.
Wilson (1957).
From a theoretical standpoint, the study of harbor oscillations
begins with the determination of the frequencies at which resonance
occurs and the wave profiles of the resonant modes. Much of the work to
date uses linearized, inviscid wave theory, in which the wave amplitude
is presumed so small that nonlinear effects are negligible. Rectangular
harbors were studied by Miles and Munk (1961), Ippen and Raichlen (1962),
Ippen and Goda (1963), Raichlen and Ippen (1965), and Mei and Unluata
- 15 -
(1976). Miles and Munk discovered the peculiar phenomenon, referred to
as the "harbor paradox," that the amplitude of monochromatic resonant
oscillations actually increases with decreasing width of the harbor
entrance. For harbors of arbitrary shape, numerical solution procedures
have been developed by Hwang and Tuck (1970), Lee (1971), and Chen and
Mei (1974).
The validity of all treatments based upon linearized wave theory is
necessarily limited by the assumption of small wve amplitudes. Fre-
quently, the resonated wave heights predicted by these theories are so
large that it is no longer justifiable to neglect nonlinear effects.
Some attempt to include nonlinear effects in the analysis of harbor
resonance was made by Gaillard (1960) and Biesel (1963), but complete
solutions were not obtained. The general properties of nonlinear waves,
and of water waves in particular, have been studied extensively by
Phillips (1960), Longuet-Higgins (1962), Benney and Luke (1964),
Whitham (1965, 1967 a,b), Benjamin and Feir (1967), Chu and Mei (1970),
Kim and Hanratty (1971), Mei and Unluata (1972), Bryant (1972), etc.
A review article by Phillips (1974) is especially useful. For the most
part, these authors restrict their attention to progressive waves travel-
ling in an infinite medium having no obstacles in the path of propaga-
tion. Finite amplitude standing waves in a basin have been treated by
Tadjbakhsh and Keller (1960), Benney and Niell (1962), and by Verhagen
and van Wijngaarden (1965), who also performed several experiments.
Laboratory investigations of harbor resonance have been hampered by
several difficulties, among them inadequate modelling of the ocean
- 16 -
domain and Reynolds number dissimilarity. Biesel (1954) discusses the
similitude of scale models for the study of harbor oscillations. Ippen
and Raichlen (1962) and Raichlen and Ippen (1965) investigate the prob-
lem of coupling between a large but finite ocean basin and a smaller
harbor basin. Ippen and Goda (1963) and Lee (1971) performed experi-
ments to corroborate their theoretical work on harbor resonance. Much
of this experimental work is indirectly predicated on the assumptions
of linearized theory. For example, it is assumed that the response of
a harbor to monochromatic incident waves is itself monochromatic so that
wave records need not be Fourier analyzed. Further, the experiments are
often conveniently performed in deep water, because the extrapolation
to shallow water is trivial in the context of linearized theory. For
nonlinear theory, this extrapolation cannot be made.
1.2 Elements of Water Wave Theory
In a wide variety of water wave problems, the fluid motion is
essentially incompressible, inviscid, and irrotational. For constant
depth, the governing equations are
V2 +ezz = 0 in fluid -h < z <_ n(x,y,t)
z=0 on z=-h
t + 4'xx + yfny = z on z = n(x,y,t)
+ 12 + (vfl2) + gn = 0 on z = n(x,y,t) (1.2.1)
where
- 17 -
n(x,y,t) = displacement of the free surface from its undisturbed
position
U(x,y,z,t) = (y ,z) = fluid velocity
(x,y,z,t) = velocity potential
g = acceleration of gravity
h = water depth
Subscripts are used to denote differentiation. (1.2.la) is the conse-
quence of continuity for an incompressible fluid and irrotationality;
(1.2.lb) states that the vertical component of velocity must vanish at
the bottom; (1.2.lc) is the kinematic free surface boundary condition
which ensures that particles on the free surface move tangentially to
it; and (1.2.ld) is the dynamic boundary condition stating that atmos-
pheric pressure is constant (e.g., zero) on the free surface and that
surface tension is negligible. Note that the two boundary conditions
on the free surface are nonlinear because both p and n are unknown.
In linearized wave theory, the wave amplitude is assumed to be
infinitesimal and (1.2.lc,d) are replaced by
1t ~1 z =0 on z = 0
t + gn = 0 on z = 0 (1.2.2)
In an infinite ocean of constant depth h, one solution to (1.2.la,b)
and (1.2.2a,b) is the monochromatic progressive wave
n(x,y,t) = Re {Aei(kx - wt)I
p(x,y,z,t) = Re {-igA cosh k(z+h) ei(kx - wt)} (1.2.3)W cosh kh
- 18 -
where
W2h/g = kh tanh kh =..1p2 (1.2.4)
is the dispersion relation. In deep water, p 2 >> 1, the dispersion
relation is approximately w2h/g = kh, and the phase velocity is twice
the group velocity. In shallow water, p 2 «<1,
2h/g = (kh) 2[1 - ]-(kh) 2 + (kh) 4  (kh) 2  (1.2.5)3]
the dispersion relation is approximately linear; the phase and group
velocities are both equal to gi (1 + O(kh)2); and the waves are said
to be non-dispersive. For intermediate values of p, w is a transcen-
dental function of kh, and the waves are highly dispersive.
When the assumption of infinitesimal wave amplitudes is removed,
the nonlinearity parameter
E_= A/h = wave height/water depth
is also needed to characterize the wave. Different wave motions are
2possible depending on the relative magnitudes of E and p2. In deep
water, the combination sp 2, which is a measure of the slope of the
free surface, is found to be of great importance in determining when a
wave will break. In shallow water, the combination E/ 2, known as
Stokes' number, is used to identify three qualitatively different
kinds of waves:
/P2 << 1 : Linear, dispersive waves (Linearized theory)
6/12 n 0(1) : Nonlinear dispersive waves, including cnoidal and
solitary waves (Boussinesq and Korteweg-de Vries
equations)
- 19 -
E/y2 >> : Nonlinear, nondispersive waves (Airy's equations)
In treating harbor resonance, we deal almost exclusively with
shallow water, or "long" waves. This is because the wavelengths of
resonant harbor oscillations are on the scale of the horizontal dimen-
sions of the harbor, which are normally much greater than the water
depth. For example, a harbor might be two kilometers long and only
twenty meters deep, so that p 2 10~4 for the lower resonant modes. It
is useful then to have a theory in which the assumption of shallow water
is built-in, but the assumption of infinitesimal wave amplitudes is not.
To derive such a theory, we introduce scaled variables so that
the relative magnitudes of terms in the equations of motion appear ex-
plicitly. For a shallow water wave of frequency w, an appropriate
choice of variables is
t' = ot
(x',y' = tO (xy)
Z' = z/h (1.2.6)
-n =/h
U' = u//gh
$' = o4/gh
where primes denote dimensionless variables. Henceforth, all variables
are dimensionless, and for convenience, the primes will be omitted.
Substituting (1.2.6) into (1.2.1), the equations of motion take the
dimensionless form:
2 1v P (- <z<n)
1y
- 20 -
(z=0 on z = -1
at +$' +cyn=ynz on z = n (1.2.7)
$0 + 2) O + 02 + fl= 0 on z=2
2w 2  -+ a- a a 2  2
where p2 wh/g, V 2 ( , ), and V s . Expanding
a ayax ay
c(x,y,z,t) = X 2n 0(n)(xyzt)(1.2.8)
n=o
we find from (1.2.7a) that
0(0)=0zz 0(1.2.9)
(n) =-V2 (n-1) n = 1,2,3...
zz
Integrating twice with respect to z and applying the bottom boundary
condition (1.2.7b), we have
(0) (x,y,z,t) = (0(x'y't) independent of z
(n) (-1)n(z+1)2n 2n (0)p (x,y,z,t) = (2n)! n() x,y,t)
n = 1,2,3... (1.2.10)
where integration constants are independent of z and have been grouped
into (O(. From (1.2.8)
Cpxyzt O (-1 )n 2n z+l)2nV 2n (0)(xyt(xnzt)= ( (x,y,t) (1.2.11)
n=o
Note that the leading term in the expansion for the vertical velocity
(z is 0(.P2
So far, we have made no assumption as to the magnitude of c(O) and
- 21 -
no use of the nonlinear boundary conditions on the free surface,
(1.2.7c,d). Denoting the magnitude of (0) by e, several approxima-
tions to the free surface boundary conditions are possible, depending
on the relative magnitudes of e and p2 and the degree of accuracy
desired. In Airy's approximation, for example, pc2 << and terms of
0( 2) are neglected. In the Boussinesq approximation, e «2 << 1 and
2 2 4
relative errors of 0(E ,pE ,4 ) are neglected. A detailed treatment
of the various types of approximations and their domains of validity is
found in Peregrine (1972).
The Boussinesq approximation has several advantages over Airy's
approximation. First, Airy's equations can be formally derived from
the Boussinesq equations by taking the limit p2 - 0. Secondly, it is
well known that dispersion tends to counteract the wave steepening
effect of nonlinearity, and thus to inhibit wave breaking. Airy's
equations contain no dispersive terms and lead to shock formation even
in physical situations where shocks are not observed experimentally.
By contrast, the Boussinesq equations, which contain the leading order
effect of dispersion, are found to agree well with experiments so long
as the ratio of peak-to-trough wave height to water depth is less than
.5. This is quite close to the onset of breaking, which occurs in
shallow water when the above ratio is roughly .7. Thus, the Boussinesq
approximation has a greater range of validity, so long as breaking
does not occur.
Substituting (1.2.11) into (1.2.7c,d), and neglecting relative
2 2 4
errors of 0(S ,EP ,s ), one arrives, after some algebra, at the follow-
ing form for the Boussinesq equations: -
- 22 -
"t + V.u + V.(nt) = 0
ut+ ÷ l+y 2) I 2tt =(1.2.12)
where t(x,y,t) is the depth-averaged horizontal velocity, defined by
iu(x,y,t) = 1: J$ dz
-l
÷(0) 2 2-+( 2
~-i() 2 ( )) (Z + 1)2dz
-l
vq(O) -112 V2(vqi 0)) (1.2.13)
In the limit p2-+-0, (1.2.12) reduces to Airy's equations.
The fluid pressure, which is found by integrating the vertical
equation over z, is given by (pgh)p where p is the density of the fluid
and
2 2 4-p = (n - z) + p2(z + z2/2)(7.ut) (1.2.14)
The term (n - z) is simply the hydrostatic pressure beneath the wave;
the additional term reflects the influence of vertical motion, which is
0(20(i 2).
One consequence of the nonlinearity of (1.2.12) is that monochroma-
tic waves are no longer possible solutions. If one begins with a wave
of the form (n1(x,y), t1(x,y))e-it, one finds that the nonlinear terms
generate second harmonics, having the time dependence e 2it. Interac-
tion of the first and second harmonics produces third harmonics, having
the time dependence e-3it. This process, if continued, leads to a
time periodic wave of the form:
- 23 -
TI(xgy gt)= Re 00 n , - int
I (x,y,t) n=0 u nxy -
R O fLn ]e-flt (1.2.15)
n=o Ln(xYJ
where, by convention, (n , it - ( n* ). Substituting (1.2.15)n -n = n n ~ siutn 1..5
into (1.2.12), the n-th harmonic is found to satisfy
-inan + i~ ~I (sC-O
1 2 ,) +
-in u +-(l P2 n2)$n +sU- s n-s=u n(1.2.16)
Combining the divergence of (1.2.16b) with (1.2.16a) and neglecting
2 2 4relative errors of O(E ,et ,p2 ),
(2+ k =)nn sn-s 2 usa n-s) (1.2.17)
s=-o
where
2 2 1 2 2
kn = n/(l - Pyn)
= n2.[ + n2V2 + 0(ny)4] (1.2.18)
is the dispersion relation in dimensionless form. Note that the dis-
persion relation may be written as
n2 = k ( - 4-k2A + O(k-P) 4)
k 4
-n- tanh k p + O(k k) (1.2.19)
P n n
which shows clearly itsconnection to the exact dispersion relation
- 24 -
(1.2.4). The error in the Boussinesq dispersion relation is 0(k n)
4
and can be quite large for the higher harmonics.
The near linearity of the dispersion curve when p2 << I leads to
the phenomenon of resonant linear interactions, which has been studied
by Kim and Hanratty (1971), Mei and Unluata (1972), and Bryant (1972).
Consider two plane waves
n1 = Re {EA1 e X(k1I - n2t) I (n1,9k1 > 0)
(1.2.20)
n2 = Re {EA2 ei(k2x - n2t)} (n2, k2 > 0)
where, from (1.2.17),
k1 = n1 + 0( 2)
(1.2.21)
k2 = n2 + o(V2
Through nonlinear interaction, a third wave is formed
T3 = Re {EA3 ei[(k1+k2)x - (n+n2)t] + 
ei[(k1-k2)x- (n1-n2)t
(1.2.22)
Since k1 +k2 = n1 -1+ n2 + 0(p2 ) n3 represents a propagating wave which
travels with approximately the same phase velocity as n and n2; thus,
the three waves maintain the same relative phases as they propagate.
The magnitude of n3 is found to increase linearly with distance so long
as the three waves are in phase and jn3j « 1'n111n21. This is known
as a resonant three-wave interaction. Eventually, the influence of
dispersion limits the distance of coherent propagation to 0(1/u2). The
net result is that nonlinear interactions of 0(s2) add constructively
-over a distance of 0(1/p2), and the magnitudes of sA3 and EB3 are
- 25 -
found to be sO(s/p2). The Stokes parameter, E/p2, therefore governs
the extent of harmonic generation. The subject of resonant nonlinear
interactions is discussed at length in Chapter 3.
When e/p2 << 1, higher harmonics are generated in profusion, and
the horizontal length scales used in (1.2.6) are no longer appropriate.
This is the situation just before breaking. The Boussinesq equations
cannot be expected to apply to such cases because the scales used in
their derivation are no longer the correct physical scales.
1.3 Linearized Theory of the Long Narrow Bay
The harbor we shall consider is a long, narrow bay of width 2a,
length L, and constant depth h. The geometry of the ocean-harbor system
is shown in Fig. 1.1. The coastline and harbor walls are taken to be
rigid, impermeable, vertical surfaces along which the normal component
of velocity must vanish, viz.
u=xgySt)-n= 0 (1.3.1)
where (x5,yS) is a point on the surface, and 'As is a vector normal to
the surface.
The presence of finite boundaries in the fluid introduces new
length scales, in addition to those which characterize the wave length
and amplitude. We shall see that the dimensionless parameters
2 L(w/gii)
(1.3.2)
6 2 vr-h
- 26 -
L'
a f & & & 1 .0 A 4
- .-e--- -- 4
Ft' F cr N ~ (7 FFUFF grr ,j UFV'
4
HARBOR
INC-tDENT
WAVE
OCSAN
Fig. 1.1. Geometry of long, narrow bay
07c.
I
determine the linear resonance properties of the harbor basin. For a
narrow bay (a/L << 1), resonant oscillations typically occur when
z = 0(1) and 6 = t-(a/L) << 1.
Because 6 << 1, the ocean-harbor system can be split into three
distinct regions: the far field of the ocean, r =/x2 + Y2 >> 6; the
near field of the entrance, r r 0(6); and the far field of the harbor,
-z < x << -6. In the far field of the ocean, the fluid motion varies
on the scale of the wavelength and (x,y) are suitable horizontal coord-
inates. In the near field of the junction, the scale of motion is 6
and the coordinates
(x*,y*) = (x,y) / 6 (1.3.3)
are most appropriate. Finally, in the far field of the harbor, the
natural choice for coordinates is (x, y*).
Before embarking on a nonlinear treatment of the resonance proper-
ties of a long, narrow bay, we summarize known results from linearized
theory. The geometry of Fig. 1.1 has been studied by Miles and Munk
(1960) and Unluata and Mei (1973). Starting from linearized water
wave theory, Unluata and Mei obtain asymptotic expansions for the wave
field in each of the three regions. By matching the asymptotic solutions
in adjacent regions, they determine the response of the harbor in terms
of the incident wave amplitude and frequency. For an incident wave of
the form
n (x,y,t) = Re(4A e i(k1x-t) (1.3.4)
where k1 is the dimensionless wave number satisfying the dispersion
relation (1.2.19), the response in the far field of the harbor is found
- 28 -
to be
H (xt) = Re (T cos k1 (x + z)e-it
uH(x,t) = Re ( T sin k 1(x + z)e-it
k11
where
T E A[cos k I- sin k k]~(1.3.6)A[o k1  k1
and
(k1 6)[l + -in k16 + zn )] (1.3.7)
zny E .577216 = Euler's constant
Z is the radiation impedance of the harbor entrance; its real and
imaginary parts are the radiation damping and mass reactance, respective-
ly.
Letting Z1 = Z R + iZ 1P, we have from (1.3.6),
Max nH(x,t) = jT1 I
x,t
- A1 I[(cos kz + Zjj sin k1k)2 + (ZiR sin k 1)21-1/2
(1.3.8)
When kI satisfies the "resonance condition"
cos k1zk + Z11 sin k1 z = 0 (1.3.9)
k
(1.3.8) yields
- 29 -
Max n(xt) = ZR IAI kT (1.3.10)
x~t I1Rsn 1;
k,
so that the incident wave amplitude is amplified by a factor proportional
to 1/ZlR, which is in turn proportional to 1/6. As the width of the
harbor entrance is made smaller, 6 decreases and the resonant amplifica-
tion increases; this is the so-called "harbor paradox.'
Since Z 1 ' u 0(6 zn) << 1, the resonance condition (1.3.9) can be
written as
cos (k1 Z - Z 1) 0 (1.3.11)
F-k1
which yields the resonance criterion
k1  = (m + 1)r +-Z m = 0,1,2... (1.3.12)
1
The lowest mode, m = 0, is often called the quarter-wave mode because
1]_ 2w7T 1kz i/2 or z~ (A)(T) = 4-, where xl is the dimensionless wave-
length.
Notice that, in (1.3.5)., nH and uH are independent of y; that is,
the far field of the long, narrow bay is essentially one-dimensional.
This allows a great simplification in the nonlinear treatment of the
harbor, as will be seen in Chapter 3.
1.4 Scope of this Investigation '
The theoretical part of this study analyzes the response of a long
- 30 -
narrow harbor to time-periodic incident waves. In Chapter 2, the wave
field of the ocean is shown to be of sufficiently small amplitude that
it may be treated as linear. Assuming the entrance width to be small,
e.g. << 1, the interaction of the ocean with the harbor is then
approximated by an impedance boundary condition applied at the harbor
entrance. The advantage of this approach lies in its mathematical
simplicity; its disadvantage is that higher harmonics have larger
values of knS and are therefore less accurately represented. Where
greater accuracy is required, a semi-numerical method, such as that of
Lee (1971) or Chen and Mei (1974), should be used. In Chapter 3, the
nonlinear response of a long, narrow harbor is examined. With the
assumption 6 << 1, the wave field far from the harbor entrance is
shown to be essentially one-dimensional, allowing a simplification of
the nonlinear equations. The field variables are Fourier analyzed into
time harmonics, and an infinite set of coupled, second-order, ordinary
differential equations is found to govern the spatial variation of the
complex harmonic amplitudes. Perturbation techniques and direct numeri-
cal methods are used to solve these equations.
Laboratory experiments were carried out using a model harbor placed
in a large shallow wave basin. Chapter 4 deals with the experimental
apparatus and procedures, and with the effects of friction and multiple
reflections from the walls of the wave basin.
Results of the theoretical and experimental work are compared in
Chapter 5. In addition, some theoretical predictions for large scale
harbors are presented. A summary of the main conclusions and suggestions
for future research are found in Chapter 6.
- 31 -
CHAPTER II
WAVE SYSTEM OF THE OCEAN
2.1 Motivation for a Linear Treatment
The wave field of the ocean consists of three parts; the incident
wave, n(x,y,t); the wave reflected from the coastline, nr(x,y,t); and
the wave radiated away from the harbor entrance, nR(x,y,t). To a first
approximation, the total wave field, nT(x,y,t), is simply a linear
superposition of these three, viz.
nT i + nr + nR (2.1.1)
+T
The total velocity field, u (x,y,t), satisfies a similar relation.
In analyzing the ocean, we begin with the linear approach of (2.1.1).
After the linear solution is obtained, the magnitude of the nonlinear
terms will be estimated, and the assumption of linearity will be reexam-
ined. We shall find that, when the harbor entrance is narrow (6 << 1),
the radiation damping is 0(6). It is then possible for an incident
wave of O(sa) to excite a resonant harbor response of 0(E). To determine
a solution that is accurate to 0(62), one need only account for nonlinear
effects in the harbor, not in the ocean. This is the motivation for a
linear treatment of the ocean.
Deleting nonlinearterms in (1.2.16) and (1.2.17), we have for the
n-th harmonic
-inpn + v-u = 0
n n
-in n + nn = 0 (2.1.2)
n
- 32 -
and
(72 + k2))n = 0 (2.1.3)nfn
2 2 1 2 2
where k2 = n (1 + n P ). Note that n (x,y) is proportional to then 3n
velocity potential for t (x,y). From (2.1.2b), $n = 0 so that n (x,y)veoiyptnilfrun 0 0
is a constant everywhere in the ocean. By properly defining the water
depth, the constant can be set to zero, to the leading order, without
loss of generality. The term t0(x,y) represents a steady current and
will be presumed to be zero, to the leading order. Thus, the zeroth
harmonic may be deleted.
2.2 The Radiated Wave
If the harbor entrance were sealed, the boundary condition at the
coastline would be
T
u= .(o,y) = < y<) (2.2.1)
k
n
and the total wave field would be given by
Tii
n (xy) = n (xy) + ny(-x,y) (2.2.2)
where n1(-x,y) is simply the reflection of the incident wave about then
y-axis. The function n4(x,y) can be any solution of the Helmholtz
equation in the domain x > 0, - co < y < o. For example, n(x,y)
might be an obliquely incident wave
-iaX' iyn ( 2 n 2)
r$(x,y) = cne nXe (n'n > 0, an +$g 2 = k2 (2.2.3)
- 33 -
in which case
T nisTn(x,y) = 2cn cos ax e (2.2.4)
represents a standing wave whose amplitude varies sinusoidally in the
y direction.
With the harbor entrance open, the boundary condition on the y-axis
is more complicated:
TT= ian T(o,y) 0 jyj > 6
ux(oy)x 1 (2.2.5)
kn Un(y) ly < 6
where Un(y) is the unknown velocity distribution in the harbor entrance.
The total wave field
nT(x,y) = n (xy) + n (-xy) + n (x,y) (2.2.6)
contains a radiated wave component which is the solution to the
following boundary value problem:
(V2 + k2)n = 0n n
an (o,y) 0 |yl > 6
n x f .A2 (2.2.7)9X ik2
nnU(Y) lyI < 6
RR
lim 1,r ('-n ik ) = 0ar n n
where r x2 + y Equation (2.2.7c), known as the radiation condi-
tion, ensures that the radiated wave will propagate away from the har-
bor entrance, as required by causality. In essence, the radiation
- 34 -
condition replaces the initial conditions that would be present in a
fully time-dependent treatment of wave scattering. For a detailed
derivation of the radiation condition, see Stoker (1957), pp. 174-181.
A formal solution for the radiated wave can be constructed using
Green's functions. The appropriate Green's function, which satisfies
(2.2.7) with Un(y) replaced by the unit singularity 6(y-y'), is found
to be
Gn(xyy') = ( n2)( H x2 + (-y'))) (2.2.8)
(1)where H is the zeroth order Hankel function of the first kind.
Since
6
Un(y) =f Un(y) 6(y-y') dy' (2.2.9)
-6
Rthe solution for qn(x,y) is given by
k2 6
(xyH(k + (-y) 2)U (y)dy' (2.2.10)n2n 0on / + __)ny
-6
Let us examine (2.2.10) in the "far field" of the ocean, that is,
2 22
where r = x + y2 >> 6. Substituting the Taylor series
(1) (1)
H)(kn x2 + (-y1 2) = H 1 (knr) + y'[Hl + 2.2.11)
(2.2.11)
into (2.2.10), we have
R k2  6 aH(1) 6
n (xy) = --- {H)(knr)f U(y')dy' + [0 ] { y'U=(y'dy'+...
-6 -6
(2.2.12)
- 35 -
The first term on the right hand side represents the field of a mono-
6
pole whose strength is proportional to the total flux, {Un(y)dy. The
second term is that of a dipole whose strength is proportional to the
6
dipole moment of the velocity distribution, f YUn(y)dy. When 6 << 1,
-6
the dipole term is smaller than the monopole term by a factor of 6, and
higher terms in the multipole expansion are still smaller. Defining
tte average velocity by
U n Un(y)dy' (2.2.13)
-6
we have
(-~k( ) kn6U )Hl) (knr)-[l + 0(6)] (2.2.14)
Note that, if U is 0(s), the radiated wave is O(skn6) in the farn n
field. We turn our attention now to the near field of the harbor
entrance.
2.3 The Impedance Boundary Condition
In the near field of the harbor entrance, r is 0(S) and, assuming
k n6 << 1,
H l)(knx2 + (y-)y) ) 1 + 2n( 2 + 7-y-)2) + O(kn6znk6)2
(2.3.1)
where any 2 Euler's constant = .57722... . Substitution of (2.3.1) into
- 36 -
(2.2.10) yields
k n yT, (x3,Y) ~- 1I+ in( /x2 + G-y_ )2)U (y')dy'* (2.3.2)
-6
Defining the average wave amplitude in the entrance by
6
Tin {R(o,y)dy (2.3.3)
-6
we have
2
= {I + n( f )(y)dy'y(2.3.4)n n n 1T 26jf j 2 ''n
-6 -6
SZ n Un
where
k2C 6 yk U (y')
Z n = ( 2 ){l J{n( 2ny-y') dy'dy} (2.3.5)n 
-2 6 -6n
is called the "radiation impedance" of the ocean-harbor junction. The
real and imaginary parts of Zn are the radiation resistance and mass
reactance of the junction, respectively. Note that only the mass
reactance depends upon the details of the velocity profile, Un(y).
Introducing the variable v = y/6, and using (2.2.13), we arrive at the
following expression for the radiation impedance:
k .
Zn ( f)(kn6) [1 + -v (Zn kn6 + cn)] (2.3.6)
where
1 1
cn =zin + < nlv-v'(Un(V')/n)dv'dv (2.3.7)
-37 -
Observe that the radiation resistance and mass reactance are O(kn6)
and O(kn6znk n6), respectively. In deriving (2.3.6), we assumed that
k 6 << 1. This assumption is progressively worse for the higher har-
monics, which have larger value of k n6. We accept this deterioration
of accuracy in the hope that the magnitudes of the higher harmonics
are progressively smaller. Where greater accuracy is required in treat-
ing the high harmonics, the impedance approach can be abandoned in
favor of a semi-numerical method such as that of Lee (1971) or Chen
and Mei (1974). To predict the essential features, we choose to avoid
these elaborate alternatives.
The form of the velocity profile, Un(y), depends upon the geometry
of the harbor entrance and upon the wave field in the harbor, which is,
of course, still undetermined. In the linearized theory of harbor
resonance, the wave field in the harbor is expressible as an integral
over the Green's function for the harbor, in exact analogy with (2.2.10).
The requirement that rn(x,y) be continuous at x=o leads to a linear
integral equation for Un(y). The integral equation can then be solved
by variational approximation or a numerical quadrature procedure. This
approach is only feasible when the harbor is treated as linear, because
the use of Green's functions relies upon the principle of superposition.
To solve the nonlinear problem, a different approach must be found.
When flow spearation at sharp corners is not significant, one
may exploit the assumption that kn << I and determine an approximate
form for the velocity profile, without prior knowledge of the wave
field of the harbor. This has been recognized long ago by Rayleigh and
is easily demonstrated as follows: In the near field of the junction,
- 38 -
the appropriate length scale is 6, and the terms in the Helmholtz
equation are in the ratio,
(k )n v2) ' 2 O(kn6)2 (2.3.8)'nan''n n
It follows that
V2 n = 0 + O(kn6)2 (2.3.9)
so that Laplace's equation governs the motion in the near field with
a relative error of O(kn&)2. This is called the "quasi-static"
approximation. The name derives from the fact that an observer in the
near field sees steady flow and is unaware that it is being driven by
oscillatory wave motion.
To solve Laplace's equation, it is only necessary to specify the
geometry of the harbor entrance; the details of the adjacent wave
fields are irrelevant. For the geometry of Fig. 1.1, the potential
flow field can be found by applying a Schwartz-Christoffel transforma-
tion. Details are given in Unluata and Mei (1973). The radiation
impedance is then found to be
Z2 . )f(kl6)[l +- nkn6 + zn ](2.3.10)
where Zny = Euler's constant = .5772.
From (2.2.6), the average total wave amplitude at x=o is given by
T= 2H1(0) + = A + Zn U (2.3.11)
n n n nn
where
- 39 -
A =2n 1(O) 2n i r(o,y)dy (2.3.12)
-6
Equation (2.3.11) will be referred to as the "impedance boundary condi-
tion." Observe that it was derived without even specifying the geometry
of the harbor basin.
2.4 On the Validity of the Linear Approximation
Having obtained the linear solution for the wave field in the
ocean, we now estimate the magnitude of the neglected nonlinear terms.
The following table shows the orders of magnitude of the n-th harmonic
incident, reflected, and radiated waves in both the near and far fields
of the ocean.
Table 2.1
i r R
n n n
Far Field A A k 6U
n n n n
Near Field A A (k 6znk 6)U
n n n n n
The magnitude of Un/An is thus far arbitrary.
In section 1.3, Eq. (1.3.6), we have seen that, when the n-th
harmonic is resonated,
Un/An 10(g16) (2.4.1)
n
This result can be anticipated if one thinks of the harbor as a damped
- 40 -
harmonic oscillator, with damping coefficient k 6, being driven by
n
an externally applied force of magnitude An. The response of the
nn
oscillator U at its natural frequency, is then given by (2.4.1).
Of course, a harbor is a continuous system and therefore has many
resonant frequencies, rather than just one.
In analyzing the harbor, we shall require that the maximum
value of U n be 0(c). Equation (2.4.1) then constrains An to be O(cekn6).
From Table A, the largest wave amplitude in the ocean occurs in the
near field, whereTi is O(ekn6znkn6). Thus, the largest nonlinearn n n
terms are O(ekn6znkn 6)2 in the ocean and 0(c2) in the harbor. It is
consistent to treat the harbor as nonlinear and the ocean as linear
only if
(skn6znkn6)2 «E2 (2.24.2)
which is the case whenever kn6 << 1. The conclusion, therefore, is
that, so long as kn << 1, it is consistent to neglect nonlinear effects
in the ocean, and to include leading order nonlinear effects in the
harbor.
- 41 -
CHAPTER III
RESONANT OSCILLATIONS OF A LONG, NARROW BAY
3.1 Reduction to One Space Dimension
The harbor domain, -k < x < 0, consists of two regions: the near
field, lx| O 0(6), and the far field, -z < x << -6. We shall show that
the waves in the far field are essentially plane waves varying only in
the x-direction. Intuitively, this is because oscillations in the y-
direction have a minimum cutoff frequency given by a = 7r/2. Since
the frequencies under investigation are far below this cutoff, the y-
dependence of the near field disturbance cannot propagate outward into
the far field of the harbor.
A more rigorous demonstration of the fact that, for low frequency
oscillations, the far field of the harbor is one-dimensional, proceeds
as follows: Let us define a velocity potential q(x,y,t) for the
depth-averaged horizontal velocity u(x,y,t) by the relation u E v.
Integrating the momentum equation (1.2.12b) over space, the Boussinesq
equations take the form
t + VY.[(1 + T)0V] = 0
$t+ n +12tt +1V 2 =C(t) (3.1.1)
The integration constant C(t) can be absorbed into the term it by re-
defining c; hence we may set C(t) equal to zero, without loss of
generality. The boundary conditions in the far field of the bay are
- 42 -
x (-z,y,t) = 0
y (x,6,t) = 0 (3.1.2.)
y)(x,-6,t) = 0
(3.1.2b,c) suggest that we rescale y, viz.
y* = y/6
From (3.1.1), with C(t) = 0, we then have
((1 +nf*)*2 = 62{t + [0 +n)( ] }
(3.1 .3)
1 2 = 2 + 1 2 + 12
7 * = -62 ft + n + } V 24tt + x
Integrating (3.1.3a) over y* and using the boundary conditions (3.1.2b,c)
we find
4Y)(xy*,t) = 0 + 0(62) (3.1.4)
so that, to 0(62), (x,y*,t) is independent of y*. From (3.1.1b), it
follows that n(x,y*,t) is also independent of y*, to 0(62).
We conclude, therefore, that to 0(62), n, 4, and u = +xare
independent of y and that the transverse velocity component v = 4Y is
zero. To be consistent with the neglect of 0(62) terms in the quasi-
static approximation of the near field (section 2.3), we shall terminate
the demonstration of the one-dimensionality of the far field at 0(62
- 43 -
3.2 Equations Governing Time-Periodic Solutions
The one-dimensionality of the far field of the harbor greatly
simplifies the equations governing harmonic interactions. In one
dimension, equations (1.2.16) and (1.2.17) reduce to
-inn + u' + (nsun-s
(3.2.1)
2 1
-inun + k(n + (Uu0
an
and
1" + k 2 
-in I snn n n 2 s n-s sl sn-s)" (3.2.2)
where primes denote differentiation with respect to x. From (3.2.1),
u'= inn*(1 + 0(c)), and u - n (1+ 0(sp)), for n 0. Ton n n nsn
within the accuracy of the Boussinesq approximation
(nsuYn-s
s su-s + n'un-s5
-(i(n 
- s)n sTin-s
s
1
(usun-s
s
(3.2.3)
i - . . )
s-s
(n f s)
=(usus + sn-s
(3.2.4)
~ n-s s n-s
- s(n - s)fsfn-s
s
(n /s)
so that (3.2.2) can be written as
- 44 -
Til 2 =1 n -5 s2)nn I1Iz(n +s),_,_ (3.2.5)
n"+k ln(2 2 sn-s ~2s n-s s n-s
(n f s)
(3.2.5) is a second-order ordinary differential equation governing
the spatial variation of the n-th harmonic.
To solve (3.2.5), we require two boundary conditions on fn. The
no flux condition at x = -Z provides one boundary condition:
n(-z) = 0 (3.2.6)
The other boundary condition is obtained by requiring that, at x = 0,
the wave field in the harbor match smoothly with that in the ocean.
This matching will only be performed in an average way. Equation
(2.3.11), which expresses the average wave amplitude at x = 0 in terms
of the incident wave amplitude and the horizontal fluid velocity at
x = 0, yields the "impedance boundary condition,"
in(0) = An + ZnUn(0) (3.2.7)
Since Zn is O(knSznkn6) and un = in- + 0( 2), we have
kn
(0)=A nk +0(sk nk) (3.2.8)
n
(3.2.5), (3.2.6), and (3.2.8) comprise the boundary value problem for
link().
The zeroth harmonic can be found directly by integrating (3.2.1),
viz.
- 45 -
1
u0(x) = - jns(x) u-s(x) +C0
s (3.2.9)
n0(x) = -1 us 2 + D0 us (x)
Since all velocities, us(x), vanish at x = -z, it follows that C = 0
and u (x) is at most of O(c2). In fact, by the following argument,
u0(x) is seen to be at most of 0(s3): Multiplying (3.2.la) by n- and
(3.2.1b) by u-n and summing n from -co to CO, we have:
I -in (nnf-n + unu-n) + I (nu- =0(s 3,52C 2) (3.2.10)
nn n
where we have used the fact that k = n2 + 0( 2). Replacing the
n
summation index n by -p, we see that
Y -in (fn n, + uu)-n =X pNGpn + u uP) = 0 (3.2.11)
n pp-
so that
S ('nUn)' = 0(63) (3.2.12)
n
Integrating over x and using the boundary condition un(-k) = 0 for all
n, it follows that
I (nnu-n) = OCs3)(3.2.13)
n
and, from (3.2.9a), with CO = 0,
u=(x) - 0(33) (3.2.14)
The constant D in (3.2.9b) is determined by the impedance condi-
tion at the entrance, (3.2.7). In section 2.1, we noted that, by
- 46 -
properly defining the mean depth in the ocean, A0 = 0. From (2.3.5),
Z= 0, so that i1(0) = 0 and
D0 = IY us(0)I2 (3.2.15)
s
From (3.2.9b), nO(x) is then
,n(x) = Jus(0)1 2 - I Ius(x) 2 (3.2.16)
s s
Thus n0 is Q(c2). At the back wall of the harbor, all velocities
us(-) vanish and no takes on its maximum value, viz.
IO(-k)= I us(0)j > 0 (3.2.17)
which shows that there is a mean set-up at the back wall.
Since n0 and u0 are strictly 0(c2) and O(c3) , they do not have
any feedback to other harmonics up to 0( 2). Thus, they can be calculat-
ed after all the other harmonics are obtained. For this reason, we
shall discuss the zeroth harmonic separately.
Equation (3.2.5) represents a system of coupled, nonlinear,
second-order differential equations for the harmonic amplitudes fl 1(x),
n2(x), '3(x),... . When -n (x) is 0(s), the coupling terms are 0(s2)
and are thus relatively small. By earlier theories, i.e., Mei and
Unluata (1972) and Bryant (1972), it is known that weak nonlinearity
does not produce first order effects within a distance of 0(1). There-
fore, for the first few resonant harbor modes, for which z = 0(1), the
effect of weak nonlinearity may be treated as a perturbation, as in
the next section.
- 47 -
3.3 Regular Perturbation Analysis for the Lowest Resonant Modes or for
Very Small Amplitude Waves
Substituting the perturbation series
nn(x) = En (x + s2rj2) (x) + .. (3.3.1)
into (3.2.5), we obtain a sequence of linear problems. At 0(E),
(1) " + -k (1 ) = 0
n n n
(1)
n '(-z) = 0
(1) (0) = An - Zn1l)'(0) (3.3.2)
n
For simplicity, we shall assume that the incident wave is monochromatic,
so that only A1 is nonzero. Then,
n1 (x) = T1 cos k1(x + z)
1 1 1
T = [cos k1zY - (i/k1 )ZI sin k k]~ A1
n () = 0 n=2,3,4... (3.3.3)n
This is just the solution of linearized theory, which exhibits resonances
whenever cos kIz + (IM(Z)/k 1 ) sin k2 = 0.
At ( 2)
"(2) + (2)= 0
T11 + k11
"(2) + k2 n(2) _ 3 (1) (1) 3 (1 ) ,(1)
2222 2 I l ~2 l (3.3.4)
=3 2 2 + 211'+24
- T [cos2 k(x+ k) - k2 sin2 k(x +01
- 48 -
l"(2) +k (2) =0
n n
n = 3,4,5...
with the homogeneous boundary conditions
(2)( 0
n
ri(2) (0) + Zn( 2)-(0) = 0
k n
n
The solution to (3.3.4) and (3.3.5) is
(2)rj =0
n
n = 1,3,4,5...
(2) T2 {f(x) - (
f(0) + 2i
. 2 2f -'O) ) cos k2x +
cos k2Z - 2- Z2 sin k2z
2
(3.3.6)
I+ k222 1-k2
f(x) 2 2 2) cos 2k (x + ) + 2
k2 .4k1 2k2
The solution for n2)(x) is comprised
term, a forced oscillation having the form
oscillation having the form cos k2(x + O)j.
of three parts: a constant
cos 2kI(x + z), and a free
Since
k - 4k2 = 4(l + Aw2) - 4(l+p 2) = 4p2
k2 -2k ~=2(1 + i2) - 2(1 + y2 =2 (3.3.7)
it follows that the constant term is 0(p2) and the forced term is
0(1/12). Thus, the magnitude of e n2 )(x) is generally of order
E a(E/2) , in which case the perturbation series of (3.3.1) is valid only
for e/p 2 <<I1
- 49 -
(3.3.5)
For a harbor which is not too long, e.g., one for which p2<< I,
the perturbation series of (3.3.1) is in fact valid even when e/p2 is
0(1). To show this, we expand
cos 2k1(x + Z) = cos (k2 ~ 2 )(x + )
cos k2(x + Z) + P2(x + k) sin k2(x+ z) (3.3.8)
and, from (3.3.6),
g(0) + 2 Z2 g
n 2) (x) 1 T {g(x -t 2 /cos k2(x + W)+} ()
cos k2k - -2- Zsin k2i
k2
g(x) =.(x + z) sink + ) (3.3.9)
so that, for a harbor having v29 an1, C 2  2)(x) is 0(s 2
In summary, the regular perturbation expansion is valid for
(i) E/p2 << 1 (small amplitudes) and any length harbor, or (ii)
E/P2 = 0(1) and "short" harbors having p2 << 1. To be consistent with
the Boussinesq approximation which neglects terms of o(33), we terminate
the perturbation series at (s 2
From (3.2.16), the zeroth harmonic first appears at O(E2) and is
given by
n(2) ) = 1 (1 (0)|2 - ( 1 2
1 (1),'(O) 2 _I(1) 2)(3.3.10)
2k1
- 50 -
(3.3.3) yields
n 2)(x) = 2IT1 |2 (sin2 k z - sin2 k1(x +
2k12
= 12 fT1j2 (cos 2k1(x + z) - cos 2kgz) (3.3.11)
4k1
Figures 3.1(a), (b), and (c) show the frequency dependence ofjn 1j,
In2I, and no evaluated at the back wall of the harbor (x = -t), for a
harbor having the fixed dimensions a = 50 m., h = 20 m., and L = 1000 m.
The incident wave is taken to be monochromatic, viz. A2 = A3 = A4... = 0.
with A1 = .03. The large peaks occurring at z = 1.41, 4.36, and 7.36
correspond to linear resonances of the first harmonic. In Figure 3.1(b),
the smaller peaks flanking each of these large peaks correspond to
linear resonances of the second harmonic. As frequency increases,
dispersion has the effect of shifting the smaller peaks to the right
relative to the central peaks. Eventually, the small peak to the left
of the central peak overlaps the central peak. When this occurs, the
first and second harmonics are resonated simultaneously, and the
amplitude of the second harmonic can be a large fraction of the first.
As a result of this phenomenon, the peaks in Figure 3.1(b) do not
diminish with z as fast as the first harmonic peaks of Figure 3.1(a).
In Figure 3.1(c), the only peaks in n0 occur at values of k for
which the first harmonic is resonated. In contrast with the second
harmonic, the zeroth harmonic has no additional peaks because it has
no linear resonances of its own. Since n0 is generated by the square of
ag , the peaks in Figure 3.1(c) diminish with k more rapidly than the
first harmonic peaks of Figure 3.1(a).
- 51 -
)'1.1
.3
C,,
.1
Fig. 3.la.
I I-
02
Regular perturbation theory:
2a = 100 m., L = 1000 m., A1
1 I
14 I
I I
H
I I
/0
Frequency response of large-scale harbor (h = 20 m.,
= .03). First harmonic amplitude at the back wall.
NAM
0 4t so
.oao -
.olo-
01
A ' 4a 6 to
Fig. 3.lb. Regular perturbation theory: Frequency response of large-scale harbor (h = 20 m.,
2a = 100 m., L = 1000 m., A1 = .03). Second harmonic amplitude at the back wall.
.10 -
.oS -
No
.o(o G
Cs,
0.02
I
'p
w
Fig. 3.1c. Regular perturbation theory:
2a = 100 m., L = 1000 m., A 1
Frequency response of large-scale harbor (h 20 m.,
= .03). Zeroth harmonic at the back wall.
I I
'A
T I
S
I I
/0
\41rw
.3
-dt 
- 8 .
Fig. 3.2a. Regular perturbation theory: First resonant mode of large-scale harbor (h 20 m,
2a = 100m., L = 1000 m., A1 = .03, k = 1.41). Spatial variation of the harmonics.
SOa mmmama am
a manWM O
-4..
oftan.6am4, m 4
Fig. 3.2b. Regular perturbation theory: First resonant mode of large-scale harbor (h = 20 m.,
2a = 100 m., L = 1000 m., A = .03, z = 1.41). Surface profiles at times t = mr/4,
Mi = 0,l,2...7.
til
L.1
12 39
Fig. 3.3a. Regular perturbation theory: Second resonant mode of large-scale harbor (h = 20 m.,
2a = 100 m., L = 1000 m., A = .03, x = 4.36). Spatial variation of the harmonics.
rz
040
03 
_
Fig. 3.3b. Regular perturbation theory: Second resonant mode of large-scale harbor (h = 20 i.,
2a= 100 m., L = 1000 i., A 1 = .03, = 4.36). Surface profiles at times t =mi/4,
in =0,1,2 ... 7.
.10
VI,
.o2.
.0(9-
% +,
Fig. 3.4a. Regular perturbation theory: Third resonant mode of large-scale harbor (h = 20 m.,
2a = 100 m., L = 1000 m., A = .03, z = 7.36). Spatial variation of the harmonics.
.2
0
r\.4
Fig. 3.4b. Regular perturbation theory: Third resonant mode of large-scale harbor (h =20 in.,
2a = 100 in., L = 1000 n., A = .03, 7.36). Surface profiles at times t =m
mn=o0 ,a,2 ...7m1
Figures 3.2(a), (b), 3.3(a), (b), and 3.4(a), (b) show the
spatial dependence of the harbor response for the first three resonant
modes, occurring at z = 1.41, 4.36, and 7.36, respectively. The (a)
figures show the spatial variation of the harmonic amplitudes n(x) I,
1n2(x)l, and no(x). (Recall that n, and n2 are complex quantities
whereas n 0 is real.) Observe that the amplitude of the second harmonic
increases with increasing distance from the back wall, in accordance
with (3.3.9), and that the zeroth harmonic is largest at the back wall
and zero at the entrance, in accordance with (3.3.11). In figures 3.3(a)
and 3.4(a), the zeroth harmonic is slightly less than zero, indicating
a mean set-down, at places where the first harmonic is smallest. This
slight set-down, however, is considerably smaller than the mean set-up
which exists over most of the harbor.
Figures 3.2(b), 3.3(b) and 3.4(b) show the total surface elevation,
n(x,t), at times t = mr/4, m = 0,1,2...7, for each of the first three
resonant modes. Note that, in each figure, maximum height of the wave
crest is greater than the maximum depth of the wave trough. Further,
in contrast with linearized theory, the free surface has no nodal points
and is at no time perfectly flat. These characteristics are also found
in the study of finite amplitude standing waves.
3.4 An Energy Theorem and Relation to the Korteweg-de Vries Equation
The regular perturbation method of section 3.3 breaks down when
s/12 and ex are 0(1) because resonant nonlinear interactions lead to
significant growth of higher harmonics. The remedy is to include higher
- 61 -
harmonics at 0(E) in the perturbation procedure. Some consequences of
physical interest can be deduced by using the technique of "multiple
scales." From section 3.3, the spatial growth rate of the second har-
monic is Q(y-2) or 0(e71), so we introduce a slow spatial variable
x = Ex and assume
nx) E -n(1) (xi) + 62n2)(xi) + ... (3.4.1)
where x and x are treated as separate independent variables. Then
d a
dx ax
d2  a2  a2  2 a2  (3.4.2)
= + 2 _2dx ax axax ax
For the sake of clarity in identifying secular terms, we shall take
p2 to be 0(E) and rewrite (3.2.5) as
2 1 24 12 2 -2 1 n +s
n + nnn2nn+s2 -s n-s2sn sn-s
(3.4.3)
Substituting (3.4.1) and (3.4.2) into (3.4.3), we have at 0(e)
(4.+ n2 )l(xx) = 0 (3.4.4)
ax
and at 0(e2)
2 + 2(2) (x-) = - I P2n4n(l) - 2 2n (l)
ax Eaxa
aCl) a(1)1 2 2 (1) (1)ln + s1 sa an-s
sj s n-s 2 ) 
-s ax ax
(3.4.5)
- 62 -
The general solution to (3.4.4) is
1) (xk) = c+(x)ef+inx + &(i)einx (3.4.6)n n n
where c (x) and c(x) are so far arbitrary. The boundary conditions onn n
n (x) will be imposed at the end of the calculation.
Substituting (3.4.6) into (3.4.5), we have
( +2(2)+ +nx + f x einx
Dx n ~xx n n
+ terms with e i(n-2s)x (3.4.7)
where
dc 2 4+ 1 2 S2) 1 I(n + ss]cc
f (x) = 12in - - n c + [y (n - ) + n- s(n - s)ccn-s
(3.4.8)
In (3.4.7), the forcing terms f+e nx and f-e-inx are secular, that is,
they lead to solutions for n(2) x) which are unbounded in x. To
eliminate these terms, we set f (j) = 0, viz.n
dc . 2
+ - n3c + (n + s) +sn-s =0 (3.4.9)
dx 6ss
Equation (3.4.9) governs the spatial variation of the coefficients
Cn(~x), which first appeared at 0(E) in the expression for n l)(x,i).n n
The boundary conditions (3.2.6) and (3.2.8) yield, at 0(E);
in(c (-p)e inz - C(-,s)ein) = 0
n n
+ 2 (3.4.10)
cn(0) + c~(0) = An + nZ(c+(0) - C(0))
kn
n
- 63 -
These two conditions, together with the two first-order differential
equations (3.4.9), completely determine c(X) and c~(i).
The system of coupled, nonlinear differential equations (3.4.9) is
no easier to solve than the original system (3.2.5). However, without
actually solving (3.4.9), we can make several important observations.
First, note that the differential equation for c+(x) does not involve
n
c (x) and vice versa. This means that, to 0(s), there is no resonantn
interaction between the right-going wave c ()ein(x-t) and the left-
n
going wave cn)ein(x+t). Further, by multiplying (3.4.9) by cn and
summing from n = -O to n = +o, we obtain
c2 
-i2 n3 c|2 + i (n + s)c csc- = 0 (3.4.11)2 - cn n6ns n sn-s
Replacing the summation indices by p = -n and q = -s, we have
n3 2 _ P3IC+|2=0
~n cnn- Z jcn n p (3.4.12)
I (n + s)c-ncsc=n-s=- (p + q)c qcpc _ = 0
n s p q p-
so that (3.4.11) leads to
SIc (i)I|2 = 2E+
n (3.4.13)
SIc~(x)I2 = 2E_
n
where E+ and E_ are constants, independent of x. From (3.4.10a), it
follows that Ic+(-sz)I = icn(-sk)I for all n, and therefore E+ = E_ = E.
(3.4.13) can be interpreted as a statement of energy conservation.
To Q( 2), the time-averaged energy density at position x is given by
- 64 -
27r
ED(x) =L (q2(x,t) + u2(x,t))dt
D~2r2
= 2 2r1{(' (lx)e-int)2 + ( (1 (x)e-int 2}dt
27 8 nn n n
E 21In ()x)I2 + |u l (x)2 + O(3) (3.4.14)
From (3.4.6)
r7 (x) = c (x)einx + c(~x)e-inx
nn n
(1) ~ 1 )
uO (x) = n = c+(i)einx - c(x)e~inx (3.4.15)nan x n n
so that
D(X) = c2 + Ic()2 2 (3.4.16)ED 4 n nn
where, from (3.4.13), E = E+ = E~ is independent of x. Thus, the total
energy, which is a sum of the energies contained in all the harmonics, is
independent of x. This means that an increase in the energy content of
the higher harmonics must necessarily decrease the energy content of
the first harmonic by the same amount.
Next, we shall establish a connection between (3.4.9) and the
Korteweg-de Vries equation. Consider the total wave amplitude, at
leading order,
nno x intn( X9t) = 1- n -i t
n=-a>
- 65 -
= c7()e in(x-t) + 0~(i)e-in(x+t)
E E F +(,x-t) + e F~(x,x+t) (3.4.17)
where
F (x,) E 0 c (R)e ine (3.4.18)
-00
Using the fact that
(n + s)cscn-s = I (n - s)csc n-s + I 2scscn-s
s S S
= 3 1 scs c(3.4.19)
s
we find from (3.4.9)
+ 2 3F 3+ aF
FF + +3F =0 (3.4.20)6s E: 3  2 as
(3.4.20), which is known as the Korteweg-de Vries (or KdV) equation,
describes the propagation of nonlinear, weakly dispersive, shallow
water waves progressing in one direction.
In the limit of no dispersion, p2 + 0, (3.4.20) becomes
+ 1 F a = 0 (3.4.21)2 as
which is a quasi-linear, first order partial differential equation
possessing the general solution
F = g(e - F X) (3.4.22)
- 66 -
where g is an arbitrary differentiable function. (3.4.22) is an
implicit relation for F. The characteristic curves, along which F(x,e)
is constant, are given by
3
C = g(e - Z- ci) (3.4.23)
where c is a constant. If the function g(e) possesses a maximum at a
finite value of e, then the characteristics of (3.4.23) can be shown to
intersect. At the point of intersection, say (x0,e0), F takes on two or
more different values. In physical terms, a shock is formed at (X0,o0).
Therefore, in the limit of no dispersion, the KdV equation leads to
the formation of shocks.
When 12 is not zero, (3.4.20) is a third-order partial differential
equation, and its solutions are radically different from those of
(3.4.21); in fact, no shocks have been found in all the known solutions,
analytical or numerical. One special class of solutions, known as perm-
anent waves, is of the form
F(x,e) = G(e - yx) (3.4.24)
where y is a constant, and, from (3.4.20), G(o) satisfies the ordinary
differential equation
dG 2 3 3 d
-Y + =4p-4+ G d o = 0 (3.4.25)
(3.4.25) can be solved exactly in terms of the Jacobi elliptic function
"cn"; hence the name "cnoidal waves." These waves typically have sharp
crests, shallow troughs, and phase velocities which depend on amplitude.
- 67 -
From (3.4.20), the equation for FC does not involve F_ and vice
versa. This means that the oppositely directed waves F+ and F do not
interact resonantly to 0(s) in our perturbation analysis. Benney and
Luke (1964) have used this fact to suggest that nonlinear standing waves
can be constructed by superposing two cnoidal waves of equal amplitude
travelling in opposite directions. Such waves differ from linear stand-
ing waves in that they do not have fixed nodal points and the free sur-
face is at no time perfectly flat.
The solution to (3.4.9), which satisfies the complicated impedance
boundary conditions at x = 0, in general will bear no simple resemblance
to cnoidal waves or nonlinear standing waves. In the next section, we
present a numerical procedure for obtaining the solution for the wave
field in the harbor.
3.5 Numerical Solution by Iteration
The infinite set of differential equations (3,2.5) is truncated at
the N-th harmonic, yielding the system
n + k n = R ns'n- +)Y S '%s
s s
n'(-z) = 0
nnnilO) = A - =3 Z -'Z~(Q)
n n k 2 n n
n = 1,2,3...N (3.5.1)
where
- 68 -
k2 En 2(l 1-P2 2kn ( 3~-
S 1 2 2
R =ns= (n -s)
Si (n + s)
ns2 n - s
Zn = k6 [ + 2i en k 6 - Zn ) ( ) (3.5.2)
(3.5.1) is as a nonlinear, two-point, boundary value problem. Unlike
an initial value problem, which can be solved by direct numerical inte-
gration, the solution of (3.5.1) requires an iterative procedure.
At the (p + l)-th iteration, we solve the linear boundary value
problem
(p+l)" + k2 (p+l) = R n(p)I(p+) V Sn'(P)n(P+l)
n n n ns s n-s + s ns s n-s
ryjp+-) (-s)4=O0
n
(p+) (0) = A - (ln)Z p+i)
n
n = 1,2,3...N (3.5.3)
using finite difference methods. The rate of convergence of successive
solutions can be increased by using a relaxation technique, which
consists of replacing n (Pl)(x) by
7-(p-I)(x) = Xn (x) + (1 - X)fn (x) (3.5.4)
With X = 1, convergence of successive iterates was found to be slow and
oscillatory; with X = .5, convergence was much faster and, in most cases,
- 69 -
monotonic. The iterative procedure is terminated when successive solu-
tions for nn(x) differ by less than .001, that is, when
Max ) (x) - 4P(x)I < .001
x
n = 1,2,3...N (3.5.5)
If the maximum value of InN(x)I exceeds .001, truncation at the N-th
harmonic is deemed inadequate and the system of differential equations is
expanded to include the (N + 1) harmonic. The expanded system is then
solved by iteration. This procedure is continued until a "penultimate"
solution is found, for which
Max In,*(x)I < .001 (3.5.6)
x
where N* is the highest harmonic in the solution.
In the finite difference solution of (3.5.3), the domain -z < x < 0
is represented by NPT equally spaced points located at positons
x = -P + (J - l)-H, J = 1,2,3...NPT, where H is the interval size
H = Z/(NPT - 1) (3.5.7)
When the solution requires many harmonics, the choice for H must be
made correspondingly small. For example, if we insist that one wave-
length of the highest harmonic shall be represented by at least ten
points, then the requirement on H is
2w
kN*H<W-o= .628 (3.5.8)
Since the number of harmonics in the solution, N*, is not known a priori,
- 70 -
it is necessary to make some initial choice for H, to solve the problem,
and then to check to see if (3.5.8) is satisfied. If it is not, the
problem must be solved again with a smaller value for H.
The finite difference representation of (3.5.3) is a system of
(N-NPT) linear algebraic equations for the complex amplitudes n+ )(x
n 1,2,3.. .N, J = 1,2,3...NPT. In principle, the solution can be
found by direct Gaussian elimination or matrix inversion, in O(N-NPT) 3
operations. A more efficient method is to first calculate N "basis
solutions" en(x,q) n = 1,2,3...N, q = 1,2,3...N, which are solutions to
the initial value problem for en(x,q) defined by
e" + k2e = Rs s$ e + yI s sn nn ns s n-s Snss n-s
e(-)0
en
1
0
if n = q
if n f q
n = ,2,3...N
q = ,2,3...N (3.5.9)
The solution to (3.5.9) can be rapidly calculated by "stepping" through
the difference mesh from x1 to xNPT. This requires 0(N2-NPT) operations
for each basis solution, using the differencing scheme of (3.5.13)
below. Altogether, 0(N3-NPT) operations are needed to calculate all
N basis solutions.
The complete solution to the boundary value problem (3.5.3) is a
linear superposition
- 71 -
(p+l ) N
n (x) = aqen(x,q) (3.5.10)
q=q
The boundary condition (3.5.3c) leads to the system of N equations
N+.
1 [en(0,q) + ! Ze(0,q)]-a = A
q=l n k nn q An
n
n = i,2,3...N (3.5.11)
which determines the N coefficients aq, q = 1,2,3...N. The solution of
(3.5.11) by Gaussian elimination requires 0(N3) operations, which is in-
3
significant as compared with the Q(N3*NPT) operations needed to calculate
en(xq).
Thus, by this method, the total number of operations required to
solve the boundary value problem for p+ Cx) is 0(N3.NPT). On the
IBM 370 computer at the M.I.T. Information Processing Center, the compu-
tation time per iteration was found to be
Time (in P sec)~ 8(NPT).N3  (3.5.12)
This compares quite favorably with the O(N-NPT) 3 operations that would
be required if direct Gaussian elimination were used.
The finite difference approximation to (3.5.9) and (3.5.11) was
chosen to be
k2
A2en(x3) + [en($x+) + en (xJ-1)] =
Rns 4 )n(x-)e s(x3) + S An (x) Ae (x)ns $ J n-s J s 5 - J
- 72 -
en(x ) = n q
0 ntq
en(x) enxl +H2 [ R nsrs (x.)en(xi) knen(xP] (3.5.13)
and
N .
q= [e(x q) +i Aen Tq) Aq = (3.5.14)
where
Ae(x3) N-[e(xJ) - e(xj) ] = - e(x)] +0(H)
A2e(x3) H [e(x3 ) - 2e(x3) + e(xJ 1)]
2
= d2  e(x)] + 2(H2 (3.5.15)
dx 2 =
This finite difference representation was chosen for stability, accuracy,
and computation efficiency. In (3.5.13a), the advantage of
1 k2 [en(xj) + en(XJ-)] in place of simply knen(XJ) is that it yields
an unconditionally stable scheme for the initial value problem, in which
the right hand side of (3.5.13a) is zero. The expression for en (x2) is
derived from the Taylor series
en(X n(x) + He'(x1 ) + H2e"(x) + H3eg(x) + 0(H4) (3.5.16)
i nwhichn 1eix e ) =
in which e'n(x ) = 0 by the boundary condition (3.5.9b), and the values
of en(x ) and e" (x 1) are found directly from the differential equation,
- 73 -
(3.5.9a). The relatively high accuracy in e nX2) is desirable for
solving initial value problems with en (x) and en(x2 as starting
values. The use of backward rather than central differences for A
increases the overall truncation error from 0(H2) to 0(H), but has the
advantage of increased efficiency in the computational algorithm. If
central differencing is used, a system of NxN linear equations for
en(XJ), n = 1,2...N must be solved at each point x3 of the finite
difference mesh. This would require an additional 0(N3.NPT) operations
for each basis solution and O(N4-NPT) operations for the N basis solu-
tions needed to solve the boundary value problem.
A printout of the computer program, with complete documentation,
is found in Appendix A.
- 74 -
CHAPTER IV
EXPERIMENTAL STUDY OF THE LONG, NARROW BAY
4.1 Equipment and Procedures
4.1.1 Wave Basin
The wave basin, schematically shown in Figure 4.1, measures
43 ft. long by 25 ft. wide by 1.5 ft. deep and is made of concrete.
Care was taken to make the basin floor horizontal. The interior of
the basin is lined with two layers of plastic to prevent leakage and
to provide a smooth inner surface.
The wave generator consists of two aluminum plates, each 12.5 ft.
long by 1.5 ft. wide by .25 in. thick. The plates are connected to
each other by a thin plexiglas strip and are hinged at the bottom to
an eight inch wide plate which is bolted to the floor. Each plate is
connected to a separate flywheel by a bar whose point of attachment
can be adjusted to give different stroke lengths. A rotary drive
shaft connects the flywheels to a 1.5 horsepower motor whose maximum
output torque is 515 in.-lbs. and whose period is continuously variable
from .6 to 3.0 sec. The motor and flywheel assembly are mounted on
the basin wall behind the wave generator. A row of closely spaced
concrete blocks is placed behind and parallel to the aluminum plates
in order to dampen waves generated behind the wave maker.
The open-sea condition is simulated by placing absorbers along
the three walls facing the wave generator. The wave absorbers consist
of a 4 in. layer of rubberized horsehair mounted on a plywood base
- 75 -
1 . *ammmim 
- -.
_- 
- - 13'I
ABSO RBER
.i
"I
ABSORBER
717wlkee1
0 Jotor
Flywkeel
Fig. 4.1. Experimental set-up.
j4
De
wn
00,
--4
cri
I
I
having a 1:5 slope. In preliminary tests, the wave absorbers were
found to have a reflection coefficient of 20% for a wave period of 1.5
sec. and a water depth of 6 in.
Two metal rails, mounted on the side walls of the basin, support
a movable aluminum carriage. Measuring instruments can be placed on
the carriage and swept along the entire length of the basin.
To inhibit the growth of algae, a dilute solution of CuSO4 was
added once every two weeks. The basin was completely emptied every
two months; the liner was washed clean; and deposits of aluminum
hydroxide were scraped from the wave generator.
4.1.2 Wave Gauges
Resistance wave gauges were used in conjunction with a Model 296
Sanborn recorder. A typical wave gauge consists of two parallel
stainless steel wires, measuring 5/64" in diameter and 8" long, separa-
ted by approximately .5". A bridge circuit is used to measure the
resistance between the parallel wires, which is proportional to their
depth of immersion in the water, with an error of roughly 3%. The
stylus of the Sanborn recorder deflects by an amount proportional to
the probe resistance and produces a permanent record on heat-sensitive
paper. Typical calibration curves are shown in Fig. 4.2.
4.1.3 On-Line Digital Computer
In order to analyze nonlinear wave profiles, the voltage output
terminal of the Sanborn recorder is connected by a long coaxial cable
- 77 -
(xa) (x5')
Lt
'3-
0
.2
.oS .1 l+0 .2o
Depth of Immersion (ft.)
Fig. 4.2. Typical calibration curves of a wave gauge, at two attenuator
settings.
- 78 -
to a Hewlett Packard 3450A Multi-Function Meter, which has a maximum
sampling rate of 14 counts per sec. The output values of the digital
voltmeter are read directly into the central processing unit of the
computer, where they are used as input to a BASIC language program.
The computer program is remotely controlled by means of a switch
inserted in the coaxial cable joining the recorder to the computer.
In one position the switch allows the voltage of the recorder, which is
always less than three volts, to pass directly to the computer. In
the other position, the voltage across the terminals of a six volt
battery is fed into the computer. In each program for data reduction
and analysis, a special subroutine directs the execution of the program
in accordance with the position of this switch. (See Appendix B.)
4.1.4 Model Harbor
The model harbor consists of two L-shaped sections placed adjacent
to each other to form a narrow channel, as shown in Fig. 4.1. Two by
four in. beams join the legs of each L for structural rigidity. The
model is made of 3/4 in. marine plywood and is coated with epoxy based
paint, to increase its durability in water. Each leg of the L's is
eight ft. long and consists of a one ft. high wall attached to a six
in. wide foundation. Concrete blocks are placed on the foundations to
ensure that the model does not move under the action of incident waves.
To reduce flow separation, the right angle corner of each L was sanded
to a radius of curvature of roughly .5 in. After the harbor was
placed in position, strips of sponge rubber were used to seal any gaps
- 79 -
between the basin floor and the harbor walls.
The back wall of the harbor consists of a four in. wide by one ft.
high by 3/8 in. thick plexiglas strip which is mounted on a wooden
foundation. The back wall is inserted between the two L's at any
desired position and clamped into place.
In the harbor resonance experiments, the wave amplitude in the
harbor was Fourier analyzed at intervals of two inches along the center-
line of the harbor channel. While the flow field in the vicinity of
the harbor entrance was observed to be two-dimensional, no lateral vari-
ation in the wave amplitude was observed beyond a distance of roughly
four inches into the harbor.
4.2 Multiple Reflections from Wavemaker
In the laboratory, the infinite ocean is modelled by a large but
finite wave basin. We must therefore consider to what extent reflec-
tions from the walls of the basin will perturb the steady-state wave
field in the "ocean." Raichlen and Ippen (1965) have found that the
effect of wall reflections is quite pronounced when all sides of the
basin are perfectly reflective. In our experimental set-up, wave
absorbers suppress reflections from three sides of the wave basin.
However, the aluminum wavemaker on the fourth side is an almost per-
fectly reflective surface.
To estimate the effect of multiple reflections from the wavemaker,
consider the geometry of Fig. 4.3. We shall treat the ocean as linear
and confine our attention to waves of a single frequency, w. The wave
- 80 -
field in the ocean satisfies the Helmholtz equation
(V2 + k2 )n(x,y) = 0 (4.2.1)
where
n(x,y,t) = Re(n(x,y)e-it)
2 1 2 (4.2.2)k = 1 +F(w h/g)
For simplicity, the harbor entrance will be replaced by a point source
located at (0,0), so that
= 2iQ6(y) at x = 0 (4.2.3)
where Q is proportional to the source strength. Denoting the depth-
averaged, horizontal velocity of the wavemaker by Um, we have
"x Urm2Um at x = km (4.2.4)
where kmZ= (w//gT) -Lm. When the harbor is shut (Q = 0), the solution
to (4.2.1) subject to (4.2.3), (4.2.4) is simply the standing wave
-i kU
ns(xy) = sin km ) cos kx (4.2.5)
m
Note that ls is theoretically infinite whenever km = (integer)e-r; in
practice, nonlinearity and viscous dissipation will limit ns to a
finite value.
When Q is nonzero, we must add to rs(x,y) the field due to the
radiating point source at (0,0). If there were no wall at x = zm, the
radiated wave would be
- 81 -
R (xy) = QH( kr) (4.2.6)
just as it is in the infinite ocean. In order to satisfy the boundary
condition
= 0 at x = (4.2.7)
ax k
we place an "image" source of strength Q at position (2m,0). But now,
to satisfy the zero flux condition at x = 0, we must place another
image source of strength Q at position (-2zm,0). Proceeding in this
way, we obtain an infinite string of image sources, all of strength Q,
located on the x-axis at x = 2zm' 4 m, 6,..., as shown in Fig. 4.3.
The total wave field due to all these image sources is
fi(x,y) = Q IH (kr - 2nzmxI) + H(l)(klit + 2nzmx|) (4.2.8)
n=l 0 m
Along the x-axis
i(x,o) = E(-x) + E(+x) (4.2.9)
where
E(x) = Q H0l(2nkzm + kx) (4.2.10)
n=l
To approximate E(x), we follow Morse and Feshbach (1953) and assume
kim >> 1, so that, in the domain 0 < x < zm all Hankel functions may
be represented by their asymptotic forms. Then
CO 2inkp,
E(x) -/fQeikx ei m (4.2.11)
n=1 v'2nkkm + kx
- 82 -
yj
qJ
0QS
K
C
4 JhA
d
0
1.
0I
d
E
0
6
A K
(*4M )
94
Fig. 4.3. Theoretical study of multiple reflections from wavemaker: Placement of images.
(0
Except in the case k'm = (integer)-w, for which the sum in (4.2.11)
diverges, we may approximate E(x) by
Co i(2kz )v
E(x) ~ Qedikx d e r (4.2.12)Q7jdv 2kkmv + kx
m
The integral in (4.2.12) can be evaluated explicitly in terms of Fresnel
integrals; but, for the purpose of estimation, we simply integrate
once by parts to obtain
2ikt
E(x) 7 Qeikx f(2ikz) e } + 0(k.m)-5/2 (4.2.13)
m 2kkM + kx
From (4.2.9), the net result of all multiple reflections is
i(2kz + kx) i(2kz - kx)
ii(x,o) ~ {L " m e m )+0k -/
/ZXo) k m /2kim + kx v2kz -kx m
(4.2.14)
The interesting feature of (4.2.14) is that, whereas the effect of a
single reflected wave is O(kmF)- 1/2 , the combined effect of all reflec-
tions is only of O(kkm)-
Our conclusion therefore is that, so long as
kk >> *I
(4.2.15)
kkm t (integer)-w
successive reflections from the wavemaker add destructively to produce
a net effect whose magnitude is only 0(kM)3/2 When (4.2.15) is
satisfied, the total wave field in the laboratory "ocean" is
- 84 -
n(x,y) = (si k )cos kx + QH(l)(kr) + 0(kzm)-3/2 (4.2.16)n~~) sin ky 0m (..6
and the effect of multiple reflections is small.
When either of the conditions in (4.2.15) is violated, the effect
of multiple reflections can be quite significant. When this occurs,
the harbor basin and ocean basin are said to be strongly coupled.
Condition (4.2.15) states that strong coupling will occur unless (a)
the ocean basin is very much larger than a wavelength and (b) the
frequency of oscillation is not close to a resonant frequency of the
ocean basin.
Figs. 4.4 (a),(b),(c)show experimental measurements of the wave
field in the model ocean, along the centerline x = 0. The experimental
values
Lm = 31'2"
o = 2f/(1.545 sec) = 4.067 sec 1
h = .5'
yield
km = 32.96 = 10.49n
so that (4.2.15) is satisfied. The graph at the top shows the standing
wave amplitude ns(x,o) in the ocean, when the harbor entrance is
closed; the solid curve is (4.2.5). The middle graph shows the total
wave field when the harbor entrance is open, and the first resonant
- 85 -
.0
Co0
2f 'itrrf Sir /Orr
Fig. 4.4a. Experiment (.) vs. theory (-) for the wave amplitude in the ocean (h = .5 ft., 2a =
.33 ft., L = 1.211 ft., w = 4.067 sec-1 , L = 31.17 ft.). Top: harbor closed.
Middle: harbor open. Bottom: the radiatei wave. A1 = .015.
I 5
CO .0
.5
2r '-r 4w trr
Fig. 4.4b. Experiment () vs. theory (-) for the wave amplitude in the ocean (h = .5 ft., 2a =
.33 ft., L = 1.211 ft., w = 4.067 sec~ , L = 31.17 ft.). Top: harbor closed.
Middle: harbor open. Bottom: the radiatd wave. A1 = .027.
ofo
O..
CO
Fig. 4.4c. Experiment () vs. theory (-) for the wave amplitude in the ocean (h = .5 ft, 2a =
.33 ft., L = 1.211 ft., w = 4.067 sec-, L = 31.17 ft.). Top: harbor closed.
Middle: harbor open. Bottom: the radiat~d wave. A1 = .040.
harbor mode (having L = 1.211 ft.) is excited; the solid curve is
(4.2.16). The bottom graph shows the component of the total wave field
which is the radiated wave, nR(x,o) = QH l)(kx). Since the agreement
is fairly good, it confirms experimentally that (i) the ocean can be
treated linearly and (ii) reflection from the wavemaker is not signi-
ficant.
4.3 Real Fluid Effects
In the following two sections, all variables and coordinates are
dimensional.
4.3.1 Viscous Dissipation in Wall Boundary Layers
In a laboratory setting, the damping of surface waves in a slightly
viscous fluid is a well-known problem. For infinitesimal waves of
frequency w, propagating in a fluid of viscousity v, the fluid motion
is essentially irrotational except near boundaries, where viscous
boundary layers of thickness
6v = /2v/w(4.3.1.1)
are formed. Viscous energy dissipation occurs in (a) the boundary
layers near solid walls, (b) the boundary layer near the free surface,
and (c) the main body of the fluid. Ursell (1952) has shown that the
solid wall boundary layers account for most of the dissipated energy.
For oscillatory boundary layers, the transition from laminar to
turbulent flow on a smooth flat surface occurs at the critical
- 89 -
Reynolds number (see Phillips (1966), p. 43)
R - 6 -~160
c v
(4.3.1.2)
where U is the particle velocity. In shallow water,
U9 = A-/iY (4.3.1.3)
where A is the wave amplitude, and
c vw (4.3.1.4)
In the experiments, h = .5 ft., gh = 16 ft2/sec2, w = 4.067 sec" 1,
v = l0-5 ft2/sec, and
R = (g)-887 (4.3.1 .5)
so that the transition to turbulence occurs when
(c) c = .18 (4.3.1.6)
In all but one of the harbor resonance experiments, the value of A/h was
below this critical value; therefore, laminar flows prevailed.
Batchelor (1976, pp. 353-358) derives the following expression for
the time-averaged power loss per unit area in an oscillatory laminar
boundary layer:
dP vP 1 2
S - p IUTI
v
(4.3.1.7)
where Re(UTe-ut ) is the local tangential velocity at the boundary
according to inviscid theory and p is the fluid density.
- 90 -
To estimate the viscous losses in the far field of the long, narrow
bay, we shall use the linear expressionsfor the wave field,
n(x,y,t) = Re(T cos k(x + L)e-iwt)
(4.3.1.8)
tI(x,y,t) = Re(i gkT sin k(x + L)e-lwt)x + Q(y2)
where 0(p2) terms include the vertical component of velocity and the
depth-dependent terms of the horizontal velocity. The expression for T
in terms of the incident wave amplitude is found in Equation (1.3.6).
From (4.3.1.7), the power loss per unit area on all surfaces parallel
to the x-axis is
dP
vv 1 PVgkT12 sin2 k(x + L) (4.3.1.9)dA 2 v6 xW Lv
Integrating over the side walls and bottom of the channel, we obtain
for the total power loss
P = (2a + 2h) f()dx
-L
=I pv |w2(2a + 2h) -L(kL - } sin 2kL)
~-jgkT|2(a + h)L (4.3.1.10)
v
Let us compare P v with the energy loss due to radiation damping at
the harbor entrance. The time-averaged radiated power is given by
27T/c
P= (2a-h).- dt pgn(ot)u(o,t)
0
= (ah)pg(g9)Tj2 sin 2kL
- 91 -
!~ah)pg(9K) Tj 2  (4.3.1.11)
and, from (4.3.1.10),
P /P L(k) a +h L
yRR a a
= (a2h)(kL)(v) (4.3.1.11)
In the harbor resonance experiments, a = 2 in., h = 6 in., W = 4.067
sec , v =10- ft2/sec and
- =/- 2.22 x 10-3ft.
kL (n + )Tr n = 0,1,2 (4.3.1.12)
so that
P RR=~(n + )-(2.79 x 102 (4.3.1.13)
Thus, the effect of laminar viscous damping is relatively unimportant
for the first few modes. It should be noted, however, that a small
amount of surface roughness can lead to the formation of turbulent
boundary layers, for which v is effectively much larger than its
laminar value. For this reason, care was taken to make the walls of
the harbor as smooth as possible.
In a large scale harbor, surface roughness and the Reynolds
number in (4.3.1.4) are much greater than in the laboratory, and so,
the boundary layers can be turbulent. From the experiments of Jonsson
and Carlsen (1976), the effective viscosity ve is typically of 0(100v).
- 92 -
Taking the following typical values for a large scale harbor
h = 20 m
2a = 100 m
L = 1 km
Ve = 100V = 10r4 m2/sec
kL(1/4 wave mode) = 1.41
w(1/4 wave mode) = .02 sec 1
we obtain
= .1 m
P RP~R(.7)(1.41)(.005) ~ .005
so that the power loss in the turbulent boundary layers is again
negligible.
4.3.2 Separation Loss at the Entrance
Except at very low velocities, flow separation occurs at the sharp
corners of the harbor entrance, and energy is expended in the production
of eddies. The mechanism of flow separation is related to the adverse
pressure gradient caused by the sharp curvature of the boundary. If
the flow in the entrance is assumed to be quasi-static, then the
pressure drop due to flow separation can be described by the hydraulic
formula for steady flow
n(0~,t) - (0+,t) = |%Iu(0,t)Iu(0,t) (4.3.2.1)
- 93 -
where 0 is on the harbor side of the entrance and 0+ is on the ocean
side. The constant f is an empirical friction factor whose value
depends on the Reynolds number and the geometry of the entrance. For
an asymmetrical entrance, f may take on different values depending on
whether u(0,t) is positive or negative.
An expression similar to (4.3.2.1) was used by Ingard and Ising
(1967) and Ingard (1970) in studying nonlinear sound transmission
through an orifice, with application to the absorption of sound at
high pressure levels. Experimental measurements of both the energy
loss and the distortion of the frequency spectrum of the transmitted
sound were found to agree with the hydraulic, steady-flow formula. In
particular, for a monochromatic wave incident on a symmetrical orifice,
only odd harmonics are generated by the nonlinear term fuful; whereas,
for an asymmetrical orifice, both even and odd harmonics are generated.
A reduction in the mass reactance, or inertia, of the orifice was
also observed by Ingard and Ising (1967), but this is not accounted for
in (4.3.2.1). The magnitude of the mass reactance depends on the
details of the velocity field, which is difficult to calculate when
the flow is separated. Since the inertia of the orifice is associated
primarily with the irrotational part of the flow, it is expected to be
smaller when flow separation occurs. At high sound pressure levels,
the mass reactance is observed to be one half its unseparated value,
which suggests that, at any given time, the flow field on one side of
the orifice is fully separated.
In water wave problems, separation loss has been studied in
connection with the construction of perforated breakwaters; see
- 94 -
Jarlan (1965); Terrett, Osorio, and Lean (1968); and Mei, Liu, and
Ippen (1974). The separation loss at a narrow harbor entrance has
been studied numerically by Ito (1970) and analytically by Unluata and
Mei (1975). The latter show that the higher harmonics generated by the
nonlinear term -!-Ju(O,t)|u(0,t) in (4.3.2.1) are not large under usual2g
circumstances; rather, the primary effect of separation is to reduce
the amplitude of the fundamental harmonic. We shall therefore present
a simplified treatment of flow separation based on a monochromatic
entrance velocity, u(0,t).
To allow for an asymmetrical harbor entrance, we shall approximate
f by
f 1 when u(0,t) > 0
f = { (4.3.2.2.)
f 2 when u(0,t) < 0
This is a mathematical simplification since f is more likely to be a
continuous function of time. We then have
f-uluj = uju + +fu2 (4.3.2.3)
where
I f +f2 (1  f2)
(4.3.2.4)
2 1
Next, we Fourier analyze each term of (4.3.2.3):
-u | = UI LY eint2g lul -2g n
- 95 -
fU2= nXf - int (4.3.2.5)2g 2g _0n
If we assume that u(0,t) has the simple form
u 'Ue-it + U*eit) = jul cos(t-t) (4.3.2.6)2
where is the phase of U, viz. U =uele then the Fourier coefficients
in (4.3.2.5) are found to be
fi f ,T u Iueint dt
0
4 sin(nir/2) fUJ2 ein
'r n(4-nt )
2 2 . (4.3.2.7)
f = -L f u2 e'nt dt
(o1/2 for n = 0
= fjUl 2 ein . 1/4 for n = +2
0 for all other n
Observe that fn is zero for even values of n; thus, we retrieve the known
result that, for a symmetric entrance (f=O), only odd harmonics are
generated by separation. For an asymmetrical entrance, the second and
zeroth harmonics are also generated; higher even harmonics would also be
generated if u(O,t) were more complicated than in (4.3.2.6).
For steady flow through the entrance of Fig. 1.1, the values of
f and f2 are found to be 1.0 and .4, respectively, according to Daily
and Harleman (1966), pp. 314-319. For oscillatory flows, the separation
zone on either side of x = 0 has only a half wave period during which
to develop. It may therefore be expected that the effective friction
factors are smaller for oscillatory flow. Moreover, in the experimental
- 96 -
harbor model, the sharp corners at x = 0 are slightly rounded, which
further reduces flow separation. Thus, the average friction factor
1(f1 + f2) in the experiments is expected to be smaller than .7.
We shall see in section 5.1 that the value of f in the experiments is
roughly half this value. Taking f1 = .5 and f2 = .2 in the experiments,
we have f = .35, f = .15 and
f = TIU|2 1 - .15|UU1 3'T .j
f3 = 4.f|Ui2e3i4 = .03 U3/JUJ3 157T
(4.3.2.8)
f 0 = 2f UI = .075jUI
f2 = If IU 2e2 ic1 = .0375 U2
Since U is 0(s), the above coefficients are all 0(62). Note that the
main effect of separation is to alter the first harmonic, through the
coefficient f1. The term fi will lead to a slight increase in the
mean set-up, n0(x).
Taking the first harmonic component of (5.3.2.1) and using the
impedance condition, n1O(0+) = A1 + Z1U, on the ocean side, we obtain
YO(0) = A1 + Z1U + T/g
(4.3.2.9)
A I+ (Z1 + c )U
where
c, = $ TU|/g (4.3.2.10)
- 97 -
The constant ce is an effective damping constant due to separation.
Eq. (4.3.2.9) is an effective harbor boundary condition which includes
(i) the forcing of the incident wave, A.; (ii) the inertia of the har-
bor entrance, Im(Z); (iii) the damping due to the radiated wave,
Re(Z1); and (iv) the damping due to separation loss at the entrance, ce.
Using (4.3.1.8) for the wave field in the harbor, we have
ri (0~) = T cos kL (4.3.2.11)
U = i T sin kL (4.3.2.12)
and, from (4.3.2.9),
T = A[cos kL - i9kS(Z + c ) sin kL]
co 1 e
= A[(cos kL + 9k Im(Z1) sin kL) - i9 (Re(Z ) + 4 k T|T|) sin kL 1-~W 1 W 1 37r w
(4.3.2.13)
Resonance occurs, as before, when cos kL + Im( Z1) sin kL = 0. At
resonance,
ITI = IAl[ (Re(Z 1) + - k K |TI) Isin kL|]~I
= ITJ|[l + aT/|Tfol] (at resonance) (4.3.2.14)
where
T| =AJ(9 t Re(Z 1)|sin kL)~
= ITI when f = 0 (no separation)
S-ky |T/Re(Z) (4.3.2.15)37JTJwRe
- 98 -
Physically, a measures the amount of energy loss due to separation as
compared with that due to radiation damping. When a >> 1, separation
losses dominate radiation losses.
Solving (4.3.2.14) for ITI/ITJ|, we find
|Tj/jTo| = 1 1 + 4cc
~ a-1/2 when a >> 1/4 (4.3.2.16)
A plot of fT/IT vs. a is shown in Fig. 4.5. Suppose we allow Re(Z)
to approach zero while maintaining f and A constant (and nonzero). This
corresponds to the well-known harbor paradox, in which Re(ZI) is made
smaller by decreasing the width of the harbor entrance and JT0j is
observed to increase without limit. From (4.3.2.16), we see that
|TI/IT0! + 0 and that, specifically,
IT! jT a-1/2
= (4 f)-/2 (Re(Z)fTf)1/2
= (I 4k -/2 (JJAl 1/2 (4.3.2.17)
37 r w kisin kLj
Instead of increasing without limit, ITI approaches a finite value
proportional to (IAI/f)l/2. Thus, the presence of separation losses
effectively removes the harbor paradox.
- 99 -
00
I /.I/-
V
So /oo
Fig. 4.5. Entrance loss: Normalized first harmonic amplitude, JT/T0j, vs. normalized friction
factor, a.
CHAPTER V
EXPERIMENTAL AND THEORETICAL RESULTS
5.1 Comparison of Theory with Experiment
The experimental set-up is shown in Fig. 4.1, where the coordinates
are defined.
Five length scales are involved in the study of the long, narrow
bay: harbor width, 2a; mean water depth, h; harbor length, L; wave-
length of the first harmonic of the incident wave, approximately
vg"hi-2ir/w; and incident wave amplitude, |A1 Jh. In the experiments,
these had the values
2a = 1/3 ft.
h = 1/2 ft.
L = 1.211, 4.173, and 7.136 ft.
/gh.(2w/w) = (4.0125 ft/sec) (1.545 sec) = 6,200 ft.
A1 -h = .0075, .0135, and .020 ft.
Three harbor lengths and three incident wave amplitudes were tested.
The corresponding dimensionless parameters are
6 E a-(w//gh) = .169
9 E L.(w/vg-h) = 1.227, 4.230, and 7.233
2 2P = W h/g = .257
A1J = .015, .027, and .040
The values of z were chosen so as to produce the first three linearly
- 101 -
resonant harbor modes. Each resonant mode was studied for three
different values of 1A1 I; so, altogether, nine resonant harbor
oscillations were examined.
The incident waves produced by the wave generator are not per-
fectly monochromatic. With the harbor entrance closed, the standing
wave amplitude in the ocean was measured at x = y = 0 and was found to
have the harmonic composition given in Table 5.1. The phases, Pn, are
defined by An E aAnt n and are given in radians.
Table 5.1
|A1 1 .015 .027 .040
A21 .001 .003 .012
A31 .000 (2) .001 .003
2 - 21 6,210 5.646 5.506
3 - 3 .671 1.795 5.314
These are used as inputs in the numerical calculations.
Table 5.2 compares the wavenumbers kn calculated by the approximate
Boussinesq dispersion relation
(knh)2 =n2 2(1 +TI n 2) + 0(np)6n3
2 = ,2
where p = W h/g = .257, with those calculated by the exact dispersion
relation
n2 2 = (knh) tanh (knh)
- 102 -
which is valid for arbitrary depth, e.g. arbitrary (np).
Table 5.2
n Exact kn (ft:.) Boussinesq kn (ft:1 )
1 1.060 1.056
2 2.448 2.350
3 4.710 4.048
4 8.228 6.244
5 12.850 8.986
Note that the error in the Boussinesq values for kn is quite large for
the higher harmonics.
The theoretical assumption that the far field of the harbor is one-
dimensional is only valid when transverse modes of oscillation, or
"cross modes," are not present. The minimum wavenumber, kc, for which
cross modes are possible, is given by
(kC).(2a) = T
so that, with 2a = 1/3 ft.,
kc = 9.42 ft:1
From Table 5.2, we see that cross modes can exist for the fifth harmonic
or higher. In the experiments, the amplitudes of these harmonics were
negligible, and the far field of the harbor was observed to be approx-
imately one-dimensional.
- 103 -
The wave amplitude in the harbor was Fourier analyzed at intervals
of two inches along the centerline of the harbor channel. The zeroth
harmonic, or mean set-up, was removed by the Fourier analysis. The vari-
ation of the first three harmonic amplitudes with distance from the back
wall of the harbor is shown, for each of the nine experimenal cases,
in Figs. 5.1(a) through 5.1(i). The horizontal dotted line indicates
the maximum amplitude of the first harmonic predicted by linearized
theory. The solid lines are theoretical curves based on the inviscid
nonlinear theory of Chapter 3. These solutions were computed by the
numerical procedure of Section 3.5, with the incident wave composition
given in Table 5.1 and the Boussinesq values of kn given in Table 5.2.
Convergence of the iteration scheme was achieved in all cases with five
harmonics or less.
For the shortest harbor (L = 1.211 ft.) the higher harmonics are
small, owing to the short distance over which nonlinear interactions
take place, and the inviscid nonlinear theory agrees closely with
linearized theory. The reduction of the first harmonic amplitude is
due primarily to separation loss at the entrance, rather than to
nonlinearity, and will be discussed later.
In the two longer harbors (L = 4.173 ft. and 7.136 ft.), there is
considerable harmonic generation; in some instances, the second harmonic
is as large as fifty percent of the first. In such cases, the second
harmonic siphons an appreciable amount of energy from the first harmonic,
and therefore, the maximum amplitude of the first harmonic is observed
to be significantly less than that predicted by linearized theory. For
the longest harbor (L = 7.136 ft.), the agreement between experiment and
- 104 -
and inviscid theory is fairly good. In particular, the significant
reduction of the first harmonic and the presence of higher harmonics
are well accounted for by the nonlinear theory. However, for the second
harbor (L = 4.173 ft.), we note the qualitative feature that the ob-
served second harmonic is higher than calculated, and that the observed
reduction of the first harmonic is larger than that calculated. These
discrepancies cannot be attributed to harmonic generation by flow separa-
tion at the entrance, as this was shown to be quite small, in Section
4.3.2.
We believe that part of the discrepancy between theoretical and
experimental values of the higher harmonics is most likely due to errors
in the Boussinesq values for kn, cf. Table 5.2. Errors in the values
of k n lead to corresponding errors in the positions of nodes and anti-
nodes of the harmonic amplitudes. These errors increase with increasing
distance from the back wall of the harbor. Furthermore, because the
harbor response is large at resonant values of knL, a small error in
k n can significantly alter the absolute magnitude of the n-th harmonic.
We suspect that this happens for the second harbor (L = 4,173 ft.). In
this case, the Boussinesq value of k2L is 9.807 which does not lie close
to a resonant mode; whereas the exact value of k2L is 10.216, which
7
coincides with the fourth resonant mode, having kL = f + Im(Z).
Because of this, the second harmonic predicted by Boussinesq theiry is
expected to be smaller than that observed in the laboratory. The
energy conservation theorem of Section 3.4 implies that the first har--
monic must also be affected.
Figures 5.2(a), (b), and (c) show the results of a numerical
- 105 -
calculation in which the exact values of kn are used instead of the
Boussinesq values, for the second harbor. Note that the second har-
monic is now larger and the first harmonic smaller than that obtained
previously (cf. Figs. 5.l(d),(e),(f)). This improves the agreement
between theory and experiment, though the quantitative agreement is still
not as good as for the long harbor (L = 7.136 ft.) case. Of course,
this ad hoc correction of the values of kn must await rigorous justifica-
tion by a nonlinear theory which is valid for arbitrary depth.
We have used the shortest harbor to determine the average friction
factor, ? = (f + f2), which appears in the entrance loss theory of2 1
Section 4.3.2. The value f = .35, which is half the upper limit expect-
ed on the basis of steady flow values for f1 and f2, was found to give
good agreement with the experimental data on the shortest harbor. In
Figs. 5.3(a) through 5.3(i), the experimental data for all nine harbor
resonances is compared with theoretical curves computed using the
entrance loss theory of Section 4.3.2, with ? = .35. Since harmonic
generation by separation has been shown to be small, only the main
effect of separation, e.g. reduction of the first harmonic, has been
included in the numerical computation.
The inclusion of entrance loss does not significantly improve the
agreement between experiment and theory for the second and third harbor
lengths. In particular, the quantitative agreement is still better
for the longest harbor than for the medium -ength harbor. By using
T = .35 and an ad hoc application of the exact dispersion relation, as
explained previously, we obtain the theoretical curves shown in Figs.
5.4(a),(b),(c) for the medium length harbor.
- 106 -
*10
.4L,
Fig. 5.la. Experiment vs. inviscid nonlinear t'heory (f = 0). (A,+,x): measured first, second,
and third harmonic amplitudes. (-): nonlinear theory. (----): linear theory
for first harmonic amplitude. L = 1.211 ft., A1 = .015.
C
00
2&
L 1/.
Fig. 5.1b. L = 1.211 ft., A = .027 (See Fig. 5.la).
.301
.151
10
Fig. 5.lc. L = 1.211 ft., A1 = .040 (See Fig. 5.la).
to. ;t0 9
'4.
'4
'4
'4
'4
Sb
&
.05 & a
&
&
U, * , hi: 'a:'. *:
L = 4.173 ft., A1 = .015 (See Fig. 5.la).Fig. 5.1d.
&4
5 & & &  &
4 
42K32:
Fig. 5.le. L = 4.173, A1 = .027 (See Fig. 5.la).
.1-1
r%33t
Fig. 5.1f. L = 4.173 ft., A1 = .040 (See Fig. 5.la).
IPM am Oft
4e 
OR
.. o -
e/m1%9I}\\1/
*os
X 10
L = 7.136 ft., A1 = .015 (See Fig. 5.la).Fi g. 5.l1g.
*1\
X+
Fig. 5.1h. L = 7.136 ft., A3 .027 (See Fig. 5.1a).
.30
/ \ C- 39
I\ /
/ \
01
- +4*.+ I\t
L = 7.136 ft., A = .040 (See Fig. 5.la).Fig. 5.1i.
.lc -
WAD
+2 3
Fig. 5.2a. Experiment vs. inviscid nonlinear theory (f = 0) with exact dispersion relation.
(A,+,x): measured first, second, and third harmonic amplitudes. (-): nonlinear
theory. (----): linear theory for first harmonic amplitude. L = 4.173 ft.,
A1 = .015.
.et
foo am\u
a x 3c
33
Fig. 5.?b. 1 = 4.173 ft., A1 = .027 (See Fig. 5.2a).
930
wo 3cIr-r+-93
itt rf
Fig. 5.2c. L = 4.173, A = .040 (See Fig. 5.2a).
)42
Fig. 5.3a. Experiment vs. nonlinear theory with separation loss ( = .35). (A,+,x): measured
first, second, and third harmonic amplitudes. (-): nonlinear theory. (---):
linear theory for first harmonic amplitude. L = 1.211 ft., A1 = .015.
sawoh
46
Fig. 5.3b. L = 1.211 ft., A1 = .027 (See Fig. 5.3a).
-4
13'4
I
Fig. 5.3c. L = 1.211 ft., A1 = .040 (See Fig. 5.3a).
,30
.a5S-
-J
A6
/.2
w6a It
.10~
a
&3
Fig. 5.3d. L = 4.173 ft. , A1 = .015 (See Fig. 5.3a).
4m I %%
&a & /
393
Fig. 5.3e. L = 4.173, A = .027 (See Fig. 5.3a).
.30 -
/9x39 31/
i i i i : : i/
.153
Fig. 5.3f. L = 4.173 ft., A1 = .040 (See Fig. 5.3a).
.0 am
J112 01P/ l o%
/1~4 & +
L = 7. 136 ft., A1 = .015 (See Fig. 5.3a).Fig. 5.3g.
e4t* ql ~ %
/-9 /
3tX 3/3
/C
L = 7.136 ft., A1 = .027 (See Fig. 5.3a).
I 4e we %* %.
Fig. 5.3h.
/ \
/ \
'*+Is+. 
. ./+/ & x
L = 7.136 ft., AI = .040 (See Fig. 5.3a).
.30 w7
Fig. 5. 3i .
\ /
\ /4
ama
Fig. 5.4a. Experiment vs. nonlinear theory with separation loss (F = .35) and exact dispersion
relation. (A,+,x): measured first, second, and third harmonic amplitudes. (-):
nonlinear theory. (----): linear theory for first harmonic amplitude. L = 4.173 ft.,
A1 = .015.
5lo W, &m/Gr
4m
&~
4D.
at s 3a ax
/ a3 9
Fig. 5.4b. L = 4.173 ft., A1 = .027 (See Fig. 5.4a).
.30
r -4.. r r
C
r
'4 / 
'4,
/ '4
/ 
"S
.13-
4
&
&
03C a /
a
'6
5.4c. L = 4.173 ft. ,A1 = .040 (See Fig. 5.4a) .
5.2 Some Theoretical Predictions for a Large Scale Harbor
Having obtained some experimental confirmation of the nonlinear
theory, we shall investigate theoretically the response of a harbor
having the practical dimensions 2a = 100 m, h = 20 m, and L = 1000 m.
The first three lineary resonant modes occur for z = wL//g-i = 1.41, 4,36,
and 7.36 corresponding to wave periods of 5.305, 1.716, and 1.016
minutes, respectively. Table 5.3 summarizes the linear results per-
taining to these harbor modes.
Table 5.3
_ p_ Amplification (.n1C-z)/A1 )
1.41 7.95 x 10-4 .0705 14.35
4.36 7.60 x 10-3 .218 4.83
7.36 2.17 x 10-2 .368 2.99
Note that each value of P 2 is much smaller than the corresponding value
of 6; this is typical for large scale harbors. In practice then, the
assumption 6 << 1 is likely to be more restrictive than the assumption
P2 << 1.
p2
Consider the response of the harbor to monochromatic incident
waves, with normalized amplitude A1 = .03, A2 = A3 = A4... = 0.
Linearized theory predicts that, for the first resonant mode (z = 1.41),
the normalized wave amplitude at the back wall of the harbor will be
,nl(-k) = .430. This corresponds to a peak-to-trough wave height of .86
times the depth, which is greater than the value at which waves break in
- 131 -
shallow water.
The predictions of inviscid nonlinear theory are shown in Figs.
5.5(a) through 5.5(d), for the case A = .03. Fig. 5.5(a) shows the
number of harmonics required in the numerical solution, N*, versus the
normalized incident wave frequency, z = wL/vgW. As expected, the
greatest number of harmonics are required at values of z corresponding
to resonant harbor modes, for which the wave amplitude in the harbor
is large.
Fig. 5.5(b) shows the variation of the first harmonic amplitude
at the back wall, (rI(-9J), with incident wave frequency, w. The solid
line is the prediction of linearized theory (and regular perturbation
theory, also); the points are nonlinear results. According to nonlinear
theory, the first resonant mode occurs at z = 1.37 and has a peak height
of .382; as compared with z = 1.41 and a peak height of .430 in linear-
ized theory. For the second and third resonant modes, the resonant
&mplification is smaller, and the discrepancy between nonlinear and
linear theory is therefore smaller.
Fig. 5.5(c) shows second harmonic generation due to nonlinearity.
The solid line is calculated from the regular perturbation theory of
section 3.3; the points are numerically-generated results. From
Table 5.3, we see that, in the range of z covering the first three
2
resonant modes, z << 1/B2. Therefore, we expect that regular
perturbation will adequately describe second harmonic generation so long
as the wave amplitudes are not too large. Comparing 5.5(c) and 5.5(a),
we see that the major discrepancies between regular perturbation and
numerical theory occur where the numerical solution required many more
- 132 -
than two harmonics. The various peaks in second harmonic generation
have been discussed in Section 3.3; here we simply recall that the
three central peaks correspond to linear resonance of the first har-
monic, whereas the two small peaks flanking each central peak correspond
to linear resonance of the second harmonic. The effect of dispersion
increases with increasing k and tends to shift the small peaks to the
right relative to the central peak. At k = 7.36, one of the small
peaks has almostoverlapped its adjacent central peak; this corresponds
to a situation in which both the first and second harmonics are reson-
ated simultaneously. As a result of this phenomenon, the peaks in
Fig. 5.5(c) do not diminish with z as fast as the first harmonic peaks
of Fig. 5.5(b).
The variation of J|l 3 (-z)J with z is shown in Fig. 5.5(d). Note
that the largest peak in Ir3(-R)| occurs at z = 1.41, corresponding to
the first resonant harbor mode, and has a magnitude of .023, which is
comparable to the magnitude of n2(-0). This is contrary to the usual
situation in which the higher harmonics are progressively smaller. The
reason for this can be seen from Table 5.3. Since the second resonant
harbor mode (z = 4.36) occurs at a frequency very nearly three times that
of the first resonant harbor mode (z = 1.41), we see that, when z = 1.41,
both the first and third harmonics are resonated simultaneously. This
is quite generally true for the long, narrow bay, because its resonant
modes appear at roughly odd multiples of '/2. This resonance mechanism
for the third harmonic is similar to the one occurring at z = 7.36 for
the second harmonic, though the underlying reasons are different for the
two cases.
- 133 -
Fig. 5.5(e) shows the mean set-up at the back wall, n(-z), vs. Z.
The solid line is the prediction of regular perturbation theory, cf.
Fig. 3.1(c); the points are calculated from the numerical solutions
using (3.2.17). Note that the peaks occur at values of z for which the
first harmonic is linearly resonated. Since no depends primarily on
the square of nl, the peaks of no diminish in height more rapidly than
thoseof T,.
The numerically-predicted spatial dependence of the harbor response
is shown in Figs. 5.6(a)-(c), 5.7(a)-(c), and 5.8(a)-(c) for each of
the first three resonant modes. These figures may be compared with the
predictions of regular perturbation theory, shown in Figs. 3.2, 3.3, and
3.4, Note that the higher harmonics are progressively smaller; this
is the a posteriori justification for truncating the numerical solution
at a finite number of harmonics. The differences between the numerical
solutions and regular perturbation theory are mainly due to the fact
that the latter overestimates the first harmonic amplitude.
For incident waves with amplitudes greater than A1 = .03, the
numerical schemefrequern1y requires more than ten harmonics to achieve
convergence. This is not entirely satisfactory because the value of 6
associated with the tenth harmonic, at say k = 4.36, is 2.18. This
exceeds the value 6 =7T/2 = 1.5708 at which cross modes first appear in
the bay. Thus, at larger amplitudes, the assumption of one-dimension-
ality in the far field of the harbor is no longer tenable. Further-
more, the peak-to-trough wave heights predicted by the nonlinear theory
exceed half the depth. At such large amplitudes, the waves are no longer
correctly described by the Boussinesq equations.
- 134 -
Thus far, we have not accounted for entrance loss due to separation,
whose primary effect is to reduce the resonant amplitude of the first
harmonic. The magnitude of this reduction depends upon the average
friction factor f, which in turn depends strongly upon the particular
shape of the harbor entrance. As argued in Section 4.3.2, we expect
? to be less than .7 for the long narrow bay. To estimate the influence
of separation, we shall use the conservative value f = .35, which
coincides with the value chosen for the laboratory experiments.
Of practical interest is the relative importance of nonlinearity
and separation in reducing the resonant amplitudes of the first harmonic.
Fig. 5.9 shows the variation of 1n1(-z| with increasing incident
wave amplitudes, A1 I) for each of the first three resonant harbor modes.
The dashed lines represent linearized theory; their slopes are just the
amplification factors of Table 5.3. The symbols "A" and "x" indicate
the values of Inj-z)I according to inviscid nonlinear theory (f = 0)
and nonlinear theory with separation (f = .35), respectively. Note
that, for the first resonance mode (t = 1.41),, the combined reduction
due to both separation and nonlinearity is roughly twice that due to
nonlinearity alone; thus both effects are of roughly equal importance.
For the second and third harbor modes (z = 4.37 and 7.36), the effect
of nonlinearity dominates the effect of separation. This is reasonable
since nonlinearity increases with increasing z, whereas separation does
not.
The spatial dependence of the harbor response according to numerical
theory with f = .35 and |A1 | = .03 is shown in Figs. 5.10(a)-(c),
5.11(a)-(c), and 5.12(a)-(c), for each of the first three resonant
- 135 -
modes. These figures may be compared with the predictions of inviscid
nonlinear theory, shown in Figs. 5.6, 5.7, and 5.8. As expected, the
inclusion of separation losses leads to a reduction in the predicted
wave heights. For the first resonant mode, this reduction is sub-
stantial; for the second and third resonant modes, the viscous theory
agrees fairly closely with the inviscid theory.
When the incident wave is not monochromatic, nonlinear interactions
generate a great many new frequencies, some of which may be resonated by
the linear mechanism. Consider, for example, an incident wave having
A1 = 0 but A2 and A3 nonzero, as in Fig. 5.13(a). The incident wave
spectrum is dominated by two frequencies, 2w and 3w, where w coincides
with the first resonant harbor mode, viz. z= wL//gh = 1.41. The
harbor response according to linearized theory is simply an amplification
of A3, corresponding to the resonant mode with k = 4.36, and is shown in
Fig. 5.13(b). However, nonlinear theory allows for the generation of
additional harmonics; among them the first and fifth harmonics, both
of which coincide with resonant harbor modes. Fig. 5.13(c) shows the
harbor response at x = -z predicted by inviscid nonlinear theory. Nine
harmonics are present in the numerical solution. This means that the
wave spectrum of an incident wave can be significantly changed by non-
linearity. This possibility, somewhat similar to subharmonic resonance
in simpler physical systems, can be important to the design of mooring
systems for ships in the harbor.
- 136 -
-I
.2
I
9.
I
Ct)
Fig. 5.5&-f Inviscid nonlinear theory
2a = 100 m., L = 1000 m.,
numerical solutions.
= 0):
A1 = .03).
Frequency response of large-scale harbor (h 20 m.,
Number of harmonics, N*, required in the
N *
s
-- )
I
10
mmim | | | iI I .I.I. I I I I I1 .& &SJ J A A~.I.4..bj.A i..~sriikvAa UI I U9 a
O
I
'pi
00 .
L Af
WLt/r
Fig. 5.5b. First harmonic amplitude at the back wall according to nonlinear theory (A) and
regular perturbation theory (-). (See Fig. 5.5a).
I
4%:
8./0
Fig. 5.5c. Second harmonic amplitude at the back wall according to nonlinear theory (+) and
regular perturbation theory (- ). (See Fig. 5.5a).
I.01-
a
it a
3sv% x x"x orLa a a ax a KW
I aaII 1 1
02 I. 2 P/0
Fig. 5.5d. Third harmonic amplitude at the back wall according to nonlinear theory (x).
(See Fig. 5.5a).
.1
'7~ (-#)
. 09
I
(41 L/v37
Fig. 5.5e. Zeroth harmonic at the back wall according to nonlinear theory (A) and regular
perturbation theory (-). (See Fig. 5.5a).
I I I I I I I
SF
I I
l0
I
.5"
.4'--
.3-
I)
.1I
Fig. 5.6a.
'ja
.I I AF
.8
-F
-. 2lia
T
Inviscid nonlinear theory (~ = 0): First resonant mode of large-scale harbor
(h = 20 m., 2a = 100 m., L = 1000 m., A = .03, z = 1.41). Spatial variation
of the lower harmonics.
*vac 3 1131
.016-~ 
*
.oIC
.0105
-I .os.
-X'S
Fig. 5.6b. Spatial variation of the higher harmonics. (See Fig. 5.6a).
.5.
-. 5
am /.
Fig. 5.6c. Surface profiles at times t = m7/4, m = 0,1,2...7.
W4e
S21
R R 3 t
Fig. 5.7a. Inviscid nonlinear theory C? = 0): Second resonant mode of large-scale harbor
(h = 20 m., 2a = 100 m., L = 1000 m., A1 = .03, z = 4.36). Spatial variation
of the lower harmonics.
C)010 I3
.oloS
I 1 t
3dt-
Fig. 5.7b. Spatial variation of the higher harmonics. (See Fig. 5.7a).
evesoo
oo 0. COEM 0'
0..*.gof f. o P
w -e
4::m
I01103600
Fig. 5.7c. Surface profiles at times t = mfr/4, m = 0,l,2...7. (See Fig. 5.7a).
4t.ob
00
C2 -7
Fig. 5.8a. Inviscid nonlinear theory (~ = 0): Third resonant mode of large-scale harbor
(h = 20 m., 2a = 100 m., L = 1000 m., A1 = .03, t =1.36). Spatial variation
of the lower harmonics.
.ol -
- 13
.ooC
/ a 3 7
Fig. 5.8b. Spatial variation of the higher harmonics. (See Fig. 5.8a).
.II
C,
00g
+aow
3 4L- 7
Fig. 5.8c. Surface profiles at times t = mir/4, m = 0,1,2...7. (See Fig. 5,8a).
040
-V0
60
#40
4e
.3 -
-
a
a
a
I
I =to+H
.00
It
044
a mom;2:3
POWe
wool,,
*40# ~.k7.3L~g- 0 *a
- S
r 5
I a
.*a
a I A
A1
First harmonic amplitude at the back wall of the large-scale harbor vs. incident
wave amplitude. (l----): linear theory. ( ): inviscid nonlinear theory U 0).
(x): nonlinear theory with separation loss (U = .35).
04e
M"
-'
I ,
-p
Fig. 5.9.
x
4w
spow
4opp-
Ono"
moomd"
, s O
5.104.. Nonlinear theory with separation loss (f .35): First resonant mode of large-scale
harbor (h = 20 m., 2a = 100 m., L = 1000 m., A1 = .03, 9, = 1.41). Spatial variation
of the lower harmonics.
."os
13t
pa(. 5l.a.
5.l0b. Spatial variation of the higher harmonics. (See Fig. 5.l0a).
x3
6'1
5.1Oc. Surface profiles at times t = m7/4, m = 0,1,2...7. (See Fig. 5.10a).
ILIL
CY,4
Fig. 5.lla. Nonlinear theory with separation loss (~ = .35): Second resonant mode of large-
scale harbor (h = 20 m., 2a = 100 m., L = 1000 m., A1 = .03, z = 4.36).
Spatial variation of the lower harmonics.
*a em
.010 -
.00
a 3 4
Spatial variation of the higher harmonics, (See Fig. 5.lla).Fig. 5.11b.
'fx)t)
U,
R
Won
%Ira 
dOo
\40.
4MUD owe con" 
WOW GPM low qpftn
W-M 4 moo am" *AND O*ft
Goo WOM IWO ago* %ft arm* 
4w wom
am, 4wb 4,m 40M wow
WOW
0#90
-4r-r-
- -
0*0
.2 3
Fig. 5.llc. Surface profiles at times t = mfr/4, m = 0,1,2...7. (See Fig. 5.lla).
010 -
.00
Fig. 5.12a. Nonlinear theory with separation loss (f = .35): Third resonant mode of large-
scale harbor (h = 20 m. , 2a = 100 m., L = 1000 m., A1 = .03, = 7.36).
Spatial variation of the lower harmonics.
.010-
113
I a 3 9- 4 7
Spatial variation of the higher harmonics. (See Fig. 5.12a).Fig. 5.1 2b.
06
ra r 0-r*ra
o2 3+7
Fig. 5.12c. Surface profiles at times t = m-f/4, m = 0,1,2...7. (See Fig. 5.12a).
(a) INCIDENT WAVE
ols
1Aql
*05
(a) HARBOR RESPONSE
(NONLINEAR, INVISCID THEORY)
.003
1 .OvaI .t0.3 ,so'
3
i
Fig. 5.13. Response of large-scale
L = 1000 m., A2 = .050,
harbor to nonmonochromatic incident wave (h = 20 m.,, 2a = 100 m.,
A3 = .015, z = 1.41).
.o-s
(b) HARBOR RESPONSE
(LINEAR THEORY)
I~
p
.o47
I
.ol7
I'
t 0?
s3 + -7 9
,Oli
96 7
CHAPTER VI
CONCLUDING REMARKS
6.1 Validity of the Various Approximations
The proposed nonlinear theory for harbor resonance hinges on the
smallness of three parameters--E, 2, and 6--corresponding to the
assumptions of weak nonlinearity, weak dispersion, and a narrow harbor
entrance. The assumption c n p2 << 1 is required by the depth-averaged
Boussinesq equations; the smallness of 6 permits a linear treatment of
the radiated wave in the ocean and a one-dimensional treatment of the
nonlinear response in the long, narrow bay.
Two major approximations are made in the present work, namely to
decompose the solution into a finite number of harmonics and to use
an impedance boundary condition at the harbor entrance. The impedance
boundary condition is, however, only valid for the lowest few harmonics,
hence the calculated higher harmonics may be quantitatively unreliable
even if numerical convergence is achieved by including many harmonics.
Nonetheless, the lowest few harmonics are expected to be substantially
correct.
In a large scale harbor, the depth is usually comparable to or
smaller than the entrance width, so that p2<s 62 << 6. Therefore the
condition 6 << 1 is generally more restrictive than the condition
P2 << 1; in other words, the assumption 6 <<1 is likely to break down
before the assumption «2 << 1. For this reason, future efforts should
2first be directed at removing the restriction on 6, rather than on p
- 162 -
More will be said about this in Section 6.3.
6.2 On the Importance of Nonlinear and Frictional Effects
For monochromatic waves incident on a long, narrow bay, the princi-
pal effect of nonlinearity is to transfer energy from the linearly
resonated first harmonic to higher harmonics. The efficiency with which
higher harmonics are generated depends upon e for a harbor having
2 2 2
z << l/2 and upon s/2 for a harbor with z 1/p 2. In some cases,
harmonic generation is enchanced by linear resonance occurring at an
integral multiple of the fundamental frequency. For example, when the
quarter-wave mode is excited by the incident wave, odd harmonics of
the incident wave frequency generally lie close to the higher harbor
modes and are therefore amplified. Similarly, when the fundamental
frequency resonates one of the higher harbor modes, the influence of
dispersion may cause the second harmonic to be linearly resonated also.
This phenomenon gives partial explanation for some observed features in
the experiments, for the second and third harbor lengths studied.
The energy relation of Section 3.4 indicates that energy is trans-
ferred to the higher harmonics at the expense of the first harmonic.
The resonant amplification of the first harmonic, fn1(-z)/A1J, is
therefore found to decrease with increasing incident wave amplitudes.
Thus, linearized theory overestimates the resonant response of a harbor
subject to monochromatic incident waves.
When the incident wave is not monochromatic, the difference tone
of two high frequency components may excite a highly resonant harbor
- 163 -
mode occurring at a low frequency. In this case, nonlinear theory
indicates that the harbor response may in fact be more severe than that
predicted by linearized theory.
In all cases, harmonic distortion is to be expected from nonlin-
earity. In particular, the existence of a mean set-up, n(x), which
is highest at the back wall of the bay, is of practical importance,
even though it is usually quite small. Further, the details of the
frequency spectrum in the harbor may be important in designing mooring
systems for ships in the harbor.
The primary effect of friction is to increase the total damping in
the harbor. Boundary-layer dissipation is shown to be small compared
to radiation damping, except in unusually long harbors. The effect
of separation loss depends largely on the shape of the harbor entrance.
For an entrance with sharp corners, separation loss may be as important
as nonlinearity in damping the fundamental harmonic. For an asymmetrical
harbor entrance, separation may contribute substantially to the mean
set-up, nr(x); however, the generation of higher harmonics is usually
much less than that due to nonlinearity in the far field of the harbor.
Finally, separation is known to reduce the mass reactance of the entrance
and this may slightly increase the characteristic frequencies of the
resonant harbor modes.
6.3 Suggestions for Future Research
In this investigation, the long, narrow bay was chosen for its
simplicity, and certain assumptions were made for convenience in the
- 164 -
mathematical formulation. Below we sketch how some of these limitations
might be removed:
(1) Extension to harbors of arbitrary shape:
For a two-dimensional harbor, there is a greater variety of possible
nonlinear interactions than in the one-dimensional, long, narrow bay.
In particular, nonlinearity will couple waves intersecting at oblique
angles. Since the density of linearly resonant modes (in frequency
space) is greater in a two-dimensional harbor, harmonics generated by
nonlinearity are more likely to be linearly resonated. Thus, nonlinear
effects are expected to be even more important in a two-dimensional
harbor than in the long, narrow bay.
At present, numerical procedures exist for solving the linear theory
for harbors of arbitrary shape, see Lee (1971) and Chen and Mei (1975).
The latter employ a hybrid analytical-finite element technique which is
based upon a variational principle; in this regard, it may be useful to
note that the Boussinesq equation can also be cast in variational form;
see Whitham (1967b).
(2) Extension to slowly varying depth:
This involves a slight modification of the constant-depth
Boussinesq equations and should present no theoretical difficulty; see
Peregrine (1972).
The extension of the present theory to water of arbitrary depth
involves the complete, three-dimensional theory sketched in Section 1.2.
This is both difficult and unnecessary for most large harbors, which
- 165 -
typically resonate at frequencies that are well covered by the shallow-
water approximation.
(3) Extension to harbor with wide entrances:
In the present theory, the matching conditions at the entrance are
satisfied in an average way using an effective radiation impedance. This
is unsuitable for wide entrances or for high frequency waves, such as
might be generated through nonlinearity. A numerical procedure could
be used to ensure exact matching at every point in the entrance.
(4) Entrance Loss:
When flow separation at sharp corners is significant, inviscid
theory cannot be used to describe the flow field in the entrance. At
present, an understanding of separation is perhaps more accessible through
experiment than theory.
- 166 -
LIST OF REFERENCES
Batchelor, G.K., An Introduction to Fluid Dynamics, Cambridge University
Press, Cambridge, England, 1974, pp. 353-358.
Benney, D.J. and Luke, J. C., "On the Interactions of Permanent Waves of
Finite Amplitude", J. of Mathematics and Physics, Vol. 43,
1964, pp. 309-313.
Benney, D. J. and Niell, A. M., "Apparent Resonances of Weakly Nonlinear
Standing Waves", J. of Mathematics and Physics, Vol. 41,
1962, pp. 254-263.
Biesel, F., "The Similitude of Scale Models for the Study of Seiches
in Harbors", Proc. of the Fifth Conference on Coastal Engin-
eering, Grenoble, France, 1954; Council Wave Research,
Berkeley, California, 1955, pp. 95-118.
Biesel, F., "Radiating Secc d-Order Phenomena in Gravity Waves",
Proc. of Tenth Congress of the International Association for
Hydraulic Research, Vol. 1, London, England, 1963, pp. 197-203.
Bryant, P. J., "Periodic Waves in Shallow Water", J. of Fluid Mechanics,
Vol. 59, Part 4, 1973, pp. 625-644.
Chen, H. S. and Mei, C. C., "Oscillations and Wave Forces in a Man-Made
Harbor in the Open Sea", Proc. of the Tenth Symposium on Naval
Hydrodynamics, June, 1974 (in press).
Chu, V. H. and Mei, C. C., "On Slowly Varying Stokes Waves"', J. of Fluid
Mechanics, Vol. 41, Part 4, 1970, pp. 873-887.
Daily, J. W. and Harleman, D. R. F., Fluid Dynamics, Addison-Wesley,
Reading, Mass., 1966, pp. 314-319.
Gaillard, P., "Des Oscillations Non Lineaires des Eaux Portuaires",
Houille Blanche, Vol. 15, 1960, pp. 164-172.
Hwang, L. S. and Tuck, E. 0., "On the Oscillations of Harbours of
Arbitrary Shape", J. of Fluid Mechanics, Vol. 42, Part 3,
1970, pp. 447-464.
Ingard, U., "Nonlinear Distortion of Sound Transmitted through an Ori-
fice", J. of the Acoustical Society of America, Vol 48, No. 1
(Part 1), 1970, pp. 32-33.
Ingard, U. and Ising, H., "Acoustical Nonlinearity of an Orifice", J. of
the Acoustical Society of America, Vol. 42, 1967, pp. 6-17.
- 167 -
Ippen, A.
Ippen, A.
T. and Goda, Y., "Wave Induced Oscillations in Harbours: The
Solution for a Rectangular Harbour Connected to the Open Sea",
Technical Report No. 59, Ralph M. Parsons Laboratory for
Water Resources and Hydrodynamics, Dept. of Civil Engineering,
M.I.T., Cambridge, Mass., July 1963.
T. and Raichlen, R., "Wave Induced Oscillations in Harbors:
The Problem of Coupling of Highly Reflective Basins",
Technical Report No. 49, Ralph M. Parsons Laboratory for
Water Resources and Hydrodynamics, Dept. of Civil Engineering,
M.I.T., Cambridge, Mass., 1962.
Ito, Y., "On the Effect of Tsunami-Breakwater", Coastal Engineering in
Japan, Vol. 13, 1970, pp. 89-102.
Jarlan, G. E., "The Application of Acoustic Theory to the Reflective
Properties of Coastal Engineering Structures", Quarterly
Bulletin, National Research Council of Canada, Ottawa,
Canada, No. 1965, 1 (Apr.), pp. 23-63.
Jonsson, I. G. and Carlsen, N. A., "Experimental and Theoretical Investi-
gation in an Oscillatory Turbulent Boundary Layer", J. of Hy-
draulic Research. Vol. 14, No. 1, 1976, pp. 45-60.
Kim, Y. Y and Hanratty, J. J., "Weak Quadratic Interations of Two-
Dimensional Waves", J. of Fluid Mechanics, Vol. 50, Part 1,
1971, pp. 107-132.
Lee, J. J., "Wave-Induced Oscillations in Harbors of Arbitrary Geometry",
J. of Fluid Mechanics, Vol. 45, Part 2, 1971, pp. 375-394.
Longuet-Higgins, M. S., "Resonant Interaction Between Two Trains of
Gravity Waves", J. of Fluid Mechanics, Vol . 12, 1962, pp.
321-332.
Mei, C. C. and Chen, H. S., "Hybrid-Element Method for Water Waves",
Proc. of the Symposium on Modeling Techniques, American
Society of Civil Engineers, 1975.
Mei, C. C., Lui, P. L-F., and Ippen, A. T., "Quadratic Loss and Scatter-
ing of Long Waves", J. of the Waterways, Harbors, and CoastaZ
Engineering Division, WW3, 1974, pp. 217-239.
Mei, C. C. and Unluata, U., "Harmonic Generation in Shallow Water Waves",
Waves on Beaches and Resulting Sediment Transport, Academic
Press, New York, 1972, pp. 181-202.
Mei, C. C. and Unluata, U., "Resonant Scattering by a Harbor with Two
Coupled Basins", J. of Engineering Mathematics, Vol. 10,
No. 4, 1976, pp. 333-353.
- 168 -
Miles, J. W. and Munk, W. H., "Harbor Paradox", J. of the Waterways
and Harbors Division, Proc. of the American Society of Civil
Engineers, Vol 87, 1961, pp. 111-130.
Morse, P. M. and Feshback, H., Methods of Theoretical Physics, Vol, I,
McGraw-Hill, New York, 1953, pp. 814-815.
Peregrine, D. H., "Equations for Water Waves and the Approximations
behind Them", Waves on Beaches and Resulting Sediment Trans-
port, Academic Press, New York, 1972, pp. 95-121.
Phillips,
Phillips,
Phillips,
0. M., "On the Dynamics of Unsteady Gravity Waves of Finite
Amplitude. Part I. The Elementary Interactions", J. of Fluid
Mechanics, Vol. 9, 1960, pp. 193-217.
0. M., The Dynamics of the Upper Ocean, Cambridge University
Press, Cambridge, 1966, pp. 42-43.
0. M., "Nonlinear Dispersive Waves", Annual Review of Fluid
Mechanics, Vol. 6, Annual Reviews, Inc., Palo Alto, California,
1974, pp. 93-110.
F. and Ippen, A. T., "Wave Induced Oscillations
J. of the Hydraulics Division, American Society
Engineers, Vol. 91, No. HY2, 1965, pp. 1-26.
in Harbours",
of Civil
Stoker, J. J., Water Waves, Interscience Publishers, Inc., New York,
1957, pp. 174-181.
Tadjbakhsh, I. and Keller, J. B., "Standing Surface Waves of Finite
Amplitude", J. of Fluid Mechanics, Vol. 8, 1960, pp. 442-451.
Terrett, F. L., Osorio, J. D. C., and Lean, G. H., "Model Studies of a
Perforated Breakwater", Proc. of the Eleventh Conference on
Coastal Engineering,.1968, pp. 1104-1109.
Unluata, U. and Mei, C. C., "Long Wave Excitation in Harbours--An
Analytical Study", Technical Report No. 171, Ralph M. Parsons
Laboratory for Water Resources and Hydrodynamics, Dept. of
Civil Engineering, M.I.T., Cambridge, Mass., 1973.
Unluata, U. and Mei, C. C., "Effects of Entrance Loss on Harbor Oscilla-
tions", J. of the Waterways, Harbors and Coastal Engineering
Division, Proc. of the American Society of Civil Engineers,
Vol. 101, No. WW2, 1975.
Ursell, F., "Edge Waves on a Sloping Beach", Proc. Royal Soc. A, Vol.
214, 1952, p. 79.
Verhagen, J. and van Wijngaarden, L., "Nonlinear Oscillations of Fluid
in a Container", J. of Fluid Mechanics, Vol. 22, 1965, pp. 737-
751.
- 169 -
Raichlen,
Whitham, G. B.,
Waves
1965,
"A General Approach to Linear and Non-Linear Dispersive
Using a Lagrangian", J. of Fluid Mechanics, Vol. 22,
pp. 273-283.
Whitham, G. B., "Nonlinear Dispersion of Water Waves", J. of Fluid
Mechanics, Vol. 27, 1967, pp. 399-412.
Whitham, G. B., "Variational Methods and Applications to Water Waves",
Proc. Roy. Soc. of London, Series A, Vol. 299, 1967, pp. 6-25.
Wil son , B. N. , XIX International Navigation Congress, Secti on II,
Communication 1, London, 1957, pp. 13-61.
- 170 -
APPENDIX A
NUMERICAL ITERATION PROCEDURE FOR SOLVING A NONLINEAR,
TWO-POINT BOUNDARY VALUE PROBLEM
The following is a FORTRAN IV computer program which solves
the boundary value problem posed in Section 3.5. The variable names
follow closely the notation used in Section 3.5. The principal variables
are:
XLONG -- harbor length (L)
XWIDE -- Harbor width (2a)
XDEEP -- harbor depth (h)
FRICT -- friction factor (?) in separation loss theory
NHMAX -- maximum number of harmonics
NH -- number of harmonics at a given point in the iteration
procedure. NH < NHMAX.
Z -- an array of size NHMAX containing the radiation
impedance for each harmonic
A -- an array of size NHMAX containing twice the incident
wave amplitude for each harmonic
R,DR -- arrays of size (NHMAX)2 containing the coupling
constants for the nonlinear equations
B an array of size (NPT)(NHMAX) containing the current
solution for in(xj), n = 1.2...NH, J = 1,2...NPT
E -- an array of size (NPT)(NHMAX)2 containing the NH
basis solutions en(xj), n = 1,2...NH, J = 1,2...NPT
On the IBM 370 computer at the M.I.T. Information Processing
Center, the total storage required by the program was found to be
Storage (in Kilobytes) = 52 + 8(NPT)(NHMAX)(NHMAX + 1)/1024
- 171 -
The execution time for an iteration with NH harmonics was found to be
approximately
Time (in psec) = 8 (NPT)(NH)3
For example, with NPT = 135 and NH = NHMAX = 10, the total storage
required is 168 kilobytes, and the execution time per iteration is
approximately 1.08 seconds.
- 172 -
C*****MAIN EPOGRAM: CCMPUTES INPUTS TO SUBROUTINE BVP
CCMPLFX A (20)jZ (20)
REAL W(2C)
CCMMCN /AREA 1/AIZ,Wr,H,NPTNHMAXFRICT
DATA XDEEP,XWIDE,XLONG/20.,100.,1000./
CALL TIMING (ISTART)
READ(5,*) ITIMF
ISTART=ISTART+6000*ITIME
NHMAX=10
WPITE(6,902) NUNAXXDEFPXWIDE,XLCNG
WFITE (7,902) NHMAXXDFEPXWIDEXLONG
902 FORMAT(' NiMAKXDEEP,XWIDE,XLCNG=',I3,3F10.1)
DO 1 T=2,NHMAX
1 A(I)= (0.,0.)
CALL COUJLE
5 READ(5,*) XL1,XL2,NXLXHIGHIPT,NH1,FRICT
IF (NXL .LT. 0) G0 TC 999
FRICT=FRICT*4./(3.*3.1415926)
A(1)=(1.,0.)*XHIGH/XDEEP
DO 20 J=1,NXL
XL = XL1+(J-1)*XL2
XA = XL*.5*XWIDE/XLCNG
XH = XL*XDEEP/XLONG
NPT=2+IPT*XL/.0628
HXL/ (NPT-i)
WFIT"El(6,901) NPT,XEIG,XLXAXH,H,FRICT
W RIT E (7, 901) NPTFXHIGH,XLXA,XHfH,FRICT
901 FORMAT (' NPT, XH IGH, XL,XAXH ,H=',I3,F6.2, 4F10.6/' FRICT='F10.6)
DO 10 I=1,NHMAX
Q=I*SQRT(1.+I*I*XH*XH/3.)
Z(i)=(0.,-1.)*XAt+.637*XA*(ALOG(Q*XA)-.875)
W (I) =1.+. 5*H*H*o*Q
Q=1./CABS (COS (Q*XL) +Q*Z (I) *SIN (Q*XL))
10 WPITF(6,900) I,Q
900 FORMAT (' AMPLIFICATICN OF HARMONIC I,12,'=',F8.4)
CALL PVP(NH1,ISTART)
MAIN0001
MAIN0002
MAIN0003
MAIN0004
MAIN0005
MAIN0006
MAIN0007
MAIN0008
MAIN0009
MAIN0010
MAIN0011
MAIN0012
MAIN0013
MAIN0014
MAIN0015
MAIN0016
MAIN0017
MAIN0018
MAIN0019
MAIN0020
MAIN0021
MAIN0022
MAIN0023
MAIN0024
MAIN0025
MAIN0026
MAIN0027
MAIN0028
MAIN0029
MAINO030
MAIN0031
MAIN0032
MAIN0033
MAIN0034
MAIN0035
MAIN0336
IF (NPT .LT. 0) GO TO 999
20 CONTINUE
GC TO 5
999 STOP
END
MAIN0037
MAIN0038
HAIN0039
MAIN0040
MAIN0041
-J:
SUBROUTINE CCUPLL
C*****COMPUTES COUPLING CCNSTANTS FOR TH'r NONLINEAR EQUATIONS
REAL R(20,20),DR(20,20)
COMMCN /AiEA2/R,DR
DC 10 1=1,20
DC 5 L=1,20
IF (L .EQ. I) GO TC 5
Q=.5-.25* (1/(1+3*(2*L-I)))
R (I,L)=Q*(I*I-2*L*L+2*I*L)
DR(I,)=-(I+L)*Q/(I-L)-(2*T-L)*Q/L
5 CCNTINUE
10 CONTINUE
fETUEN
END
COUP0001
COUP0002
COUP0003
COUP0004
COUP0005
COUP0006
COUP0007
COUPOO08
COUPOO09
COUPO010
Coupooll
COUP0012
COUP0013
COUP0014
C,,
SUBROUTINE BVP(NH1,ISTART)
C*****SOLVES A BOUNDARY VALUL PROBLEN CONSISTING OF N COUPLED,
C*****NONLINEAF, SFCONE-CREEP LIFFERENTIAL EQUATIONS SUBJECT TO
C*****TWO-POINT BCUNDARY CONDITIONS
COMPLEX B(1350),E(13500)
COMPLEX S(20,20),C(20),Z(20),CCNJGCQ,A(20)
REAL W (20) , P(20 ,20) ,DR(20,20)
REAL BPRINT(10)
COMMON /ALEA1/A,Z,W,H,NPTNHMAX,FRICT/AREA2/RDR/AREA3/SC
IF (NH1 .GT. 0) G TO 5
Dc 10 J=1,NPT
JX= (J-1)*NHIIAX
DO 10 I=1,NHMAX
10 B(I+JX)=(O.,0.)
S NP=NPT-1
NHP=NHMAX*NPT
DO 80 NH=1,NIIMAX
QRELAX=1.
ITER=O
15 ITER=ITEE+1
IF (ITEB .GT. 50) GO TO 62
CALL TIMING(ISTCP)
ILEFT=ISIART-ISTOP
WRITE(6,899) ILEfT,NH
899 FORMAT(115X,I10,2X,13)
IF (ILEFT .LT. 1000) GC TO 81
DC 50 K=1,NH
KX=NHAP* (K-1)
DO 20 I=1,NH
20 E(T+KX)=1/(1+3* (I-K))
DO 35 J=1,NP
JX= (J-1)*t NHMAX
NX=NHMAX*(1-1/J)
10 35 I=1,NH
IX=I+JX+KX
CQ= (0.,O.)
BVP 0001
BVP 0002
BVP 0003
BVP 0004
BVP 0005
BVP 0006
BVP 0007
BVP 0008
BVP 0009
BVP 0010
BVP 0011
BVP 0012
BVP 0013
BVP 0014
BVP 0015
BVP 0016
BVP 0017
BVP 0018
BVP 0019
BVP 0020
BVP 0021
BVP 0022
BVP 0023
BVP 0024
BVP 0025
BVP 0026
BVP 0027
BVP 0028
BVP 0029
BVP 0030
BVP 0031
BVP 0032
BVP 0033
BVP 0034
BVP 0035
BVP 0036
12=1/2
DC 300 L=1,NH
IF ((I2/L)-(L/(I+l))) 200,303,100
100 CQ=CQ+R (IL)*13(L+JX)*E(IX-L)+DR(IL)*(B(L+Jh)-B(L+JX-NX))*(E(TX-L)
* -E (IX-L-NX))/(H*H)
GO TO 300
203 LX=T+JX
CQ=CQ+P(1,L)*CONJG(B(LX-T))*E(LX+KX)
* +D(IL)*CONJG(B(LX-I) -B(LX-I-NX))*(E(LX+KX)-E (LX+KX-NX))/(H*H)
300 CONTINUE
CQ=(.5*H*H*CQ+L (IX))/W(I)
35 E (IX+N'iMAX)=CQ+(NX/NH iMAX)*(CQ-r (IX-NX))
C (K)=h(K)
DO 50 =1,NH
JX=I+ NBP-NHIAX
IX=JX+(K-1)*NHP
CQ=Z (I)+FRICT*9. ,-1.) *CABS (B(JX) -B(JX-NHMAX) )/(H*T*I)
50 S(IK)=E(IX)-CQ*(E(IX)-E(IX-NH4AX))/H
CALL SYSN(NH)
Q=QRELAX
ICVRG=O
IGCOD=O
DO 70 J=1,NPT
JX=(J-1)* NHMAX
DC 65 T=1,NH
IX=T+JX
CQ=-Q*t B(IX) +Q'rC(1)*E(IX)
IF (NH .LT. 2) GO TC 61
DC 60 K=2,NH
60 CQ=CQ+Q*C(K)*F(IX+ (K-1) *NfHP)
61 ICVRG=CVRG+CABS(CQ)/.001
65 B(IX)=B(IX)+CQ
70 IGOOD=IGCCD+CABS (B (NH+JX))/.001
WFITE(6,900) (B (I),I=1,NH)
900 FORYMAT(1X,3P10F12.6)
QRELAX=.5
BVP 0037
BVP 0038
BVP 0039
BVP 0040
BVP 0041
BVP 0042
BVP 0043
BVP 0044
BVP 0045
BVP 0046
BVP 0047
BVP 0048
BVP 0049
BVP 0050
BVP 0051
BVP 0052
BVP 0053
BVP 0054
BVP 0055
BVP 0056
BVP 0057
BVP 0058
BVP 0059
BVP 0060
BVP 0061
BVP 0062
BVP 0063
BVP 0064
BVP 0065
BVP 0066
BVP 0067
BVP 0068
BVP 0069
BVP 0070
BVP 0071
BVP 0072
IF (ICVRG .GT. 0) GC TO 15 BVP 0073
IF (IGCOD .EQ. 0) GC TO d5 BVP 0074
30 CONTINUF BVP 0075
NH=NMAX BVP 0076
WRITE(6,%80) BVP 0077
WRITE(7,980) BVP 0078
980 FORMAT(' MOrE HARMONICS UFPDED') BVP 0079
GO TO 85 BVP 0080
81 WRITE(6,961) BVP 0081
WRITE(7,981) BVP 0082
981 FOtPMAT(' TIME RAN CUI') BVP 0083
NPT=-1 BVP 0084
GC TO 85 BVP 0085
82 WRITE(6,$82) BVP 0086
WRITE (7,982) BVP 0087
982 FORMAT (' NOT CONVERGFNT') BVP 0088
85 WRITE(6,901) Nil BVP 0089
WRITE(7,901) NH BVP 0090
901 FCJRMAT(' NH=',13) BVP 0091
o 9C J=1,NPT BVP 0092
DO 91 I=1,NH BVP 0093
91 BPRINT(I)=CAES(B(I+(J-1)*NHMAX)) BVP 0094
WRITF(6,902) (BPINT(I),B(I+(J-1)*NHMAX), 1=1,NP) BVP 0095
902 FORMAT(1X,3P15F8.3) BVP 0096
90 CONTINUE BVP 0097
WRITE(6,903) BVP 0098
903 FCRMAT(1HI1) BVP 0099
RETURN BVP 0100
END BVP 0101
00
SUBROUTINE SYSN (NH)
C*****SOLVES A SYSTEM OF N LINEAR ALGEERAIC EQUATIONS
COMPLEX S (20,20) ,C (2C) ,PVT
COMMON /AREA3/S,C
D0 50 K=1,NH
PVT=S(K,K)
IF (CABS (PVT) .LI. 1.E-5) WPTTE(1,90)
C (K)=C (K) /PVT
D0 20 J=K,NH
20 S(K,J)=S(K,J)/PVT
DO 40 1=1,NH
IF (I .EQ. K) GO TO 40
PVT=S (1,4)
C(I)=C(I)-PVT*C(K)
DO 30 J=K,N]
30 S(I,J)=s(IJ)-PVT*S(K,J)
40 CONTINUE
50 CONTINUE
90 FCRMAT (' **SYSN SINGULAR**')
RETURN
END
SYSN0001
SYSNO002
SYSNO003
SYSNO004
SYSN0005
SYSNO006
SYSNO007
SYSNO008
SYSNO009
SYSNO010
SYSN0011
SYSN0012
SYSNO013
SYSNO014
SYSNO015
SYSNO016
SYSNO017
SYSNO018
SYSNO019
SYSNO020
SYSNO021
to
C*****SAMPLE INPUT DAT
20
.25 .25 4 .6 1 -1 *35
1.21 .04 11 .6 1 -1 .35
1.75 .25 8 .6 1 -1 .35
3.72 .12 11 .6 1 -1 .35
5.25 .25 5 .6 1 -1 .35
6.36 .20 11 .6 1 -1 .35
1.41 1. 1 .6 1 -1 0.
1.41 1. 1 .5 1 -1 0.
1.41 1. 1 .4 1 -1 0.
1.41 1. 1 .3 1 - 1 0.
1.41 1. 1 .2 1 -1 0.
1.41 1. 1 .1 1 -1 0.
1.41 1. 1 .6 1 -1 .35
1.41 1. 1 .5 1 -1 .35
1.41 1. 1 .4 1 -1 .35
1.41 1. 1 .3 1 -1 .35
1.41 1. 1 .2 1 -1 .35
1.41 1. 1 .1 1 -1 .35
4.36 1. 1 .6 1 -1 0.
4.36 1. 1 .5 1 -1 0.
4.36 1. 1 .4 1 -1 0.
4.36 1. 1 .3 1 -1 0.
4.36 1. 1 .2 1 -1 0.
4.36 1. 1 .1 1 -1 0.
4.36 1. 1 .6 1 -1 .35
4.36 1. 1 .5 1 -1 .35
4.36 1. 1 .4 1 -1 .35
4.36 1. 1 .3 1 -1 .35
4.36 1. 1 .2 1 -1 .35
4.36 1. 1 .1 1 -1 .35
7.36 1. 1 .6 1 -1 0.
7.36 1. 1 .5 1 -1 0.
7.36 1. 1 .4 1 -1 0.
7.36 1. 1 .3 1 -1 0.
(output for this case on p. 181)
DATA0001
DATA0002
DATA0003
DATA0004
DATA0005
DATA0006
DATA0007
DATA0008
DATA0009
DATA0010
DATA0011
DATA0 012
DATA0013
DATA0014
DATA0015
DATA0016
DATA0017
DATAO18
DATA0019
DATA0020
DATA0021
DATA0022
DATA0023
DATA0024
DATA0025
DATA0026
DATA0027
DATA0028
DATA0029
DATA0030
DATA0031
DATA0032
DATA0033
DATA0034
DATA0035
DATA0036
NFT,XHIGIXLXAXH,H= 24
FRICT= 0.0
AMPLIFICATION OF HARMONIC
AMPLIFICATION OF HARMONIC
AMPLIFICATION OF HARMONIC
AMPLIFICATION OF HARMONIC
AMPLIFICATION OF HARMONIC
AMPLIFICATION OF HARMONIC
AMPLIFICATION OF HARMONIC
AMPLIFICATION OF HARMONIC
AMPLIFICATION OF HARMONIC
AMPLIFICATION OF HARMONIC
18.729173 215.369761
18.729173 215.369761
-20.925488 209.342420
-19.937202 209.556282
-19.463349 209.658980
-19.661296 208.032846
-19.409839 208.172917
-19.073937 208.113015
-19.061219 208.125353
-19.050982 208.125830
NH= 5
208.996 -19.051 208.126
208.602 -19.016 207.734
207.423 -18.912 206.559
205.462 -18.739 204.606
202.724 -18.496 201.879
199.218 -18.184 198.387
194.955 -17.804 194.140
189.948 -17.356 189.153
184.212 -16.840 183.440
177.767 -16.258 177.022
170.634 -15.612 169.919
162.840 -14.901 162.157
154.412 -14.130 153.764
145.382 -13.300 144.772
135.785 -12.413 135.217
125.661 -11.474 125.136
115.050 -10.486 114.571
103.998 -9.453 103.567
92.551 -8.379 92.171
80.760 -7.269 80.432
68.677 -6.127 68.403
56.355 -4.960 56.137
43.850 -3.772 43.687
31.217 -2.570 31.111
0.30 1.410000 0.070500 0.028200 0.061304
1=
2=
3=
4=
5=
6=
7=
8=
9=
10=
14.3488
0.9706
3.5589
0.9543
2.1271
0.9675
1.4810
0.9950
1.1839
1.0105
-3.293194
-3.307885
SAMPLE OUTPUT
108635 1
108634 1
108634 2
.108631 2
108629 2
2.618542
2.323011
-3.315073 2.182728
-3.134884 1.597423 3.473077 -1.074277
-3.127156 1.602835 3.448527 -1.062287
108626 3
108621 3
108614
-3.032238 1.585308 3.693797 -0.573599 0.079876 0.300598
108601
-3.032259 1.588185 3.690494 -0.579309 0.080550 0.300239
108588
-3.027946 1.591547 3.683764 -0.564491 0.058331 0.301515 -0.156799 -0.128566
3.421
3.492
3.714
4.098
4.654
5.374
6.240
7.220
8.276
9.366
10.444
11.464
12.378
13.142
13.711
14.044
14.105
13.865
13.299
-3.028
-3.125
-3.414
-3.886
-4.526
-5.314
-6.222
-7.219
-8.270
-9.335
-10.373
-11.339
-12. 189
-12.881
-13.373
-13.626
-13.606
-13.284
-12.640
12.394 -11.658
11.144 -10.332
9.559 -8.666
7.667 -6.669
5.540 -4.363
1.592
1.559
1.462
1.301
1.080
0.803
0.474
0.099
-0.314
-0.756
-1.218
-1.688
-2.155
-2.606
-3.026
-3.402
-3.721
-3.970
-4.135
-4.207
-4.176
-4.035
-3.782
-3.414
3.727
3.663
3.473
3.166
2.755
2.262
1.716
1.175
0.793
0.885
1.341
1.875
2.384
2.828
3.185
3.440
3.586
3.617
3.532
3.332
3.018
2.595
2.070
1.458
3.684
3.619
3.428
3.117
2.699
2.190
1.610
0.983
0.334
-0.310
-0.1924
-1.481
-1.961
-2.342
-2.610
-2.756
-2.776
-2.672
-2.454
-2.136
-1.737
-1.282
-0.797
-0.309
-0.564
-0.563
-0.559
-0.556
-0.556
-0.567
-0.593
-0.642
-0.719
-0.829
-0.973
-1.150
-1.357
-1.586
-1.825
-2.059
-2.270
-2.438
-2.540
-2.557
-2.468
-2.256
-1.911
-1.425
0.307
0.314
0.332
0.360
0.390
0.414
0.425
0.414,
0.373
0.299
0.192
0.092
0.189
0.382
0.593
0.800
0.985
1.131
1.223
1.247
1.195
1.064
0.859
0.593
0.058
0.058
0.055
0.051
0.047
0.042
0.037
0.035
0.036
0.043
0.059
0.085
0.125
0.178
0.245
0.325
0.412
0.501
0.583
0.648
0.683
0.678
0.622
0.506
0.302
0.308
0.328
0.356
0. 397
0.412
0. 424
0.413
0.372
0.296
0.183
0.035
-0.142
-0.338
-0.540
-0.731
-0.895
-1.014
-1.075
-1.065
-0.980
-0.820
-0.592
-0.310
0.203
0. 195
0. 173
0. 139
0.099
0.063
0.053
0.075
0.101
0. 120
0.129
0.127
0.118
0.106
0.103
0.118
0.149
0. 188
0.227
0.257
0.269
0.258
0.217
0.148
-0.157
-0.152
-0. 138
-0.117
-0.091
-0.063
-0.037
-0.015
0.002
0.013
0.021
0.028
0.037
0.052
0.074
0.106
0.144
0.186
0.225
0.253
0.262
0.244
0.196
0.114
4
4
5
-0.129
-0.122
-0. 104
-0.076
-0.040
-0.001
0.038
0.073
0.101
0.120
0.127
0.124
0.112
0.093
0.072
0.052
0.037
0.031
0.034
0.046
0.064
0.082
0.095
0.094
-D
00
,--I
APPENDIX B
FOURIER ANALYSIS OF EXPERIMENTAL DATA
The following is a BASIC computer program which Fourier
analyzes the experimental data, while the experiment is in progress.
The execution of the program is remotely controlled by a switch, as
explained in section 4.1.3. The principal variables in the program are:
T -- number of'voltage readings taken in one complete
wave period
KI -- number of harmonics analyzed
C -- a calibration factor
F -- an array containing (T+l) voltage readings
A -- a discrete Fourier cosine transform
B -- a discrete Fourier sine transform
The instruction CALL(1,D,F) reads the external Hewlett Packard 3450A
Multi-Function Meter and stores the voltage reading under the variable
name D.
- 182 -
10 DIM F[30J
20 READ T,KI,C
30 DATA 22,4,1
40 LET T1=T+1
50 LET I=-2
55 LET 1=1+2
60 GOSUB 1000
70 FOR J=1 TO TI
80 CALL (1,D,F)
90 LET F[J]:D
100 NEXT J
110 GOSUB 3000
120 GOTO 55
130 STOP
1000 REM ****************** SWITCHING ROUTINE *************
1010 CALL (1,D,F)
1020 IF D >= 5 THEN 1040
1030 GOTO 1010
1040 CALL (I,D,F)
1050 IF D <= 5 THEN 1070
1060 GOTO 1040
1070 WAIT (500)
1080 RETURN
3000 REM ******************** FOURIER ANALYSIS *********
3010 PRINT I
3020 LET F[IJ1=(F(1J1+F[ TI1) /2
3030 LET H=6.28318/T
3040 FOR K= TO K!
3050 LET A=B:0
3060 FOR J=1 TO TI-1
3070 LET A=A+F[J]*COS((J-1)*K*H)
3080 LET B=B+F[JI*SIN(CJ-1)*K*H)
3090 NEXT J
3100 PRINT 2*SQR(A*A+B*B)/T/C,
3110 LET P[K]=ATN(B/A)+(1-SGN(A))*1.5708
3120 NEXT K
3130 PRINT
3140 PRINT " ",
3150 FOR K=2 TO K1
3160 LET P=P[K]-K*P[1]
3170 PRINT P-6.28318*INT(P/6.28318),
3180 NEXT K
3190 PRINT
3200 RETURN
9999 END
- 183 -
BIOGRAPHY
Steven R. Rogers was born on February 28, 1952, in Washington,
D. C. . His primary education was received at the Hebrew Academy of
Washington and his secondary education at the Talmudical Academy of
Baltimore and Montgomery Blair High School of Silver Spring.
He was awarded a Bachelor of Science degree in physics by
M. I. T. in June, 1972. His bachelor's thesis, "Ising Chain with
Random Magnetic Moments", appears in Studie's in Apptied Maathemaica,
Vol. 52, No. 2, 1973, pp. 163-173. A fellowship from the Fannie and
John Hertz Foundation of Los Angeles, California, has enabled him to
pursue graduate studies at M. I. T..
His employment experience includes work at the Technion--
Israel Institute of Technology (1974), the Weizmann Institute of
Science (1972), and the X-ray Exposure Control Laboratory (1970 and
1971), where he served as a commissionned officer in the U. S.
Public Health Service.
The author is a member of the American Physical Society and
the American Geophysical Union, and is a licensed amateur radio
operator. He is also a member of Phi Beta Kappa and Sigma Xi honorary
associations.
- 184 -