|
UFL/COEL-97/005
A NUMERICAL HYDRODYNAMIC MODEL FOR
INLET-RIVER SYSTEMS
by
Peter N. Seidle
Thesis
1997
A NUMERICAL HYDRODYNAMIC MODEL FOR INLET-RIVER SYSTEMS
By
PETER N. SEIDLE
A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
UNIVERSITY OF FLORIDA
1997
ACKNOWLEDGMENTS
I would like to thank my advisor, Dr. Wang, for his
guidance and time I have used in my research. I would also
like to thank Professors Dean and Thieke for serving on my
committee and reviewing this thesis.
I would like to thank the Coastal and Oceanographic
Engineering (COE) Department for providing me with a research
assistant position and funding for my studies. With the
opportunity and resources the COE Department has provided me,
I am able to fulfill my personal goals and can better
contribute to the field of coastal engineering. I would like
to thank each and every professor in the COE Department for
piqued my interest and curiosity in their areas of research.
I would like to especially thank Dr. Lihwa Lin for his
patience, help, and guidance in my studies and research. I
truly appreciate the help from Professor Sheng, Dr. Eduardo
Yasuda, and Justin Davis for helping me with numerical
modeling. Gratitude is owed to Wendy Smith for her help as
well. I also thank each and every fellow student who has
helped me with class work.
Many thanks to Becky, Helen, John, Lucy, Sandra, and
Sidney for the administrative and computer support.
And saving the best for last, I would like to thank my
parents for their love, support, and encouragement, as well as
Micheal O'Shea and Pete's Wicked.
TABLE OF CONTENTS
ACKNOWLEDGMENTS ................
LIST OF FIGURES ................
.......................... ii
.......................... vi
ABSTRACT ................................................. ix
1. INTRODUCTION ................
1.1 Background ............
1.2 Objective .............
2. GOVERNING EQUATIONS .........
3. NUMERICAL METHODS ...........
3.1 Linearized System ....
3.2 Grid System ...........
3.3 Non-Linear Terms ......
3.3.1 Friction Term ..
3.3.2 Convective Term
3.3.3 Diffusion Term
3.4 Numerical Stability ...
...........................1
...........................1
..................... .... 3
................... ......10
... ..... ............ .... .10
o ........................ 15
............... ..........17
.......................... 17
.......................... 19
..........................23
.......... ............... 24
4. BOUNDARY CONDITIONS ...................................27
4.1 Open Boundary Conditions ........................29
4.1.1 Tide .....................................29
4.1.2 Neumann ..................................30
4.1.3 River and Canal ..........................31
4.2 Solid Boudaries .................................35
5. MODEL
5.1
5.2
5.3
VERIFICATION ....................................36
Wind ............................................36
Seiching ........................................38
Massachisetts Bay ...............................41
6. MODEL VALIDATION ......................................49
6.1 Field Data ......................................49
6.1.1 Bathymetry ...............................50
6.1.2 Tide .....................................51
6.1.3 Current ..................................54
6.1.4 Drogues ..................................57
6.1.5 Wind .....................................60
6.2 Spectral Analysis ...............................60
6.3 Grid ............................................62
6.4 Open Boundary Conditions ........................64
6.4.1 Tidal OBC ................................64
6.4.2 Neumann OBC ..............................65
6.4.3 Canal OBC ................................66
6.5 Results .........................................68
7. MODEL
7.1
7.2
APPLICATION .....................................71
Single Jetty ....................................72
Pair of Jetties .................................73
8. CONCLUSIONS AND RECOMMENDATIONS .......................76
8.1 Conclusions .....................................76
8.2 Recommmendations ................................78
LIST OF REFERENCES .......................................81
BIOGRAPHICAL SKETCH ......................................83
LIST OF FIGURES
Figure Daage
2.1 Relationship between drag coefficient (shear stress
coefficient) and wind speed. ......................9
3.1 Rectilinear cell. .................................16
3.2 NumericalCells applied to friction term. .........18
3.3 Numerical cells applied to convective term. ......20
3.4 Numerical cells applied to convective term. ......22
4.1 Example of grid showing solid and open boundaries.
..................................................28
4.2 Example tidal and Neumann open boundaries. .......33
4.3 Example of a grid for a river or canal open boundary
condition. ........................................35
5.1 Comparison of analytical and numerically modeled
set-up. ...........................................38
5.2 Numerically modeled seich with 1 node. ..........40
5.3 Numerically modeled seich with 2 nodes. ..........41
5.4 Massachusetts Bay. ................................42
5.5 Geometry of 2D analytical solution. .............. 44
5.6 Analytical surface elevations at high tide and
maximim ebb velocity vectors. ....................45
5.7 Maximum velocity vectors from numerically modeling
Massachusetts Bay. ................................47
5.8 Numerically modeled surface elevations at high tide
in Massachusetts Bay. ............................48
6.1 Bathymetric survey of Ponce de Leon Inlet. .......51
6.2 Location of stations where tide and current data
were collected. ...................................52
6.3 Continuously observed tidal record and discrete
current measurements. .............................53
6.4 Flood velocity data taken inside the canal. ......55
6.5 Ebb velocity data taken in the canal. ............56
6.6 Flood current patterns studied using drogues. ....58
6.7 Ebb flow patterns studied using drogues. .........59
6.8 One tidal ccyle extracted from tide data. ........61
6.9 Discrete spectrum of the tidaal cycle. ........... 62
6.10 Area of Ponce de Leon Inlet which the grid is
applied to. .......................................63
6.11 Simulated tidal OBC signal compared to the observed
tidal signal. .....................................66
6.12 Maximum flood velocity vectors for numerical model
validation. .......................................68
6.13 Maximum ebb velocity vetcors showing numerical model
validation. .......................................69
7.1 Flood velocity vecors of a single jetty applied to
Punta Gorda, FL. ..................................72
7.2 Ebb velocity vecors of a single jetty applied to
Punta Gorda, FL. ..................................73
7.3 Flood velocity vecors of a pair of jetties applied
to Punta Gorda, FL. ...............................74
7.4 Ebb velocity vecors of a pair of jetties applied to
Punta Gorda, FL. ..................................75
8.1 Bathymetric smoothing. ...........................68
8.2 An idealized channel/bay system. .................80
Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science
A NUMERICAL HYDRODYNAMIC MODEL FOR
INLET-RIVER SYSTEMS
By
Peter N. Seidle
May, 1997
Chairperson: Dr. Hsiang Wang
Major Department: Coastal and Oceanographic Engineering
A numerical model is developed for predicting tidal
induced near shore hydrodynamics of an inlet-river system or
inlet canal network. The model couples a two-dimensional
model for the inlet region and a one-dimensional model
representing the river and/or canal network. The model is
based on depth averaged equations of motion including the
convective acceleration, friction, turbulent diffusion, and
wind stress terms. The numerical method used to solve the
equations is the factorization scheme. The model can be
driven using tides and river discharges. A canal system is
treated as a river with zero discharge.
The model is verified by comparing known analytical
solutions to the numerical results for selected cases. The
numerical model is then applied to Ponce de Leon Inlet which
serves a canal system for the City of Punta Gorda, FL. The
model is used to study the effects on the hydrodynamic
processes when jetties are added to Ponce de Leon Inlet.
CHAPTER 1
INTRODUCTION
Coastal regions are important to communities for their
leisure purposes and resources. Access to these resources may
be through coastal waterways, residential canal system,
harbor, etc. It is important to keep these resources
accessible not only for leisure purposes, but also for
commercial, fish and wildlife, and ecological reasons. There
have been numerous studies and projects to learn how keep
these resources open, accessible, and stable in all aspects.
One area of study is numerical modeling. This thesis will
develop a numerical model such that hydrodynamic studies of
inlet improvement measures such as jetties, rip-rap fill,
breakwaters, groins, etc., can be evaluated.
1.1 Background
It is important to understand the hydrodynamics of water
bodies, large and small, so that the benefits gained from the
coastal resources can grow, or at least remain stable.
Technology has given us the tools so that these hydrodynamic
2
studies can be conducted in a faster, more efficient manner.
Numerical models, confirmed by analytical models and augmented
by physical models, can be developed to learn more about the
micro- and macro-scale processes which govern the
hydrodynamics of coastal waterways, harbors, canals, and
inlets.
As an aid to the design of inlet and waterway
improvements, the effects of the tidal circulation, a
hydrodynamic process, near the inlet improvement measures can
be studied using numerical modeling. Assuming that the tidal
flow is not highly stratified and that the inlet width is
large compared to its depth, 2D numerical modeling should be
adequate to predict the tidal flows and velocities.
As previously mentioned, researchers have conducted
innumerable studies and published many papers about
hydrodynamics of estuaries and numerically modeling tidal
areas. Arcilla (1989) developed a general algorithm for
generating numerical models. Estuaries and bays are a focus
of study in two-dimensional circulation: Cheng (1982) studied
Tampa Bay, FL for circulation; Schomaker (1983) developed a
circulation model of Monterey Bay, CA; Reidel and Wilkinson
(1976) studied Cockburn Sound, Western Australia; Reeve and
Hiley (1992) studied the southern North Sea; Bathen (1976)
developed a model to study Keehi Lagoon, HI; etc. But these
3
models were developed for larger basins than a navigational
inlet and are only applicable to the water body for which they
were designed, so a more generalized two-dimensional numerical
model is needed which can be applied to different types and
sizes of inlets.
1.2 Objective
The objective of this project is to develop an implicit
two-dimensional numerical model which can easily be applied to
a multitude of inlets with different geometrical
configurations. The model will also be designed to aid in the
placement of jetties, groins, rip-rap fill, or other such
inlet improvements by simulating the tidal flows. Analysis of
this tidal process will allow for the optimum placement of
inlet improvement measures to ensure maximum benefits from the
changes.
This report will (1) present the equations of continuity
and motion used; (2) derive the numerical method used; (3)
show the differencing techniques applied to the governing
equations; (3) verify the numerical models efficacy; (4)
validate the model using data obtained from Punta Gorda, FL;
and (5) show how the model is applied as an aid to inlet
improvement design.
CHAPTER 2
GOVERNING EQUATIONS
The governing equations of motion used for this model
were derived from the Navier-Stokes equations of fluid motion
for a non-compressible, Newtonian fluid. The simplifying
approximations applied were as follows:
1)By order-of-magnitude analysis, the vertical
accelerations are negligible compared to gravity.
Hence, a hydrostatic approximation is applied.
2)Since the Rossby number, the ratio of Coriolis
force to inertial forces, is small for inlets this
model is intended for, the Coriolis force is
neglected.
3)Assume the turbulent Reynolds stresses are
proportional to velocity gradients by an eddy
viscosity coefficient and that the viscous
stresses are negligible compared to the turbulent
5
stresses. Also, for simplicity, assume the eddy
viscosity coefficient is constant for the inlet.
4)Assume atmospheric pressure is constant and
uniform for the region being modeled.
5)Assuming the inlet is well mixed and that there
are no density driven currents, the baroclinic
pressure gradient can be neglected.
6)Boussinesq approximation is applied to the
barotropic pressure gradient. This assumes that
the density gradient is small enough such that the
density is assumed to be uniform.
The continuity and resulting equations of motion are as
follows:
a +---uD+vD-= (2.1)
at ax ay
au au au a a2u a2u nd (2
-+u-a+v-=-g-- (-+-)+ (2.2)
at ax ay ax 8D ax2 ay2 pD
9v 9v av arn f 2 ua2v cv)y-wnd "
-+u- +v- = -g-i f v+A(-+- )+- (2.3)
at dx ay ay 8D x2 ay pD
where
t=time
u=depth averaged velocity in the x-direction,
v=depth averaged velocity in the y-direction,
u=depth averaged velocity vector,
D=total depth,
ri=water surface elevation,
g=gravity,
x,y=right handed Cartesian coordinates,
f=Darcy-Weisbach friction coefficient,
A,=horizontal diffusion coefficient,
Tx=x-direction wind shear,
ty=y-direction wind shear, and
p=density of sea water.
There are two empirical coefficients in the above
equations, the Darcy-Weisbach friction coefficient and the
7
eddy viscosity coefficient. In numerical modeling they are
often adjusted to tune the model based upon field data. There
are numerous studies concerning the friction coefficient for
open channel flows. In engineering applications, the most
common methods of developing a friction factor use the Manning
coefficient or Chezy coefficient. The friction coefficient is
a function of many types of environmental conditions, such as
type and amount of vegetation on both the bottom and banks,
and some of these conditions have been quantified in terms of
the Manning coefficient (Henderson, 1966). For a sandy bottom
along an open coast, a Darcy-Weisbach friction coefficient of
approximately 0.01 is commonly used.
The value, as well as the range of variability, of the
eddy viscosity coefficient is far less well established. This
is because eddy viscosity is a flow property which is related
to turbulent mixing. Since flow turbulence is a complicated
phenomenon and the scale of mixing varies greatly from case to
case, it is not a simple matter to model its mixing effects
using only one representative coefficient. In an oceanic
environment the eddy viscosity coefficient is known to vary
enormously. The order of magnitude for the coefficient
suggested by Pedolsky (1979) ranges from 10 to 10000. For
each scenario numerically modeled, a stability condition
(discussed in chapter 3.4) must be considered in selecting the
8
eddy viscosity coefficient. The stability condition relates
the grid size of the numerical model and time step of the
numerical method to the coefficient.
The wind shear, tnd is somewhat empirical as well and
may be necessary to validate the numerical model. If field
studies were conducted on a windy day, then the wind must be
accounted for in order to validate the model. Once the
numerical model has been validated, the wind speed and
direction may be changed or set to zero depending on the long
term predominant conditions.
The wind shear forcing term can be written as:
x-wind P. C
air'd- 2 WsinO (2.4)
pD pD
y-Wfod P. dW cosO (2.5)
pD pD
where
Cr=drag coefficient,
W1o=wind speed at 10m, and
O=wind bearing.
Although the drag coefficient is empirical, researchers have
uLuuI
C
~ 0.002
0001
I
fooi A
O .
-----------------------------------------'----o-
N K
..5 0 15 20 25 .. 30
0 10
20 30
Windsioeed 17.o
40 50
Figure 2.1
Relationship between drag coefficient (shear
stress coefficient) and wind speed. Figure was
taken from Roll (1965).
found that it is on the order of 10-3. The drag coefficient
could be found graphically from figure 2.1. Also, Wilson
(1960) concluded that for light winds (typically less that 10
m/sec), Cd=1.49 x 10-3 and for stronger winds, Cd=2.37 x 10-3.
For validating this numerical model, Wilson's findings were
used.
The governing equations (2.1, 2.1, and 2.3) are solved
numerically for the problems posed in this thesis. The
numerical scheme and the applicable boundary conditions are
discussed in chapters 3 and 4.
-, ,. .. ,.
I _
:
.. .
I
CHAPTER 3
NUMERICAL METHODS
The numerical method adopted for the model is the so-
called factorizationn" scheme. This method of implicitly
solving a set of linearized equations was developed to study
tidal currents in the vicinity of the Mississippi Sound in the
Gulf of Mexico (Sheng, 1983).
3.1 Linearized System
The governing equations (2.1, 2.2, and 2.3) are
linearized in order to develop the factorizationn" numerical
method. The governing equations now become:
By au av
L+D- +D-=O (3.1)
at ax ay
u--g (3.2)
at ax
--g (3.3)
at ay
The above equations 3.1, 3.2, and 3.3 can be combined as
a system of equations by letting w=[r, u, v]T The system of
equations can be written as:
Iwt+Aw,+Bw =O (3.4)
where I is an identity matrix and A and B are coefficient
matrices as follows:
0 D 0
A=g 0 0
000
0 0 D
B= 0 0 0
g0 0
The subscripts t, x, and y indicate the first partial
derivative respective to the subscript. Applying a forward
finite difference approximation to the time derivative yields:
n+1 n
w n+ w +! n 1 nI I
At +Aw. +BWy =0 (3.5)
At x y
Multiplying by At and factoring out w"'l results in:
w"n I(+AtAx+AtB6 )=w (3.6)
where 6 is the finite difference respective to the subscript
and the superscripts indicate the time step. The left hand
side (LHS) of the equation can be expanded to:
w"'(l +AtA6x+AtB8 +(At)2ABx y)-02=W (3.7)
where 02 denotes a 2nd order error term. This term will be
resolved and removed in later steps of this derivation. A
term of AtB6 is added and subtracted to the right hand side
(RHS) of the equation with no net change. The LHS of Equation
3.7 is then factored into two terms to give:
wn"(1 +AtA6x)(1+AtB6)-02=(1 -AtB8)w "+AtB6w (3.8)
An intermediate time step may be introduced, which is
also referred to as time splitting. Letting w* represent the
1-AItB8
intermediate time step, n+ w*=-A w" and equation 3.8
1 +AtA86
can be rewritten as:
w"n'(1 +AtA86)(1 +AtB8y)- 2=(1 +AtA6)w *+AtB8yw "
Moving the error term to the RHS of the equation and re-
expanding it the above equation then becomes:
w"+(1 +AtAWx)( +AtB6y)=(l +tA8x)w +AtB86w "+(At)2ABx yw "
(3.10)
Factoring out AtB6y from the two last terms on the RHS yields:
wn'(1 +AtA6 )(1 +AtB8y)=(1 +AtA8)w +AtB (1 +AtA8x)w"
(3.11)
Lastly, removing the common term, 1+AtA6x from the equation
results in:
(1 +AtB86)wn+l = +AtB8yW n
(3.12)
In summary, the time splitting equations used to solved
the linearized system of governing equations (3.1, 3.2, and
3.3) are written as follows:
(3.13)
(3.14)
(1 +AtAx)w *=(1 -AtBS y)
(1 +AtB y)w"n+=w +AtB8 "
Expanding equation 3.13 gives the following:
Tq +AtD8xu* ="-AtD6v "
(3.15)
(3.9)
u +Atg6xr*=" (3.16)
v*=v "-A tgyl" (3.17)
and expanding equation 3.14 gives:
rl' n+AtD8 vn"'+=i*+AtD8 v" (3.18)
unI=u* (3.19)
v "ln+Atg6 n+1 =v* +Atg6T"n (3.20)
Note that equations 3.15 and 3.16 are a system of
simultaneous equations where *' and u* are the two
unknowns, as well as equations 3.18 and 3.20,
where 1"n' and v"nI are the unknowns. Since equation 3.17 is
independent of q* and u* it can be substituted directly
into equation 3.20. Likewise, equation 3.19 states that
u" =u* so equations 3.15, 3.16, 3.18, and 3.20 can be
adjusted accordingly. The four resulting equations are as
follows:
S* +AtD6xu ''n n-AtD6v (3.21)
x f"fDyn
Atgx6* +un"+ =u"
(3.22)
11"l+AtD8 v"nl=*+AtD8 v n (3.23)
Atg68 n+n+n=V" (3.24)
The above factorizationn" method is similar to the
Alternate Direction Implicit (ADI) method in that equations
3.21 and 3.25 are solved in the x-direction (x-sweep) and that
equations 3.23 and 3.24 are solved in the y-direction (y-
sweep). It is also important to note that equations 3.22 and
3.24 are essentially the momentum equations in the x- and y-
directions, respectively. This will allow the non-linear
terms to be added. Moreover, when equation 3.21 is
substituted into equation 3.23 for i* ,the original continuity
equation (Equation 3.1) is recovered.
3.2 Grid System
A rectilinear staggered grid system is used for this
numerical model, shown in Figure 3.1. In applying finite
differences to equations 3.21 through 3.24 the staggered grid
system is examined more closely. A forward difference is
vij+1
ui+lj
Figure 3.1
Rectilinear cell
applied to the velocities for the continuity equations (2.21
and 3.23) and may be written as:
(3.25)
(3.26)
*A / i n+ .
Atg( -1_)+uJ =Uij
A backward difference is applied to the surface elevation for
the momentum equations (3.22 and 3.24) and are as follows:
n+il n++ n+l-t D vA,,- n n
4j +alvtvi,,l-v )=I Vi+j)(viji vi.
A -g(, nl -n+1r_ n+1 n
4gr i -r l ),i + =v
(3.27)
(3.28)
3.3 Non-Linear Terms
The non-linear terms in the governing x- and y- momentum
equations (2.2 and 2.3) which have been neglected thus far can
now be added to the numerical x- and y- momentum equations
(3.26 and 3.28), respectively. It is necessary to examine the
adjacent elements of the grid shown in figure 3.1 in order to
properly examine each non-linear term of the momentum
equations individually.
3.3.1 Friction Term:
The frictional forcing term in the x-direction momentum
equation is f-uu This term is non-linear in that it is
8D
quadratic in terms of flow velocity vector. Here the velocity
vector term, lul is expressed numerically in terms of the
known velocity components in the previous time step, such as u"
and v" Furthermore, since the x-direction momentum
n+]
equation, eq. 3.26, is solving for u,, at the grid center
depicted in Figure 3.2, lul must be calculated at the proper
location. The numerical equation used to calculate the
Figure 3.2
Numerical cells applied to friction
term.
magnitude of the velocity vector for the x-direction momentum
equation is as follows:
I|I= .ci2 +i 4 ) (3.29)
The numerical expression of the frictional forcing term which
is added to the LHS of x-direction momentum equation, eq.
3.26, is given as follows:
n n n n
f ( i-+(--- ij+ i-Iij )2 n+l
(ui+2 +
8D1 4
Similarly, the frictional forcing term added to the LHS
of the numerical y-direction momentum equation, eq. 3.28, is:
n n n n
S+ + +Ui+- 2 n+1
(( -) + vi
8D1 4
3.3.2 Convective Terms:
The convective terms have directional characteristics
which must be appropriately accounted for in the numerical
scheme. That is, the down stream behavior is affected by the
upstream convection phenomenon (Hirsch, 1994). A preferred
differencing technique is the upwind method, which is similar
to a backward difference written with respect to the direction
of the velocity vector. Because of this directional
phenomenon, the direction of the velocity must be
predetermined in the numerical expression of these convective
terms.
For the x-direction momentum equation, the convective
au au
acceleration terms are u-+v- taken from equation 2.2. In
ax ay
order to be included with the numerical equation of motion
(equation 3.26), finite differencing must be applied to these
convective acceleration terms such that they are congruous
Figure 3.3
Numerical cells applied to convective
term.
with the upstream convection phenomenon. Both convective
au au
terms, u- and v- need to be examined separately.
ax ay
au
To analyze the term u- first the velocity direction
ax
needs to be determined. An average of the velocities between
the two cells (figure 3.3), u_-1, and u,+12, are used to
identify the direction of the velocity so that the proper
n
upwind differencing method may be applied. If uin2j is
positive, then the upwind difference is u "-u"_ The
numerical equation used to evaluate the first convective term
is:
S ui ui+
n n
au n uijUi-j~ (3.30)
x Ax
However, if ui+2, is negative, then the upwind deference
is u-i-ul and the equivalent numerical difference is as
follows:
n n
u -" ui1 -uii (3.31)
U
ax /2Ax
au
The second convective term, v- can now be examined
ay
with the grid shown in figure 3.4 using the same steps in
analyzing the first convective term. If vi"n/2j is in the
positive direction, then the numerical convection term is as
follows:
n n
v vn _i ij-1 (3.32)
ay v'-1/ Ay
And, if vil"/2,+1 is negative, the numerical differencing is:
n n
-Vu- n i ij+l (3.33)
-V = -V i-1/2j +1 1
ay .Ay
At every time step, the direction of the velocity is
determined for each element so that either equations 3.30 or
Uij+1l
A
Vi-1/2j+1
uu
U..
A
Vi-1/2j
Uij-1
Figure 3.4
Numerical cells applied to
the convective term.
3.31 and either equations 3.32 and 3.33 are added to the
numerical x-direction momentum equation, eq. 3.26.
The same process used to derive the numerical equations
au au
for u-+v- can be used to find the numerical equations for
ax ay
the convection terms in the y-direction momentum equation,
av av
u-+v- For a positive velocity direction, the two
ax ay
convective terms are as follows:
n n
Uav V U (3.34)
ax i-1/2J Ax
and
n n
a yn Vij v-1
V ay Vij-/2y
(3.35)
Likewise, if the velocity is negative then the terms are as:
n n
av n l ij
-u-a u
ax12 Ax/2h
(3.36)
and
n n
SV n Vij -Vij+l
- y- i+1/2
ay i+I Ay
(3.37)
3.3.3 Diffusion Term:
The diffusion term is linear provided that the diffusion
coefficient is not a function of position, time, velocity, or
water surface fluctuation. The numerical expression for the
x-direction diffusive term:
a2u 82u
A,(-+-)
ax2 ay2
is obtained by simply applying a center difference to yield:
Aun i-2u u' -2u+,
n n n n n n
A i(.lj ij-z --lj +u Ij+ 4 Iu-1)
Ax2 Ay2
which is included on the RHS of the numerical x-direction
momentum equation, eq. 3.26. The terms added to the RHS of
the y-direction momentum equation, eq. 3.28 are:
A nI-2V n n-2v" n
Ax2 Ay2
3.4 Numerical Stability
Numerical stability criteria for the non-linear system
presented here are not established at present. However, the
stability criteria of a linear system are used instead. There
are two criteria which need to be met for numerical stability:
one related to tide propagation and another to diffusion.
For propagation stability, the well-known Courant
criterion is invoked. The Courant Number used in the present
model is defined as:
Cr= t ( +lul+vl) (3.38)
FA7+Ay2
where
At=time step,
Ax=x-direction grid spacing,
Ay=y-direction grid spacing,
h=maximum depth,
g=gravity,
lul=maximum x-direction velocity superimposed on
the tidal velocities, and
Ivi=maximum y-direction velocity superimposed on the
tidal velocities.
For a chosen grid size the time step, At, is then constrained
by the above equation for a critical Courant number, Cr.
Generally, for implicit numerical techniques the critical
Courant number can be greater than unity, whereas for explicit
methods the Courant number must be less than unity to maintain
numerical stability. For the (ADI) method, the Courant number
is suggested to be less than 5 to flows that are not
consistent with observed nature (Weiyan, 1992). The advantage
of using the factorization method is that stability is
maintained for Courant numbers of about 10 (Sheng, 1983).
For diffusive stability the eddy viscosity coefficient is
examined. The following condition has been suggested by
Hirsch (1994):
A < 1 (3.39)
At 2
This criterion is used in the present model to guide the
selection of the diffusion coefficient depending on the
spatial grid sizes and the time step size.
In summary, the grid size, Ax and Ay, are first selected.
Next, the time step, At, and the eddy viscosity coefficient
are then selected using the Courant number relationship, eq.
3.38, and the aforementioned diffusion criterion by Hirsch
(1994), eq. 3.39.
CHAPTER 4
BOUNDARY CONDITIONS
This chapter addresses "geometrical" boundary conditions
(BCs) used to define, drive, and control the numerical model.
Since the rectilinear grid is used, the numerical model must
be able to differentiate land and water cells. The land and
water cells are what define the model's geometry, so that an
inlet, lagoon, or harbor can be numerically simulated using a
rectilinear grid.
The bottom and surface BCs are addressed by the governing
equations as friction and wind stress forcing terms,
respectively, which were defined in chapter 2. The following
sections in this chapter will address the details of applying
the land and water boundary conditions.
The "geometric" boundary conditions (BCs) used in the
numerical model can be broken into two types, open and solid
An open boundary (OB) is a cell which is located on the edge
of the grid that covers the modeling area. The input to the
model, such as tide, river discharge and length, or canal
length, are applied to the OB cells. Inside the modeling
j=J
1,5
1,4 .6,4
2,3 .6,3
2,2 .6,2
j=2
j=1
i=l i=2 i=I
Figure 4.1
Example of grid showing solid and open
boundaries.
area, some cells are designated either as land or water to
define the model's geometry. The land/water interface creates
a solid boundary. Figure 4.1 shows an arbitrary example.
Cells (6,2), (6,3), and (6,4) are defined as OBs, so input
must be applied to these cells to drive the model. Whereas
29
cells (1,4), (1,5), (2,2), and (1,5) are defined as land or
solid boundary cells.
4.1 Open Boundary Conditions
Four types of open boundary conditions (OBCs) are
incorporated in the model; tidal, Neumann, canal, and river.
Though a Neumann OBC does not drive the numerical model as a
tidal OBC would, it is used to control the hydrodynamic
stability of the numerical model in certain tidal OBCs.
There are some considerations to be made when selecting
the model's grid and designating the OB cells as either a
tide, river, or canal BC. In practice, one edge of a cell may
be designated as a tidal OB and a different edge may be
defined as a canal or river OB. These methods of application
will be discussed in the following sections.
4.1.1 Tide
The tidal OB is defined by water level changes such as a
tidal fluctuation of the form ==A cos(at+() where A is the
amplitude, a is the period, and q is the phase angle. Though
it is not necessary for all scenarios, the numerical model can
accept a real tide as input using up to four Fourier series
30
harmonics. The Fourier series representing a mixed tide would
be as follows:
4
l= C (A,cos(ot+(,)) (4.1)
i=1
When applying the tidal OB, the velocities remain as unknowns,
otherwise, the model would be over-constrained and could not
be validated. The edge of the tidal OB cells should be of
sufficient distance offshore from the region of interest.
4.1.2 Neumann
If the area that is going to be modeled requires placing
tide OB cells on the corners, then the a Neumann BC must be
incorporated. Figure 4.2 shows the grid which would be
applied to such a situation. Where the cells marked with a T
indicate where a tidal forcing condition is imposed. The
cells marked with an N are where a Neumann BC is imposed. The
Neumann BC requires that the surface elevation gradient equal
zero between the OB cells and the cells in the J-1 row or the
2nd row of cells.
JI IN I IN IN N I
J-1 T
T
T
T
2 T
1 N N N N T
Figure 4.2
Example of tidal and Neumann open boundary
conditions.
4.1.3 River and Canal
In order to apply the river or canal OB, the velocities
must be perpendicular to the edge of the grid. If there is
circulation or lateral currents, with respect to the river or
canal, then the model area must be extended out to where the
velocity vectors become perpendicular to the OB. Also, in
32
applying a river OB, the head of the river must not be
influenced by the tide. The discharge at the head of the
river is used to drive the model and should be constant,
though the model can be easily modified to allow for a time-
dependent river discharge. The model becomes unstable or may
produce flows that are not indicative of nature if the tide
propagates to the head of the river. A canal is treated as a
river, except that the discharge at the head is zero.
Moreover, since the discharge at the head of a canal is zero,
the model will be stable and valid if the tide propagates to
the head, or end, of the canal.
The length, depth, and discharge at the head, as well as
grid spacing and friction coefficient, are needed as input at
the river/canal OB (figure 4.3). In addition, the user needs
to define a grid, or cell, length for the river or canal
boundaries. The cell can be significantly larger compared to
that of cells in the model area. It is recommended that the
river or canal be divided into 15 to 30 cells. If any more
cells are used, then the computing time required to run the
model is compromised. The governing equations used by the
numerical code in the river or canal cells are as follows:
S+-uD=0 (4.2)
at ax
Discharge
at head
N' Cell
River or Canal OB
Cell Size
Figure 4.3
Example of a grid for a river or canal open boundary
condition.
du+ u = l- f
at ax ax 8D
(4.3)
where
u=depth averaged velocity in the direction of the
river or canal,
t=time,
D=total depth,
l=water surface elevation,
g=gravity,
x=direction of the river or canal, and
f=Darcy-Weisbach friction coefficient.
In the previous chapter (chapter 3.3), statements
regarding the governing equations dealing with the friction
coefficient and convective term also apply to the above
equations. The model uses the surface elevation from the
previous time step and the river/canal discharge as boundary
conditions to calculate the new velocity at the model's OB
cells. The computed velocities at the OB cells are then used
as a boundary condition to drive the model. The new surface
elevations calculated for the cells at the OB are then used as
input into the river/canal OB subroutine for the next time
step. Since the flow at the OB is assumed to be perpendicular
to the boundary, the transverse velocity at the OB is set to
zero. For example, since the velocity vector at the OB in
figure 4.3 should be in the x-direction, the v-velocity is set
to zero.
At all OB cells, the model becomes unstable if diffusive
or convective terms of the momentum equations are included.
This prevents eddies from occurring at the boundaries which
35
would cause instability of subsequent computations. If field
data shows that there is circulation at the edges of the
numerical model's grid, the grid area should be enlarged to an
area where there is no circulation found by the field studies.
4.2 Solid Boundaries
Finally, at land boundaries, the flow condition is
satisfied by setting the normal velocity component to zero.
There is also a non-slip condition being imposed so that the
velocity parallel to the boundary is set to zero as well.
CHAPTER 5
MODEL VERIFICATION
The numerical code must be compared to analytical cases
to verify the efficacy. The model's efficacy was tested using
wind set-up and seiching in a closed basin, and by replicating
tidal currents in Massachusetts Bay. Once the model is
verified, it can be confidently applied to various inlets.
5.1 Wind
The analytical set-up due to wind stress is (Choi, 1992):
(x)= wnd(x- ) (5.1)
pgD 2
where
uwnd =Shear stress by wind,
p=density of water,
x=distance from middle of basin,
g=gravity,
D=depth, and
L=length of basin.
A closed, flat, rectangular basin was used to develop the
analytical solution. The length of the test basin which the
wind was blowing across, was 450 m long, with a grid size of
50 m. The depth of the basin was constant at 2 M and the wind
speed was constant at 5 m/s. The wind shearing force,
calculated using the respective drag coefficient discussed
previously in the governing equations in chapter 2, equals
0.0448 Pa. Figure 5.1 shows the surface elevation of a cross
section from the numerical model compared to the analytical
wind set-up. Only one cross section is presented since, as
expected, every cross section has the same wind set-up
profile.
The reason there is a difference between the numerical
and analytical set-ups is that the model's output is rounded
to 0.1 mm, where the analytical line does not include any
round off. For example, at x=25 m, the numerical set-up is
-0.0004 m and the analytical set-up is -0.0004386 m. If the
numerical model's output included more significant digits for
this test, then a closer comparison could be made. Moreover,
if the analytical solution was rounded off to the 1/10 mm,
then the two lines would correspond exactly.
x10- Wind Set-up
3
2-
-1 -- Analytical
1 --- Numerical
-2 -
-3-
-4 *
0 50 100 150 200 250 300 350 400 450
Distance, (m)
Figure 5.1
Comparison of analytical and numerically modeled
set-up.
5.2 Seichina
The model was also tested for its ability to replicate a
standing wave in a closed basin. The basin length required
for seiching is:
1Ln (5.2)
2
where
l=basin length,
n=number of nodes, and
L=wavelength of long wave.
For a channel which is 2 m deep and has a tide period, T, of
12 hrs, the wavelength is 191,352 m. The basin lengths which
would form a seiche with 1 and 2 nodes are 95,676 m and
191,352 m, respectively.
The rectangular test basin used for both cases is 3 m
deep and 450 m wide. Also, for both cases (1 and 2 nodes),
the number of cells used for the width and length were 9 and
50, respectively. Which means that Ax for the two different
cases are 1913.52 m and 3827.04 m, and Ay was 50 m for both
cases as well. The Courant number for both cases was 2. The
model was driven by specifying te water surface elevation at
the right hand side of the basin as a function of
time, i=Acos 21 where A, the amplitude, was set to 0.1 m
and T, period, was set to 12 hrs.
Figures 5.2 and 5.3 present the surface elevation for a
cross section of the basin for 1 and 2 nodal standing wave at
times 1/4T, 1/2T, 3/4T, and T. Note that the numerical
results are reasonably good, however, they did exceed the
I I I I I I I I I
0 10 20 30 40 50 60 70 80 90 100
0.2 -- i -- i -- i -- i -- i -- i -- i --- i --- i -
0.2
0 10 20 30 40 50 60 70 80 90 100
0 10 20 30 40 50 60 70 80 90 100
0.1
00.1 .. ..-
0 10 20 30 40 50 60 70 80 90 100
Distance, (km)
Figure 5.2
Numerically modeled seiche with 1 node at time,
from top to bottom; (a) =1/4T, (b) = 1/2T, (c)
=3/4T, and (d) =T.
input amplitude because friction was set to 0.0001 in the
numerical model. The small friction coefficient did not allow
enough energy to be damped from the system. Though it is not
evident for times 1/4T and 3/4T, there were surface
fluctuations which were less than 0.05 mm. The initial
surface conditions were evaluated using a cosine function to
represent a seiche, so the model quickly reached equilibrium
and became accurate.
I I 0 I
10 20 30 40 50
60 70 80 90 100
-1
0 10 20 30 40 50 60 70 80 90 100
0.2
-0.2 LI
0 10 20 30 40 50 60 70 80 90 100
Distance, (km)
Figure 5.3
Numerically modeled seiche with 2 nodes at time,
from top to bottom; (a) = 1/4T, (b) = 1/2T, (c) =
3/4T, and (d) = T.
5.3 Massachusetts Bay
The final verification of the model's efficacy was to
reproduce the results obtained from an analytical solution of
the 2D equations of continuity and momentum with geometry
representative of Massachusetts Bay (shown in figure 5.4).
The analytical solution for tide, derived by Briggs and Madsen
0
0.2 -
0-
I
I I i I
:. .... t . .AL
Plymouth
Cape Cod
Canal,-0."
0.0 :.5, .. .10.0 .
Scae (Nautical milss)
V . . .. ... .. .
. 7Ioston
SLigh, h p
. : :
(+0.2) .
Iace
(0.0)
:.M ::. .I c
S I "..... . .
n .x .: ...", yi.. :. : .j
: .. .
: I, ..
-. r. -
Figure 5.4
Massachusetts Bay. Figure taken from Briggs and Madsen
(1973)
E.: ;:
.............................................
C. C
..t.r 4
N.ME,
W" 6 )D,:: ,
'"
/ / / 7 / / / /
43
(1973), to the linearized continuity and momentum equations
(eqs. 3.1, 3.2, and 3.3) is:
i-=Aoe t[osmY _-E 2mosinmyo(sinkx2 sinkjx) coshmyosk] (5.3)
n=1 mk, (x2 -x) sinhmyo
where
n=
Xo
o2 2
M = ---k.
n gh "
and Xo, x,, x2, and y, are the geometric dimensions of the bay
as defined in figure 5.5. For simplicity, the depth, h, was
kept constant. The tide at a single location anywhere in the
region must be known and is used to solve for Ao. Briggs and
Madsen used the tide at the Boston Light House as the driving
force in their analysis. However, Boston Light House is
inside the bay, not at an open boundary, so the tide
measurement there could not be used to drive this numerical
model. From examining the tide (figure 5.6) at the open
boundary, found analytically by Briggs and Madsen, a value of
44
4ri
Figure 5.5
Geometry of 2D analytical solution. Figure taken from
Briggs and Madsen (1973).
4.30 ft can be used to drive the numerical model. The depth
that Briggs and Madsen used was 120 ft. They also used a
sinusoidal tide as input at Boston Light House, so it will be
assumed that the tide at the entrance to the bay will be a
sinusoidal tide as well. Figure 5.6 also shows the
analytically derived maximum ebb currents for Massachusetts
Bay.
Bay.
U I-
*
aS
o e
o 0
co
:CO
o "
analytical surface elevations at high
and Madsen (1973)
0
-4t-
4.
a In0
-s0
c.n
-tru
and Madsen (1973).
For the numerical model, the bay was divided into 30
increments in the x-direction and 10 increments in the y-
direction. The Courant number used for the model was 2.2, the
Darcy-Weisbach friction coefficient was 0.05, and the
diffusion coefficient was set to zero since the analytical
solutions neglected the diffusive terms. Note that the
analytical solutions do not include the friction or convective
terms. Briggs and Madsen used order-of-magnitude comparison
to neglect the friction term so that the governing equations
could be solved easier. Friction had to be include in the
numerical model so that energy could be dissipated and the
model would remain stable. A high friction coefficient was
used because of the deepness of the Massachusetts Bay. This
caused a discrepancy comparing the numerical results to the
analytical solutions.
Figure 5.7 shows the maximum ebb velocity vectors
developed by the numerical model. By visual comparison the
numerical model is in good agreement with the analytical
results. Figure 5.8 shows the numerically modeled surface
elevations at high tide in the bay. Again, the discrepancies
could be due to the friction and convective terms which the
numerical model includes, and the analytical solution
neglects.
47
Velocity Vectors of Massachusetts Bay
II I I I I
0 10 20 30 40 50 6(
Nautical Miles
Figure 5.7
Maximum velocity vectors from numerically
modeling Massachusetts Bay.
, ttt it1 \\\.... .....
. l4 %t l \ \ \ \ \ \\ ... .
I t t I\ .- . .
- - - - - .. .
l l l ~\1 ...sec
* .. ...... =1 ftlsec
Surface Elavation (ft) of Massachusetts Bay
0 10 20 30 40 50
Nautical Miles
Figure 5.8
Nummerically modeled surface elevations at high
tide in Massachusetts Bay.
CHAPTER 6
MODEL VALIDATION
Since the model's efficacy was verified, the numerical
model can now be applied to various inlets, but must be
validated for each inlet. Validation can be an iterative
process; the model's coefficients (friction, diffusion, and
wind drag) need to be tuned until the numerical results
simulate the actual hydrodynamic process (i.e. circulation and
eddies) and the currents of the inlet being modeled. Once the
model is validated, any inlet improvement measure, such as
jetties, groins, or rip-rap fill, can be confidently
simulated.
6.1 Field Data
The Coastal and Oceanographic Engineering Department of
the University of Florida was commissioned by the City of
Punta Gorda, FL to examine the potential for shoreline erosion
due to the navigational channel at Ponce de Leon Inlet in the
vicinity of Punta Gorda City Park (Wang, et al., 1996). Ponce
de Leon Inlet is located on the east side of Charlotte Harbor
50
at the mouth of the Peace River. The channel serves as the
public boat ramp at Punta Gorda City Park as well as a
residential canal system. The data which was collected and
presented in a report to the City of Punta Gorda will be used
to show the validation and application (chapter 7) of this
numerical model.
Data were collected from Ponce de Leon Inlet in November
and December 1995 and March 1996. The data necessary to
validate the numerical model are: a bathymetric survey; a tide
record; current measurements; tidal flow pattern studies; and
wind speed and direction measurements. However, other field
data may be necessary for coastal structure design purposes.
6.1.1 Bathymetry
The bathymetric data collected should be complete enough
such that bottom contours can be drawn accurately. If
available, navigational charts may be used in lieu of, or to
supplement the field data. No navigational charts were
available for this inlet so a bathymetric field survey was
conducted. The bottom contours obtained from this data are
shown in figure 6.1. The depth for each cell in the numerical
grid can be interpolated from these contours.
1 -.51 0) Contour in m (NGVD)
Figure 6.1
Bathymetric survey of Ponce de Leon Inlet, figure
taken from Wang, et al. (1996).
6.1.2 Tide
Tide gages were placed at channel markers 8 and 18
(figure 6.2) to measure tide height for the period of data
collection, these data are shown in figure 6.3. The time lag
of the tide between channel markers 8 and 18 is on the order
of 5 to 10 minutes. The maximum difference in amplitude
between the two tidal records is approximately 1.5 cm, though
most of the tidal peaks show no difference in amplitude.
Figure 6.2
Location of stations where tide and current data
were collected, figure taken from Wang, et al.
(1996).
0.5
0.4
0.3
A '
SOddmeter data
-0.2-
0 I I I I
03 2
Figure 6.3
Continuously observed- Tide at Channel Marker 8
- Tide at Channel Marker 18
-0.4current measurements, Current at inlet entrance
al. Oddmeter data
December 1995
Figure 6.3
Continuously observed tidal record and discrete
current measurements, figure taken from Wang, et
al. (1996) .
6.1.3 Current
Two instruments were used to collect current data: an
electro-magnetic meter, and a vane-type meter. The electro-
magnetic meter manufactured by Marsh-McBurney, was installed
on a docking pile just inside the inlet (station 1 in figure
6.2). The data collected from this meter are also given in
figure 6.3. Comparing the current data to the tide data, it
may be concluded that the maximum currents occur when the tide
is near zero. It is also easily seen that there is slack
water when the tide changes either from flood to ebb or ebb to
flood.
Other current data were collected from an Ottmeter, a
hand-held device which uses vanes driven by the flow to
measure the current. One minute averages were taken at about
0.5 m below the surface and 0.5 m above the bottom. Figure
6.4 shows the data collected during the flood tide and figure
6.5 shows the ebb velocity data measured at station 1, 2, and
3 (see figure 6.2) in the inlet.
Figure 6.3 shows the Ottmeter data taken at the same
location as the electro-magnetic meter. The flood velocities
between the two different methods are in agreement, though
there is a large discrepancy at the ebb tide. This
discrepancy can be attributed to the improper installation of
miT e=16:43
00.41
00.21
0.s 9 0.28
-0. s cu1ren)95
Flood current(m/s)
Station 2
0.42
e 0.51
0 0.54
0.e0.36
0039 4.36
Station 3
M 0.74 0 0.50 0.26
0- 0.54
00.59
Time=1 6:55
Time=17:09
Horizontal Distance (m)
Figure 6.4
Flood velocity data taken inside the canal using
the Ottmeter, figure taken from Wang, et al.
(1996).
Station I
V*
-2-
e0.44
0
t-21
LU
0.
-2-
S t I... ...... .. .
a 0.25
8 0.32
Time=1 4:07
0 0.31
^ 0022 a0.26 -'
-S 0 -0.27 -- 1/30/95
Ebb current(m/s)
Station 2 Time=14:34
0 0.27
I-Z
E
.2
w
w>
0 0.30 a 0.30
o 0.23 1-
Time=1 5:02
Horizontal Distance (m)
Figure 6.5
Ebb velocity data taken in the canal using the
Ottmeter, figure taken from Wang, et al. (1996).
tantinn I
02
-2 -
Station 3
YIUIIV
,
57
the electro-magnetic current meter. The meter should have
been installed in deeper water, instead the meter was
installed such that it was less than 0.2 m of water during low
tide therefore giving the low, inaccurate measurements. Also,
the current meter was about 0.75 m from the pile, there could
have been some adverse affects on the current caused by the
pile.
6.1.4 Drogues
The ebb and flood current patterns were studied using
drogues. As the drogues were placed in the water, their
motions were video taped. From a video tape the velocities
were estimated from using estimated distances and time elapsed
on the video tape. When extracting the distances the drogues
traveled from relation to channel markers of a known distance
a problem is introduced with the perspective of the video
camera, because the video camera was at fixed, stationary
point on land. This is a crude method to obtain the current
velocities; therefore, the values may vary greatly from those
obtained through more non-Lagrangian methods. As a result,
the quantitative velocity values should only be used as a
relative reference between individual drogue path lines. The
main purpose of the drogues is to qualitatively study the flow
0 50m
Ponce De Leo
n Inlet
O.5m/sec
CHARLOTTE
HARBOR
Figure 6.6
Flood current patterns estimated using drogues,
Wang, et al. (1996).
IL
Figure 6.7
Ebb flow patterns studied using drogues, Wang, et
al. (1996).
60
patterns. Figures 6.6 and 6.7 show the flow patterns from the
drogue study for the flood and ebb tides. This is one of the
more important field studies conducted for validation of the
numerical model.
6.1.5 Wind
The days that the tide and current measurements were
taken, the atmospheric conditions were ideal, there was a high
pressure system, a slight breeze, and no synoptic events had
occurred in the previous days or were predicted. Hence, wind
effects were negligible and not included for validating the
numerical model.
6.2 Spectral Analysis
Spectral analysis is performed on the real tide, so that
the predominate frequencies and amplitudes can be used to
simulate the real tide using Fourier series (eq. 4.1). The
tide record observed at Punta Gorda (figure 6.3) shows that
the tide is mixed diurnal and semi-diurnal. Moreover, if the
record was taken for 30 days, a spring and neap tide should be
evident. For the model validation, spectral analysis will be
S0
-0.05
-0.1
-0.15
-0.2
61
Observed Tidal Cycle
Hours
Figure 6.8
One tidal cycle extracted from the tide data.
performed on only one period of the recorded tide in figure
6.8. The rational for selecting only one tidal period is:
(1) excluding the scattered current measurements at
the beginning of he current record (figure
6.3), the current data reflects only one tidal
cycle; and
(2) since the model is being used as only an aid to
inlet morphological improvement measures,
62
residual currents and tidal flushing are not
of interest.
The continuous tide data taken were digitized to 1 hour
increments for a cycle during the time that current
measurements were being recorded. Twenty four data points
were taken from the observed data record. Figure 6.9 shows
the resulting spectrum of the tidal cycle.
6.3 Grid
The grid chosen for the Ponce de Leon Inlet as shown in
figure 6.10. The grid size Ax and Ay, is 6 m x 6 m, which
results in 59 rows and 29 columns. This resolution will be
sufficient to allow jetties to be inserted and show the
currents close to the shore. A higher resolution grid could
be used, but the computation time required will increase. The
information gained by using a higher resolution grid does not
justify the excess time needed to run the numerical model.
Both the northern and southern boundaries were extended,
so that the eddies which were studied using the drogues
(figure 6.7) do not cause a numerical instability at the open
boundaries. The drogue studies also show that the tidal flow
!2
E
150
100-
50-
0
0 50 100 150
meters
Figure 6.10
Area of Ponce de Leon Inlet
which a grid will be applied
to.
becomes one-dimensional in the canal; henceforth, the grid
chosen should allow the model to be validated against the data
collected in the field studies.
Using the stability criteria discussed in chapter 3.4,
the time step and eddy viscosity coefficient can be
determined. The Courant number used was 10, which resulted in
64
a time step of 18 sec and an eddy viscosity coefficient of
0.12. The Darcy-Weisbach friction coefficient was set to
0.01, and, again, wind was not used as a parameter for this
model, since field studies were conducted on a calm day.
6.4 Open Boundary Conditions
The open boundary conditions (OBCs) needed to drive the
model for the case of Punta Gorda are tidal, Neumann, and
canal open boundaries (OBs). The tidal OBs will be placed on
the western edge of the grid, the Neumann OBs will be placed
on the northern and southern edges of the grid which extend
into Charlotte Harbor, and canal OBs will be placed on the
cells on the eastern edge of the grid. The following sections
will discuss the parameters used for these OBCs.
6.4.1 Tidal OBC
The amplitudes, circular frequencies, and phase shifts
for the most dominant harmonics, determined by spectral
analysis, are shown in table 6.1. These are the parameters
which will be used to drive the numerical model at the open
boundary using equation 4.1. The simulated tidal signal is
n-th Amplitude, Ai Circular Phase Angle, 4
Harmonic (m) Frequency, a0 (rad)
__(rad/sec)
1 0.0847 7.2722 x 10-5 0.8178
2 0.1654 1.4544 x 10-4 2.7241
3 0.0119 2.1817 x 10-4 0.9671
6 0.0079 4.3633 x 10-4 2.5830
Table 6.1
Parameters used for the tidal OBC.
compared against the measured signal from figure 6.8 in figure
6.11.
6.4.2 Neumann OBC
Since the grid extends into Charlotte Harbor where the
tidal OBC is, the appropriate measures discussed in chapter
4.1 must be applied here. In comparing the area the numerical
model's grid will cover (figure 6.10) to the bathymetric data
(figure 6.1), it is shown that the boundaries were extended
roughly 60 m to the north and 60 m to the south. This adds a
buffer zone so that any circulation in the area of interest
does not cause an instability in the OBC and so that any
effects caused by the Neumann OBC do not adversely affect the
flow patterns in the area of interest. During the model
validation or its application, if the eddies occur near the
Simulated Tidal Cycle
Hours
Figure 6.11
Simulated tidal OBC signal compared to the
observed tidal signal in fig. 6.8.
Neumann OB cells, the model's grid must be extended farther to
the north and south.
6.4.3 Canal OBC
The parameters necessary for this OBC are depth at the
end of the canal, the Darcy-Weisbach friction coefficient, the
size of the increment (Axcana), and the number of increments
(Icanal) These parameters are somewhat empirical. The
67
friction coefficient may be greater in the canal than in the
inlet due to the different bottom sediments characteristics,
excess submerged vegetation (e.g. sea grass), or overgrown
foliage on the banks (e.g. shrubs and mangroves). In the
interest of computational time, as previously discussed, the
total number of increments should not exceed 15. The measured
or estimated length of the canal can be divided evenly into a
number of increments. The right side of figure 6.2 shows the
beginnings of a complex finger canal system. This may add
damping, or friction, to the long wave, so the friction
coefficient should be increased. The finger canal system also
makes it very difficult to obtain an "effective" length if the
canal were straight.
For the Punta Gorda Isles finger canal system, the
surface area of the system was estimated from maps. Since
there are five canal OBCs used on the numerical model and 10
increments were used for each OBC, the surface area was
divided by 50 to obtain the size of increment, Axcanai, which
was 1500 km. The friction factor used in the numerical model
was 0.05. The canal system was assumed have uniform depth,
hence the depth at the end of the canal was set to the depth
at the OBC, or beginning of the canal.
- - - -
y) UU - - - - - - - -
E 0 . '--- -- -- - - -. ' - .' ' ---
E 80
40-
0 0.3 mlsec
20 . . . . . . .
-20 0 20 40 60 80 100 120 140 160 180
meters
Figure 6.12
Maximum flood velocity vectors for numerical
model validation.
6.5 Results
The numerical model was ran for 5 tidal cycles, a
sufficient time for the model to reach steady state from the
initial conditions.
Figure 6.12 shows the maximum flood
velocities and figure 6.13 shows that maximum ebb velocities.
The maximum velocities, both flood and ebb, are about 0.3 m/s,
this value is a little low compared to the electro-magnetic
current meter readings of about 0.4 m/s and the estimated
100 - - - - - -
......................
40- .
20
-20 0 20 40 60 80 100 120
meters
Figure 6.13
Maximum ebb velocity vectors showing
model validation.
drogue velocities of 0.5 m/s.
140 160 180
numerical
The difference in velocities
may be as a result of the canal OBC. One could adjust the
parameters used in the canal OBC to produce better results,
though that could take many iterations.
Comparing figures 6.13 and 6.7, the numerical model does
not show the circulation patterns that the drogue studies
show. The eddy viscosity coefficient was increased to the
point that the stability criteria (eq 3.44) was slightly
violated, though no instabilities became apparent. A possible
70
causes for this difference are that the bathymetric data may
have errors. Another possible cause is that the bathymetric
relief near the channel may be too great for the model's
resolution. This will be addressed later in chapter 8.
Though a more likely reason that the numerical model was not
able to produce the circulation is because the circulation may
be initiated by superimposing the wave conditions on the ebb
tidal flow. The predominant waves were found to come from the
north which could cause a circulation pattern on the northern
side of the channel during the ebb cycle. Also, the Peace
River empties into Charlotte Harbor just north of Punta Gorda.
The current from the river, a north to south current, or an
ebb current throughout Charlotte Harbor can initiate a
circulation pattern which was observed.
CHAPTER 7
MODEL APPLICATION
Recall, the objective of this research is to develop a
numerical model to aid in the placement of coastal structures
as an inlet improvement. When applying the numerical model,
any differences or inconsistencies between the field data and
the model validation results must be kept in mind. The
effects of these discrepancies on the numerical output must be
realized when scrutinizing the inlet improvement measures,
such as where a jetty or groin is placed and its length.
Ponce de Leon inlet in Punta Gorda, FL will be used to
show how the numerical model can be used. Two cases will be
demonstrated where a single jetty is installed and a pair of
jetties are installed. The results presented in this section
are a "zoomed in" portion of the model's grid. This allows
the currents in the area of interest to be shown in better
detail.
-20 0 20 40 60 80 100 120 140 160 180
meters
Figure 7.1
Flood velocity vectors of a single jetty applied
to Punta Gorda, FL.
7.1 Single Jetty
Figures 7.1 and 7.2 show the maximum flood and ebb
velocities, respectively, for a single jetty case.
From
comparing figures 7.1 and 6.12, the ebb currents become
greater along the southern bank, the bank with no jetty. The
jetty also causes return flow during the ebb current on the
southern shoreline, though the velocities may not be large
enough to initiate shoreline erosion.
The improvement
measures recommended to the City of Punta Gorda were to instal
100 -
E80 -
60-
40-
20
0--'-----'--
-20 0
Figure 7.2
Ebb velocity
Punta Gorda,
vectors of a single jetty applied to
FL.
the single jetty as shown and to add rip-rap fill to the
southern shoreline, Wang et al. (1996).
7.2 Pair of Jetties
Figures 7.3 and 7.4 show the maximum flood and ebb velocities,
respectively, for the case with a pair of jetties. The pair
of jetties virtually eliminates the erosional currents in the
proximity of the shoreline. Again, the discrepancies between
-20 0 20 40 60 80 100 120 140
meters
Figure 7.3
Flood velocity vectors of a pair
applied to Punta Gorda, FL.
of jetties
the field data and the validation results, namely the drogue
study and circulations eddies, must be kept in mind.
It is
possible that the pair of jetties may cause significant eddies
outside of the channel near the shoreline which may result in
higher, erosional velocities.
ion
160
140
100-
. . . . . . .. 0.3 m/sec
II I I I I I
160 180
ioo -- -- i-- i --i --- i --- '--- -- i i -- -
160. . . . . . . .
140. .
E 80- ...... .... . .
6 0 . . . . . . . . . .
0. . . . .
0- 0.3 m/sec
60
4o0 : - iii. -
-20 0 20 40 60 80 100 120 140 160 180
meters
Figure 7.4
Ebb velocity vectors of a pair of jetties applied
to Punta Gorda, FL.
CHAPTER 8
CONCLUSIONS AND RECOMMENDATIONS
A two-dimensional numerical model was developed using a
factorization method for numerical solution. The model
developed accounts for constant friction, convective
acceleration, turbulent diffusion, and wind stress. The
numerical model can be adopted to a variety of inlet-river and
-canal systems to study the tidal currents. This model can be
applied to coastal engineering design projects involving the
installation of jetties, breakwaters, groins, or rip-rap fill
where it is necessary to study the feasibility of such inlet
improvement measures.
8.1 Conclusions
The numerical model's efficacy was verified reasonably
well. There is a discrepancy when comparing the numerical
model to the analytical solutions of Massachusetts Bay. A
possible cause of the difference could be the friction factor.
The analytical solutions did not account for friction, where
as if friction is not included in the numerical model, then
77
the model may become unstable or produce erroneous results.
The analytical solutions also neglected the inertial terms.
Perhaps, if the convective terms were set to zero in the
numerical model, then the results might change.
Overall, the numerical model that was developed can be
applied to a multitude of tidal inlets and areas. The model
was shown that it could produce reasonable results. The
numerical model developed was not able to produce ebb tidal
patterns as measured in the field when applied to Ponce de
Leon inlet at Punta Gorda, FL. Possible reasons for the
differences between the numerical model and the field data
could be bathymmetry and the canal OBC.
Most navigable inlets consist of a deep channel dredged
across a relatively shallow area. When numerically modeling
these inlets, the model may produce better results if the
channel is "numerically smoothed". Figure 8.1a shows the
channel at Punta Gorda and figure 8.1b how one might "smooth"
the relief of the bathymmetry. The sudden relief in
bathymmetry may cause a large gradient in the friction and
convective terms with profound affects on the numerical model.
78
a
0
-0.5
-2
40 -30 -20 -10 0 10 20 30 40
b
0
-0.5
S-2
-1.5
-40 -30 -20 -10 0 10 20 30 40
Distance (m)
Figures 8.1 a & b
(a) shows the original bathymmetry of channel and
(b) shows the "smoothed" bathymmetry of channel.
8.2 Recommendations
The following studies and improvements to the numerical
model are recommended:
1)The affects of "numerically smoothing" a channel
should be closely studied. From testing the
numerical model as it was being developed, it is
reasonable that a "smoothed" channel will produce
79
more realistic results. Information was not
presented in this thesis on the grounds that no
analytical solutions or data were available at
the time to support this hypothesis.
2)Develop anther OBC which allows for a lagoon
connected by a small canal as shown in figure
8.2. The governing equation presented by Bruun
(1978), but developed by earlier researchers, for
this case are as follows:
d'lB F AB dfB dflB gA gA
t2 LcAC dt dt LA B (8.1)
where
lB=surface elevation of bay,
t=time,
F=frictional losses,
Lc=length of channel,
AB=surface area of bay,
Ac=cross sectional area of channel,
g=gravity, and
io=surface elevation at entrance.
a AREA-A-, -. BAY FRESH
EPTH WATER
OCEAN \ OLUME v0 INFLOW
ELEVATION SURFACE AREA AS O
E oN ELEVATION
Figure 8.2
An idealized channel/bay system. Figure was
taken from Bruun (1978).
This lagoon/canal system can be applied to the
numerical model as an OBC the same as the canal
and river are discussed in chapter 4.1, except
that the ocean elevation shown in figure 8.2 can
be the surface elevation can be used to drive the
model as at an OB.
3)Morphological changes and sediment transport
should be added to the numerical model. This
will allow the long term sediment transport,
shoaling, and bank erosion/accretion to be
evaluated for a given inlet improvement measure.
LIST OF REFERENCES
Arcilla, A.S. "An Integrated Numerical Approach for Coastal
Engineering Problems." Journal of Coastal Research.
Vol. 5, 1989: 603-616.
Bathan, Karl H. "A Tidal and Wind Driven Numerical Model of
a Lagoon Before and After a Major Modification." J.K.K.
Look Laboratory of Oceanographic Engineering,
University of Hawaii, Vol. 6:2, 1976.
Briggs, Douglas A. and Madsen, Ole S. "Analytical Models for
One- and Two- Layer Systems in Rectangular Basins."
Report No. 172, School of Engineering, Massachusetts
Institute of Technology, Cambridge, MA, 1973.
Bruun, Per. Stability of Tidal Inlets. Elsevier Scientific
Publishing Complany, Amsterdam, 1978.
Cheng, Ralph T. "Modeling of Tidal Residual Circulation in
San Francisco Bay, California." Two-Dimensional Flow
Modeling. U.S. Army Corps of Engineers, The
Hydrographic Engineering Center, Davis, California,
1982.
Choi, Jei-Kook. "Three-Dimensional Curvilinear-Grid Modeling
of Baroclinic Circulation and Mixing in a Partially
Mixed Estuary." Dissertation, University of Florida,
Gainesville, FL, 1992.
Henderson, F. M. Open Channel Flow. Macmillan Publishing
Co., Inc., New York, 1966.
Hirsch, Charles. Numerical Computation of Internal and
External Flows. Volume 1, John Wiley & Sons Ltd.,
Chichester, 1994.
Pedlosky, Joseph. Geophysical Fluid Dynamics. Springer-
Verlag, New York, 1979.
Reeve, D. E. and Hiley, R. A. "Numerical Prediction of Tidal
Flow in Shallow Water." Computer Modelling of Seas and
Coastal Regions. Elsevier Science Publishers Ltd., New
York, 1992: 197-208.
Reidel, H. P. and Wilkinson, F. L. "Numerical Modeling An
Aid to Assessing Field Data." Coastal Engineering.
1976: 3243-3262.
Roll, H. U. Physics of the Marine Atmosphere. Academic
Press, New York, 1965.
Schomaker, Christine W. "A Model for Tidal Circulation
Adapted to Monterey Bay, California." Thesis, Naval
Postgraduate School, Monterey, California, 1983.
Sheng, Y. P. Mathematical Modeling of Three-Dimensional
Coastal Currents and Sediment Dispersion: Model
Development and Application. Technical Report CERC-83-
2, Aeronautical Research Associates of Princeton,
Princeton, New Jersey, 1983.
Wang, Hsiang; Lin, Lihwa; and Wang Xu. "Ponce de Leon
Channel Improvement Studies, Punta Gorda, Florida."
UFL/COEL-96/004, Coastal and Oceanographic Engineering
Department, University of Florida, Gainesville,
Florida, 1996.
Weiyan, Tan. Shallow Water Hydrodynamics. Elsevier Science
Publishing Company, Inc., New York, 1992.
Wilson, Basil W. "Note on Surface Wind Stress over Water at
Low and High Wind Speeds." Journal of Geophysical
Research. Vol. 65, No. 10, 1960: 3377-3382.
BIOGRAPHICAL SKETCH
Peter N. Seidle was born in Binghamton, NY on 28
January, 1972. His family moved to Maryland in suburban
Washington DC where he started kindergarten and graduated
from Thomas S. Wootton High School, Rockville, MD, in 1990.
In 1994 he graduated from Widener University, Chester, PA,
with a Bachelor of Science in Mechanical Engineering.
Shortly after graduating from college, he became employed by
the ARIST Corporation to work as a program analyst in the
Ballistic Missile Defense Organization, which is part of the
Department of Defense. In August 1995 he began his graduate
studies at the University of Florida in the Coastal and
Oceanographic Engineering Department. In May 1997 he
graduated with a Master of Science in the field of study he
has and will always enjoy and be fascinated with.
|
RSS
HIDE
