NUMERICAL MODELING OF WAVE-INDUCED CURRENTS USING A PARABOLIC
WAVE EQUATION
By
HARLEY STANFORD WINER
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1988
[UNIVERSITY OF FLORIDA LIBRARIES
ACKNOWLEDGEMENTS
The writing of a dissertation is a long process involving interaction with and support
from many people. The author would like to thank all those who in one way or another
offered encouragement, support, advice and helped in the problem solving process. The
contributions of those who are not mentioned by name are in no way diminished by the lack
of individual acknowledgement.
First, the author wishes to express deep appreciation to the members of his committee:
the chairman, Dr. Hsiang Wang, who enticed him to come to Florida and enroll in the
doctoral program, and the cochairman, Dr. James Kirby, who shared his knowledge of wave
equations and numerical modeling, along with Dr. Robert Dean, Dr. Bent Christensen and
Dr. Joseph Hammack. Each of these men has been both a teacher and a friend. The author
wishes to also thank Dr. Peter Sheng who at one time served on the supervisory committee
but due to a prior commitment was unable to be at the defense and officially be on the
committee.
The author would like to thank the staff of the Coastal and Oceanographic Engineering
Laboratory at the University of Florida and the staff of the department of Coastal and
Oceanographic Engineering who helped and aided the author in innumerable ways. Thanks
are also extended to Ms. Lillean Pieter who helped prepare the figures and Dr. Claudio
Neves, a fellow student, who manipulated the I^TgX software bo as to produce a document
in thesis mode.
A special acknowledgement is given to the author’s parents who encouraged and pushed
him to complete this doctorate. The author would also like to thank his teachers, friends
and fellow students who offered encouragement and interesting discussions.
Lastly the author thanks his soul-mate and partner, Esther DeJong, who managed to
be patient with each extension of the final date. Baby cakes, it’s PAinally Done. I love you.
m
TABLE OF CONTENTS
ACKNOWLEDGEMENTS Ü
LIST OF FIGURES vi
ABSTRACT ix
CHAPTERS
1 INTRODUCTION 1
1.1 Literature 4
1.2 General Description of the Model 8
2 CIRCULATION MODEL 16
2.1 Governing Equations 16
2.2 Radiation Stresses 24
3 WAVE MODEL 31
3.1 Governing Wave Equation 31
3.2 A Parabolic Approximation 36
3.3 Energy Dissipation Due to Wave Breaking 45
4 FINITE DIFFERENCE SCHEMES AND BOUNDARY CONDITIONS 49
4.1 Numerical Solution of the Governing Equations 49
4.2 Finite Differencing of the Parabolic Wave Equation 52
4.3 Boundary Conditions in the Circulation Model 56
4.4 Boundary Conditions in the Wave Model 59
5 RESULTS 63
5.1 Wave Set-up 64
5.2 Longshore Currents 68
iv
70
5.3 Waves On a Rip-Channel . . .
5.4 Shore-Perpendicular Breakwater 76
5.5 Shore Parallel Emerged Breakwater 88
5.6 Comparison with Experimental Results of Gourlay 92
6 CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER STUDY ... 103
6.1 Improved Wave Formulations 104
6.2 Including Wave Reflection 105
6.3 Extending the Model to Include Sediment Transport 105
APPENDICES
A DERIVATION OF THE DEPTH AVERAGED EQUATIONS 107
B DERIVATION OF THE RADIATION STRESS TERMS 113
C LAG RANG LAN DERIVATION OF THE GOVERNING WAVE EQUATION . . 117
D OBTAINING EQUATION (3.79) 120
BIBLIOGRAPHY 124
BIOGRAPHICAL SKETCH 129
v
LIST OF FIGURES
1.1 Computational domain of the model 8
1.2 Flow chart of computer program 11
1.3 Flow chart of circulation model 13
1.4 Flow chart of wave model 14
2.1 Longshore current profiles for different friction formulations for an 11
second, 1.028 meter wave at 17.28 degrees to a 1:25 plane beach 22
2.2 Experimental results for wave set-up and set-down. Reprinted from
Bowen, Inman, and Simmons (1968), Journal of Geophysical Research,
vol. 73(8) page 2573, copyright by the American Geophysical Union. . . 30
4.1 Definition sketch of grid nomenclature and coordinate system 56
4.2 Start-up function according to Eq. 4.43 with Ci=30 and C2=3 60
4.3 Waves passing through a lateral boundary unaffected by the boundary. 61
5.1 Comparison of the numerical solution of Vemulakonda et al. (1985) for
set-up with experimental data of Bowen (1968). Source: Vemulakonda
et al. (1985) 65
5.2 Comparison of numerical results from the model with experimental re¬
sults of Hansen and Svendsen (1979). Wave set-up and wave heights for
a 2 second wave on a beach slope of 1:34 and incident wave height of 7
centimeters 67
5.3 Non-dimensional longshore velocities vresus the non-dimensional dis¬
tance offshore from the analytical solution of Longuet-Higgins. Reprinted
from Longuet-Higgins (1970b), Journal of Geophysical Research, vol.
75(33) page 6793, copyright by the American Geophysical Union 68
5.4 Comparison of experimental results of Visser (1982) and numerical result
for longshore currents due to a 1.1 second wave on a slope of 1:20 for
the deep water conditions H0 = .065m and 0O = 20 degrees 69
5.5 Top: Depth contours, in meters, using Eq. 5.7 with the coordinates as
shown. Bottom: Depth contours as used in the present model. Grid
spacing is 5 meters 71
vi
5.6 Velocity vectors from the results of Ebersole and Dalrymple for a wave
angle of 30 degrees. Scale for velocity vectors is 1 inch = 2.87 m/sec.
Source: Ebersole and Dalrymple (1979) 72
5.7 Velocity vectors obtained from the present model using the same input
conditions as for Ebersole and Dalrymple 73
5.8 Mean water surface elevation contours showing that the water surface is
lower above the trough. Values indicated are in millimeters 74
5.9 Depth contours in meters for the case of a bar with a gap 75
5.10 Velocity vectors for the case of waves over a bar with a gap. Dashed
lines indicate depth contours. Grid spacing equal to 10 meters 75
5.11 Depth contours for the case of normally incident waves upon a shoal.
Grid spacing is .05 meters. Contour interval is .02 meters 77
5.12 Velocity vectors for the case of waves breaking upon a shoal 78
5.13 Maximum velocity across the top of a shoal versus incident wave height
for waves over a shoal on a plane beach. Incident wave height less than
9.1 centimeters are non-breaking waves. Incident wave heights greater
than 9.1 centimeters produced waves breaking over the shoal 79
5.14 Position of the groin in relation to the grid rows and columns 80
5.15 Results for the numerical model of Liu and Mei (1975). Deepwater
wave height equal to 1 meter, deepwater wave angle equal to 45 degrees.
Source: Liu and Mei (1975) 83
5.16 Results from the present numerical model using the same conditions as
were used by Liu and Mei. Offshore boundary condition of constant wa¬
ter surface level which precludes tangential flow at the offshore boundary. 84
5.17 Results from the present numerical model using the same conditions as
were used by Liu and Mei. Offshore boundary condition of no on-offshore
flow. Grid spacing equal to 20 meters 85
5.18 Experimental set-up for the tests conducted by Hales (1980). Source:
Hales (1980) 86
5.19 Experimental measurements of the mid-depth averaged velocities obtain
by Hales. Measurements at one foot intervals. Source: Hales (1980). . . 87
5.20 Rip-current velocities adjacent to the down wave side of the jetty. Com¬
parison of experimental and numerical results 88
5.21 Velocity vectors of the currents obtained numerically for the same con¬
ditions as for the Hales experiment. Grid spacing equal to one foot. . . 89
5.22 Definition sketch of grids for the case of a shore parallel breakwater. . . 90
vu
5.23 Wave crests in the lee of a semi-infinite breakwater in a region of constant
depth 91
5.24 Velocity vectors obtained from the model for a 100 meter long breakwater
100 meters offshore for a 1 meter, 6 second wave for increasing angles of
incidence 93
5.25 Stream function contour lines obtained from the Liu model for a 10
second 1 meter wave normally approaching a 700 meter long breakwater
350 meters offshore on a 1 on 50 slope. Source: Liu and Mei (1975). . . 94
5.26 Flow lines obtained from the present model for a 10 second 1 meter wave
normally approaching a 700 meter long breakwater 350 meters offshore
on a 1 on 50 slope 95
5.27 Velocity vectors obtained from the present model for a 10 second 1 meter
wave normally approaching a 700 meter long breakwater 350 meters
offshore on a 1 on 50 slope 95
5.28 Set-up contour lines obtained from the present model for a 10 second
1 meter wave normally approaching a 700 meter long breakwater 350
meters offshore on a 1 on 50 slope 96
5.29 Physical layout for the experiment of Gourlay. Source: Gourlay (1974),
Proceedings 14*^ Coastal Engineering Conference, copyright by the Amer¬
ican Society of Civil Engineers, reprinted with permission 97
5.30 Computational domain for the experiment of Gourlay. Grid spacing
equal to .1 meter 98
5.31 Comparison of experimental and numerical results for the experiment of
Gourlay: Wave Set-up Contours. Top figure is from Gourlay (1974) Pro¬
ceedings 14^ Coastal Engineering Conference, copyright by the Amer¬
ican Society of Civil Engineers, reprinted with permission. The bottom
figure is the computer results 100
5.32 Results from the experiment of Gourlay: Contours of current veloci¬
ties and streamlines. Source: Gourlay (1974), Proceedings 14^ Coastal
Engineering Conference, copyright by the American Society of Civil En¬
gineers, reprinted with permission 101
5.33 Results from the numerical model: Streamlines of the flow 101
5.34 Results from the numerical model: Contour lines of the current veloci¬
ties 102
viii
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
NUMERICAL MODELING OF WAVE-INDUCED CURRENTS USING A PARABOLIC
WAVE EQUATION
By
HARLEY STANFORD WINER
August 1988
Chairman: Dr. Hsiang Wang
Major Department: Aerospace Engineering, Mechanics and Engineering Science
Prediction of waves and currents in the nearshore region in the presence of proposed
breakwaters is needed in the design process in order to optimize the placement of the
structures. A numerical model is developed which solves for the wave field, the depth
averaged mean currents and the mean water surface elevation within a given computa¬
tional domain with the presence of breakwaters. The model iterates between a wave model
and a circulation model. The wave model uses a parabolic approximation to a combined
refraction-diffraction mild slope equation which includes currents. The circulation model
uses an alternating direction implicit method solution to the depth averaged equations of
momentum using the radiation stress terms obtained through solution of the wave model as
the forcing terms. Each of the terms in the depth averaged equations of momentum involve
assumptions that are discussed in the report. The model is constructed in modular form so
that improvements can easily be incorporated into the model.
The model was employed for several test cases such as rip channels, bar gaps, and
both shore perpendicular and shore parallel breakwaters. Good agreement was obtained
when the model was used to simulate a laboratory experiment of Michael R. Gourlay
(“Wave set-up and wave generated currents in the lee of a breakwater or headland,”
Proc. 14*^ Coastal Eng. Conf., Copenhagen, 1974, pp. 1976-1995).
IX
Suggestions are made for improving the model and for including sediment transport in
X
the model.
CHAPTER 1
INTRODUCTION
As waves propagate they are transformed by variations in currents and depths, along
with interaction with obstructions. The major processes of wave transformation include
shoaling, reflection, refraction, diffraction and dissipation. Shoaling is the change of wave
energy, propagation speed and wave height due to changes in the water depth and/or current
field. Refraction is the change in the direction of wave propagation due to the spatial
gradient of the bottom contours and/or current field. Diffraction is the phenomenon of
wave energy being transmitted laterally along a wave crest. Dissipation is the process of
energy loss due to a multitude of possible causes such as bottom percolation, damping
due to viscous effects, turbulence, and wave breaking. As waves approach a shoreline and
break, momentum exchanges can produce currents and mean water level variations. These
currents can be in the form of longshore currents, rip currents, or circulation cells in the lee
of obstructions. They in turn modulate the incident wave field. The process can be quite
complicated depending upon the complexity of the bathymetry, the ambient current field,
and the presence of obstructions.
Coastal engineers and planners have many reasons for wanting to obtain accurate pre¬
dictions of waves and currents and water elevations in the coastal rone. For instance, pre¬
dictions of extreme conditions are needed for site planning and erosion mitigation studies.
Harbor authorities need wave predictions for structure design and for operational decisions.
This dissertation develops a numerical model for obtaining predictions of waves, currents
and mean surface elevations in a given domain with the offshore wave conditions and the
bathymetry as the input. The model uses a linear wave equation and depth integrated
equations of momentum and mass conservation. Comparison with experimental results
1
2
indicates that the model produces useful and valid predictions of the waves and currents in
the nearshore region in the presence of obstructions such as groins and breakwaters. The
predictions obtained through use of the computer model herein reported are not intended
to be the sole resource for the engineer. A good engineer or planner will seek answers from
multiple sources and cross check them, recognizing the limitations of each method used to
obtain the answer.
Predictions of waves and currents can be obtained for different cases analytically,
through experiments or by numerical solutions. Analytical solutions are obtained only
for simple cases involving idealized beaches using simplified assumptions which are not al¬
ways satisfied in nature. Longuet-Higgins and Stewart (1963) give an analytic solution for
wave set-up for steady state conditions on a plane beach with straight and parallel con¬
tours, assuming linear waves approaching the beach normally, and spilling breakers where
the wave height shoreward of the breaker line is a constant times the total depth of the
water column. Longuet-Higgins (1970a) gives a cross shore profile of the longshore current
for a plane beach with straight and parallel contours, for linear waves approaching the beach
at an angle, spilling breakers, and linearized bottom friction, while neglecting the influence
wave set-up has upon the total depth of the water column. In a companion paper, Longuet-
Higgins (1970b) extended the solution to include lateral diffusion of the current momentum.
There are other analytical solutions for waves and currents but all are limited to very simple
cases.
The literature is replete with laboratory investigations of waves and currents. Physical
models have advantages and disadvantages. The advantages of experimentation are that
results can be obtained for complex situations where there is uncertainty about the exact
form of the governing equations or about the precise form of any of the component terms
of the equation. Any reasonable amount of detail can be included. Visualization of the
physical processes provides an opportunity to identify problems which otherwise might not
be evident.
3
However there are limitations to the use of laboratory experiments. There is the problem
of scale effects which means that all scaling requirements cannot be satisfied simultaneously;
thus, all features may not be interpreted properly to the prototype. For example, the
time scale dictated by Froude modelling laws in which gravitational forces dominate is
different from that obtained from Reynolds modelling laws in which viscous forces are most
important. Another disadvantage of physical models is the cost of labor to construct and
run large scale experiments. Also, because of the cost of construction, it is not always very
easy to change details in the model. Thus, it is not always possible to isolate one parameter
for study. For example, an investigation of the influence of the size of the gap between
segmented breakwaters might require many reconstructions of the rubble mound structure.
The large size of a model means that it is not easily stored once the initial experiment
is completed. The options are to either dismantle the model upon completion of testing
or to incur some maintenance costs. The significance of this is that the model cannot
always be easily retested if, for example, additional input information becomes available or
man or nature makes a substantial alteration of the prototype after the model has been
dismantled. Also the results of a model for one prototype are not always transferable to
other prototype situations. Often the available space for the physical model will require
the use of a distorted model in which the vertical scale is different from the horizontal
scale. This presents the additional problem that for a distorted model, when modelling
a situation where both refraction and diffraction are equally important, the laws of scale
modelling dictate different time scales for refraction and diffraction. For refraction the
time ratio is Tr — \fLy, where the subscript r is for the ratio of the value in the model
to the value in the prototype and T and Lv are time and vertical length respectively. For
diffraction Tr = ~n-1" in shallow water and Tr = \/Ljjt in deep water, where Ljjr represents
y/Lvr
the horizontal length scale.
Numerical models do not have many of the disadvantages found with physical models.
Any mechanism for which the governing equations are known can be modeled numerically.
4
There are no limits in terms of scale effects. Numerical models, once developed, can be
adapted easily to many different prototype situations. There is no problem in the storage or
duplication of the model. Today computer facilities are far more available than laboratory
facilities. Computer models easily allow for sensitivity testing of the importance of an
individual variable and any specific feature of the model can be monitored.
Numerical models also have their disadvantages. In fact some of their advantages
are also their disadvantages. For example, being able to monitor any feature also means
that the computer will deliver only the information that is explicitly specified. Thus an
investigation of one aspect will not indicate problems with other aspects. Though any
mechanism for which the governing equation is known can be modeled, the model is only as
good as the assumptions involved in the equations. For example, the forces exerted on the
water column by the complex interaction between the fluid motion and the bottom are not
totally understood and are often greatly simplified. Another example is the use of depth
integrated two-dimensional models of fluid flows when the reality is three-dimensional. Both
of the above are incorporated in the present study. Knowing the correct equation does not
always make for successful computer modelling since the complexity of the equations may
make for difficulties (for example solving for the complete three-dimensional flow field using
the Navier-Stokes equations). Computer time and cost for large programs can often be
significant.
1.1 Literature
This section is a brief history and literature review of the attempts of the coastal
engineering profession to obtain predictions of waves and currents in the nearshore region
given the offshore wave conditions.
The extension of Snell’s law which governs optical wave refraction to the analogous
problem of water waves over straight and parallel contours has been known for some time.
Johnson, O’Brien, and Isaacs (1948) reported graphical methods for constructing wave
refraction diagrams. Two methods were reported, one producing wave fronts and the second
5
producing wave rays. Pierson (1951) showed that the method using wave fronts introduced
more error than the wave ray method. These methods proceeded upon the assumption that
the contours are locally straight. The wave rays obtained were important in that the wave
height was determined to vary inversely as the square root of the spacing between the wave
rays. Using these ray construction techniques O’Brien (1950) was able to demonstrate that
the April 1930 destruction of the Long Beach Harbor breakwater on a seemingly calm day
most probably resulted from the focusing of long wave energy from a South Pacific storm.
Munk and Arthur (1952) showed a method to obtain the wave height directly from the
curvature of the ray and the rate of change of the depth contours. The wave ray tracing
techniques did not include the effects of currents and had the disadvantage (from a numerical
modeling point of view) that wave heights and directions were obtained at positions along
a ray and not at pre-specified grid positions. Noda et al. (1974) solved this problem by
devising a numerical scheme that solved a wave energy equation and a statement for the
irrotationality of the wave number to obtain wave height and direction at grid points. This
scheme also included the effects of current refraction but neglected diffraction.
Numerical models of wave induced currents were developed by Noda et al. (1974) in
which time dependent and advective acceleration terms were ignored. Birkemeier and Dai¬
ry mple (1976) wrote a time-marching explicit model which neglected lateral mixing and
advective acceleration. Ebersole and Dalrymple (1979) added lateral mixing and the non¬
linear advective acceleration terms. Vemulakonda (1984) wrote an implicit model which
included lateral mixing and advective acceleration terms. All of the above used the numer¬
ical scheme of Noda et al. (1974) to obtain wave heights and wave angles. Each of these
models neglected diffraction.
Diffraction of water waves represents some very difficult mathematical problems. His¬
torically most solutions were based upon the assumption of a flat bottom, since to quote
Newman (1972),
The restriction to constant depth h is important.... Problems involving wave
propagation in a fluid of variable depth are of great interest, and the resulting
6
diffraction effects are among the most important of all ... but these problems
are less tractable than those involving a fluid of constant depth, (page 2)
Penney and Price (1944) showed that Sommerfeld’s solution of the optical diffraction prob¬
lem as described by Bateman (1944) is also a solution of the constant depth water wave
diffraction problem. Penney and Price (1952) solved the constant depth diffraction problem
for breakwaters assuming a semi-infinitely long, infinitesimally thin barrier. The solution
is in terms of Fresnel integrals which need to be evaluated numerically. The Shore Protec¬
tion Manual (SPM,1977) produced many graphical solutions based upon Penney and Price
(1952) that give the wave height and direction in the lee of breakwaters.
The problem of combined refraction and diffraction according to the SPM could be
handled by assuming a constant depth for a few wavelengths so as to obtain the diffraction
effects and then continue with the refraction effects. Perlin and Dean (1983) used essentially
the same procedure in a numerical model which calculated the longshore thrust in terms of
the breaking wave height and angle.
Liu and Mei (1975) solved for the wave field around shore parallel and shore perpendicu¬
lar breakwaters by using a curvilinear coordinate system where one coordinate corresponded
to the wave ray and the other to a curve of constant phase. Their solution was for “an am¬
plitude modulation factor D(£,f) which accounts for diffraction” (page 72), where £ and f
are the two coordinates. An equation for .D(£,f) was obtained and a boundary value prob¬
lem was established by matching conditions in the incident zone and the shadow zones. For
constant depth the solution for £>(£,f) is in terms of Fresnel integrals. Upon obtaining the
wave field the radiations stresses were obtained and used to force a circulation model that
solved for the stream function and mean water surface level while ignoring the non-linear
advective acceleration terms and the lateral diffusion of current momentum. Although Liu
and Mei produced results for the wave field, currents, and set-up resulting from the diffrac¬
tion of waves by breakwaters, their study is not a wave-current interaction model since
there are no current terms in the wave equations and no feedback of current terms from the
circulation model.
7
Berkhoff (1972) derived an equation governing combined refraction and diffraction with
the assumption of a mild slope. Radder (1979) applied the splitting method of Tappert
(1977) to obtain a parabolic equation governing the propagation of the forward propagat¬
ing component of the wave field. Booij (1981) extended the work of Berkhoff to include
the effects of currents and also obtained a parabolic approximation. Kirby (1984) made
corrections to the equation of Booij.
Parabolic wave equations have recently been used in wave-induced circulation models.
Liu and Mei (1975) solved for the wave field by matching conditions between the shadow
zones and the illuminated zone in the lee of a structure. The wave information was then
used to force a circulation model which formulated the governing depth integrated time
averaged momentum equations in terms of a stream function. However this model did not
have current feedback in the wave equation. Yan (1987) and the present study both include
a current feedback for the parabolic wave equation and solve directly for the mean currents.
Prior to the use of parabolic wave equations to obtain the wave forcing in circulation
models, refraction wave models using a numerical scheme derived by Noda et al. (1974)
were used. Birkemeier and Dairymple (1976) produced an explicit numerical solution of
the depth averaged equations of momentum and continuity using the wave model of Noda
et al. (1974). Their model did not include the nonlinear advective acceleration terms
or any mixing. Ebersole and Dalrymple (1980) included the advective acceleration terms
and lateral mixing in an explicit model which also used the refraction scheme of Noda
et al. Vemulakonda (1984) Produced an implicit model which also was based upon the
wave refraction scheme of Noda et al. Yan (1987) used a parabolic approximation to the
mild slope equation which accounts for combined refraction and diffraction. The present
model Improves upon Yan’s model in numerical efficiency, by introducing structures and
most importantly by obtaining the radiation stress terms directly from the complex wave
amplitude.
8
Figure 1.1: Computational domain of the model.
1.2 General Description of the Model
This section is a description of the numerical model that has been constructed to calcu¬
late the wave and current field within a given domain. This domain usually is bordered by
a shoreline along one boundary and an offshore region along the opposite boundary. The
other boundaries are sufficiently far from the region of interest (i.e. a structure or shoal,
etc.) so as to not affect the region of interest. This will be discussed further in the section
on boundary conditions. The computational domain for the purposes of this investigation
consists of a rectilinear coordinate system. Figure 1.1 shows the computational domain.
The model consists of two parts, a circulation model and a wave model. The circulation
model computes the mean water surface level and the depth averaged mean currents using
the radiation stresses from the waves as the forcing terms and includes bottom friction,
advective acceleration terms and the lateral diffusion of current momentum. The wave
model determines the radiation stress terms by first solving for the complex amplitude as a
9
function of the total depth and the depth averaged mean currents, and then computing the
radiation stresses from the complex amplitudes. The model starts with initial conditions
of zero mean currents and zero mean water surface elevation. The offshore incident wave
height is initially zero and then is smoothly increased up to its given value using a hyperbolic
tangent curve. This is done so as to minimize any transients due to a shock start-up.
The model iterates between the wave model and the circulation model until steady state
convergence is obtained. Each of the two parts of the model will be described separately.
The circulation model uses an alternating direction implicit (ADI) method of solution.
Three equations are solved for three unknowns. The unknowns are fj the mean water surface
level, and U and V, the two components of the depth averaged mean current. The equations
are the continuity equation and the x- and y-direction momentum equations. Each iteration
of the circulation model solves for the three unknowns on the next time level. In other words
each iteration advances the values of fj, U and V by a time increment, At. For the purposes
of this investigation the model solves for steady state conditions using the time marching
character of the circulation model. The model can be used to solve for time dependent
wave fields, but this would require the substitution of a time dependent wave equation in
the wave model. The ADI method first solves, grid row by grid row, for U at the next time
level and an intermediate value of fj, by solving the x direction equation of motion and the
continuity equation using an implicit numerical scheme. This solution is in terms of the
known values of 17,V,and fj at the present time level. Once all the grid rows are solved in
the x direction, the y equation of motion and the continuity equation are solved, again by
an implicit scheme, grid column by grid column, for V and fj at the next time level in terms
of U at the new time level, fj at the intermediate level and V at the present level. For each
of the implicit schemes in the two directions, the governing equations yield a tridiagonal
matrix. The tridiagonal matrix and the method of solution will be discussed later.
The wave model also yields a tridiagonal matrix equation. The governing equation is
a parabolic equation for the complex amplitude A. The complex amplitude contains phase
10
information meaning among other things that y variations of the wave field are contained
within the complex amplitude. The term “parabolic equation” indicates that there are
first order derivatives in the x direction and second order derivatives in the y direction.
The method of solution is to solve for the amplitude on the 1+1 grid row in terms of the
known amplitude values on the I grid row. The solution scheme starts on the first grid row
using the prescribed offshore incident wave and then advances grid row by grid row to the
last grid row. To advance to the next grid row the amplitude values on the 1+1 grid row
are implicitly formulated using a Crank-Nicolson scheme. This will be further discussed in
section 4.2 which describes the finite-differencing scheme for the parabolic wave equation.
Once the complex amplitudes are known the radiation stress terms are computed in terms
of the absolute value of the amplitude and the gradients of the amplitude. This is described
in section 2.2.
The circulation model and the wave model are the heart of the program; however there
is much more. In all there are some 27 subroutines in the program. These subroutines cover
everything from input and output to convergence checks and flooding of the shoreline. The
main program is compartmentalized into many subroutines so as to facilitate any changes
or future upgrading. For example, the bottom shear stress, a term of importance in the
equations of motion, currently uses a crude approximation. If a better approximation is
developed that is also computationally efficient, it will be an easy task to replace the present
subroutine that computes the bottom shear with the newer subroutine. There are many
places where the model can be improved with better assumptions and better approxima¬
tions. These will be discussed in detail in the final chapter outlining recommendations for
further study. The modular construction of the model also allows for ease in substituting
one wave model for another.
Figure 1.2 shows a flow chart of the computer program. The subroutine INPUT is called
to initialize the values of the unknowns and to input the water depth, the offshore wave
height and direction and the position of any structures. The three subroutines WAVE,
11
Figure 1.2: Flow chart of computer program.
12
SHEAR, and LATMDC compute the radiation stresses, the bottom shear stress, and the
lateral mixing coefficients. These terms are the forcing terms in the equations of motion
and are used in the subroutine CIRC, which is the circulation model, to solve for fj, U and V.
There is a decision tree which directs the program to the three subroutines WAVE, SHEAR,
and LATMIX at selected intervals. This is purely arbitrary and is done in the interest of
computer speed. It is assumed that the forcing terms do not change at a great rate and can
thus be updated periodically. This updating is more frequent during the start-up time while
the offshore incident wave height is being increased. The CHECK subroutine is to check the
model for convergence. Since steady state conditions are assumed, this subroutine indicates
convergence when the absolute value of the percentage difference between the updated values
of U,V, and fj and the previous values is less than an arbitrarily small specified value. The
UPDATE subroutine simply updates the values of the unknowns for the next iteration. The
FLOOD subroutine allows for the addition or elimination of grid rows at the shoreline to
allow for beach flooding due to set-up.
Figure 1.3 is a flow chart of the circulation model. The subroutines XCOEF and
YCOEF determine the coefficients of the tridiagonal matrix equation for the general cases
in the x and y directions. The subroutine TRIDA solves the tridiagonal matrix equation
using the double-elimination scheme described by Carnahan, Luther, and Wilkes (1969).
Subroutines XCOJ1 and XCOJN are adaptations of XCOEF to determine the coefficients
of the tridiagonal matrix for the special case along the lateral boundary where J=1 or
J=N. Likewise YCOIM is an adaptation of YCOEF for the I=M grid row which borders
the shoreline. These adaptations are necessary because boundary conditions necessitate
different numerical formulation of some of the derivative terms. Similarly with the presence
of an emerged impermeable barrier bordering the shore side of the JJJ grid row, special
subroutines YCJJJ and YCJJJ1 are used on the grid rows I=JJJ and I=JJJ+1 so as to
properly include the boundary conditions for these grid rows, and for shore perpendicular
groins subroutines XCJGR and XCJGR1 incorporate the boundary conditions for a groin
13
Figure 1.3: Flow chart of circulation model.
14
Figure 1.4: Flow chart of wave model.
located between the J = JGR and the J = JGR+1 grid columns.
Figure 1.4 is a flow chart of the wave model. Subroutine WAVNUM determines the wave
number k, the group velocity Ct and the intrinsic frequency a at each of the grid centers.
These terms are used in the governing parabolic equation. Subroutine AMPCO determines
the coefficients of the tridiagonal matrix and CTRIDA solves the matrix equation for the
complex amplitude on the next grid row. CTRIDA is essentially the same as TRIDA except
that it has been modified to handle complex numbers. The subroutine WW is called by the
AMPCO subroutine in order to determine the dissipation coefficient following the work of
Dally (1980,1987).
Chapter 2 describes the circulation model in greater detail. Individual sections will
cover the governing equations and the formulation of each of the component terms. Chap¬
ter 3 describes the wave model, the derivation of the linear wave equation, the parabolic
15
approximation, the manner of including wave energy dissipation due to wave breaking and
the method of obtaining the radiation stress terms directly from the complex amplitude.
The fourth chapter gives details of the numerical methods used to obtain solution of the
governing equations and the boundary conditions employed. The fifth chapter presents the
results of applications of the model and comparisons with other studies both numerical and
experimental, as well as some analytical solutions to some idealized cases. The final chapter
summarizes the study and makes recommendations for future research.
CHAPTER 2
CIRCULATION MODEL
2.1 Governing Equations
The governing equations for the circulation model are the depth integrated time aver¬
aged horizontal equations of momentum
dU TTdU „dU dfj 1 1
aT + ual + + 9al+ jDr“ ~ jDr-
+
J_ (dstx asxv\ 1 ¿>r,
pD \ dx dy ) p dy
dV TrdV T rdV dfj 1 1
aT + 0 ai + v a7 + 3ai + Jd'1" ~ Jdt-
dx
+
_i_ f^xv 1 = Q
pD \ dx dy ) p dx
(2.1)
(2.2)
and the continuity equation
w + é'£'D)+¿ _ ■•,A(J»V) c°sh k[h0 + z) iut
9 a cosh kh
a
(2.18)
21
where A is the complex amplitude, u/ is the absolute frequency, a is the intrinsic or relative
frequency following the current, and h0t h and k are as previously explained. Equations
governing $ and a parabolic approximation solution for A are presented in the next chapter.
The total velocity vector tT< is expressed as
ut =
= £,+s(ff)]r+[v+*(!;)]
where SR indicates the real part, so that its magnitude is given by
«i - f2+v!+(£)+(£)+[«(£)]’+[* (£)F
The x and y components of the time averaged bottom friction are thus
“ 'rj£[a + *{i)]-''51*
Tty - >4 Jo \v+*(%)}■ w
(2.19)
(2.20)
(2.21)
(2.22)
As in Ebersole and Dalrymple, the assumption is made that only the gradients of A are
retained and the gradients of the depth and spatially varying quantities such as k and a
are ignored. This is essentially the assumption of a flat bottom and very slowly varying
currents.
Figure 2.1 plots the magnitude of the longshore current as a function of the distance
offshore for the various friction formulations discussed above. These velocity profiles were
obtained numerically for an eleven second wave with a deep water wave height of 1.028
meters (1 meter at offshore grid of the computational domain, using a AX of 10 meters
and 30 grid rows) and a deep water wave angle of 17.28 degrees (10 degrees at the model
boundary) approaching a 1:25 plane beach. A friction factor, F, of .01 was used for each case.
Using a friction formulation that ignores the contribution of the waves results in a very weak
friction and thus a large longshore current. Using the Longuet- Higgins formulation gives a
friction slightly less than the more proper time averaged formulation, while a combination of
the Longuet-Higgins formulation and the strong current model (no wave contribution) yields
a slightly stronger friction. For practical numerical considerations the latter formulation is
VELOCITY (M/SEC)
22
in
Figure 2.1: Longshore current profiles for different friction formulations for an 11 second,
1.028 meter wave at 17.28 degrees to a 1:25 plane beach.
most often used since it requires less computation time and more importantly the Gaussian
integration involving the reed part of the gradient of complex amplitudes at times produces
stability problems.
The lateral mixing term in Eqs. (2.1) and (2.2) also involve several simplifying assump¬
tions. The first assumption is found in the derivation of the two depth integrated equations
of motion where it is assumed that the turbulent shear stress term -pu'v1 is independent
of depth. The other assumption is that the process of momentum transfer due to turbulent
fluctuations can be represented by a product of a mixing length parameter and derivatives
of the mean current. This is what is done in the present study and also in previous studies.
Ebersole and Dalrymple (1979), Vemulakonda (1984), and Yan (1987) all used the formu¬
lation given by Longuet-Higgins (1970b) in which he developed an analytical expression for
the cross shore distribution of the longshore current velocities using a linearized friction
formulation as described above, linear waves at an angle to the normal of a plane beach,
23
and a spilling breaker assumption. Longuet-Higgins wrote
/ dU dV\ ..
" = a¿- + t*al) (2'23)
in which ex and ev Eire mixing length coefficients.
Longuet-Higgins reasoned that the mixing length parameter ex must tend to zero as
the shoreline was approached. He assumed that ex is proportional to the distance from the
shoreline, X, multiplied by a velocity scale which he chose as the shallow water celerity
y/gh. Thus ex is written as
ex = NX \/gh (2-24)
where N is a dimensionless constant for which Longuet-Higgins chose the limits as
0 < N < 0.016
(2.25)
Obviously, for physical reasons there needs to be some limit on the scale of these eddies.
Ebersole and Dalrymple arbitrarily chose the value of ex near the breaker line to be the
maximum value and for ex to be constant offshore from the breaker line. The same limit
upon ex is used in this study. The coefficient ev is assumed to be a constant and it is
arbitrarily assigned the maximum value of ex.
In the present study the further assumption is made in computing that the cross
derivatives are negligible (a common assumption in fluid mechanics). Therefore, the lateral
mixing terms become
(*»f)
(2.26)
/ sv\
\£xdx)
(2.27)
in the x direction and
in the y direction. This means that the lateral mixing term is represented by a purely lateral
diffusion of momentum. That the lateral mixing term is a lateral diffusion term has some
significance in that the diffusion term is expressed explicitly and thus there results a time
step limitation to an implicit model. This will be discussed in section 4.1 dealing with the
finite differencing of the equations governing the circulation model.
24
2.2 Radiation Stresses
In a series of papers in the early 1960s Longuet-Higgins and Stewart (1962, 1963, 1964)
introduced and developed the mathematics of the radiation stresses. These stresses which
can be conceptually thought of as the “flux of excess momentum” can be obtained through
integration of the momentum equations over the water column as is shown in Appendix A.
Using linear wave theory the radiation stress components are found to be
Szz = E n(cos2 6 + 1) — ^
At
Svv = E n(sin2 0 + 1) — ^
At m
Szv = En cos 0 sin 6
(2.28)
(2.29)
(2.30)
where E is the wave energy equal to |pgH2, n is the ratio of the group velocity to the wave
celerity and 6 is the angle the wave rays make with the x axis.
If the wave height and wave angle are known at each grid point, the radiation stresses
and their gradients are readily determined for use in the governing equations. To determine
the wave height and angle at each grid point Birkemeier and Dalrymple (1976), Ebersole
and Dalrymple (1979) and Vemulakonda (1984) each used the following scheme which is
based upon the irrotationality of the wave number and an equation for the conservation of
wave energy. The irrotationality of the wave number determines changes in the wave angle
due to refraction while the energy equation determines changes in the wave height due to
shoaling and refraction. This scheme was first introduced by Noda et al. (1974).
The irrotationality of the wave number
V x í = 0 (2.31)
where
k = \k\ cos Oi + |fc| sin Oj
(2.32)
is expanded in Cartesian coordinates to obtain
ae dO
cos 6-—I- sin 0— + sin 0
ox ay
|fc| dx
— cos 9
J_d\k\
1*1 dy
(2.33)
25
For waves on a current the dispersion relation is
oj = a + k
Ü = tanh \k\h + |fc| cos 9 + \k\sin 9
(2.34)
Equation (2.34) is differentiated to eliminate and from Eq. (2.33). Using a
forward difference derivative in x and a backward difference derivative in y, 9{j is obtained
as a function of cosí, sin 0, U,V,k,h and their derivatives as well as 6 at the adjoining grids.
The sin 6 and cos 0 terms are expressed in terms of the surrounding four grids using a 2nd
order Taylor expansion. The solution is an iterative process alternating solution of the 6
equation with solution of the dispersion relationship (2.34) for k = |Jt| which is accomplished
by a Newton-Raphson method until acceptable convergence is obtained for k and 9.
Next an energy equation is solved for the wave height. A conservation of energy equation
is given by Phillips (1969 ed. Eq. 3.6.21) as
_+_{m + CJ> + Sj,._i = 0 (2.35)
Assuming steady state conditions and expanding in Cartesian coordinates gives
(V + C,cos t)1|| + (V + C, sin «)i+ A(U + C, cos S)
dU
dU
+ aZ.(y + cs sin 9) + gzx~q¿ + Tvz-^[ + f¡
dy
dy
'XV
dV _ dV
dx °vv dy
0 (2.36)
where
a
XX
avv =
'TV
— Tyl —
n(cos2 9 + 1) — i
2
n(sin2 9 + 1) — -
2
n cos 9 sin 9
Defining
Qi j = \{^U + C>cose) + ^V + Co8in + 9 }
(2.37)
where
dU dU dV _ dV
zxd¿ + Tvx~dí + Tzv~d¿(Tvv'd¿.
(2.38)
a =
26
and noting that U,V,Cly and a are known quantities; taking forward difference derivatives
in x and backward difference derivatives in y and solving for the wave height, Hij, at the
*,j grid yields
Hij =
(V+C,»in g)g,,y_l _ (U+C, cot e)Hi+1j
AY ~ AX
V+C4»ini U+C, coei r,
ÍY" ¿X
(2.39)
AY Af LJjj
The computation is a row by row iteration proceeding from large i (offshore in the
coordinate system of the above cited authors) to small i. During the computation of each
Hij the breaking height is also calculated and if Hij exceeds this value, then the breaking
wave height is used. Several breaking height indices are possible. Miche (1944) gave as the
theoretical limiting wave steepness the condition
^r- = .142tanh(fcj,h¡,)
Lb
(2.40)
where the b subscript indicates breaking conditions. Le Méhauté and Koh (1967) indicated
from experimental data that the limiting steepness is
— = .12tanh(fc¡,hí,)
Lb
(2.41)
or
J5T* = ^.12tanh(Jfe*fc*) (2.42)
Kb
This is the breaking height criterion used by Noda et al. (1974). Perhaps the simplest
breaking height criterion is the spilling breaker assumption
Hh = .78ht (2.43)
which is used by Vemulakonda (1984) and Vemulakonda et al. (1985).
The scheme described above of determining the wave angle 6 from the irrotationality
of the wave number (2.33) and the dispersion relation (2.34) and the wave height from
the energy equation (2.35) has its limitations. Principally this scheme does not take into
account the effects of diffraction. This limitation is surmounted by use of a parabolic wave
equation. However this presents another problem in that the wave angle is not determined
27
in a parabolic equation solution. Only the complex amplitude, which contains the mag¬
nitude and the phase of the wave, is obtained. Yan’s (1987) solution to this problem is
to obtain 0 directly from the complex amplitude by assuming that the slope of the water
surface resulting from the waves indicates the direction of wave propagation. This value
for 0 and twice the absolute value of the complex amplitude for the wave height are then
used in Eqs. (2.28-2.30) to obtain the radiation stresses. Hie value thus obtained for 0 is
then also used in the dispersion relationship. For plane waves this works fine; however, in
diffraction zones where short-crested wave are encountered problems may arise, since the
instantaneous slope of the wave at a particular location can be other than in the direction
of wave propagation.
In this report the radiation stresses are obtained directly from the complex amplitude
and its gradients using equations developed by Mei (1972) and repeated in his book (1982).
The basis of the derivation is to use the definition of the velocity potential to substitute for
the wave induced velocities in the definition for the radiation stresses.
For an inviscid incompressible fluid the instantaneous equations of motion are
+ = — Vp — pgez (2.44)
V-f = 0 (2.45)
where q represents the velocity components, p is the pressure and ex is the z or vertical unit
vector. Denoting the horizontal components of q as ua,a = 1,2 and the vertical component
of q as w, the boundary conditions are
= 0 (2.46)
= tv (2.47)
on z = rj, where r¡ is the instantaneous free surface, and the bottom boundary condition
dh
u, — = -«, (2.48)
q is assumed to be composed of a time averaged part q and a fluctuating part q'\ i.e.
dr)
dr¡
dt +Up'dxp
f= 9t*,y ,z) + q,(x,y,z,t)
(2.49)
28
Integrating a horizontal component of (2.44) vertically over the depth of the water column
from z = -h to z = r¡, using the boundary conditions (2.46-2.48) and then taking time
averages (the same procedure as in appendix A) the following definition of the radiation
stresses is obtained:
ro
Sap = p u'au'pdz + Safi
J — h
[i
Pdz~ y(»i + *0J
h ¿
+ yW)1
(2.50)
where the overbar is used to indicate the time averaged values and 6ap is the Kronicker
delta function
tap
for a = P
for a ^ f3
(2.51)
Mei (1972) was able to express Eq. (2.50) in terms of the general complex amplitude
function B by relating it to the velocity potential and the free surface height by
_ ig cosh k[h0 + z) iut
a cosh kh
t) = B e~iut
(2.52)
(2.53)
An expression for p is obtained by integrating the vertical component of Eq. (2.44)
dw dw
fai+in‘f3rt+,m
dw
dz
(2.54)
Integrating from any depth z to the surface, using Leibnitz’s rule for interchanging the
order of differentiation and integration, employing the two surface boundary conditions,
Eqs. (2.46) and (2.47) and the continuity Eq. (2.45), Eq. (2.54) yields:
d fi d H
= pgin ~ z) +p—J wdz + p—J upwdz-pw\ (2.55)
Taking time averages gives
rfl du'-w'
p{x,y,*) = P9{r) ~ z) + Pjt dx dz ~ p(w'm)2 (2.56)
Substitution of Eqs. (2.52),(2.53) and (2.56) into (2.50) and making the assumption that
horizontal derivatives of k, h and a can be neglected (this is essentially a mild slope as¬
sumption), the following is obtained after straightforward but tedious calculations which
29
are shown in Appendix B:
c P9 L (dBdB*) 1 f 2kh \
4 | ydxa dxp J k2 V sinh2fch/
°ap
IBI
2 kh
sinh 2 kh
+
dxa dip
(|VhB|2 - *2|£|2)¿(2*hcoth2*h - 1)] } (2.57)
where B* denotes the complex conjugate of B. In obtaining Eq. (2.57) use was made of the
flat bottom Helmholtz equation
Vj;£ + Jfc2fi = 0 (2.58)
The general complex amplitude function B is related to the complex amplitude A which
is solved for in the wave model by
B = Ae^flcolidz^
(2.59)
There are some problems with Eq. (2.57). One may question the assumption of a
flat bottom; however the inclusion of the slope and gradients of the wave number and the
intrinsic frequency leads to intractable mathematics. The use of the localized flat bottom
is also implicit in the derivation of Longuet-Higgins and Stewart to obtain Eqs. (2.28- 2.30)
since in carrying out the integrations involved it is assumed that horizontal derivatives of
the velocity potential are limited to derivatives of the phase function. For plane waves
using a spilling breaker assumption, results using Eq. (2.57) differ little from those using
Eqs. (2.28-2.30); however in diffraction zones Eq. (2.57) is superior for reasons mentioned
previously.
The major problem with the use of Eq. (2.57) is also a major problem associated with
the use of Longuet-Higgins’s formulation. This problem is that in both formulations the
initiation of wave breaking implies the dissipation of energy and thus the forcing of set¬
up and longshore currents starting at the breaker line. The physical reality is somewhat
different in that in the initial stages of breaking, energy is not immediately dissipated but
rather is transformed from organized wave energy to turbulent energy and to a surface roller
(Svendsen 1984). There is a space lag between the initiation of breaking and the dissipation
30
30
23
20
1.3
l 0
03
0
03
6
4
2
>
0
•2
-4-
eat.
Figure 2.2: Experimental results for wave set-up and set-down. Reprinted from Bowen, In¬
man, and Simmons (1968), Journal of Geophysical Research, vol. 73(8) page 2573, copyright
by the American Geophysical Union.
of energy. This can be seen in figure 2.2 which gives the results for the experiments of Bowen
et al. (1968). Analytical and numerical solutions for the mean water surface elevation will
give a set-up starting at the break point while the solid fine in the lower graph shows
that set-up measured in the experiment begins shoreward of the break point. This will be
discussed further in Chapter 5.
CHAPTER 3
WAVE MODEL
3.1 Governing Wave Equation
Severed methods have been used to derive the governing equation for the propagation
of waves over varying topography and varying currents. Among the methods used are per¬
turbation techniques, multiple scales, a Lagrangian and the use of Green’s second identity.
With the assumption of a mild slope and weak current conditions each of these methods
should yield the same equation. For illustrative purposes the method of using Green’s sec¬
ond identity will be used here to derive a linear equation. A Lagrangian derivation of the
same equation is presented in Appendix C. The linear equation is used since second order
waves exhibit singularities as the water depth goes to zero.
A velocity potential t is assumed such that the water particle velocities are given by
V$, where V is the three-dimensional gradient operator
d - d * d -
(31)
and *,j, and k are the unit vectors in the x,y, and z directions, respectively. The velocities
are a superposition of steady currents and wave induced motion. The velocity potential is
given by
t = to + eft (3.2)
where to is defined such that its gradient yields U,V, and W which are the steady current
terms, and ft is the wave component of the velocity potential. An undefined scaling e
is used so as to separate the wave and steady current part of the velocity potential. The
wave part of the velocity potential is composed of two parts: /, which contains the depth
dependency and which is a function of the horizontal coordinates and time. The depth
31
32
dependency / is given by linear theory as
j _ cosh k(h0 + z)
cosh kh
(3.3)
where h0 is the depth of the bottom referenced to the still water line, and h is the total
water depth and k is the wave number vector. It should be noted that / is also a function
of the horizontal coordinates since k and h0 vary with horizontal position.
For incompressible, irrotational flow the governing equation for the velocity potential
is the Laplace equation
VV = 0 h0zz — 0
h0 < z < T]
(3.5)
(3.6)
The governing equation is complemented by the following boundary conditions. On the free
surface the kinematic free surface boundary condition is given by
• VhV = 4>x
z = r¡
The dynamic free surface boundary condition with zero atmospheric pressure
1,
4>t + ^(V¿)2 + gz = 0
z = r¡
(3.7)
(3.8)
and the bottom boundary condition
- Vh • Vfcho = 4>t
z = -h0
(3.9)
The derivation technique is to write Green’s second identity for and / as defined by
Eqs. (3.2) and (3.3) respectively.
rfi rfi
í fzzdz - Í dz = [ft - fz]lt
W — ft0 J ~ ItQ
(3.10)
33
Then Eq. (3.6) is used to substitute for zz in the first term of Eq. (3.10) and the boundary
conditions are used to substitute for z at the mean free surface and at the bottom in the
third term of Green’s identity. Equation (3.2) is used and terms of order e only are retained
in order to obtain the governing equation for .
In order to obtain z at fj it is assumed that
tj = tj + eq.
(3.11)
The kinematic free surface boundary condition, Eq. (3.7), is then expanded in a Taylor
series about z = fj,
_ d
+ e"a¡
dn
- z
= 0
J*=fJ
(3.12)
solving for z using Eqs. 3.2 and 3.11 while retaining only terms of order e yields
dfj
4>* U—
A total horizontal derivative
dt
+ U • Vhfj + Vh • Vfcfj - fjWz
*=tf
(3.13)
D 9 ft yn
Di = ai + u’Vh
(3.14)
is introduced so that the first two terms on the right hand side of Eq. (3.13) can be written
as the total horizontal derivative of fj. The final manipulation of Eq. (3.13) is to recognize
that — ^ at the surface is the horizontal divergence of the horizontal current (= V>, • Ü).
Thus
* |fl= + fjVh • Ü
(3.15)
rj is eliminated from Eq. (3.15) by expanding the dynamic free surface boundary condition,
Eq. (3.8), in a Taylor series about z = fj.
4>t + + gz
*=T)
= 9*1 +
+cñ
M=f¡
¿« + ¿(V*)2
, 1 3
9 + *t+ 2^(V^)
= 0
J*= the order t terms yield a solution for rj
(3.17)
34
Substituting into Eq. (3.15) gives
The first term of the Green’s identity using Eq. (3.6) becomes
r f.,dz=-r fvi+dz=-r w^p^dz (3.19)
* “ ho * ” ho * ho
Now
Vh = U + efVh4> + e4>'Vhf (3.20)
and
^h • = Vfc • Ü + eVfc • (/ V/,^) + • V/,/ + zVhf dz (3.22)
J ~ho * “/lo
These two terms are handled using Leibnitz’s rule of integration.
vaT /(/v„¿)dz = r vhf(fvj)dz+r fvh(fvht)dz
* ho * /lo * hp
+^hVf2 U + |_a„ V/,<£ (3.23)
Using this relationship and recognizing that / |r fid.—*™* (3.26)
J-h0 J-ho 9
The third term in Eq. (3.10) is
H. \-H = Hu In -Hu Uo (3.27)
The first term on the right hand side is equal to t |q since / 1^= 1 and has already been
determined. The second term on the right hand side of Eq. (3.27) is easily determined using
the bottom boundary condition Eq. (3.9)
- Hm |-fc„=/Vfc^*Vh/i0 |_fco (3.28)
Using Eq. (3.20) this becomes
- H. U„= fU • Vkh0 |_fco +e/2 |_fco Vhj> ■ Vhh0 + efivhf ■ Vhh0 |_fco (3.29)
The first term is dropped since it not order e and the third term is clearly higher order in
the slope. Thus, again dropping the t
-Mk=/!|-lVkfVA (3.30)
This term cancels the bottom boundary term resulting from the Leibnitz integration above.
The fourth term of Eq. (3.10) is also easily determined.
“ 4>fu |-A0= ~^81DcMh °kh ^ *"fco= ~k tanh kh ” Y* ^ (331)
Evaluating at the mean water surface and retaining only terms of order e yields
-Mu th=~+ (3-32)
Collecting terms from (3.24),(3.26),(3.18),(3.30) and(3.32) yields the following hnear
governing equation for .
+ (V* • ü)££ - Vfc(CCfVfc* + (a2 - k'CCM = 0 (3.33)
This equation was first derived by Kirby (1984).
36
3.2 A Parabolic Approximation
Equation (3.33) is a hyperbolic wave equation. It can be altered to an elliptic equation
by assuming that any time dependency is solely in the wave phase. Before doing this an
energy dissipation term DD is introduced as follows. It is assumed that Eq. (3.33) holds
for the case of no energy dissipation; then for the case with dissipation the left hand side is
balanced by any energy that is dissipated. The governing equation is then written as
+ (V* • - Vk(CCfV*¿) + (
DD=- WD( (3'35)
where w is an undefined coefficient indicating the strength of the dissipation. Eventually
the w term will be related to the energy dissipation due to wave breaking following the work
of Dally, Dean and Dalrymple (1984).
Using Eq. (3.35) in Eq. (3.34) and making the assumption that the only time dependency
is in the phase, i.e.
dd>
= -iu)* (3.36)
reduces the hyperbolic equation (3.34) to the elliptic form
-2iuU ■ + 0 • Vh(U • Vj) + (Vfc • ff)(Ü ■ Vh¿) - V„ • (CCtVh¿)
+{ = úrw<£ (3.37)
where only the phase contribution to the horizontal derivative of ^ is retained in obtaining
the term on the right hand side of (3.37).
Before deriving a parabolic approximation to Eq. (3.37), the need for such an approxi¬
mation should be discussed. There are two major computational drawbacks to numerically
37
solving an elliptic equation. The first is that solution is required simultaneously for each
grid. This means that for a large computational domain, say m by n grids, it is necessary
to invert anmxnbymxn matrix. If m and n are sufficiently large, problems may be
encountered in terms of the interned memory storage, and even if the available storage is
not a problem, the number of computations required to invert the matrix may make for
an exceedingly slow solution. In some situations this may be acceptable; however for the
purposes of a wave induced circulation model in which repeated solution for the wave field
is necessary the program would take a very long time to run. The other disadvantage in¬
volved in the numerical solution of an elliptic equation is that boundary conditions must be
specified at all of the boundaries. In many applications the only known boundary condition
is the offshore, or incident wave amplitude. For example, the domain of interest may be the
wave field in the lee of an obstruction where the only down wave condition is that the waves
freely pass through the down wave boundary. Or in the present study in which the wave
amplitude is known to go to zero at the shoreline a problem ensues in that the shoreline is
(for purposes of the circulation model) at the grid edge and not at the grid point.
A parabolic equation eliminates both of these problems. Since the solution proceeds grid
row by grid row where each solution uses the results of the previous grid row, no down wave
information is required. Since only one grid row is solved at a time, the solution requires
only that a tridiagonal matrix equation be solved to obtain values for the grid row. The only
required boundary conditions are the conditions on the first grid row (usually the offshore
boundary) and lateral boundary conditions which are usually a derivative condition that
allows for the wave to pass through the boundary without any reflection at the boundary.
There are several ways to transform the elliptic Eq. (3.37) to a parabolic equation. The
most direct is to make certain assumptions concerning the length scale of variations in the x
direction. At the heart of the parabolic approximation is the assumption that the direction
of wave propagation is essentially along the x axis. For waves propagating at an angle to
38
the x axis the proper form of 4> is
¿ = -||g y) e*(/ * co" * dz+f * ,in * «*»> (3.38)
and the proper form of the dispersion relationship is
u> = a + kcosOU + fcsin 6 V (3.39)
where ff is the angle of the wave propagation relative to the x axis. The essence of a parabolic
approximation is to assume that the angle that the wave rays make with the x axis is small
so that the ksin0 term can be neglected in Eqs. (3.38) and (3.39) and cos 6 is assumed to
be unity. The manner of doing this varies with different authors. Yan (1987) retains the
complete dispersion relation (Eq. (3.39)) but deletes the y dependency in Eq. (3.38) in the
derivation of his governing wave equation. He solves for 6 by determining the slope of the
instantaneous water surface. For long crested waves this method is satisfactory although
requiring additional computations. However, this method may produce erroneous results in
diffraction zones where the presence of short crested waves will produce a horizontal com¬
ponent of the Blope normal vector that is different from the direction of wave propagation.
Kirby’s (1983) method is to make the assumption that
¿=-t-gA(?lVle<(ftdz) (340)
o
so that the e,fktin6dv Gf ¿he phase function is absorbed into the amplitude function,
A.
Using Eq. (3.40) directly in Eq. (3.37) and further assuming that derivatives of Ax are
small compared to derivatives of the phase function (i.e. that « ikAz ) will yield a
parabolic equation. For illustration the case of linear diffraction in the absence of currents
and in a region of constant depth is examined. Substitution of Eq. (3.2) into the Laplace
equation (3.6) yields the Helmholtz equation
V£¿ + *2¿ = 0
(3.41)
This equation is also obtained directly from Eq. (3.37) for the stated assumptions of no
currents, no dissipation, and constant depth. Substitution of (3.40) into the Helmholtz
equation gives
Axx + 2ikAx + AVV = 0 (3.42)
and employing the above mentioned assumption concerning derivatives of Ax yields the
linear parabolic equation
2 ikAx + AVV = 0 (3.43)
Another method of obtaining a parabolic equation is to use a splitting method in which
it is assumed that 4> is composed of forward moving and backward moving components that
are uncoupled.
4> = + ¡j>~ (3.44)
Uncoupled equations for + and are developed and since the equations are uncoupled
only the equation for + is retained. In the following the work of Kirby (1983) and Booij
(1981) is closely followed.
The second order elliptic equation
d_ /1££\
dx Vt dx )
+ = 0
(3.45)
can be split exactly into equations governing forward and backward disturbances $+ an $
which are given by
d$+
dx
d$~
dx
+¿7$+
—»7**~
Booij chose the relation
(3.46)
(3.47)
(3.48)
40
where £ is an as yet unknown operator on . Substituting into Eq. (3.45) and neglecting
derivatives of ( alone leads to the result
4>zz +
<¡>z + = o
(3.50)
The elliptic equation (3.37) is now put into the form of Eq. (3.50). Letting p = CCg
for notational convenience and expanding some of the terms of Eq. (3.37) in their spatial
coordinates gives
[p4>z)z + (p4>y)v + k2p + (u>2 - a2 + iw(V/, • U))4> + iaw4>
~{U2k)z ~ {V2iv)v - (UV]>z)v - (UV4>v)z + 2 tuff -Vj = 0 (3.51)
Booij chose to neglect terms involving the square of the mean current components assuming
that they are smaller than terms involving p and then arranged the above equation as
4>zz + P l[Pz + 2iuU)4>x + k2 ^1 + 4> — Q (3.52)
where
M4> = (w2 - a2 + tw(Vfc • Ü) + «0w)<¿ + 2iuV4>y + [p4>v)v (3.53)
Kirby points out that this choice by Booij leads to his inability to derive the correct wave
action equation. Kirby instead does not drop any of the squares of the mean current
components until after the splitting is complete. By assuming that the waves are oriented
in the x direction the dispersion relation is written as
a = oj — kU
(3.54)
leading to the result
(jj2 — (j2 — 2ukU - (kU)2
(3.55)
Using this relationship Kirby then arranged Eq. (3.51) as
h, + (p - V')-'(p - U*)J, + (1 + ‘^ ) i = 0
(3.56)
41
where
Ml = {2ukU + iüj{Vh-Ü) + i = (3-60)
Í - *»<,-*)»(!+ 35^)* (3.61)
Booij points out that the splitting matrix approach employed by Radder (1979) is equivalent
to the choice
1R = *(1 + 2A:2(p - Z72))
£r = ki(p-U2)2
Using the results (3.60) and (3.61) in Eq. (3.46) along with the definition (3.49) gives the
parabolic equation
^ (1 +
M
k2(p- U2)
By expanding the pseudo-operators in the form
M
) *i+ } = )) ** (3'65)
42
where
•{\
Radder
i Bóoij
(3.66)
and where P2 = Pi + j for both approximations. Kirby chose the approximation with
Pi = 0 since for his derivation which was for weakly nonlinear waves it should correspond
to the information contained in the nonlinear Schrodinger equation. With Pi = 0 Eq. (3.65)
becomes
1
k(p - U2)fc + - { fc(p - u2)}x 4>+ = ik\p - U2)j>+ + '-M4>+
(3.67)
where M+ is defined by Eq. (3.57). An equation for the complex amplitude A is then
obtained by making the substitution
(3.68)
which after dropping squares of the components of the mean current results in
(C, + V)A, + | >v \ (7) t A
4[cc<(p)J/iA=0 (3'69)
Kirby then further modified the equation by a shift to a reference phase function using the
relation
A = A'ei(lx-fkdx) (3.70)
where k can be taken as an averaged wave number, resulting in
(C, + U)A'X + i{k - k)(C¡ + U)A' + \ ^
The shift to a reference phase function is intended to incorporate the phase information
in the complex amplitude so that the instantaneous water surface can be easily recovered
in the form
q(x,y) = 3i{AV/£d1}
(3.72)
43
where 9? designates the real part. This is an important consideration especially when ob¬
taining the radiation stresses.
It is important that in obtaining Eq. (3.69) that squares of the products of the current
components be retained until the last step. The reason for this is seen in that terms such
as kUV and kU2 are combined with terms of ojV and uU using the dispersion relationship
to give terms of aU and aV.
Equation (3.71) has been successfully used by Kirby in the form given and with modifi¬
cations. The modifications are the addition of a non-linear term which makes the equation
applicable for weakly non-linear waves. However for present purposes Eq. (3.71) presents a
certain problem which will shortly be corrected. The problem is that for waves propagating
at an angle to the x axis it is not possible to obtain a correct form of the conservation
of wave action. The case of no current on a beach with straight and parallel contours is
given as an illustration. The equation for the conservation of wave action reduces to the
conservation of wave energy flux. In mathematical terms this is
^Cscos0a2J =0 (3.73)
where a is the wave amplitude. For waves approaching the beach in a direction normal to
the beach Eq. (3.71) reduces to
+¿(C,M' = 0 (3.74)
Multiplying Eq. (3.74) by the complex conjugate of Aand add to it the complex conjugate
equation of (3.74) multiplied by A', yields
[C,|>l|2]s = 0 (3.75)
which is Eq. (3.73) for the special case of cos# = 1. However it is not possible to obtain
Eq. (3.73) from Eq. (3.71) for the general case of waves at angle 6 on a beach with straight
and parallel contours.
The importance of this inability to give the correct form of Eq. (3.73) is that if energy
flux outside the breaker line is not conserved then gradients of the Sxv radiation stress
44
will generate currents off-shore of the breaker line in the circulation model. Although
such currents in general will be of a small magnitude, they can cause some numerical
computational problems. The problem though is easily corrected by using as the dispersion
relationship
u) = a + k cos 9 U (3.76)
and the substitution
_*.A(x,y) j(f kcottdx)
(3.77)
Rather than using a splitting method, the procedure will be to go directly from Eq. (3.51)
using Eq. (3.77) and writing
u>2 — a1 = 2uik cos 9U — (k cos 6 U)"1
(3.78)
in place of Eq. (3.55). Dropping the (j)zx term and all terms containing products of the
current components that cannot be reduced as explained above results in the following
parabolic equation.
Cg cos 6 + U'
(Cgcos 9 + U) Ax + -
At
A + V Av +
5(7),
-cos29)A- i CCg(£)„j +^A = 0
Now adding a phase shift
icoiídz—f tcostdz)
so that the free surface can be recovered by
rj(i,y) = AV'(/tco,fdz)
(3.79)
(3.80)
(3.81)
where 6 is the wave angle obtained by Snell’s law in the absence of any obstructions or y
variation of the depth or current. The following governing equation is thus obtained.
CgCosO + U'
(Ct cos 6 + U)A'X -|- t(k cos 0 - k coe 0)(Ct cos 9 + U)A' + —
At
+VK + \ (£) A' - - coS! H)A' - j [cc.(^)
A'
+ ^A' = 0 (3.82)
45
Details of the algebra and assumptions involved in obtaining Eq. (3.79) are given in Ap¬
pendix D. For computational purposes it will be assumed that cos# = cos#.
From Eq. (3.82) the correct form of the conservation of wave action can be obtained.
For the case of no y variation of the wave number Eq. (3.82) reduces to
(CtcoB6 + U)A'x + -
Cg COS # + U
A' + VA[
1+5©/
■i*C,(l-cos!9)A'-Í[cC((^)„] +|A = 0
(3.83)
Multiplying Eq. (3.83) by the complex conjugate of A and adding to it A times the complex
conjugate equation of (3.83) yields
(C, cos# + U)
+
\V\A\2]
a
X
a
= 0
(3.84)
which is a closer approximation to the exact expression which should have a (C, sin 0 + V)
term instead of V in the y derivative term.
3.3 Energy Dissipation Due to Wave Breaking
Kirby (1983) found a method of relating the dissipation coefficient w to the work of
Dally (1980) who proposed that the decay of energy flux in the surf zone was proportional
to the amount of excess energy flux. The excess energy flux is the difference between the
actual flux and a stable energy flux. In mathematical terms this relationship is expressed
as
(ECt)x = -j[ECg-(ECt).] (3.85)
where
EC,
E
C„
K
h
the energy flux
wave energy = ^pgH2
the rate of energy propagation which is —
ok
constant
the water depth
46
(ECt)t = the ‘stable’ energy flux.
H = the local wave height
The stable energy flux is obtained by assuming it to be a function of the height of a wave
propagating over a flat bottom after breaking ceases, and that this height is a constant
times the water depth, (y) = 7. The constants K and 7 were determined by analysis of
experimental data obtained from Horikawa and Kuo (1966).
Kirby related Daily’s work to the dissipation term w by deriving a conservation equation
for the wave action Garrett (1967) and Bretherton and Garrett (1968) showed that for
waves propagating in a moving and variable medium, the quantity that is conserved is the
wave action (7) > which is the wave energy divided by the intrinsic frequency. Kirby derived
a wave action conservation equation from the linear hyperbolic wave Eq. (3.34) using the
definition of W as given by Eq. (3.35):
^ + (Vfc • u)^ - Vk(CC9Vkj) + ( (!)],=»(!)
(3.99)
48
Setting the right hand sides of Eqs. (3.98) and (3.99) equal to each other results in Eq. (3.97).
At first glance it may appear unusual that the same expression for the dissipation
coefficient results for the case of currents as for the case of no currents. This can be
explained in part by the fact that Daily’s work was with waves where breaking was depth
limited. For the case of current limited breaking (for example a stopping current), another
expression is necessary. In this study Eq. (3.97) will be used, since in general the current
field will be moderately weak. Such an assumption of weak currents was invoked in the
derivation of the governing parabolic equation.
CHAPTER 4
FINITE DIFFERENCE SCHEMES AND BOUNDARY CONDITIONS
4.1 Numerical Solution of the Governing Equations
The finite differencing of the governing equations in the circulation model is derived
by a matrix analysis as suggested by Sheng and Butler (1982). A rational analysis of the
complete Eqs. (2.1-2.3) including all the nonlinear terms is not possible. However guidance
can be obtained through an analysis of the linearized equations. For this purpose the linear
equations for the special case of a flat bottom and without any of the forcing terms will be
examined.
fíñ an f)v
(4.1)
(4.2)
(4.3)
Let W be a vector of fj,U and V
dfj au dV n
-7T + D— + D— = 0
dt dx dy
au dfj
ar + sai = 0
d-L + gdA = 0
at ydy
w =
n
u
V
(4.4)
Using Eq. (4.4), the governing equations are written as
Wt + AWz + BWV = 0
where A and B are the following matrices
(4.5)
’ 0
D
0 ‘
' 0
0
D '
A =
9
0
0
and
B =
0
0
0
0
0
0
. 9
0
0
(4.6)
To find an implicit finite difference scheme for the governing equations, Eq. (4.5) is
finite differenced as
jyn+i _ jyn
AT
+ AW?+1 + BWf+1 = 0
(4.7)
49
50
This equation is solved for Wn+1 in terms of Wn
(1 + 2AZ + 2A y)Wn+l = Wn
(4.8)
where
1 AT
= 2AAX{-
and
A„ = -B
1 „AT
2 AY
(4.9)
4Aj;AvWn+i is added to the left hand side of the equation and 2AyWn is added and sub¬
tracted to the right hand side of Eq. (4.8) so that the equation can be factored to get
(1 + 2A.)(1 + 2Xv)Wn+1 = (1 - 2Av)Wn + 2XvWn (4.10)
An intermediate value W * is introduced so that the above equation can be written as
two equations
(l + 2AI)W* = (l-2A„)Wn
(4.11)
(1 + 2Xv)Wn+1 = W* + 2X vWn
(4.12)
Note that by eliminating W* from the two Eqs. (4.11) and (4.12), Eq. (4.8) is recovered
with the addition of an error of 4A2AyWn+l which has been introduced.
Writing out the component equations of Eq. (4.11) gives
AT AT
V* = yn-9^vfjn
(4.13)
(4.14)
(4.15)
(4.16)
While writing out the component equations of Eq. (4.12) gives
Vn+1 + = fj* + D^6yVn
un+1 = u*
yn+1 , c ¿¡n+1 _ y» . c =n
v +9AY6yV ~v +gAY6^
(4.17)
(4.18)
(4.19)
51
Using Eq. (4.15) to eliminate V* in Eq. (4.19) and (4.18) to eliminate U* in Eq. (4.14) the
following are obtained. For the x direction
AT AT
«* + daxSiU = “ days,v’'
AT
Un+1 + g—6xf¡* = Un
and for the y direction
AT AT
r+1 + v'"+1 = «• + D-^e,vn
AT
yn+1 + g^8vfjn+l = Vn
(4.20)
(4.21)
(4.22)
(4.23)
The above equations are numerically stable for any combination of AT, AX, and AY.
However as a circulation model they are incomplete in that they do not contain any of the
forcing terms that drive the circulation. In particular in this model wave forcing represented
by gradients in the radiation stresses and bottom friction are of concern. Bottom friction is
especially important because at steady state conditions gradients in the radiation stresses
must be balanced by either gradients in the water surface elevation or by bottom friction
forces. In the absence of friction forces the model could yield continuously accelerating
currents. The method of adding the forcing terms is somewhat arbitrary in that the rational
derivation of the finite differencing scheme of Sheng and Butler cannot be expanded to
include these essentially non- linear terms. The addition of these terms to the model also
limits the latitude of the ratio of jn a manner that is not analyzable. That both the
friction and the radiation stresses are nonlinear should be recognized. The non-linearity
of the velocity friction term is quite evident since it is represented by the quadratic term
pF\U\U. The radiation stress terms, which represent excess momentum flux induced by
waves, are functions of the wave height squared, based on linear wave theory. In shallow
water the wave height H is a function of the total water depth, h0 + fj so that H2 is a
function of fj2.
All the forcing terms and the nonlinear acceleration terms are added at the present or
time level with the exception of the bottom friction term. The quadratic friction term
52
is finite differenced as pF\Un\Un+1 and the friction term as suggested by Longuet-Higgins is
differenced as pF |uofj| Un+1, where |uorj| is in terms of values at the n time level. However,
when using the integral quadratic expression, it is evaluated entirely in terms of the present
or known values.
The nonlinear advective acceleration terms are added at the nth time level using central
difference derivatives. In the x direction momentum equation these terms are
TT + Ui—lJ ^ AT/" I V i TA i T/- T + l — l
U,J 2AX h ^y
(4.24)
and in the y direction equation they are
jWj + u<-u + fu» + '— (4'25)
In the model of Vemulakonda (1984) which is based upon the WIFM (Waterways Ex¬
periment Station Implicit Flooding Model) (Butler 1978) model, a three time level leap-frog
scheme is used. This gives the advantage that all of the forcing terms (and the nonlinear
acceleration terms) are added at the time centered level. However this has some disad¬
vantages in that additional computer memory and computational time are required and
leap-frog schemes can exhibit instabilities. In the present study a two time level scheme
is used and the wave forcing terms as well as the non-linear acceleration terms are added
at the present time level. Since, for all of the applications in this report, the model is run
until steady state conditions are obtained, it is reasoned that there is no great necessity to
properly time center all of the terms. If the use of the model is extended to time dependent
applications, this may need to be corrected.
4.2 Finite Differencing of the Parabolic Wave Equation
Dividing Eq. (3.82) by (Cycosí + U) yields
f C, cos B+U j
A!z + i(k cos Ó — k cos 6) A1 + ^ ^
;A' +
+
2 (C,cosí+ 17) (C,cosí + 17) v
, i kCa(l — cos2 Í) A,
v 2 (Ct cosí + U)
(Y.) LA
2(Cycosi + i7) \a )
2(Cg cosí + U)
CC,
4A
+
w
2(CjCosí + U)
A' = 0
(4.26)
53
The following notation is used to describe the finite differencing of Eq. (4.26): Fj
represents the value of F[i,j) where F is any term and i and j designate the grid row
and grid column. Because of the staggered grid system used for the circulation model as
previously explained, in which the velocity components are specified at the grid edges, the
velocity terms U and V are averaged across the grid so as to produce a grid centered term
for use in the wave model. Thus 17) in the wave model is actually \[Uj + £^’+1) and similarly
for V and other grid designations.
A Crank-Nicholson finite differencing scheme is used for Eq. (4.26). All terms not
involving an x derivative are an average of the value of the term on the * grid row and on
the i + l grid row. Terms having y derivatives are centered at», j on the i row and í +1, j
(X'+1—X'.)
on the « + 1 grid row. The one term with an x derivative is written as —* ^ 1 . Thus all
terms are centered at » + |,j. The finite differencing of Eq. (4.26) is written as
1
+ 2
£lzA
AX
+ '- [(jfc’+1 cos 0,+1 - fc‘+1 cos 9i+1)A$+1 + (£’ cos 0* - k) cos 0‘)A']
1 f [(C7g cos g + C/)/o-]«->-1 ~ [(Cg cos g +- f7)/cx]«-
h 2AX \ \[Ca cos 0 + U)/(r]fl + \[Ca cos 0 + U)/o
Vcr
(Ce cos 6 + U)
r
- A;-ii
2AY
l}(4+‘+4)
f Vo V |a;+i - a;.,]
\[Cacos6 + U)J. 2AY
1
+ 2
2(Cycos0 + U)
[Ca COS0 + U)
<+' r(V/»);tt - (V/a))t\
i
+
2(C¡ cos 0 + 17)
)!(
2AY
w*)Ui - PH-1
2AY
4
t
4
kCg{l -cosJ0)l*+1 A<+1 +
(Cg cos 0 + U)
kCg{ 1 — cos2 0)
■(#),
[Cg cos 0 + U)
•+l
P
«'+!
V]
(«),
+
4+1
V)
[Cg cos 0 + i/);+1 (C, cos 0 + £7)‘
/
1 Í w*+1 A,+1 w*- A\ )
+ 4 { (Cycos0 + £0*+» + [Cg cos0 + U)) J ° ^'27)
54
where p = CCg for notational convenience and the y derivative of p ^ j is written as
‘ (rUi + pj) [(£)’+, ~ (£))] - (Pj + p;~i) [(f)’ - ($)'„
2(AY)2
(4.28)
In order to evaluate w)+1 by the method to be explained in section 3.3 an intermediate
value of A)+1 is obtained by the explicit formulation
A,+1 — A*-
3 -Í- + i(k* cos O' - k'j cos Ol)A'j
AX
1 Í [(Ct cos 0 + V)/<,])+' - \{C¡ cos 0 + U)/a]
1 [ic„
' AX 1 [(C, cos 0 + U)/afAi + [(Cy cos e + )/a]
(
!H
(Cg cos 0 + U)/a
j+i
2AY
+
( - V
(YMUi -
A* *
kCt(l - cos2 0)
\2[Cg cos 0 + U) J .
2AY
1 2
(Cpcosi + i/)
AÍ
n
V)
T +
2 (Cg cos 0 + U);. _r 2(Cg cos 0 + U)) °
When Eq. (4.27) is written out using Eq. (4.28) and w*+1 obtained through use of
(4.29) there results for any particular values of i and j an equation for the three unknown
values A)+\, A)+1, and Alt.\. The equation is put into the form
ai aYA + hi 4+1 + Cj A’ti = dj (4.30)
where the coefficients a¡, bj, and c; are known quantities involving the properties of the
environment (depth, current, celerity, etc.) and the dj involves the environmental properties
and the known values of the amplitude on the » grid row. When Eq. (4.30) is written for
all the j values, j equal to 1 through N, for a specified i value, the ensemble of N equations
takes the form of a tridiagonal matrix equation. For illustrative purposes let N be equal to
four. The resulting matrix equation is
r^i
aj ci
A’+1
‘ d\ ‘
a2 &2 c2
4+1
d2
03 b$ Cj
A*3+1
dz
£
*-
*
i
4+1
. d* .
L 4+1J
(4.31)
55
Equation (4.31) contains four simultaneous equations with six unknowns. Two bound¬
ary conditions should be specified. At each boundary a functional relationship is established
between the grid inside the domain and the grid outside of the domain. For the wave prop¬
agation model there are two types of boundary conditions that are readily available for use
in a computer model. These are a totally reflective boundary condition which would be
used to model a laboratory experiment in a wave basin, and non-reflecting non-interfering
boundary which will allow waves to freely pass through the boundary. For the reflective
boundary, the mathematical statement of the boundary condition is ^ = 0. At the lower
boundary this is finite differenced as = 0 which gives the functional relationship
A\ = A0. (4.32)
At the upper boundary the functional relationship for a reflective boundary is
An+i = An.
(4.33)
For the non-reflective boundary the mathematical statement is = —imA where m is
the longshore component of the wavenumber. Since the lateral boundaries are sufficiently
far from any diffraction effects, only refraction and shoaling effects are present. Thus
using Snell’s Law m = ksin# = kosin0o. The non- reflective boundary condition is finite
differenced at the lower boundary as
—£y~~ =Y[Al + A°] (434)
which gives
Aq =
fj.únAY)
1 +
tmAYA
2 )
Ai
(4.35)
Similarly at the upper boundary the functional relationship used to eliminate Ajv+i is
(l+ imAY)
AN+1 = (! ,mAY) An (4-36)
The use of either pair of boundary conditions, (4.32) and (4.33) or (4.35) and (4.36) to
eliminate Aq and A;v+i will then yield an N by N tridiagonal matrix equation comprising N
56
Figure 4.1: Definition sketch of grid nomenclature and coordinate system.
equations for N unknowns. Returning to the illustrative example where N is equal to four,
the matrix Eq. (4.31) becomes
6i + 7\d\ ci
02 f>2 C%
O3 63
CS
' A\+1 '
‘ di '
A'2+1
¿2
4+1
¿3
[ 4+1J
. d< .
(4.37)
where 7\ and 7n represent the functional relationships as described above on the j= 1 and
j—N boundaries.
Equation (4.37) is solved easily for the unknown AJ+1 using a double sweep method.
In the computer program the subroutine CTRIDA uses an algorithm of Carnahan, Luther,
and Wilkes (1969) to efficiently accomplish the double sweep solution of eq. (4.37).
4.3 Boundary Conditions in the Circulation Model
As shown in figure 4.3 a staggered grid system is used so that velocities are expressed
at the grid edge and fj and the wave properties are defined at the grid center. So that the
wave model will have full grids it is desirable to place the boundaries of the model at the
grid edges on all sides, with the exception of the offshore grid row where the wave conditions
are an input or boundary condition for the wave model. Figure 1.1 can also be referred to.
The offshore boundary condition that is used in the circulation model is that the mean
57
water surface level is not changed by any of the circulation within the model, mathematically
ETA(1,J) = 0 J = 1, N
(4.38)
This condition necessitates that the offshore boundary be sufficiently far from any set-up or
set-down. Other possible offshore conditions are a no-flow condition or a radiation condition.
The no flow condition would conserve water within the computational domain so that the
volume of water involved in set- up would be compensated by an equal volume of set-down
offshore. In an open ocean environment the open ocean can be conceived of as an infinite
reservoir. However in modeling laboratory experiments a no flow condition as the offshore
boundary condition would be the most appropriate. A radiation boundary condition as
proposed by Butler and Sheng (1984) imposes the condition
dr) dr)
—- + C —- = 0
dt dx
(4.39)
This condition allows long waves to freely propagate out of the computational domain
without any reflection. This tends to eliminate any seiche within the computational domain
which is induced by the start-up transients. This condition was applied at one time in the
present model but it did not preserve the offshore zero set-up condition and was changed in
favor of the more rigid zero water surface elevation condition. If the model is adapted for
use with time varying wave trains (an adaptation that the model can easily accommodate)
then the radiation condition will become imperative.
At the shoreline the boundary condition is that there is no flow into or out of the beach.
Since the model allows for flooding of the beach the index of the shoreline changes. The
model allows for differential flooding along the shoreline, for example, due to differences of
wave intensity in the lee of structures. For each grid column in the y direction the last grid
is the IWET(J) grid, where IWET(J) is allowed to change with changes of the mean water
level. The numerical statement of the no flow boundary condition at the shoreline is
t/(IWET(J)+l,J) = 0
J=1,N
(4.40)
58
This boundary condition is used in other models, though not all models allow for differential
flooding of the shoreline. Yan (1987), whose model was used to investigate waves and
currents around a tidal or fluvial discharge, used the no flow condition away from the
discharge and defined £/(IWET(J)+l,J) as equal to the discharge per grid divided by the
water depth at the grids where discharge occurred.
Other conditions are also imposed at the shoreline boundary. These are that the wave
height and the longshore velocity are both zero. These conditions are imposed at the
shoreward grid edge of the IWET(J) grids even though the wave height and longshore
velocities are specified a half grid space form the shoreline. This is done by writing x
derivatives of the radiation stress terms and of the V velocities that need to be positioned
at the center of the IWET(J) grid in terms of three points. The three points are at the
IWET(J)-1 and IWET(J) grids and the zero value at the shoreline at the edge of the
IWET(J) grid.
For the lateral boundary conditions several alternatives are available. Birkemeier and
Dalrymple (1976) and Ebersole and Dairymple (1979) both used a periodic or repeating
beach. This meant that computations were performed between the second and grid,
and values outside of the computational domain are determined according to
Q(I,1) = Q{I,N) (4.41)
Q{I,N + 1) = Q{I,2) for all I (4.42)
for any quantity Q.
At one point this was used in the present model until the problem of what should have
been transient edge waves generated by wave obstructions during the start-up were found
to be trapped in the model. The periodic boundary condition allowed for disturbances
propagating out of the domain at one lateral boundary to enter the domain at the other
boundary. For straight and parallel contours, the periodic beach condition worked well in
the present model but with the presence of breakwaters another condition had to be used.
For the lateral boundary conditions there are two choices, either closed or open bound-
59
aries. Closed boundaries are easily applied by the condition of no flow across the boundary.
Such a choice of boundary conditions would accurately duplicate the physical conditions ex¬
isting in a laboratory basin (or for prototype conditions within a harbor). Open boundaries
are obtained by specifying derivative conditions. In the present study an open boundary
is obtained by specifying that there are no gradients of the depth and the wave proper¬
ties across the lateral boundaries. This necessitates that the lateral boundaries are placed
sufficiently far from any obstructions or any y gradients caused by the obstructions. An
alternative open boundary condition would be that there is no gradient of the flow. This
would, however, put the position of the boundary condition at the center of the bordering
grid rather than at the edge of the grid and could possibly introduce error at the boundary.
4.4 Boundary Conditions in the Wave Model
In the wave model only lateral boundary conditions are needed for computational pur¬
poses. The offshore condition is merely the specified wave amplitude and wave angle. No
condition is needed at the shoreline since the parabolic wave model marches grid row by
grid grow and needs only the lateral conditions.
The offshore condition is that each time the wave model is called the wave amplitude
on the first grid row is specified as
A(1,J) = C(t) ~~A0(J) J=1,N (4.43)
where A0(J) is defined in the input subroutine according to
A0(J) = e‘*osin 0oJaY (444)
where k0 is the wave number along the first grid row and 90 is the angle the wave ray
makes with the x axis. This offshore condition necessitates that the water depth be constant
along the offshore grid row. This is another reason for choosing the offshore condition in
the circulation model to be a constant zero. With the presence of wave obstructions and
current deflectors such as jetties it is possible that with other offshore conditions the mean
water surface elevation could vary along the offshore grid row. The use of one condition in
60
o
Figure 4.2: Start-up function according to Eq. 4.43 with Ci=30 and C2—3.
the wave model which is contrary to the condition in the circulation model can introduce
instabilities in the overall model. The interaction of the two models is more sensitive than
either of the models individually.
The C(t) term in equation 4.43 is a start-up function which smoothly goes from zero
to unity so as to bring the incident wave height from zero to H0 and minimize any shock
loading of the model. Mathematically it is given by
C(t) = .5
t
1 + tanh(— C2)
w
(4.45)
where C\ and C2 are arbitrary constants to adjust the slope and offset of the curve. Fig¬
ure 4.2 shows equation 4.45 for C\ and C2 equal to 30 and 3 respectively.
The lateral boundary conditions are either open or reflective. The reflective condition
is expressed as
(4.46)
while the open condition is expressed as
= ■'•"A (4.47)
where m is the longshore component of the wavenumber which in the absence of diffraction
effects remains constant
m = k0 sin 60 = k sin 9
(4.48)
61
Direction of
Wave Propagation
x
Figure 4.3: Waves passing through a lateral boundary unaffected by the boundary.
according to Snell’s law. Equation 4.46 or 4.47 are used to eliminate the amplitude outside
of the domain from the tridiagonal matrix equation. To illustrate equation 4.47 is finite
differenced at the lower lateral boundary as
A(I,1) - A(I,0) _ . A(I,1) + A(I,0)
AY im 2
Solving for A(I,0)
i _ imAY
A(I,0) = tVa(I,1) (4.50)
1 + im^-
Thus A(I,0) is eliminated in terms of A(I,1). similarly at the upper boundary A(I,N+1) is
eliminated in terms of A(I,N). The same procedure is used with equation 4.46 which yields
the simple relationships
A(I,0) = A(I,1) (4.51)
and
A(I,N+1) = A(I,N)
(4.52)
62
It is important that the boundary conditions 4.46 and 4.47 be used as prescribed
above. The parabolic equation gives N equations containing the N+2 unknowns A(I+1,J)
for J=0,N+1. The two boundary conditions are thus used to reduce the number of un¬
knowns to a number equal to the number of equations. If equation 4.47 is written in terms
of A(I+1,1) and A(I+1,2) at the lower boundary and in terms of A(I+1,N-1) and A(I+1,N)
at the upper boundary and the governing wave equation is written for J=2,N-1, then errors
may be introduced at the boundaries. This is especially true for the case of a reflective
boundary. Figure 4.3 shows waves freely passing through the lateral boundaries without
any noticeable reflection for the case of a plane wave with a period of 10 seconds propagating
over constant depth at an angle of 10 degrees to the x axis.
CHAPTER 5
RESULTS
The numerical model combining the iterative interaction of the circulation model and
the wave model as described in the previous chapters was tested for several cases in which
either analytical or experimental results are available. The model was also run for cases
where previous numerical results are available.
The two cases involving analytical solutions are for set-up induced by normally incident
waves on a plane beach and longshore currents generated by incident plane waves at an angle
on a plane beach. The first section of this chapter covers set-up, comparing the results
with the analytical solutions of Longuet-Higgins and Stewart (1964) and the experimental
results of Bowen et al. (1968) and Hansen and Svendsen (1979). The second section covers
longshore currents and compares results with the analytical solution of Longuet-Higgins
(1970a,b) and experimental results of Visser (1982).
The model was run for the case of a rip-channel on a plane beach. Comparisons with the
numerical results of Birkemeier and Dairymple (1976) and Ebersole and Dalrymple (1979)
are given in the third section. A rip current was also generated for the case of normal waves
breaking on a shore parallel bar with a gap. No comparisons are available for this case.
The model was run for the cases of shore perpendicular and shore parallel breakwaters
represented in the model as infinitesimally thin impermeable barriers. For the shore per¬
pendicular barrier, the results are compared with the experimental results obtained at the
Waterways Experiment Station as reported by Hales (1980). For the shore parallel barrier
the results are compared with the numerical results of Liu and Mei (1975) who did not
include a wave current feedback mechanism.
In the final section the model was used to simulate the laboratory experiments of
63
64
Gourlay (1974) in which the shoreline in the lee of a breakwater was made to curve and
meet the breakwater as occurs with tombolo formation. The results of the model were in
close agreement with Gourlay’s experimental values.
5.1 Wave Set-up
Longuet-Higgins and Stewart (1964) gave a simple but elegant derivation for the wave
induced set-up on a plane beach assuming linear shallow water theory. For the steady state
case of normally incident waves on a plane beach without any y variations and without any
currents the depth integrated equations of motions reduce to
gdi + ±as~=0
(5.1)
dx pD dx
In shallow water for normally incident waves the radiation stress term from linear wave
theory, Eq. (2.28), becomes
S„ = E(n+l) = ^=^Wff2
(5.2)
A spilling breaker assumption is made so that the wave height if is a constant ratio k to
the water depth D = h0 + fj,
H = k (h0 + fj) (5.3)
Substituting equations (5.2) and (5.3) into Eq. (5.1) gives
5rj _ _ 3 dK2(h0 + fj)1 _ 3 2 fdhQ dfj\
dx 16(h0 -I- fj) dx 8* V dx + dx)
Solving for gives the result
dfj _ |/c2 dh0
dx 1 + |/c2 dx
that the slope of the water surface inside the breaker line is a constant times the slope of
the beach.
Previous models are able to produce the results of the analytical expression but to
date none, including the present, have matched the experimental results. In figure 5.1 it
is seen that the slope of the mean water surface produced by the wave set-up is nearly
65
DISTANCE FROM STILL-WATER LINE
IN BEACH, x CMS
3.0
2.5
2.0
1.5
CMS
1.0
IS?
0.5
0
-0.5
Figure 5.1: Comparison of the numerical solution of Vemulakonda et al. (1985) for set-up
with experimental data of Bowen (1968). Source: Vemulakonda et al. (1985)
a constant. For the numerical results the set-up starts at the breaking point, whereas in
the experimental results the water surface elevation is nearly constant for some distance
shoreward of the breaking point and then there is a change in the slope of the water surface
elevation. The numerical results follow the analytical solution. This is an area that needs
further investigation. The reason for the difference between the analytical and experimental
results is that the analytical formulation is predicated upon the radiation stress term being
a function of the wave height, so that as the wave height is reduced after the breaking point
so is the radiation stress term. As pointed out by Svendsen (1984) the radiation stress term
remains nearly constant after breaking for a distance roughly equal to five or six times the
breaking depth. This distance is termed by Svendsen as the outer or transition region in
which there is a rapid change of wave shape. In the inner region the wave shape according
to Svendsen remains nearly constant and is characterized by a surface roller or bore. By
incorporating the surface roller Svendsen is able to obtain close agreement for wave height
and set-up to the experimental results of Hansen and Svendsen (1979). However this close
66
agreement is just for the inner region. For the outer region there is no viable solution other
than that the radiation stress term and the water surface elevation remain nearly constant.
Further there is no analytical separation of the outer and inner regions other than the a
po8teri definition of where there is a sharp change in the slope of the water surface elevation.
This definition is not useful in a model that seeks to determine the water surface elevations.
Svendsen’s work shows promise and should be pursued; it is not however at this point
viable for inclusion in a numerical circulation model. For one it needs to be expanded for the
case of waves at a angle (this need will be discussed further in the next section). Another
limitation is that Svendsen’s model is for monotonically decreasing depths, which is a severe
restriction (for example excluding barred profiles).
Svendsen conceptually addresses the question concerning the spacial difference between
the initiation of the set-up in the analytical and experimental results. The analytical model
assumes that once breaking commences energy dissipation also begins. Svendsen states that
the organized potential energy of the wave prior to breaking is transformed into forward
momentum in the outer region. In keeping with the definition of the radiation stresses as
the flux of excess momentum, i.e. momentum that would not be present in the absence of
the wave, it is easy to conceptualize that the radiation stress term should remain nearly
constant just after breaking. The problem is in assuming that the wave form in the breaking
zone is characterized the same as outside the breaking zone. The author, at present, offers
no solution other than to recognize the existence of the problem.
In the present study results for the set-up on a plane beach with normally incident
waves differ somewhat slightly from the analytical solution as given by Longuet-Higgins and
Stewart. This is because the present model does not use the spilling breaker assumption.
Instead wave energy is dissipated according to a model suggested by Dally (1980) which is
explained in section 3.3. Thus the slope of the mean water surface is not a constant times
the bottom slope. A result from the experiments of Hansen and Svendsen (1979) show this
clearly in figure 5.2. The figure shows a slight curvature in the wave height after breaking
67
mm
Figure 5.2: Comparison of numerical results from the model with experimental results of
Hansen and Svendsen (1979). Wave set-up and wave heights for a 2 second wave on a beach
slope of 1:34 and incident wave height of 7 centimeters.
and in the slope of the set-up. Although the problem of where the set-up begins is still
present, the curvature of the mean water surface is the same for both cases. It is also quite
clear that the diminution of the wave height is not constant as is assumed by the spilling
breaker assumption.
In figure 5.2 the breaking wave height in the experiment is much greater than in the
model. This is because linear wave theory cannot adequately predict wave shoaling. Ven-
gayil and Kirby (1986) were able to duplicate the shoaling recorded in the experiments of
Hansen and Svendsen (1979) by numerically solving coupled equations for the many har¬
monic components. Their results showed that in the shoaling process prior to wave breaking,
wave energy was transmitted to the higher frequency components. Their work however is
not adaptable for use in the current model since it is not known how to dissipate energy
68
Figure 5.3: Non-dimensional longshore velocities vresus the non-dimensional distance off¬
shore from the analytical solution of Longuet-Higgins. Reprinted from Longuet-Higgins
(1970b), Journal of Geophysical Research, vol. 75(33) page 6793, copyright by the Ameri¬
can Geophysical Union.
from the various frequency components once breaking is initiated.
5.2 Longshore Currents
Longuet-Higgins (1970b) gives an analytical solution for the longshore current generated
by obliquely incident waves on a plane beach. The solution involves a spilling breaker
assumption and assumptions concerning the form of the bottom friction and the lateral
mixing. For steady state conditions the y equation of motion reduces to a balance between
the gradient of the Sxy term and the bottom friction with the lateral mixing serving to
distribute the current. Details of the solution can be found in the cited paper.
Figure 5.3 gives the results of the analytical solution of Longuet-Higgins. In the figure
the ordinate V is the non-dimensional velocity v/v0, where t; is the actual velocity and v0
is the velocity at the breaking point obtained for the case of no lateral mixing.
57t. , .isinfl .
v0 =-¿-(fffis)i (5.6)
6 C
The abscissa X is the non-dimensional distance z/zj, where zj is the distance from the
69
Wave Set-Up Line Breaker Line
Figure 5.4: Comparison of experimental results of Visser (1982) and numerical result for
longshore currents due to a 1.1 second wave on a slope of 1:20 for the deep water conditions
H0 — .065m and 60 = 20 degrees.
shore to the breaker line.
Previous models have reproduced the results of the Longuet- Higgins solution by incor¬
porating all of the same assumptions (i.e. spilling breakers and bottom friction). Since the
present model does not contain the same assumptions comparison will not be made with the
analytical results of Longuet-Higgins. Instead comparison is made with the experimental
results of Visser (1982). For both the analytical and the numerical results it should be noted
that the same inadequacy is present as was discussed in the previous section. In both cases
it is implied that wave energy dissipation commences as soon as the wave begins breaking.
Thus both analytical and numerical solutions initiate a longshore current starting at the
breaker line. It is reasonable to expect that the same holds for the y component of the
onshore flux of excess momentum as for the x component. That is, that the Sxv term also
remains nearly constant over a certain unspecified distance as the organized momentum flux
of the non- breaking wave becomes turbulent momentum upon breaking. This is illustrated
in figure 5.4 which shows a comparison of the experimental results of Visser (1982) and the
model results for the same test conditions. The maximum current is nearly the same in
both the model and the experiment but the maximum is further offshore in the model. This
can be attributed to the initiation of a longshore current in the model at the breaker line.
70
5.3 Waves On a Rip-Channel
Both Birkemeier and Dalrymple (1976) and Ebersole and Dalrymple (1979) examined
the case of waves incident upon a shore with the presence of a rip-channel. The bottom
topography was given according to a formula first used by Noda et al. (1974)
h(x,y) = mx jl + Ae-s^)1/S sin10 ^(y - ztan^)| (5.7)
where ft is the angle the rip-channel makes with the x axis which is, in their coordinate
system, directed offshore from the still water line, m is the beach slope, A is an amplitude
of the bottom variation and A is the length of the repeating beach. The origin of the
coordinates is at the intersection of the still water line and the trough as is shown in
figure 5.5 which gives the depth contours for the case where the slope is .025 and ft is equal
to 30 degrees.
Due to the restriction of the present model that there be no y gradients at the lateral
boundaries the above formulation is modified to be
h(x ) = f mx{n-Ae-s(^)1/S8in10f(y-ztan^)} | J(y - xtan£)| < x
V 1 [mi | j(y — x tan ft)\ > x V ’ '
Also shown in figure 5.5 are the contours obtained through use of the modified Eq. 5.8
for ft equal to 30 degrees and a .025 slope. The two formulas differ only in that the repeating
rip channels in the Noda formulation are removed. For application in the present model the
presence of the trough in the farthest offshore grid row is removed so that the depth in the
offshore grid row remains constant in compliance with the requirements of the boundary
conditions as discussed in sections 4.3 and 4.4.
Figure 5.6 shows the results obtained by Ebersole and Dalrymple for the case of a wave
angle of 30 degrees a wave height of 1.0 meters and a period of 4 seconds. Figure 5.7 gives
the results for the present model. Even though the scales of the two figures are different it
is evident that there is agreement between the two. The length scale of the eddy is about
30 meters for both cases. In the caption to the figure of his results, Ebersole indicated that
the lateral mixing coefficient may have been large. The present author does not believe
71
Figure 5.5: Top: Depth contours, in meters, using Eq. 5.7 with the coordinates as shown.
Bottom: Depth contours as used in the present model. Grid spacing is 5 meters.
72
Figure 5.6: Velocity vectors from the results of Ebersole and Dalrymple for a wave angle
of 30 degrees. Scale for velocity vectors is 1 inch = 2.87 m/sec. Source: Ebersole and
Dalrymple (1979).
73
Figure 5.7: Velocity vectors obtained from the present model using the same input condi¬
tions as for Ebersole and Dalrymple.
74
Figure 5.8: Mean water surface elevation contours showing that the water surface is lower
above the trough. Values indicated are in millimeters.
this to be the case. In the present study the lateral mixing coefficient was reduced to as
low a value as 2 meter2/second without substantially changing the magnitude of the eddy
adjacent to the rip trough. However Ebersole did find a significant difference between the
case of no lateral mixing and with lateral mixing. Results for the present model for the
case of no lateral mixing are not available since the model sometimes exhibits an instability
without the stabilizing influence of the lateral mixing terms. Figure 5.8 shows the mean
water surface elevation contours obtained from the model which show that the water mean
surface is lower above the trough.
The agreement between the results of the Ebersole model and the present model for
the case of waves incident upon a rip-channel should not be surprising since diffraction in
this case is small compared to current and depth refraction.
Another test result is the case of normally incident waves breaking on a bar in which
75
Figure 5.9: Depth contours in meters for the case of a bar with a gap.
Figure 5.10: Velocity vectors for the case of waves over a bar with a gap. Dashed lines
indicate depth contours. Grid spacing equal to 10 meters.
76
there is a gap in the bar. Figure 5.9 gives the contour lines for this test case. The topography
was created by adding an inverted parabola on to a píeme beach with a slope of 1:40. The gap
was created by suppressing the bar with a symmetric smooth tangent function. Figure 5.10
shows the velocity vectors resulting from a 6 second wave with an incident wave height of 1
meter. The contours of the bottom topography are represented by dashed lines. The driving
force of the circulation is the return flow through the gap in the bar of the water that is
carried over the bar by the breaking process. These results may not be totally correct given
the arguments of section 5.1. If there were no gap in the bar the waves breaking on the bar
would push water shoreward until sufficient set-up was produced. However with a gap in
the bar lateral flow is allowed and set-up cannot be maintained.
Kirby (conversation with the author) has remarked that he has observed the initiation
of a very strong current in the direction of the waves produced by waves breaking on a sand
mound in a small wave basin. This was tested in the model using normally incident waves
on a slope of .02 with a shoal described by the following equation which is the shoal used
in the Delft experiments of Berkhoff et al. (1982).
Í.45 i < —5.82m
.45 — .02(5.82 + x) (Ój)1 + (Í)1 > > 0 (5.9)
.45 - .02(5.82 + x) - .s/l - (*%)’ - (*)! + .3 (rfg)1 + (I)1 < 1.0
Figure 5.11 shows the depth contours for this test case. Figure 5.12 shows the velocity
vectors obtained for the case of a barely breaking wave. For the case of a non- breaking
wave the velocity vectors plotted to the same scale are not discernable. Figure 5.13 plots
the max velocity obtained in the shoal region versus the incident wave height. At a wave
height of about 10 centimeters there is a sharp change in the slope of the curve indicating
increased velocities once breaking commences.
5.4 Shore-Perpendicular Breakwater
The model is tested here for the case of waves interacting with a shore-perpendicular
breakwater. The breakwater is represented as an infinitesimally thin impermeable barrier
located in between the J = Jbrk and J = Jbrk+i grid columns and extending from
77
Figure 5.11: Depth contours for the case of normally incident waves upon a shoal. Grid
spacing is .05 meters. Contour interval is .02 meters.
78
Figure 5.12: Velocity vectors for the case of waves breaking upon a shoal.
79
Figure 5.13: Maximum velocity across the top of a shoal versus incident wave height for
waves over a shoal on a plane beach. Incident wave height less than 9.1 centimeters are
non-breaking waves. Incident wave heights greater than 9.1 centimeters produced waves
breaking over the shoal.
80
Figure 5.14: Position of the groin in relation to the grid rows and columns.
the shore boundary to the I = Ibrk grid row. Figure 5.4 illustrates the position of the
breakwater. Mathematically, this thin impermeable barrier behaves as an internal solid
boundary that prevents any activity on one side of the barrier from being transmitted to
the other side. This requires modifications to both wave and current models.
In the wave model an internal boundary condition is imposed so that at the breakwa¬
ter face wave energy is totally reflected. This is accomplished by imposing the reflection
condition
f* = o
dy
(5.10)
on both sides of the breakwater from grid row I = I brk to the shoreline, or
A[I,Jbrk) = A[I,Jbrk+\) for all I > Ibrk (5-11)
When the wave equation is expressed in finite difference for J = Jbrk, A[I, Jbrk+i)i which
is on the other side of the barrier, is eliminated by virtue Eq. 5.11. Likewise when the wave
equation is applied to the other side of the barrier A[I, Jbrk) is also eliminated. Further
81
reflective conditions are imposed at the barrier by assuming that the derivatives of any
property perpendicular to the breakwater has a zero value at the breakwater.
In the circulation model the interned boundary condition of impermeability at the bar¬
rier is imposed by setting the velocity component normal to the barrier equal to zero.
V(I,Jbrk+i) = 0 for all I > Ibrk (5.12)
To insure that current-induced momentum is not diffused across the barrier, subroutines
were introduced in the computation for the x sweep for the two grid columns adjacent to
the barrier.
Two sample cases were run for comparison purposes. The first case is a breakwater
extending 400 meters from the shoreline on a beach with a composite slope of 1 on 10 for
the first one hundred meters and a slope of 1 on 100 beyond that point. This is a geometry
selected by Liu and Mei (1975). A deep water wave height of 1 meter and a deep water
wave angle of 45 degrees was used by Liu and Mei. Using linear shoaling and refraction
this gives a wave height of .84 meters and a wave angle of 31.15 degrees as the wave input
conditions at the offshore grid row 800 meters from shore in a depth of 16.95 meters.
The solution method of Liu and Mei was to solve for the wave field in the absence of
currents and mean water level changes and then use the wave field as a constant forcing in the
solution of the circulation. By ignoring the lateral mixing and advective acceleration terms
and assuming steady state conditions, and by expressing the two horizontal components of
the depth integrated flow as derivatives of a stream function ip, the two horizontal equations
of motion were reduced to one Poisson equation
VV = /(*,y) (5.13)
where the / term contains the forcing due to water level gradients, friction and wave forcing.
This equation represents a boundary value problem with two unknowns ip and fj. Upon
specifying boundary values an initial guess for the water surface fj is made and then Eq. 5.13
is solved for ip. With this new solution for ip the value of rj is updated by solving steady
82
state equations of motion. This process is repeated until convergence.
Figures 5.15 and 5.16 show the results for Liu and Mei and the present study respec¬
tively. Figure 5.15 shows the contours of the stream function whereas figure 5.16 shows
flow lines which give the direction of the flow but not the intensity of the flow. The present
model shows a separation eddy at the tip of the breakwater which is absent in Liu’s results.
Both models show a series of rip cells in the up wave side of the jetty which result from
the interaction of the incident and reflected waves. There are not as many rip cells in the
present model as the weaker ones farther upstream are buried by the stronger longshore
current. Neither the Liu model nor the present one show a rip current at the downwave side
of the jetty. In figure 5.16 the flow lines in the offshore region are directed in the onshore
- offshore direction which is different from the stream function contour lines in figure 5.15.
This is because of a difference in the offshore boundary condition in the two models. In
the Liu model the offshore condition is a constant xp which means that the flow must be
tangential to the offshore boundary. In the example used in the present model the offshore
condition is a constant water surface (fj = 0) which forces the tangential velocity to be zero
at the offshore boundary. An attempt to change the offshore condition to a zero normal
flow condition did not accurately conserve water within the computational domain, even
though convergence was achieved. Figure 5.17 shows the flow lines resulting from using no
onshore flow at the offshore boundary.
The second case is to simulate the experiments conducted at the Waterways Experiment
Station and reported by Hales (1980). A jetty perpendicular to the shoreline extended 15
feet out from the still water line on a plane beach with a slope of 1 on 20. The plane beach
extended beyond the jetty tip until a depth of one foot was reached and further offshore
the laboratory basin was a constant depth of one foot. Waves were sent at an angle by
use of a moveable wave paddle and wave guides which were set up to follow the wave rays
at the ends of the wave paddle. Figure 5.18 shows the experimental set-up for these tests.
These wave guides in essence made the experiments to be in a closed basin and greatly
83
Figure 5.15: Results for the numerical model of Liu and Mei (1975). Deepwater wave height
equal to 1 meter, deepwater wave angle equal to 45 degrees. Source: Liu and Mei (1975).
100 200 300 400
_j 1 1 1
Scale (meters)
84
Figure 5.16: Results from the present numerical model using the same conditions as were
used by Liu and Mei. Offshore boundary condition of constant water surface level which
precludes tangential flow at the offshore boundary.
85
Figure 5.17: Results from the present numerical model using the same conditions as were
used by Liu and Mei. Offshore boundary condition of no on-offshore flow. Grid spacing
equal to 20 meters.
86
Figure 5.18: Experimented set-up for the tests conducted by Hales (1980). Source: Hales
(1980).
87
WOVE PERIOD = 1.0 SECONDS
WOVE HEIGHT NEQR GENERATOR = 0.139 PEFT
Figure 5.19: Experimental measurements of the mid-depth averaged velocities obtain by
Hales. Measurements at one foot intervals. Source: Hides (1980).
affected some of the circulation patterns away from the jetty. Figure 5.19 shows the average
mid-depth velocities recorded on the down wave side of the jetty. Note the strong offshore
current next to the boundary. Since the boundaries are less than two jetty lengths from the
jetty, one must ask whether the circulation near the jetty was also affected. The numerical
simulation of the Hales experiments was performed using open lateral boundaries, since
to fully represent the experimental conditions in which the lateral boundaries are neither
straight (i.e. curving wave rays) nor perpendicular to the shoreline boundary would require
mapping techniques that are not yet included in the capability of the model.
Figure 5.20 shows a comparison of the experimental and numerical results obtained for
the magnitude of the rip current adjacent to the down wave side of the jetty for the case of
a 1 second wave with an angle of incidence of 30 degrees and a one foot depth wave height
88
Figure 5.20: Rip-current velocities adjacent to the down wave side of the jetty. Comparison
of experimental and numerical results.
of .139 ft. One would expect the depth averaged velocity to be somewhere between the
mid-depth velocity and the bottom velocity. The model appears to slightly underpredict
the rip current. In figure 5.19 a circulation cell in the lee of the groin is barely discernable.
In figure 5.21 the circulation cell is very obvious.
5.5 Shore Parallel Emerged Breakwater
In order to solve for the wave and current field resulting from the interaction of waves
and an emerged breakwater several boundary conditions at the breakwater are incorporated
into the model in addition to the boundary conditions already discussed in sections 4.3
and 4.4. The breakwater is represented as being infinitesimally thin and impermeable.
Figure 5.22 shows the position of the breakwater on the shoreward side of the I = Ibrk grid
row. The impermeability condition implies that neither waves nor currents pass through
the breakwater. The no current condition is easily imposed in the circulation model by
making Cf(lBjyr+i, J) equal to zero for the entire length of the breakwater. In addition the
circulation model is modified so that momentum diffusion and derivatives do not take place
89
Figure 5.21: Velocity vectors of the currents obtained numerically for the same conditions
as for the Hales experiment. Grid spacing equal to one foot.
90
> X " "
y
y
Figure 5.22: Definition sketch of grids for the case of a shore parallel breakwater,
across the barrier.
In the wave model an adaptation is made so as to prescribe the condition that the
wave height is zero on the lee side of the breakwater. This condition on the lee side of the
breakwater is physically not realistic since common observation reveals that waves actually
propagate along the lee side of a breakwater. The parabolic model however marches waves
in but one direction. This assumption however will give good results several grids away
from the breakwater. Figure 5.23 shows the crests of the diffracted waves in the lee of a
semi-infinite breakwater for a region of constant depth. This figure shows that the zero
wave height condition on the lee side of the breakwater does not adequately represent the
waves along the breakwater or allow wave energy to properly propagate in the lee of the
breakwater; however, important information can be obtained through use of the marching
parabolic wave equation method. To set the wave height to zero on the lee side of the
breakwater the ACOEF subroutine is modified so as to advance the wave one half delta x
from the center of the I = BRK grid row to the position of the breakwater at the edge of
91
Figure 5.23: Wave crests in the lee of a semi-infinite breakwater in a region of constant
depth.
the grid row. The amplitude is then set to zero at the breakwater and remains the same
elsewhere. The wave is then advanced another one half delta x so as to be back at the grid
centers and then is marched grid row by grid row to the shoreline. With the exception of
the two one half delta x steps the wave model is otherwise unchanged.
A test case was arbitrarily chosen of a breakwater 100 meters long 100 meters from
the shoreline on a plane beach with a 1 on 40 slope. The wave conditions again were
arbitrarily chosen to be a six second period and a wave height of 1 meter on the offshore
grid row. Figure 5.24 shows the velocity vectors obtained for the case of normal incidence
and for the cases of increasing angle of incidence. For normal incidence two circulation cells
and a well defined rip current in the lee of the breakwater result. The rip current is fed
by the longshore currents which result from the difference in the set-up in the lee of the
breakwater and the region exposed to direct waves. As the angle of incidence increases a
longshore current develops and the circulation pattern becomes skewed. The rip current
is increasingly skewed in the direction of the longshore current with increasing strength of
92
the longshore current and the upstream circulation cell becomes increasingly weaker. The
other cell remains but is shifted down current.
The emerged breakwater program was also run for the case of a shore parallel breakwater
700 meters in length, 350 meters offshore on a 1 on 50 slope. This is a case tested by Liu
and Mei and represents the condition for the Santa Monica offshore breakwater. The wave
conditions were a 10 second wave of normal incidence with a deep water wave height of 1
meter. Figure 5.25 shows the symmetrical stream function contour lines obtained from the
Liu model for half of the region in the lee of the breakwater. Figure 5.26 shows the flow lines
obtained from the present model and figure 5.27 shows the velocity vectors for this case.
Figure 5.28 shows the contour lines of the mean water surface or set-up. There is obviously
greater set-up in the exposed region than in the lee of the structure but no well defined rip
current is seen in figure 5.27 or in the results of Liu shown in figure 5.25. Both models show
a strong longshore current at the shoreline directly behind the end of the breakwater, but
this current does not extend to the center of the breakwater since there is a distance over
300 meters in which the intense longshore current can be diffused and directed offshore.
5.6 Comparison with Experimental Results of Gourlay
Gourlay (1974) reported results of laboratory experiments in which currents were gen¬
erated by the diffracted waves in the lee of a breakwater. His experiment was designed to
demonstrate that longshore variation of the wave height could generate longshore currents.
The experimental wave basin was thus constructed so that the shoreline was everywhere par¬
allel with the diffracted wave crests. Figure 5.29 shows the laboratory set- up for Gourlay’s
experiments. A 1 on 10 concrete beach is parallel to the incoming wave crests in the ex¬
posed zone and in the shadow zone the slope is also 1 to 10 but the beach is curved with
a constant radius centered on the breakwater tip. Gourlay’s reasoning in choosing this
arrangement was that “since the diffracted wave crests have theoretical curves which can
be approximated by a circle the influence of waves breaking at an angle to the beach was
almost completely eliminated.” (page 1978)
93
5o
Figure 5.24: Velocity vectors obtained from the model for a 100 meter long breakwater 100
meters offshore for a 1 meter, 6 second wave for increasing angles of incidence.
94
Figure 5.25: Stream function contour lines obtained from the Liu model for a 10 second 1
meter wave normally approaching a 700 meter long breakwater 350 meters offshore on a 1
on 50 slope. Source: Liu and Mei (1975).
95
Figure 5.26: Flow lines obtained from the present model for a 10 second 1 meter wave
normally approaching a 700 meter long breakwater 350 meters offshore on a 1 on 50 slope.
Figure 5.27: Velocity vectors obtained from the present model for a 10 second 1 meter wave
normally approaching a 700 meter long breakwater 350 meters offshore on a 1 on 50 slope.
96
Figure 5.28: Set-up contour lines obtained from the present model for a 10 second 1 meter
wave normally approaching a 700 meter long breakwater 350 meters offshore on a 1 on 50
slope.
An experiment of Gourlay was numerically duplicated by making several adaptations
to the overall scheme. Figure 5.30 shows the numerical computational grids used in the pro¬
gram. Rather than model the entire wave basin it was decided for computational economy
to assume the wave paddle location to be situated 3 meters offshore of the breakwater and
that the wave paddle produced waves only in half of the domain seaward of the breakwater.
The offshore boundary condition in the circulation model was a no-flow condition through
the wave paddle. This condition may slightly influence the results obtained for the mean
water levels since the computer model will not have as great an off shore reservoir from
which the wave set up must come as in the physical experiment. An adaptation to the
flooding routine was made to allow for variable flooding. This meant that the last grid
of each column could vary independently of the neighboring grid columns. This change
necessitated the incorporation of many decision statements within the circulation model so
that there would be no flow between wet grids and dry grids unless the water elevation was
high enough to initiate flooding. At present flooding is allowed only in the x or on-off shore
direction.
97
Figure 5.29: Physical layout for the experiment of Gourlay. Source: Gourlay (1974), Pro¬
ceedings 14^ Coastal Engineering Conference, copyright by the American Society of Civil
Engineers, reprinted with permission.
98
Figure 5.30: Computation sil domain for the experiment of Gourlay. Grid spacing equal to
.1 meter.
99
In the wave model part of the program, the problem of wet and dry grids is handled
using a technique proposed by Kirby and Dalrymple (1986) by treating the dry grids as
though they have a very small (say one millimeter) depth of water. This technique allows for
keeping a uniform computational domain for the wave model while avoiding any numerical
problems such as dividing by zero depths.
Gourlay’s experiment used a wave height of 9.1 centimeters with a wave period of 1.5
seconds. The depth in the offshore constant depth region was 20 centimeters. For the
numerical computation a grid size of .1 meters was used. This created a 60 by 63 grid
starting domain which grew to a 64 by 63 grid domain once set up was established.
The numerical results show close agreement with the experiment results. Figure 5.31
shows the experimental results for the mean water surface contours resulting from the wave
set-up and the corresponding results obtained from the numerical model. The agreement
appears to be quite good. There is however an area of slight differences where the shoreline
connects to the breakwater. The computer model fails to show any set-up in contrast to the
experimental results which show a slight set-up in this region. An explanation for this is
that the parabolic approximation used in the wave model does not propagate waves along
the breakwater as occurs in the physical experiment. This in turn may also help explain
the absence of the small stagnation eddy shown in figure 5.32.
Figure 5.32 shows the experimental results for the circulation. The figure shows both
the contours for the magnitude of the velocity and some streamlines indicating the direction
of the flow. Figure 5.33 shows the stream lines obtained from the numerical values while
figure 5.34 shows the contour lines for the magnitude of the velocities obtained numerically.
Again good overall agreement is in evidence. The rip current in the exposed area in par¬
ticular is very well shown in the numerical results. The differences are that, as mentioned
earlier, there is no stagnation eddy shown in figure 5.33 and the center of the primary eddy
is slightly farther behind the breakwater in the physical experiment than in the model re¬
sults. Similarly as seen in figure 5.34 the maximum currents obtained in the model are not
100
1.1 i 1 / /JJ ///>/;?/ 7-7-71
Figure 5.31: Comparison of experimental and numerical results for the experiment of
Gourlay: Wave Set-up Contours. Top figure is from Gourlay (1974) Proceedings 14^
Coastal Engineering Conference, copyright by the American Society of Civil Engineers,
reprinted with permission. The bottom figure is the computer results.
101
Figure 5.32: Results from the experiment of Gourlay: Contours of current velocities and
streamlines. Source: Gourlay (1974), Proceedings 14^ Coastal Engineering Conference,
copyright by the American Society of Civil Engineers, reprinted with permission.
Figure 5.33: Results from the numerical model: Streamlines of the flow.
102
Figure 5.34: Results from the numerical model: Contour lines of the current velocities.
as far behind the breakwater as they were in the experiment. Both of these spacial dispar¬
ities between numerical and physical results may result from the failure of the parabolic
approximation to adequately propagate wave energy fully into the lee of the breakwater.
In comparing the velocities obtained in the model to the experimentally measured values it
should be noted that Gourlay’s results are for surface velocities whereas the model results
are for depth-averaged currents. This may be a factor in explaining that the maximum in
figure 5.32 for the experimental results is larger than the maximum in figure 5.34 for the
numerical results.
The differences in the experimental and numerical results which were mentioned above
should be kept in perspective. Overall the agreement is good both qualitatively and quan¬
titatively.
CHAPTER 6
CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER STUDY
A numerical solution for the circulation and waves in coastal regions has been presented.
The solution method uses an alternating direction implicit method of solving the depth
integrated equations of momentum and continuity. Wave forcing terms are obtained through
a parabolic approximation to the linear mild slope equation for wave propagation over
varying depths and currents. The numerical solution involves an iteration between the
parabolic wave solution and the ADI circulation solution.
The model was run for several test cases. These include cases with varying topogra¬
phy, such as rip channels and bar gaps. The model was also tested with shore parallel and
shore perpendicular breakwaters represented as infinitesimally thin impermeable barriers.
Comparisons were made with analytical and experimental results. These comparisons were
favorable. In particular comparison of the results from the model in simulating the experi¬
ments of Gourlay (1974) show the predictive capability of the model. The model was also
run for an idealized shore parallel breakwater varying the aspect ratio of the length of the
breakwater to the distance offshore of the breakwater to determine the relative importance
of this parameter.
The primary conclusion of this study is that it is possible to get good predictive re¬
sults using a numerical model that incorporates some simplifications of the physical reality.
In the present study the three-dimensional flows that are encountered in nature are rep¬
resented by two-dimensional depth averaged velocities. Simplifying assumptions are also
made concerning bottom friction, momentum transfer or mixing, wave breaking and wave
energy dissipation. For each of these topics it is possible to argue that nature is greatly dif¬
ferent and more complex than the assumptions embodied in this report. Yet the results (in
103
104
particular the comparison to the laboratory experiments of Gourlay) confirm the predictive
capability of the numerical model. The simplifications mentioned above make this possi¬
ble, for there must be a mathematical statement of the physical process before numerical
modelling of the process can even be attempted.
In making recommendations for further study we should distinguish between short range
objectives and the possibilities that may become open to numerical modelers by advances
in the computer hardware field. It is not difficult to envision within a few decades or less,
super computers with thousands of parallel processors. Such advances will make many of the
present computational methods obsolete. For example the parabolic approximation is made
solely to put the wave equation in a solvable form with respect to the capabilities of present
day computers. With sufficient parallel processors it is conceivable that it may be practical
to solve the full elliptic equation explicitly. We shall not dwell further on speculation about
future generations of computers. In the near term, the numerical model presented in this
report can be refined and the applications can be extended. The model is constructed in
modular form so that any new formulation for any of the component processes can easily
be incorporated into the model. The model is thus open for continual improvement.
6.1 Improved Wave Formulations
The model can be improved by better wave modeling and by a more correct formulation
of the radiation stress terms once breaking has commenced. The formulation used in this
model and previous models for the radiation stresses initiates set-up and longshore currents
at the breaker fine whereas as explained in section 5.1 experiments show otherwise.
An improved wave model will assign a non-zero value of the amplitude in the lee of
a shore parallel breakwater. Research should be done to establish a method of doing this
so as to enable more wave energy to propagate into the lee of the structure. Dalrymple
and Kirby (1988) have found a method to determine a non-zero wave height in the lee of a
breakwater. However their solution is in the form of a Fourier integral which is not suited
to the complex amplitude form required in the present model.
105
6.2 Including Wave Reflection
The author has made unsuccessful attempts to include wave reflection in the model.
Kirby (1986) used coupled equations governing the interaction of a forward and a backward
(reflected) propagating wave that are coupled through a slope interaction. He was able to
obtain the reflection and transmission coefficients over submerged breakwaters that agree
well with the experiments of Seelig (1980). Once the arrays of the complex amplitudes of
the two waves is determined, the radiation stress terms can be determined from equations
that have been derived following the same procedure as used by Mei (1972) but defining
as being comprised of two components. In fact this procedure can also be used to obtain
the radiation stress terms for any number of components. The only limit being the amount
of computation required as the number of components increases.
Smith (1987) following the work of Miles (1967) for wave reflection and transmission
over a step derived equations for the complex reflection and transmission from a submerged
plane barrier. An attempt was made using this formulation to obtain the circulation and
wave pattern from a submerged breakwater represented as a plane impermeable barrier.
However problems of numerical stability were encountered. It is suspected that the sharp
gradients between the region with the presence of the reflected wave and regions of its
absence adjacent to the barrier produced unrealistic gradients in the radiation stress terms.
An attempt was made to include the reflected wave for the case of an emerged plane barrier
with the same unstable result. This is a topic for further investigation.
6.3 Extending the Model to Include Sediment Transport
The model as it currently stands can be extended to include sediment transport. The
wave model can be slightly amended to retain wave energy dissipation as an array. If a
formulation of sediment entrainment can be made in terms of the wave energy dissipation
and the bottom shear stress, everything else is in place to accomplish the sediment transport
model. A simple numerical mechanism would be to assume steady state conditions for
certain time intervals and solve for the wave and current field and the sediment concentration
106
as a function of known variables and then solve a separate model for sediment balance and
then update the depths and proceed iterating between the sediment balance model and the
present wave-current model marching in time. Since it has already been demonstrated that
the wave- current model gives adequate results in the lee of a shore parallel breakwater, a
good test for such a sediment transport model wold be to see if it could yield a tombolo
formation in the lee of a shore parallel breakwater.
APPENDIX A
DERIVATION OF THE DEPTH AVERAGED EQUATIONS
In this appendix a derivation of the governing equations for the circulation model is
presented. This appendix closely follows the work of Ebersole and Dalrymple (1979).
Before proceeding with the derivations of the depth integrated time averaged equations
of motion and the continuity equation the Leibnitz Rule for taking the derivative of an
integral is stated since it is invoked throughout the derivation.
Ty C /(I'V) = C dJ^Tix + my)'y)dJ£- - '^y)’y)dJJT
The Leibnitz Rule can also be used in reverse to take derivatives outside of an integral.
The kinematic boundary conditions are also reproduced since the use of the Leibnitz Rule
will yield bottom and surface terms that will be eliminated through use of these boundary
conditions. On the free surface the boundary condition is
dr) dr)
dr)
Wr' = Tt+'x-aiJrv"Ty
(A.2)
and at the bottom the boundary condition
(A.3)
In the circulation model it is assumed that the still water depth does not vary with time so
the first term of Eq. (A.3) is dropped (actually for time varying depth essentially the same
equations are obtained).
The continuity equation (or the conservation of mass equation as it is abo called) in
three dimensions is
, d(/>u) d(pv) d(pw)
dt dx dy dz
(A.4)
107
108
Integrating over the depth from z — —h0 to z = r¡ and using the Leibnitz Rule to take
derivatives outside the integrad, the continuity equation becomes
d f* . dr] d [i
atLpi’-p'Tt+TzL'n,i'
h0
f \ dr> r \ dh0 d n
Si V—■
/ s
(A.5)
The terms with an overbrace are p times the free surface kinematic boundary condition
and the terms with an underbrace are p times the bottom boundary condition. Eliminating
both of these from Eq. (A.5) gives
d n . dr* . d rn J
/ Pdz+*~ pudz + — pvdz- 0
dt J-h0 ox J—k0 dvJ-h0
(A.6)
Equation (A.6) is now time averaged. Time averaging means averaging over a wave
period. The time average of a term is denoted by an overbar. Thus
^ 1 fT
■-H
Fdt
(A.7)
where F represents any term and T is the wave period. Letting u and v be comprised of a
mean flow, a wave induced flow, and an arbitrary or random turbulence component
u = ü + ú + u'
V = V + V + v1
(A.8)
(A.9)
The turbulent fluctuations are assumed to be of a high frequency relative to the wave
frequency so that their time averages over the wave period are zero. By substituting the
above expressions for u and v into Eq. (A.6) and taking time averages the continuity equation
is
^ {p{K + ’?)}+ ^ iPüiho + »?)} + J¿f_k Púdz+ (A.10)
A {««(*.+ fl)} + = o
(AH)
109
The time averages of the vertically integrated wave induced velocities ü and v are not zero.
These quantities are known as the wave induced mass flux.
Defining
U = u + u
and
V = v + 5
(A.12)
where
u =
1 Jn
r / pudz and
o + fj)J-k0
p{h0 + *?)
the time averaged depth integrated continuity equation is
%í+-kuQlo+t,)+‘§¡¡viho+tl)
1 ri
v — ~Ti rr / pi> dz
p{h0 + fi)J-ho
(A.13)
(A.14)
The assumption that the density is constant in both time and space has used in obtaining
both Eq. (A.11) and Eq. (A.14). The assumption of a constant bottom has also been used
in obtaining (A.14). If the bottom is not constant in time a term would be included in
Eq. (A.14).
The x and y momentum equations are handled in the same manner as the continuity
equation above. The x and y direction momentum equations are
du du2 duv duw 1 dp 1 ,drxx drvx drxx.
Tt + + + IT = } Ú + + if + -sT>
dv t duv t dv2 t dvw _ 1 dp 1 fdrxv drvv drZ]/
dt dx dy dz
= zir + 1A-
+
}
(A.15)
(A.16)
p dy p dx dy dz
The same procedure is followed for both the x and y equations. Only Eq. (A.15) will be
dealt with here. Integrating the left hand side over the depth and applying the Leibnitz
Rule to take the derivative outside the integral gives
it Tip\ d /*'J , dr] dh0
(LHS> = ai/-»." ~aT +
if +
dxj-h0 ndx ° dx
dy J-
, , x dr) . . dh0
uvdz - (uv)„ — - (uv).fco — +
h0 dy dy
(UW)l ~ (UW)-ho
(A.17)
110
The terms in Eq. (A.17) that are designated by an overbrace are tif, times the kinematic
free surface condition and those designated by an underbrace are u_j,0 times the bottom
boundary condition. These are both dropped leaving
, „ d r dr, dr
(LHS) = * U UdZ+d~x U UdZ+d~y U UV dZ ^
Substituting the definitions given for u and v in Eqn. (A.8) and (A.9) and recognizing that
the turbulent fluctuating components are of a high frequency relative to the wave frequency
so that the time averages of integrals involving the product of a turbulent component and
a wave induced component are zero gives
d~r d~r d ~r d~r
LHS = — / udz+ — údz + — u?dz+— tfdz (A.19)
dtJ-ko dt]-K0 dxJ-ho dxj-h0 ’
d ~r~T~Z d ~r d ~r d ~r
+ -£- {u')2dz + — / 2UÚ dz + — uv dz + — / uiidz (A.20)
dxJ-ho dxj-h0 dyJ-h0 dyJ-ho
d ~r d ~r d ~r
+ -X- üvdz+— Üi)dz + — u'v'dz (A.21)
oyj-ho oyJ-h0 dy J-K0
using the definitions of Eqs. (A.12) and (A.13) this becomes
+aiuvv,°+«)+TyS'h' \w^)Lrd‘Lriz
(A.22)
where the two terms designated by underbraces result from being added to complete the
U2 and UV terms.
The right hand side after integrating over the depth assuming a constant density and
assuming an inviscid fluid such that no horizontal viscous stresses exist, i.e. rtx and Tyx are
zero, is
RHS =
if r ^
p J-ho dx p J—h0 dz
(A.23)
Using the Leibnitz Rule and averaging over a wave period the right hand side becomes
1 d r . 1 dn 1 dhn 1 1
RHS = —— / pdz+ + -P-Ko-z— + -Tzz~~ -rnir (A.24)
pdxj-h0 p dx p dx p ” p v '
Ill
TZXq and rMX_ko are the time averaged surface and bottom shear stresses. Assuming that
the free surface pressure pn is zero (the dynamic free surface boundary condition) the right
hand side is
___ 1 d 7* 7 1 dh~0 11 ,
RHS = —— / pdz + -p.h—± + -T¿--Ti¿ (A.25)
pdxj-h0 p dx p p ’
The mean pressure at the bottom can be expressed as the sum of the wave induced or
dynamic pressure and the hydrostatic pressure
P-fco = PdVn-ko + P9[h0 + r>)
(A.26)
so that the product P-h0fa can be manipulated to read
dh0 1 dh0 13.,. _.2, . .dfj
'’-*•17 = + 2áí*~~' l(h°+
(A.27)
Terms from both the right and left hand sides are combined as the radiation stresses.
Defining
s„ = f'jp+,07) iz - \P9(h„ + «)* - rM'Y
~r> i ~ñ r>
S,^j ^vdz - ^ ,0 i, j^dz
Sn = /_ J» + ,07) dz - 1,,(*. + «)> - ^7) {/I P5 ^
(A.28)
(A.29)
(A. 30)
and making the assumptions that the product of the dynamic pressure at the bottom and
the bottom slope is negligible, that the Jif?h0(u')2dz term can be ignored and that the
lateral friction tj caused by momentum fluxes due to turbulent fluctuations is independent
of depth and that rj is defined by r¡ = —pu'v', the x direction equation of momentum is
written as
3 »r2,l , ^ , d rnr/L . . ldSxx . 1 dS,
+ --
XV
5^ + « + a-rU '<*• + *» + TyUV^ + + ,- 8x ' , 8„
+,(». +fl|5+(ÍA±«í^-ir„ + In, = o
OX P Oy p p
(A.31)
By subtracting the depth integrated continuity equation and dividing by the total depth,
the final form of the x momentum equation is
8U dU_ dU_ 1
dt + dx dy + p{h0 + fj)
dST
dS.
xv
dx
dy .
(A.32)
112
dfj 1 dr¡
+ 9TZ + -T1-
: +
—rn* = 0
(A.33)
dx pdy p{h0 +p{h0 + r¡)
The same procedure gives as the depth integrated time averaged y momentum equation
dV TTdV rrdV 1 dSt
1- U h V 1-
dt dx dy p(h0 + fj) dx
xy
dS,
vv
dy .
dr¡ 1 dri 1
+ 9dyJr pdx p[h0 + f}) T'v
p(K + v)
rr Hv = 0
(A.34)
(A.35)
APPENDIX B
DERIVATION OF THE RADIATION STRESS TERMS
The basis of the derivation is to use the definition of the velocity potential to substitute
for the wave induced velocities in the definition for the radiation stresses.
d
, _ d
w =
zp dz
where is expressed in terms of the general complex amplitude function B as
_ 1 (_ig „coshk(h0 + z) _iut
2 V a
B cosh kh
+ cc)
where c.c. represents the complex conjugate term. Also r¡' is expressed as
r,' = l(Be-iut + c.c)
where the prime is used to denote wave induced quantities.
For convenience Eq. (2.50) of Section 2.2 is reproduced below.
S<*p = P u'au'pdz + 6ap ^ pdz - y {h0 + »j)2 + y(r?')2
An expression for p is obtained in Section 2.2 and is reproduced below
d IL .W1
dz - p{w'y
Substituting Eq. (B.5) into Eq. (B.4) gives
Sap = P I Ku'pdz + SaP
* ” /l o
(B.l)
(B.2)
(B.3)
(B.4)
(B.5)
ri „ r* ( r* du'w' \
J ^ P9Í1 ~ Z)dz + J ^ ( pj -^—dz\dz
' - .. ' _ ° ' * * '
-/
p[w'Y dz - Y{K + r))2 + ^-(»?')2
-«O L
(B.6)
113
114
The terms of Eq. (B.6) are now evaluated term by term.
Term 1:
, _ 1 / ig dB cosh k[h0 + z) _iut ig dB' cosh k(h0 + z) iu)t
u'au'p —
, = 1 (_
U“ 2 \ a dxa cosh kh
ff2 cosh2 k(h0 + 2)
a dxD
cosh kh
(B.7)
55 51T dB^dB_ _ dB_dB_e-2iut _ dB' dB* Xut
4a2 cosh2 kh dxa dxp dxa dip dxa dip
dxa dip
—r—r pg cosh2 k(h0 + z)
pU°Ut = 2Vsinh 2kh
pg cosh2 k(h0 + 2)
Jfcsinh 2 kh
dB dB' dB' dB
+
dxa dip dxa dip
‘ (dB_W\
dip )
(B.8)
(B.9)
Í pu'au’pdz =
* “* o
P9
•1
( dB 55* y
k sinh 2kh
^5xa dip )
P9
*1
( dB 55* \1
ksmh2kh
\dxa dxp )\
£»(****:) i(i+-
4 dip J k2 \ si
/ cosh2 k(h0 + 2) dz
• hp
^ (sinh 2k(h0 + fj) + 2k[h + fj))
2kh \
sinh 2khJ
The second term:
J_h pg{v - z)dz = pg ^r)z -
-h«
... „(* + *)*
-P9^~
The third term:
f If *ff dB LUÍ.U ".yo T “J _lu,f , • ># I.UOU T *J
‘ ~ 2 { a dxt 6 + " “ 6
cosh k(h0 + 2) _iult ig dB' cosh k(h0 + z)
cosh kh a dx¡ cosh kh
igk psinh k(h0 + z) _lult ^ iff* D, sinh k[h0 + 2) j
cosh kh
g sinh 2A:(h0 + 2)
• »it' — L. ^
4sinh2fch 5x, ' dx(
tv
, _ 1 / «ffk sinhk[h0 + 2) -lult iff* .sinh fc(h0 + z) and 2fc(h+r?)coth2fc(/i+fj)-l 0. In deepwater, however, coth2Jfc(h+»j) “ 1
so that 2k(h + fj) coth2k(h + fj) - 1 =£ 2kh. If the derivative terms are non-zero, the Szz
116
term will be proportional to the depth. This would mean that deep water non-shoaling and
non-dissipative waves in varying depths would produce gradients in the radiation stresses.
For plane waves this is not a problem, since the complex amplitude terms are easily shown
to be zero. This is seen by substitution of B = ae,kx and B* = ae~tkx where a = |B|.
The fourth term is straight forward
igk psinh k(h0 + z) _iut igk sinh k(h0 + z) ,•
cosh kh
+ —B'
+z) _ B.^, _
4 a2 cosh2 kh ' >
7-7T2 _ 1 g2k2 sinh2 k(h + z) n gk sinh2 k{h0 + z)
( } "4 a2 cosh2 k[h + rj) ~ sinh2Jkh BB
(B.21)
(B.22)
(B.23)
- í p(w')2dz = -pgk\B\2 Í
J—h0 "—h0
2 sinh2 k(ha + z)
sinh2¿h ~ 4 1 1 V* sinh 2kh,
The fifth term needs no evaluation and is seen to cancel exactly the second term.
The sixth and last term of Eq. (B.6) is r¡')2 and is evaluated easily
r¡' = 1 (£e-‘wi + B*eiut)
(r/)2 = 1 (Be-W + 2 BB* + B*e2,wt)
SI
W=2^*= 2
fOT=fl*|2
which cancels part of the result from the fourth term.
Adding the above results for all the terms of Eq. (B.6) yields
P91
(B.25)
(B.26)
(B.27)
(B.28)
c -P9Í~(9BdB*\ 1 ^ , 2kh \
4 I \5za dip y A:2 V sinh2kh)
+
^af)
I |2—2kh
1 1 sinh 2kh
dxa dip
(|V*£|2 - k2\B\2)^(2kh coth 2kh - 1)] |
(B.29)
APPENDIX C
LAGRANGIAN DERIVATION OF THE GOVERNING WAVE EQUATION
The governing equation for wave motion can be derived through use of a variational
principle where a Lagragian is the function that is varied. Luke (1967) showed that a
Lagrangian for waves is given by
L — Jh [^‘ + 3Z
dz
(C.l)
Kirby (1984) used a time averaged Lagrangian to derive a wave equation. For the following
derivation, equation (C.l) will be used. The third term in Eq. (C.l) is determined straight
forward, while the first two terms need to be expanded about fj the mean water line in a
Taylor series.
L=f-’4+i:
4>t +
(W)J
dz + {r¡ - r¡)4>t |„ +(f7 - fj) ^
(C.2)
We assume that rj and can be expressed as
*) = V + c»?i
= o + efi (C.3)
where e is an ordering parameter and / is the depth dependency as defined in Eq. (3.3). For
the sake of simplicity we will assume that the gradient of fo represents the steady horizontal
currents. This means that fo is not a function of time and that W will not be included as
it was in the derivation in Chapter 3. This is because the W term is not expanded about
z = f¡ as it was in the Green’s function derivation and is not needed in order to obtain the
(Vfc • U)j£ term which as will be seen is obtained more directly in the present derivation.
Including W here will only produce terms that are of higher order in the mean current.
117
118
Introducing (C.3) into (C.2) yields
¿gvi gh20
L = ?Y + Wrh + ^Y±-9-£ + «t>itJ'h fdz+^-(h0 + fj) +
cVfc^oVfc<^i [ f dz + eV/,^0^1 ( Í Vnfdz+f ftdz\ +
J-h0 U-h„ J-Kc )
(yit + crti ( -—h£°^ + cVhfoVh4>i + eV^^o^i/* |*j +0(e'
! 2
The Lagrangian L is assumed to be of the form
>)}
(C.4)
It then follows that
L — Lq + fLi + e* L2 + fSZ>s + 1
¿0-i£-iS+M(*.+,)
and
(C.5)
(C.6)
¿1 = gf>m + —¿1, + —v^oV^! { r v*/ dz + r fzdz) + m (V^o)2 (c.7)
0 fl U-fco •/-*„ J 2
ti'Vhti Í [f^hf + / /*] oifx |#J (c.8)
* “/lo
To obtain the desired linear governing equation is varied with respect to i and then Li
is varied with respect to rji in order to eliminate r¡\. The variational principle holds that the
Lagrangian is constant under small variations of the unknown parameters. In mathematical
terms
and
DL
D\ is expressed in terms of its derivatives. This requires integration by parts to
119
obtain the variation with respect to \. Upon carrying this out it is found that
P± = dLi_ d dU _ dU x
D\ is arbitrary yields after
carrying out the integrations in ¿2
CC — k2CC
- V*(—*v*^) + M '-) - nu - v*(v*^o9i) = o (C.12)
9 9
which is rewritten as
+ ^h{CCgVhi) - i( i« + Vh\ = 0
(C.14)
which can be written as
lD!
Vl = pr-
g Dt
(C.15)
which upon substitution into (C.13) yields the governing equation
^ + (V^o)^ - Vfc(CC,V»*) + (a2 - k'CC'fa = 0
which upon making the substitution that 0 = V/,^o becomes
+ (v* * " v*(cc,v*^) + (
- -ig e«'(/ * CO. i dz)
and then make a phase shift by using the substitution
í co»f dx-J k cos 6 dx)
so that in the final form
_ ig A'(X>V)e*(f Icoiidz)
(D.l)
(D.2)
(D.3)
(D.4)
and all y dependency is solely in the complex amplitude A'. For reference equation (3.51)
is reproduced below.
(P4>z)z + (pv)v + k2p(¡> + [w2 - a2 + i'w(Vfc • U)j -I- i(7w<£
-{U24>z)z-{V24>V)V-{UV4>X)V-{UVÍV)X + 2Í = o (d.5)
Writing out all the components of equation (D.5) and rearranging terms gives
{p ~ U2)xx + (p - U2)4>zx + [(p - V2)4>v] v + k2p4> + (a;2 - o2)*
+iojUx4> + iuVv4> + icrw(j> - (UVx)v - [UV~~*v)x + 2iuU**x + 2iuV^>v = 0 (D.6)
The derivatives of defined by equation D.2 are as follows
4>x = —ig [(^) + ik COS 9^ c‘(/ k cot 0dz)
120
121
ZZ —
+ 2tfccos0^ + i(Jbcosi)a
— k2 cos2 9 —1 fc co* *dx)
a a J
{4>z)v —
"'I
{(a) e*(^*co,#x)v
The above forms of the derivatives recognize that the phase function e‘(/ico»edz) y
derivatives since k = k(x, y). However the y derivatives of the phase function integral are
not evaluated since the phase shift (equation D.3) will eliminate any y dependance from
the phase function and place all y dependance in the complex amplitude. The ( —)
v ° / XX
term is neglected following the assumption explained in section 3.2. Substituting the above
derivatives into equation (D.6) gives
(p-U2)x (-) e*'(/tco,#dl) + ikcos9(p - U*)g-ei
-2 UV
! 6 dx)
-tJfccos 9UUV —e'(f kco,e dl)
—ik cos 9UVv —*co*9dz^ - i{kcos 9)vUVkcot 9 dl>
a a
-2ikcos9UV (áei(fk“»edx)'j — i(kcos 0)vUV—kcot9dz)
r ^ S
~(UV)z (^e‘(/4co,í>
+ [P(ae,^fcCO,i<,l)) + 2s — Jbcos17)] ^áeHftcctdx)J = ^D 8)
It is easily shown that
kcos0{p - U2) +- kcosOU
it is also easily shown that
ui — a = 2o>k cos 6U — {k cos 6 Uy
From equation D.2 it is apparent that for a constant o/,
(fccos0)„ = -ov
Thus equation (D.8), becomes
{• [a{Ct cosí + U)]x j + 2ia(Ct cosí + E7) ^
+ [fc2p(l - cos2 Í)] J + iwA + IVvA + «7vV^} e‘(/*“•«<**)
+2iVa ^-ei(fkcot8dx)'j + |p tco. tdx)^
where a i(k cos 0)vUV£ term has been dropped.
At this point a phase shift is introduced
v.
= 0
(D.9)
(DIO)
(DU)
(D.12)
(D.13)
^ _ ¿IgUf ico»f dx-f k co« 6 dz)
(D.14)
123
where k and 0 are a reference wave number and wave angle. For computational purposes k
will be an averaged wave number across a grid row and 6 will be taken as the wave angle in
the absence of currents using Snell’s law. Further it will be assumed solely for computational
purposes that 0 is equal to 8. This gives
»[
as being comprised of two components. In fact this procedure can also be used to obtain
the radiation stress terms for any number of components. The only limit being the amount
of computation required as the number of components increases.
Smith (1987) following the work of Miles (1967) for wave reflection and transmission
over a step derived equations for the complex reflection and transmission from a submerged
plane barrier. An attempt was made using this formulation to obtain the circulation and
wave pattern from a submerged breakwater represented as a plane impermeable barrier.
However problems of numerical stability were encountered. It is suspected that the sharp
gradients between the region with the presence of the reflected wave and regions of its
absence adjacent to the barrier produced unrealistic gradients in the radiation stress terms.
An attempt was made to include the reflected wave for the case of an emerged plane barrier
with the same unstable result. This is a topic for further investigation.
6.3 Extending the Model to Include Sediment Transport
The model as it currently stands can be extended to include sediment transport. The
wave model can be slightly amended to retain wave energy dissipation as an array. If a
formulation of sediment entrainment can be made in terms of the wave energy dissipation
and the bottom shear stress, everything else is in place to accomplish the sediment transport
model. A simple numerical mechanism would be to assume steady state conditions for
certain time intervals and solve for the wave and current field and the sediment concentration
106
as a function of known variables and then solve a separate model for sediment balance and
then update the depths and proceed iterating between the sediment balance model and the
present wave-current model marching in time. Since it has already been demonstrated that
the wave- current model gives adequate results in the lee of a shore parallel breakwater, a
good test for such a sediment transport model wold be to see if it could yield a tombolo
formation in the lee of a shore parallel breakwater.
APPENDIX A
DERIVATION OF THE DEPTH AVERAGED EQUATIONS
In this appendix a derivation of the governing equations for the circulation model is
presented. This appendix closely follows the work of Ebersole and Dalrymple (1979).
Before proceeding with the derivations of the depth integrated time averaged equations
of motion and the continuity equation the Leibnitz Rule for taking the derivative of an
integral is stated since it is invoked throughout the derivation.
Ty C /(I'V) = C dJ^Tix + my)'y)dJ£- - '^y)’y)dJJT
The Leibnitz Rule can also be used in reverse to take derivatives outside of an integral.
The kinematic boundary conditions are also reproduced since the use of the Leibnitz Rule
will yield bottom and surface terms that will be eliminated through use of these boundary
conditions. On the free surface the boundary condition is
dr) dr)
dr)
Wr' = Tt+'x-aiJrv"Ty
(A.2)
and at the bottom the boundary condition
(A.3)
In the circulation model it is assumed that the still water depth does not vary with time so
the first term of Eq. (A.3) is dropped (actually for time varying depth essentially the same
equations are obtained).
The continuity equation (or the conservation of mass equation as it is abo called) in
three dimensions is
, d(/>u) d(pv) d(pw)
dt dx dy dz
(A.4)
107
108
Integrating over the depth from z — —h0 to z = r¡ and using the Leibnitz Rule to take
derivatives outside the integrad, the continuity equation becomes
d f* . dr] d [i
atLpi’-p'Tt+TzL'n,i'
h0
f \ dr> r \ dh0 d n
Si V—■
/ s
(A.5)
The terms with an overbrace are p times the free surface kinematic boundary condition
and the terms with an underbrace are p times the bottom boundary condition. Eliminating
both of these from Eq. (A.5) gives
d n . dr* . d rn J
/ Pdz+*~ pudz + — pvdz- 0
dt J-h0 ox J—k0 dvJ-h0
(A.6)
Equation (A.6) is now time averaged. Time averaging means averaging over a wave
period. The time average of a term is denoted by an overbar. Thus
^ 1 fT
■-H
Fdt
(A.7)
where F represents any term and T is the wave period. Letting u and v be comprised of a
mean flow, a wave induced flow, and an arbitrary or random turbulence component
u = ü + ú + u'
V = V + V + v1
(A.8)
(A.9)
The turbulent fluctuations are assumed to be of a high frequency relative to the wave
frequency so that their time averages over the wave period are zero. By substituting the
above expressions for u and v into Eq. (A.6) and taking time averages the continuity equation
is
^ {p{K + ’?)}+ ^ iPüiho + »?)} + J¿f_k Púdz+ (A.10)
A {««(*.+ fl)} + = o
(AH)
109
The time averages of the vertically integrated wave induced velocities ü and v are not zero.
These quantities are known as the wave induced mass flux.
Defining
U = u + u
and
V = v + 5
(A.12)
where
u =
1 Jn
r / pudz and
o + fj)J-k0
p{h0 + *?)
the time averaged depth integrated continuity equation is
%í+-kuQlo+t,)+‘§¡¡viho+tl)
1 ri
v — ~Ti rr / pi> dz
p{h0 + fi)J-ho
(A.13)
(A.14)
The assumption that the density is constant in both time and space has used in obtaining
both Eq. (A.11) and Eq. (A.14). The assumption of a constant bottom has also been used
in obtaining (A.14). If the bottom is not constant in time a term would be included in
Eq. (A.14).
The x and y momentum equations are handled in the same manner as the continuity
equation above. The x and y direction momentum equations are
du du2 duv duw 1 dp 1 ,drxx drvx drxx.
Tt + + + IT = } Ú + + if + -sT>
dv t duv t dv2 t dvw _ 1 dp 1 fdrxv drvv drZ]/
dt dx dy dz
= zir + 1A-
+
}
(A.15)
(A.16)
p dy p dx dy dz
The same procedure is followed for both the x and y equations. Only Eq. (A.15) will be
dealt with here. Integrating the left hand side over the depth and applying the Leibnitz
Rule to take the derivative outside the integral gives
it Tip\ d /*'J , dr] dh0
(LHS> = ai/-»." ~aT +
if +
dxj-h0 ndx ° dx
dy J-
, , x dr) . . dh0
uvdz - (uv)„ — - (uv).fco — +
h0 dy dy
(UW)l ~ (UW)-ho
(A.17)
110
The terms in Eq. (A.17) that are designated by an overbrace are tif, times the kinematic
free surface condition and those designated by an underbrace are u_j,0 times the bottom
boundary condition. These are both dropped leaving
, „ d r dr, dr
(LHS) = * U UdZ+d~x U UdZ+d~y U UV dZ ^
Substituting the definitions given for u and v in Eqn. (A.8) and (A.9) and recognizing that
the turbulent fluctuating components are of a high frequency relative to the wave frequency
so that the time averages of integrals involving the product of a turbulent component and
a wave induced component are zero gives
d~r d~r d ~r d~r
LHS = — / udz+ — údz + — u?dz+— tfdz (A.19)
dtJ-ko dt]-K0 dxJ-ho dxj-h0 ’
d ~r~T~Z d ~r d ~r d ~r
+ -£- {u')2dz + — / 2UÚ dz + — uv dz + — / uiidz (A.20)
dxJ-ho dxj-h0 dyJ-h0 dyJ-ho
d ~r d ~r d ~r
+ -X- üvdz+— Üi)dz + — u'v'dz (A.21)
oyj-ho oyJ-h0 dy J-K0
using the definitions of Eqs. (A.12) and (A.13) this becomes
+aiuvv,°+«)+TyS'h' \w^)Lrd‘Lriz
(A.22)
where the two terms designated by underbraces result from being added to complete the
U2 and UV terms.
The right hand side after integrating over the depth assuming a constant density and
assuming an inviscid fluid such that no horizontal viscous stresses exist, i.e. rtx and Tyx are
zero, is
RHS =
if r ^
p J-ho dx p J—h0 dz
(A.23)
Using the Leibnitz Rule and averaging over a wave period the right hand side becomes
1 d r . 1 dn 1 dhn 1 1
RHS = —— / pdz+ + -P-Ko-z— + -Tzz~~ -rnir (A.24)
pdxj-h0 p dx p dx p ” p v '
Ill
TZXq and rMX_ko are the time averaged surface and bottom shear stresses. Assuming that
the free surface pressure pn is zero (the dynamic free surface boundary condition) the right
hand side is
___ 1 d 7* 7 1 dh~0 11 ,
RHS = —— / pdz + -p.h—± + -T¿--Ti¿ (A.25)
pdxj-h0 p dx p p ’
The mean pressure at the bottom can be expressed as the sum of the wave induced or
dynamic pressure and the hydrostatic pressure
P-fco = PdVn-ko + P9[h0 + r>)
(A.26)
so that the product P-h0fa can be manipulated to read
dh0 1 dh0 13.,. _.2, . .dfj
'’-*•17 = + 2áí*~~' l(h°+
(A.27)
Terms from both the right and left hand sides are combined as the radiation stresses.
Defining
s„ = f'jp+,07) iz - \P9(h„ + «)* - rM'Y
~r> i ~ñ r>
S,^j ^vdz - ^ ,0 i, j^dz
Sn = /_ J» + ,07) dz - 1,,(*. + «)> - ^7) {/I P5 ^
(A.28)
(A.29)
(A. 30)
and making the assumptions that the product of the dynamic pressure at the bottom and
the bottom slope is negligible, that the Jif?h0(u')2dz term can be ignored and that the
lateral friction tj caused by momentum fluxes due to turbulent fluctuations is independent
of depth and that rj is defined by r¡ = —pu'v', the x direction equation of momentum is
written as
3 »r2,l , ^ , d rnr/L . . ldSxx . 1 dS,
+ --
XV
5^ + « + a-rU '<*• + *» + TyUV^ + + ,- 8x ' , 8„
+,(». +fl|5+(ÍA±«í^-ir„ + In, = o
OX P Oy p p
(A.31)
By subtracting the depth integrated continuity equation and dividing by the total depth,
the final form of the x momentum equation is
8U dU_ dU_ 1
dt + dx dy + p{h0 + fj)
dST
dS.
xv
dx
dy .
(A.32)
112
dfj 1 dr¡
+ 9TZ + -T1-
: +
—rn* = 0
(A.33)
dx pdy p{h0 +p{h0 + r¡)
The same procedure gives as the depth integrated time averaged y momentum equation
dV TTdV rrdV 1 dSt
1- U h V 1-
dt dx dy p(h0 + fj) dx
xy
dS,
vv
dy .
dr¡ 1 dri 1
+ 9dyJr pdx p[h0 + f}) T'v
p(K + v)
rr Hv = 0
(A.34)
(A.35)
APPENDIX B
DERIVATION OF THE RADIATION STRESS TERMS
The basis of the derivation is to use the definition of the velocity potential to substitute
for the wave induced velocities in the definition for the radiation stresses.
d
, _ d
w =
zp dz
where is expressed in terms of the general complex amplitude function B as
_ 1 (_ig „coshk(h0 + z) _iut
2 V a
B cosh kh
+ cc)
where c.c. represents the complex conjugate term. Also r¡' is expressed as
r,' = l(Be-iut + c.c)
where the prime is used to denote wave induced quantities.
For convenience Eq. (2.50) of Section 2.2 is reproduced below.
S<*p = P u'au'pdz + 6ap ^ pdz - y {h0 + »j)2 + y(r?')2
An expression for p is obtained in Section 2.2 and is reproduced below
d IL .W1
dz - p{w'y
Substituting Eq. (B.5) into Eq. (B.4) gives
Sap = P I Ku'pdz + SaP
* ” /l o
(B.l)
(B.2)
(B.3)
(B.4)
(B.5)
ri „ r* ( r* du'w' \
J ^ P9Í1 ~ Z)dz + J ^ ( pj -^—dz\dz
' - .. ' _ ° ' * * '
-/
p[w'Y dz - Y{K + r))2 + ^-(»?')2
-«O L
(B.6)
113
114
The terms of Eq. (B.6) are now evaluated term by term.
Term 1:
, _ 1 / ig dB cosh k[h0 + z) _iut ig dB' cosh k(h0 + z) iu)t
u'au'p —
, = 1 (_
U“ 2 \ a dxa cosh kh
ff2 cosh2 k(h0 + 2)
a dxD
cosh kh
(B.7)
55 51T dB^dB_ _ dB_dB_e-2iut _ dB' dB* Xut
4a2 cosh2 kh dxa dxp dxa dip dxa dip
dxa dip
—r—r pg cosh2 k(h0 + z)
pU°Ut = 2Vsinh 2kh
pg cosh2 k(h0 + 2)
Jfcsinh 2 kh
dB dB' dB' dB
+
dxa dip dxa dip
‘ (dB_W\
dip )
(B.8)
(B.9)
Í pu'au’pdz =
* “* o
P9
•1
( dB 55* y
k sinh 2kh
^5xa dip )
P9
*1
( dB 55* \1
ksmh2kh
\dxa dxp )\
£»(****:) i(i+-
4 dip J k2 \ si
/ cosh2 k(h0 + 2) dz
• hp
^ (sinh 2k(h0 + fj) + 2k[h + fj))
2kh \
sinh 2khJ
The second term:
J_h pg{v - z)dz = pg ^r)z -
-h«
... „(* + *)*
-P9^~
The third term:
f If *ff dB LUÍ.U ".yo T “J _lu,f , • ># I.UOU T *J
‘ ~ 2 { a dxt 6 + " “ 6
cosh k(h0 + 2) _iult ig dB' cosh k(h0 + z)
cosh kh a dx¡ cosh kh
igk psinh k(h0 + z) _lult ^ iff* D, sinh k[h0 + 2) j
cosh kh
g sinh 2A:(h0 + 2)
• »it' — L. ^
4sinh2fch 5x, ' dx(
tv
, _ 1 / «ffk sinhk[h0 + 2) -lult iff* .sinh fc(h0 + z) and 2fc(h+r?)coth2fc(/i+fj)-l 0. In deepwater, however, coth2Jfc(h+»j) “ 1
so that 2k(h + fj) coth2k(h + fj) - 1 =£ 2kh. If the derivative terms are non-zero, the Szz
116
term will be proportional to the depth. This would mean that deep water non-shoaling and
non-dissipative waves in varying depths would produce gradients in the radiation stresses.
For plane waves this is not a problem, since the complex amplitude terms are easily shown
to be zero. This is seen by substitution of B = ae,kx and B* = ae~tkx where a = |B|.
The fourth term is straight forward
igk psinh k(h0 + z) _iut igk sinh k(h0 + z) ,•
cosh kh
+ —B'
+z) _ B.^, _
4 a2 cosh2 kh ' >
7-7T2 _ 1 g2k2 sinh2 k(h + z) n gk sinh2 k{h0 + z)
( } "4 a2 cosh2 k[h + rj) ~ sinh2Jkh BB
(B.21)
(B.22)
(B.23)
- í p(w')2dz = -pgk\B\2 Í
J—h0 "—h0
2 sinh2 k(ha + z)
sinh2¿h ~ 4 1 1 V* sinh 2kh,
The fifth term needs no evaluation and is seen to cancel exactly the second term.
The sixth and last term of Eq. (B.6) is r¡')2 and is evaluated easily
r¡' = 1 (£e-‘wi + B*eiut)
(r/)2 = 1 (Be-W + 2 BB* + B*e2,wt)
SI
W=2^*= 2
fOT=fl*|2
which cancels part of the result from the fourth term.
Adding the above results for all the terms of Eq. (B.6) yields
P91
(B.25)
(B.26)
(B.27)
(B.28)
c -P9Í~(9BdB*\ 1 ^ , 2kh \
4 I \5za dip y A:2 V sinh2kh)
+
^af)
I |2—2kh
1 1 sinh 2kh
dxa dip
(|V*£|2 - k2\B\2)^(2kh coth 2kh - 1)] |
(B.29)
APPENDIX C
LAGRANGIAN DERIVATION OF THE GOVERNING WAVE EQUATION
The governing equation for wave motion can be derived through use of a variational
principle where a Lagragian is the function that is varied. Luke (1967) showed that a
Lagrangian for waves is given by
L — Jh [^‘ + 3Z
dz
(C.l)
Kirby (1984) used a time averaged Lagrangian to derive a wave equation. For the following
derivation, equation (C.l) will be used. The third term in Eq. (C.l) is determined straight
forward, while the first two terms need to be expanded about fj the mean water line in a
Taylor series.
L=f-’4+i:
4>t +
(W)J
dz + {r¡ - r¡)4>t |„ +(f7 - fj) ^
(C.2)
We assume that rj and can be expressed as
*) = V + c»?i
= o + efi (C.3)
where e is an ordering parameter and / is the depth dependency as defined in Eq. (3.3). For
the sake of simplicity we will assume that the gradient of fo represents the steady horizontal
currents. This means that fo is not a function of time and that W will not be included as
it was in the derivation in Chapter 3. This is because the W term is not expanded about
z = f¡ as it was in the Green’s function derivation and is not needed in order to obtain the
(Vfc • U)j£ term which as will be seen is obtained more directly in the present derivation.
Including W here will only produce terms that are of higher order in the mean current.
117
118
Introducing (C.3) into (C.2) yields
¿gvi gh20
L = ?Y + Wrh + ^Y±-9-£ + «t>itJ'h fdz+^-(h0 + fj) +
cVfc^oVfc<^i [ f dz + eV/,^0^1 ( Í Vnfdz+f ftdz\ +
J-h0 U-h„ J-Kc )
(yit + crti ( -—h£°^ + cVhfoVh4>i + eV^^o^i/* |*j +0(e'
! 2
The Lagrangian L is assumed to be of the form
>)}
(C.4)
It then follows that
L — Lq + fLi + e* L2 + fSZ>s + 1
¿0-i£-iS+M(*.+,)
and
(C.5)
(C.6)
¿1 = gf>m + —¿1, + —v^oV^! { r v*/ dz + r fzdz) + m (V^o)2 (c.7)
0 fl U-fco •/-*„ J 2
ti'Vhti Í [f^hf + / /*] oifx |#J (c.8)
* “/lo
To obtain the desired linear governing equation is varied with respect to i and then Li
is varied with respect to rji in order to eliminate r¡\. The variational principle holds that the
Lagrangian is constant under small variations of the unknown parameters. In mathematical
terms
and
DL
D\ is expressed in terms of its derivatives. This requires integration by parts to
119
obtain the variation with respect to \. Upon carrying this out it is found that
P± = dLi_ d dU _ dU x
D\ is arbitrary yields after
carrying out the integrations in ¿2
CC — k2CC
- V*(—*v*^) + M '-) - nu - v*(v*^o9i) = o (C.12)
9 9
which is rewritten as
+ ^h{CCgVhi) - i( i« + Vh\ = 0
(C.14)
which can be written as
lD!
Vl = pr-
g Dt
(C.15)
which upon substitution into (C.13) yields the governing equation
^ + (V^o)^ - Vfc(CC,V»*) + (a2 - k'CC'fa = 0
which upon making the substitution that 0 = V/,^o becomes
+ (v* * " v*(cc,v*^) + (
- -ig e«'(/ * CO. i dz)
and then make a phase shift by using the substitution
í co»f dx-J k cos 6 dx)
so that in the final form
_ ig A'(X>V)e*(f Icoiidz)
(D.l)
(D.2)
(D.3)
(D.4)
and all y dependency is solely in the complex amplitude A'. For reference equation (3.51)
is reproduced below.
(P4>z)z + (pv)v + k2p(¡> + [w2 - a2 + i'w(Vfc • U)j -I- i(7w<£
-{U24>z)z-{V24>V)V-{UV4>X)V-{UVÍV)X + 2Í = o (d.5)
Writing out all the components of equation (D.5) and rearranging terms gives
{p ~ U2)xx + (p - U2)4>zx + [(p - V2)4>v] v + k2p4> + (a;2 - o2)*
+iojUx4> + iuVv4> + icrw(j> - (UVx)v - [UV~~*v)x + 2iuU**x + 2iuV^>v = 0 (D.6)
The derivatives of defined by equation D.2 are as follows
4>x = —ig [(^) + ik COS 9^ c‘(/ k cot 0dz)
120
121
ZZ —
+ 2tfccos0^ + i(Jbcosi)a
— k2 cos2 9 —1 fc co* *dx)
a a J
{4>z)v —
"'I
{(a) e*(^*co,#x)v
The above forms of the derivatives recognize that the phase function e‘(/ico»edz) y
derivatives since k = k(x, y). However the y derivatives of the phase function integral are
not evaluated since the phase shift (equation D.3) will eliminate any y dependance from
the phase function and place all y dependance in the complex amplitude. The ( —)
v ° / XX
term is neglected following the assumption explained in section 3.2. Substituting the above
derivatives into equation (D.6) gives
(p-U2)x (-) e*'(/tco,#dl) + ikcos9(p - U*)g-ei
-2 UV
! 6 dx)
-tJfccos 9UUV —e'(f kco,e dl)
—ik cos 9UVv —*co*9dz^ - i{kcos 9)vUVkcot 9 dl>
a a
-2ikcos9UV (áei(fk“»edx)'j — i(kcos 0)vUV—kcot9dz)
r ^ S
~(UV)z (^e‘(/4co,í>
+ [P(ae,^fcCO,i<,l)) + 2s — Jbcos17)] ^áeHftcctdx)J = ^D 8)
It is easily shown that
kcos0{p - U2) +- kcosOU
it is also easily shown that
ui — a = 2o>k cos 6U — {k cos 6 Uy
From equation D.2 it is apparent that for a constant o/,
(fccos0)„ = -ov
Thus equation (D.8), becomes
{• [a{Ct cosí + U)]x j + 2ia(Ct cosí + E7) ^
+ [fc2p(l - cos2 Í)] J + iwA + IVvA + «7vV^} e‘(/*“•«<**)
+2iVa ^-ei(fkcot8dx)'j + |p tco. tdx)^
where a i(k cos 0)vUV£ term has been dropped.
At this point a phase shift is introduced
v.
= 0
(D.9)
(DIO)
(DU)
(D.12)
(D.13)
^ _ ¿IgUf ico»f dx-f k co« 6 dz)
(D.14)
123
where k and 0 are a reference wave number and wave angle. For computational purposes k
will be an averaged wave number across a grid row and 6 will be taken as the wave angle in
the absence of currents using Snell’s law. Further it will be assumed solely for computational
purposes that 0 is equal to 8. This gives
»[*