UFL/COEL89/023
A TWODIMENSIONAL HYDRODYNAMIC MODEL
USING A FINITEVOLUME APPROACH
by
Joaquim Jose Areias Capitio
Thesis
1989
A TWODIMENSIONAL HYDRODYNAMIC MODEL USING A
FINITEVOLUME APPROACH
By
JOAQUIM JOSE AREAS CAPITJAO
A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF MASTER OF SCIENCE
UNIVERSITY OF FLORIDA
1989
ACKNOWLEDGEMENTS
I would like to express my gratitude to Dr. Peter Sheng, my advisor and com
mittee chairman, for his guidance and support during the time I was studying at the
University of Florida, and to Dr. Robert Dean and Dr. Max Sheppard for serving
on my committee and for their comments about this work.
The support of all my fellow students and the several postdoctoral associates
from whom I have learned so much during my stay in Gainesville was also deeply
appreciated. Dr. PeiFang Wang was particularly helpful, not only for his effort in
providing me all the elements I needed during this work but also for his unwavering
support at the most difficult times. The help of Mr. Yuming Liu and Mr. Hyekeun
Lee was also essential to the development of this work. I would also like to thank
Mr. Subarna Malakar for his support when dealing with the strongwilled computers
used in this work.
I must thank Dr. Ant6nio Melo Baptista for all the great advice and timely
encouragement I have received from him in the last seven and onehalf years. His
friendship is greatly appreciated.
I would also like to express my appreciation for the support from the Luso
American Educational Commission in Lisbon (administering the Fulbright Program
in Portugal) and from the Institute of International Education (offices in New York,
Atlanta and Houston) who, as agents of the Fulbright Program, promptly gave me
all the support I needed before and during my stay in the United States.
My gratitude must also be extended to my employer in Portugal, Laborat6rio
Nacional de Engenharia Civil, for allowing me to come to the United States for an
academic program.
Last but not least, I would like to thank my family and friends in Portugal and
in the United States for the support I received from them, without which none of
this would have been possible.
TABLE OF CONTENTS
ACKNOWLEDGEMENTS .................... .......... ii
LIST OF FIGURES ..... .... .. .................... vi
ABSTRACT ................ ....... ............. vii
CHAPTERS
1 INTRODUCTION ...... ........... ...... ........ 1
2 BRIEF REVIEW ................ ................ 4
3 THE NUMERICAL MODEL ...... ............ ....... 7
3.1 Governing Equations .... ........... ..... ....... 7
3.2 The Computational Grid ..... ........... ........ 9
3.3 The Fractional Step Method .................. .... 13
3.4 The Advection Step ........................... 15
3.5 The Diffusion Step .... ...... .................... 19
3.6 The Coriolis Step .... .. .......... ............. 25
3.7 The Propagation Step ........................... 26
3.8 The Boundary Conditions ... ............ ......... 37
3.9 Numerical Stability ............................ 38
4 APPLICATIONS ................. ................. 40
4.1 Comparison with Analytical Solution ................. 40
4.2 Square Basin with Constant Slope .................. 41
4.3 Square Basin with VShaped Bottom ........ ... .... 53
4.4 Lake Okeechobee with Constant Wind ................ 53
4.5 Lake Okeechobee with Sinusoidal Wind ............... 75
5 CONCLUSIONS ............ .... .............. 79
BIBLIOGRAPHY ................................. 81
BIOGRAPHICAL SKETCH ................. .......... .83
LIST OF FIGURES
3.1 Physical domain vs. computational domain . . . 9
3.2 Contravariant and covariant directions . . .... 10
3.3 Basic computational cell ...................... 13
4.1 Skewed grid for a square basin with uniform depth ...... ..42
4.2 Comparison with analytical solution . . . ... 43
4.3 Skewed grid for a square basin with uniform bottom slope 44
4.4 Velocity results with constant slope present model ...... ..45
4.5 Velocity results with constant slope CH3D . . ... 46
4.6 Surface elevation with constant slope present model . 47
4.7 Surface elevation with constant slope CH3D . ... 48
4.8 Time series results at points A and B present model . 49
4.9 Time series results at points C and D present model . 50
4.10 Time series results at points A and B CH3D . ... 51
4.11 Time series results at points C and D CH3D . ... 52
4.12 Skewed grid for a square basin with Vshaped bottom . 54
4.13 Velocity results with Vshaped bottom present model . 55
4.14 Velocity results with Vshaped bottom CH3D . ... 56
4.15 Surface elevation with Vshaped bottom present model .. 57
4.16 Surface elevation with Vshaped bottom CH3D . . 58
4.17 Time series results at points A and B present model . 59
4.18 Time series results at points C and D present model . 60
4.19
4.20
4.21
4.22
4.23
4.24
4.25
4.26
4.27
4.28
4.29
4.30
4.31
4.32
4.33
4.34
4.35
Time series results at points A and B CH3D
Time series results at points C and D CH3D
Curvilinear grid for Lake Okeechobee . .
Bottom contours for Lake Okeechobee . .
Velocity results for Lake Okeechobee present
Velocity results for Lake Okeechobee CH3D .
Surface elevation for Lake Okeechobee present
Surface elevation for Lake Okeechobee CH3D .
Time series results at points A
Time series results at points C
Time series results at points E
Time series results at points A
Time series results at points C
Time series results at points E
Longterm simulation results 
Longterm simulation results 
Longterm simulation results 
and B present model . .
and D present model . .
and F present model . .
and B CH3D . . .
and D CH3D . . .
and F CH3D . . .
points A and B . . .
points C and D . . .
points E and F . . .
model .
t model .
S61
S 62
S 63
. 64
. 65
S66
. 67
S68
. 69
S 70
S 71
. 72
S 73
S74
S76
S77
S78
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 TWODIMENSIONAL HYDRODYNAMIC MODEL USING A
FINITEVOLUME APPROACH
By
JOAQUIM JOSE AREAS CAPITAL
December 1989
Chairman: Dr. Y. Peter Sheng
Major Department: Coastal and Oceanographic Engineering
A twodimensional model of winddriven circulation in a closed basin was developed
using a finitevolume technique for generalized curvilinear grids and applying some
of the recent developments in hydrodynamic modeling. The terms in the Navier
Stokes equations were treated separately according to the fractional step method
and the propagation step, including the continuity equation and the pressure and
stress terms in the momentum equations, was solved using a conjugate gradient
method.
The model was then applied to a number of test cases to examine the feasibility
of the approach used by comparing with results obtained with the twodimensional
version of the threedimensional model CH3D. These included a square basin with
constant slope and with a Vshaped bottom and Lake Okeechobee, in South Florida.
To evaluate the longterm numerical stability of the model, a tenday model simu
lation with varying wind was also run for Lake Okeechobee.
CHAPTER 1
INTRODUCTION
Water is as essential to any kind of human activity as the air we breathe. Ever
since the first human settlements, their location has been dependent on the presence
of water. Big cities have naturally developed around important rivers or lakes, or by
the ocean, where water for drinking, irrigation and, since the industrial revolution,
cooling of industrial plants, is easily accessible.
Unfortunately, water bodies have been considered as infinite sources of fresh
water and, at the same time, infinite dump sites for urban, agricultural and in
dustrial waste. Only in the last few years have we started to become aware of the
limits to this source and the effect of using it as a dump site. Entire seas, like
the Mediterranean, have reached dangerous levels of pollution. Rivers and lakes all
over the world are so polluted, chemically and bacteriologically, that their water is
totally unusable without previous expensive treatment.
All this has happened over a period of time and the human mind seems to
be able to adapt to new conditions without fully realizing their implications. It
took, therefore, wellpublicized accidents like oil spills or garbage being washed
ashore in popular beaches to make the general public aware of the seriousness of
the problem. That public awareness has put pressure on the political establishment
to find solutions and this has led to increased pressure on the scientific community
to develop the knowledge of the physical, chemical and biological processes involved.
One of the basic subjects involved is the knowledge of hydrodynamic circulation
in shallow water, where the problem is more pressing. The research in this field
has traditionally used three different tools to reach the same objective, a better
2
knowledge of circulation patterns in any given area:
1. field measurements;
2. scale models; and
3. numerical models.
Field measurements are in a class by themselves and even the strongest sup
porters of scale or numerical models will accept the need for some amount of field
data to support their modeling effort. The continuous development of new instru
mentation and the increase in computer speed when processing the data has made
field measurements an important part of any study concerned with the knowledge
of circulation patterns in any shallow water area. However, field measurements are
still expensive and must be used sparingly.
Huge scale models for whole estuaries were common in the past, and some
are still being used. However, the cost of keeping and operating such a model has
increased steadily, while, at the same time, the use of ever more powerful computers
has made numerical models cheaper and cheaper. Therefore, scale models seem to
be, at present, more adequate for studying specific physical processes, instead of
general circulation patterns for vast areas.
The evolution of computers, in terms of speed and memory, has made numerical
modeling relatively cheap. But numerical models are not without problems of their
own. Although the NavierStokes equations are accepted as a good representation
of the physics involved, for instantaneous, threedimensional circulation, from a
practical point of view they must be integrated over a finite period of time, giving
rise to problems like turbulence closure.
This thesis presents a new numerical model of hydrodynamic circulation, using
a finite volume technique to solve the twodimensional, vertically integrated, Navier
3
Stokes equations in a curvilinear grid, and using a fractional step method to solve
the different terms in the equations separately.
Chapter 2 consists of a brief review of work done in the past on numerical
modeling of hydrodynamics, mainly those models which, in some way, influenced
the way the present model was developed.
In Chapter 3 the finite volume equations are derived from the twodimensional
vertically averaged NavierStokes equations.
In Chapter 4 a number of test cases are presented and Chapter 5 contains the
conclusions of the present work and some ideas about possible future developments
on the present numerical model.
CHAPTER 2
BRIEF REVIEW
A large number of numerical models of hydrodynamic circulation have been
developed over the last few decades, and most of them use finite difference or finite
element methods in the solution of the NavierStokes equations. The finiteelement
methods have traditionally been considered to have a better ability to deal with
complex geometries, since they can use nonrectangular grid cells. The need for
rectangular grid cells was, therefore, the main limitation of finitedifference methods
like the one developed by Leendertse (1967), which, on the other hand, have always
been much more intuitive in their development and, therefore, easier to modify.
Although the developers of finite difference models were not ready to move into
finite element modeling, they would recognize the limitations associated with their
own technique, which was a first step into overcoming those limitations. One ap
proach was the coupling of models using a sparse grid far from the boundary and a
much denser grid near the boundary, as in Sheng (1976). Another possible approach
utilizes the grid generation techniques presented in Thompson et al. (1985) to re
solve the physical domain with a boundaryfitted curvilinear grid. Spaulding (1984)
developed the equations of motion in generalized curvilinear coordinates, in terms
of the cartesian components of the dependent variables (the velocity components),
while Sheng (1986) derived the equations of motion in terms of the contravariant
components of the dependent variables.
The approach used in this thesis is the same followed by Sheng (1986) and Sheng
et al. (1988) which consists of transforming not only the independent variables but
also the dependent variables and, therefore, working with equations in terms of
contravariant components of velocity.
Another problem associated with finite difference methods is the nonconserva
tive nature of the discretized equations in curvilinear grids due to the use of differ
ential equations. This problem is usually taken care of by treating the geometric
terms more precisely. An alternative to this is to develop the numerical model equa
tions starting with the integral equations instead of the differential equations. This
is done in finitevolume techniques, and a comparison between the two approaches
is presented by Vinokur (1986). A more detailed explanation of a model using a
finite volume method is presented in Rosenfeld et al. (1988).
When numerical schemes are developed to solve complicated equations, special
care must be taken to guarantee numerical stability and consistency. This usually
involves an option between more elaborate numerical schemes and very small spatial
and temporal resolutions. The fractional step method, presented in Yanenko (1971)
is a convenient way to get around these limitations by breaking each time step
into a series of intermediate steps, with a number of terms in the equations being
solved at each time step. Each intermediate step can then be solved using the more
convenient numerical scheme for the specific terms involved, and having only the
limitations in terms of temporal and spatial resolution imposed by those terms. In
this model, the fractional step method is used, breaking the NavierStokes equations
into four separate equations, each one to be solved at a different intermediate step.
In the advection step, only the nonlinear terms in the momentum equations are
used. The diffusion step solves for the horizontal turbulent diffusion terms and the
Coriolis terms are solved in the Coriolis step. Finally, a propagation step solves
the remaining terms in the momentum equations, pressure terms, wind stress and
bottom friction, together with the continuity equation.
6
The propagation step was solved using a conjugate gradient procedure detailed
in Hauguel (1979) and used before, among others, in the cartesian grid, finite dif
ference circulation model developed by Liu (1988).
CHAPTER 3
THE NUMERICAL MODEL
3.1 Governing Equations
The equations describing the twodimensional, verticallyintegrated flow in shal
low water can be obtained from the threedimensional NavierStokes equations, as
suming a hydrostatic vertical pressure distribution and, in this case, a constant and
uniform density. The turbulence closure problem was solved using a constant and
uniform eddy viscosity coefficient, and the final twodimensional equations are, in
a cartesian coordinate system,
of auc aVc
at = x y (3.1)
at ax ay
aUc a UcUc\ UcVc\ 2Uc aUc
t ax U )  dy A X) +AH a +AH
at 9x H y H 8 ay2
aB 1
+fVc gH + (r,w T,) (3.2)
ax p
avc a U (Y A a2vc a 2v
at = C C y + AHX2 a2 + An a2V
at H H A H ay2
of 1
fUc gH + (r, Tbu) (3.3)
By p
where t is time, x and y are the cartesian spatial coordinates, Uc(x,y,t) and
Vc (x, y, t) are the cartesian components of the verticallyintegrated velocity in the
x and y directions, H(x,y,t) is the total water depth, AH is the eddy viscosity
coefficient, f is the Coriolis parameter, g is the acceleration of gravity, p is the
8
water density, r,(x, y,t) and rv(x,y,t) are the x and y components of the surface
wind stress and rb,(x,y, t) and (x, y, t) are the x and y components of the bottom
stress.
The bottom stress was modeled, as in Liu (1988) by
Tbz = CU +V(3.4)
C2H2
pgVc +V (35)
Tv = C2H2 (3.5)
where C is the Chezy bottom friction coefficient given, in the C.G.S. system
(cm /sec) by
C = 8.21 (3.6)
n
and n is the Manning coefficient. The hydraulic radius, R is given by the ratio
between the water crosssection and the wetted perimeter,
Hb
R H b (3.7)
2H + b
In lakes and estuaries, the width of the basin is usually much larger than its
depth, and the hydraulic radius is approximately equal to the depth H.
X
x 5
Figure 3.1: Physical domain vs. computational domain
3.2 The Computational Grid
Because of the extreme difficulty in generating an orthogonal curvilinear grid,
the equations above have to be solved in a generic curvilinear grid, not necessarily
orthogonal (figure 3.1). For that purpose, a different coordinate system, (C, r7) is
defined for the computational domain such that each curvilinear cell in the phys
ical grid is mapped into a rectangular cell in the computational grid. For sake of
simplicity, all sides of the computational cell have length AC = A77 = 1.
Changing from the (x, y) cartesian coordinate system to the (e, 17) curvilinear
coordinate system involves defining a number of geometric quantities. First of all, we
need to define contravariant and covariant directions (figure 3.2). A contravariant
base vector is perpendicular to a line along which its coordinate remains unchanged.
A covariant vector is tangent to a line along which only the other coordinate changes.
If the curvilinear grid is orthogonal, the two coordinate axis will be perpendic
ular to each other at every point and, therefore, the contravariant and covariant
directions coincide, since a vector perpendicular to one of the coordinate axis is also
tangent to the other coordinate axis.
Figure 3.2: Contravariant and covariant directions
Contravariant length vectors can be defined at the center of each side of a cell,
parallel to the contravariant base vectors defined above, with magnitude equal to
the length of the side:
L = y,7 x, (3.8)
Lf = yJ+ z~f (3.9)
where i and f are the base vectors in the cartesian coordinate system and the
subscripts denote partial derivatives.
Since the equations will be written in terms of contravariant fluxes, these must
be defined with relation to the cartesian vertically integrated velocities used in
equations 3.1 to 3.3. The contravariant flux components are obtained by perform
ing the dot product of the contravariant length vectors and the verticallyintegrated
velocity vector U = Uci+ Vcf:
Ue= ytUcxlVc
U,7 = yc + XVc
(3.10)
(3.11)
Following the procedure used by Rosenfeld et al. (1988) for their threedimensio
nal model, the velocity vector can also be written in terms of the covariant directions
as
U = LUC + LU"
(3.12)
To ensure the invariance of the velocity vector, the relationship between the
contravariant length vectors L and L'1 and the covariant "length" vectors Le and
L, must be given by
Le.LE = L7.L7 = 1
(3.13)
Lf.L, = L".L = 0 (3.14)
which allows for the simple formulation used for the advection terms. Substituting
equations 3.8 and 3.9 into equations 3.13 and 3.14 leads to
L, = + 
A A
(3.15)
12
iL l+ A' (3.16)
where A is the area of the cell surrounding the point where the vectors are defined,
and also the Jacobian of the transformation, given by
A = x y, ,ye (3.17)
To get the results in a more easily understandable form, the reverse conversion
also needs to be done, from contravariant flux components to cartesian vertically
integrated velocity components, using the following relationships:
Uc = U +L7 U" (3.18)
A A
V Y U+ + U" (3.19)
A A
The basic computational cell is defined in figure 3.3. The computed surface
elevation is an average value over the cell. In terms of physical representation, it
is assumed to be the value at the center of the cell, and is denoted by (i,j). The
contravariant flux components are solved at the right face of the cell, (i + 1,j)
for Ue and at the top face, (i,j + 1) for U". Once again, when it comes to a
physical representation, these flux components through the cell faces are converted
into cartesian verticallyintegrated velocity components at the center point of each
face of the computational cell.
U1j
i,j+1/2
ij i+1/2,j
Uij+1/2,j
Figure 3.3: Basic computational cell
3.3 The Fractional Step Method
Following the procedure detailed in Yanenko (1971), the equations of motion
(eqs. 3.1 to 3.3) were broken into a series of intermediate steps, for advection,
diffusion, Coriolis and propagation. A formal consistency analysis for this approach
has not yet been made for the full nonlinear NavierStokes equations. However, it
has been used successfully in a number of hydrodynamic models, like Liu (1988) or
Benque et al. (1982) and, depending on the validation of this particular model, is
worth trying.
Denoting the intermediate results with one, two and three stars, after the ad
vection step, the diffusion step and the Coriolis step, respectively, the equations for
the advection step can be written as
U_ U_ U UcUc 9 Ucvc
A t ~ HC) (3.20)
At 'z H ay H
Vj V 8 UcV cdVcVc1
At x H ) H ) (3.21)
At '9x H 5y H
For the diffusion step, only the diffusion terms are retained,
U* U
At
Vc Vt
At
Sa2 uc
A A
x2
8 2 Uc
+ AH
By2
a2vc a2vc
=AH +AH
8x2 ay2
and the Coriolis step is
UA** U*
At
fVc
The propagation step includes the continuity equation and the remaining terms
in the momentum equations, pressure, wind stress and bottom friction terms:
n+A t"
At
UC+1 _
a1+l U*
At
auc
Ox
avc
ay
(3.26)
(3.27)
(3.28)
as
= gH + Tx Tb
= gH + ruy Tby
ay
The contravariant equations for each of these intermediate steps will be the
object of each of the next four sections.
(3.22)
(3.23)
V^* V* *
SAV = f Uc
At
(3.24)
(3.25)
15
3.4 The Advection Step
The contravariant form of the equations for each of the intermediate steps will
always be written in the form
AU
A = L .F (3.29)
At
AU'7
A = L'.F (3.30)
At
where F depends on the specific step being solved.
The advection terms in the momentum equations deal with the advection of
the verticallyintegrated velocity vector through a computational cell. Writing the
verticallyintegrated velocity vector as a function of the covariant "length" vectors,
as in equation 3.12, F is defined as
a (UUMe UU'7" .\ Ut'7
H H By H
U UH ) (3.31)
Substituting F in equations 3.29 and 3.30 we get the contravariant equations
for the advection step
U'* U" Y 7, a ze UtUE y,7 d z, UtU'
At A a A H ) A a(A H
y, a (Xz UU'7 y, a z, U'7U
A a]\A H A A \A H H
+, ( yA UCU x, a y, UaUtM
A9(A H ) A9a(A H
S y, UtUt
+ a 'uI
A a A H
A9B(A H )
A a (UA HU
AaA yH
A8 (AH )
A8, UUA H )
+yf a
Ye a
x, a
Aa
zx a
A a?7
~Aa)^
y, UCU?7
kAH
(A H
(yr~u^)
(wAU"Uf)
Since the e and r7 derivatives are totally separable, each of these equations
can still be broken into a esweep equation and a trsweep equation to be solved
sequentially. If UeS and U'0 denote the intermediate flux components after the
esweep,
Sy,7 0
A bC
A a9
A+
eA a
A ae
y= y
A a?
(zx UEUC
[A H
AH
(x UtUE
A H )
zA UtUI
Y,7 a
A a9
+ A7
AB a
Yn( a
A 9r]
Aan
X, a
+ A7
X, UUt
(xy, UtU
A H
A H
( fUEU
\A H )
(x1Uf b
(y U t7U
aX a
~~A~9r
(3.37)
\A H
YtUbP7)
U,* U""
at
(3.32)
(3.33)
U& U"
At
U"# U'"
At
(3.34)
ax (uu A H
A AH
Ue U&
At
(3.35)
x a (ye UtU7
A (A H )
U"' U'1#
At
y a re UtU )
A A7H 
(3.36)
ae 9
A ar
+ a
Ae
Ae a
y17 UIUt
A H)
These four equations are solved explicitly, with the spatial derivatives being
approximated by centered differences around the right side of each cell, where the
Ut fluxes are computed, or around the top side of each cell, where the U' fluxes are
computed.
Denoting the center of a generic cell with the subscripts i and j, the right side
face of the cell will be denoted by i + 1,j and the top side of the cell by i,j + 1
(figure 3.3).
The discrete form of the esweep equations is therefore written as
U& 
i+2
U j
U ,, 7,+ At x,+,i i+lV i+l,
2+ Ai+I Ai+1, ,j Ai, H"
At n,;+ U t" u:tu ,i Y u ,i7
A.._+_ Yt!+i,. U A,,,;_
yii+IJ S+yu
Ai+ Ai+l, HA, Ai, Hi
xV t,," /U
+_ __+_ Y i+l, +l j +j i 7+ U ,,i 1
Ai+ j Ai+,, HH, A.,i ,
Ai i H. Ui," i
,, A+t U+,n+ V++
A i+ Ai+ ,Ii+ ji ..,
U7. +" 7" 2
xz . ,1+ I.iI
Ai j+, Ai+., HH.
2 2j 2 S+ 21 2
(3.38)
U +, n ++
, At y,+, + n"+',I +
Aj+j Ai+,j+ H 1,j ,
& i+ H ,i+ U+
Aj.~ H i 1
2 l+ 1 11 j 2
(3.39)
and the equations for the trsweep become
At# +#
=' Ui +,At si+ ,i+ i+ ,,i+ +1. ,i+ 1
+21, Ai+ Ai+ 1j+ H"
+ i+ 2U1, U ,, !+i+
Ai+ ,i HI + ,"
*,J! 1
1i+ ,iAt
&+At
'2
2 +22 2 2+
Ai+,i+ iH +,+
A H"+$+(
Ai+ ,i i+,i )+ ,J
xAi+ Hr ,i ,
Ai+d A H+"+i+,+ I i+rl', l .+,
dA+ i Ai + ,.
7 At Y& e7#
X,+ AI UE U"+ U+
 ____+_,_ Y'+* J+, t+,++ 2.v3+ 2
Ai+. i+I", +i
A, 2 2 + 1, 2
2 2 2 ly 2
(3.40)
U7* .! I,+ .^,i + 1 1
2^ 1 2e Ajr IAjj+j Hi
U. a U .a2
Ai, Hi
,i+ i,
19
AA_,i +1 i+l A, H! ,
A,+ Ai,i+1 H"i+ Ai, H",
+ t yA,+t "* + + jY,, Ui i
A,,^ ^A.,, i (. A, ^
(3.41)
Since the UE and U" equations are solved at different positions, interpolated
values of each component must be computed at the position where the other com
ponent is solved. These interpolations, like all those made in the other intermediate
steps, are made by fourpoint averaging.
U, =
U +t +U.+ U +U.
i+ ,j i+ ,j+1 t ,f+1 + +
2 2 211
(3.42)
U.+ + UL. + 7 U 4 + U177
2+13 4
(3.43)
3.5 The Diffusion Step
In the diffusion step, the turbulence closure problem is solved assuming a con
stant and uniform horizontal eddy viscosity. This makes it possible to take the
eddy viscosity out of the derivatives leading to a simple definition of the F vector
in equations 3.29 and 3.30:
F = AHV'2U
(3.44)
In the contravariant coordinate system, the Laplacian of the verticallyintegra
ted velocity vector can be written as
la
 A<9
" 119[a
(YU) 5 (Y)]
A [ x,) + ()
1a (y,) a (y,)
() + )]}
(3.45)
Writing the verticallyintegrated velocity vector U as a function of the con
travariant flux components, and performing the dot product with the length vectors,
the equation for UE, equation 3.29 becomes
ai \ A A A
* X17Y7 Un)
A
A Te
( A
+ U) + (xUe + XXU)]}
A ) A {a A(A
+Ang i 8d 1 EyU, + zy U4
A' ByQr A [9^ \AA/
U + xye ( U + xe )
{[ A ( 'A [u 9
9 A A )
Y017 yUe +Y 27U
l A [9 \ A A )
U
)
+ Y~"U' n "A7aA yUe
A )A 9\ A
+ ey Ue + eyU7 U17
197 A A )I
Ut" UV
At
xex UA
( A
X 2
A ]i
8, O
A "
A2a
9 aA
A
AH 
A2 0T
a KC U +
+x"U") +
A
Equation 3.30, for U", becomes
Ur*" U"'*
At
 AH y
A2 ae
y, a (xy, + zy, U)
A 81 A A
_Q (xey + xty^
a XG U + xnYE U"7
9r7 ( A A
X a a ra
+xA U + A xU +
A atr A
HA2 at7 A [9 A
57 A U
+x U, + a
A 577
aE a
+ xlU
+AH
A29U
ae 9(
By A
+AH
A /
A
A
A a a
A /JJ]
( A
UY + U)
A
SA
(+YE U+'7 + ] a
A A I 5
A A
[a '7 2l i
8( A A]
+ YU LU
A A A A
\ A[9\{ A A j
a 8 yqyx UE + IL
A ae A A 9q\A
C+7 Ur + A UJ
A 8A \ a( A
+ BY7 A A
4 ae y AU
X~Ti )]}
(YAYU A+ 7A U
SA A )
A )\ +
9Bi A
(3.46)
(3.47)
AC a
A Te
To make things simpler, auxiliary variables B, C, D and E were defined as
follows:
B Y[ (xe__ +w _
B = Y17 Uf + X? Y UI7
A a A A
xAl
A
9 (xl_ .U +
51 A9
( X2 U +
) Y+ A
X ny
wA
A C"
"sA '? UI
a(XT" U + U'
aA A
YeC XyUE + XYlY U7)
+ A[_ (XXr Ue + U)
A 8C A7A)
a (xye + xye U)
,a ( Ex2 u x
+ ~ U + X U"7 )
+ AAJ+
D L7[ a" ( UEU +
A [e A .
9\ A
XA b AA
E = e
A
+A
A
A )
r + zXy' U"i
A
+ay A
a (UE +
at~\A
a7 (XE
9rj\ A
+ y" U1]
A
(3.50)
A '7)]
iUe + x0 yUn)] (3.51)
A
Equations 3.46 and 3.47 then become
UE" UE
At
y, BB y., C Az, BD
AH yoi+ AH AH
A2 8( A2 on A2 ae
U"" U' ye oaB
At A2 a9
ye ac
 AH A
A2 d77
zxe D
+ AH 
A2 9
z, BE
 AH oE
A28a4
xe aE
+ AH 
A29ri
The diffusion equations are then solved explicitly, with the spatial derivatives
being approximated by centered differences. Equations 3.52 and 3.53 become
(3.48)
(3.49)
UUe + 0qU'?
577 A A
a (YEYus y+ U
9^ A A
[9 (XyEU + X A
e A A
(3.52)
(3.53)
U) +
L UC. Y7i+A i At
,i +. AH w (B,+,, B,,)
i+5Ji
YAt aAt.At
i+4, i+i,j
at. At
A At
+AH ? (E?.+ E, ) (3.55)
The discrete form of the auxiliary variables is
Si+'I ,i
_ _+ ___+ = 'Y B"' _
Aij i r 2z A,i ji^ ^
SA (C Ci,) + A ( 1 At,,+ D )
*0 1 $4 71 U 1 S' + $4 i 
AA, (,+ E,,) (3.55)
2 i+i
A,3+ A+1
S'Bi,7 Y = Ue .,+ i1 U. 1
A i,i Ai+ ,i++, Ai + ,,i '
j 2 2' 2
A,j + i A ii,_ A i j+,
, A7, ,i i ,i+ U
A,2 t^i A_
AA+i A_ i A1+,
Z2 2 2
+ ,+1 xr .+ U"7. +
+ , A+
Ai,j+ 2i 4*
A,_ Ui_ )7 (3.56)
iAd3j
+  U2
CA,=,,i Ai+,j A _ i i 1iY,
Aij Ai+, +52 A_ I
Aij 2
 n 1 ,i Yi U .
A 'i'3 ',,+
X z,i+ Y,,i+ Un +
A ,i+ +2
2 
+ 1, 7i+ Ut1 .* ,.i. i,'i U 2 i ,
, A x 2 A+ uX" X+uI ,
Ai+ A2 1,i Ai+I ,
_X,,j+1 ye,,j+ 1 +
Aid+ i
"zi Y i U:2
Ay~y_ i
X2 X2
A,2 ,i V
+ 1i ,Ai* U + U.
XCij+IXr17,i+ 17 U'i X Ci. XZ,7
A+ x i U.
AiJ++2 Ai*1 '
Uni,j 4,i Y "i+ i Y* I q i *
U A, ,* 4, A ,_,U V : +
A I_ I A A Ai+ ~
2 2 2
,i U YZ) n' U ( + a + U
Aid i d+S 2 A Aij
+ U.r 2 .U,_
2. z+,i ( + i U + z2 U
SxA, ,i_ + a rii i+ 7+, 1 + _
+ AL_ "+ A i+
tpj .21, '1 I XC XC
~iIj+ ,j
A.. 1 f
',, /_,, _
E ji = +i /?i+Ii f* Yf
Ei (Y+U. 
A,,_ Ai+ Y+
22 2i
Yl"i i+1 rr
A_ i A.,i s.j+
A,_f~J 12
+ U+ ,
A"d 2
2
i+ ,i Ti
Aidj+. *
+ i,
1 .Y ) *
4iA 1 i+
(3.58)
2
'v ,i"'  U + u .
2
+, Y,,, "U ,j+ I Ys,A,+Y i+
A _idj A.+, ,.+
A 2*'+,
(2
x2
Aij 2 
A 5 *
(3.57)
Ai,i
25
Ai Uirp (3.59)At
8t2 1? I. Y7i+,,i
2^ t u (3.Y0)
+ Ai.'_ 24 Un 1 + 7 +U71 1
The F vector in the Coriolis step is given by
F = Af(Ux ) (3.60)
which, after replacing U with the contravariant flux components and taking the dot
product with the contravariant length vectors leads to
UC'" U'" fzXC7 + YY fx + Y'U" (3.61)
At A f (3.6
Equation 3.61 is solved first, implicitly in terms of UC, and using, in the right
hand side, the last known values for the other component, U'1".
+ Y"i+J'j Y'7i Y +i YC ..
= + fAt ', ,U +,
*+ia *+fa AW 2+iu
x2 + y2
Ai+J ,i + ,i
+ At 2(3.63)
Since the two flux components are being computed at different positions, the
computed values for U" after the diffusion steps must be interpolated in order to
be used in the Coriolis step. As in the previous cases, a fourpoint interpolation
scheme was used.
After getting new values for U7, they are also interpolated to be used in equation
3.62, which, once again, is solved implicitly in terms of U1.
2 Ai,j+U
xe'j+i iJ + ii y.
f At +U. .(3.64)
Ai,,+ +2
3.7 The Propagation Step
The propagation step includes not only the remaining terms in the momentum
equations but also the continuity equation. Equations for U"+' and U'"+' are first
derived from the momentum equations and substituted in the continuity equation
to get an equation on than can be solved using a conjugate gradient method. The
new values for are the used to compute UC and U".
Going back again to equations 3.29 and 3.30, the vector F is given, in the
propagation step, by
F = H ) h + r
The continuity equation, reflecting the mass balance over each computational
cell is given, in the contravariant coordinate system, by
t BUe aU"
A au arl
At a8 ar
(3.66)
The bottom friction formulation used is outlined in section 3.1. Note that, for
a contravariant coordinate system, the flux magnitude term, which, in the cartesian
coordinate system is given by
S = U + Vj
(3.67)
becomes
S(2+y) u2+(2 + 'X ) U7 + 2 (xzX, + y,) Ue U
(3.68)
The bottom friction components to be used in the momentum equations can,
therefore, be written as
pgUCS
S= C2 (3.69)
C2H2
(3.65)
pgU S
r p H2 (3.70)
rl C2H2
All three equations in this step are solved implicitly. However, a number of
problems arise from that. First, the bottom friction terms are nonlinear, and, to
be solved implicitly, must be linearized. This is done by always using the flux mag
nitude value given by equation 3.68 after the Coriolis step, and solving implicitly
only for the other occurrence of one of the flux components. The total depth is also
taken after the last time step (since the depth is changed only at the propagation
step, these are the last known depth values):
T n+ pgU"+IS***(
=C C,_ HC 2 (3.71)
+z pgU^"+ S***
+ Cn H"2 (3.72)
The other problem is that, to use the conjugate gradient as presented by Hauguel
(1979), the r derivative terms in the U equation and the e derivative terms in the
U'7 equation must be treated explicitly, so that the e derivative terms and the r7
derivative terms in the final < equation can be solved separately.
Defining a coefficient P given by
gSAt
= 1 +A H2 (3.73)
the following equations are obtained for Ue and U":
n U' gH_"At gH"At o8 n+
*** P***A aY ) ***A (
gH"At ) gH"At a At
+ Yta {y^ )+ J AYOn)+ y7'Y
0***A P B***A Xn (xf + P"'. .
At
At x,7 (3.74)
U+ U"7 gH"At a gH"At a
^Af ((Y^ ) + J^A a, (r7n)
U" + ~#***A U) + 3***A (x,")
gH"At a gH"At 9 X \X1)
0***A 86 ***A z0977
At At
Yetw + exTw (3.75)
Substituting these into the continuity equation 3.66 and replacing all occur
rences of .n+1 with + A, a single equation was derived, where A is the only
unknown, that can be solved using the conjugate gradient method:
At gAt a Hn a Xt a ( X
At A [ ***A aE (yA) + xsa (,A)
gAt a H" a a L11
A ar2 ***A [ye (5YA) + zx (XA)j I
1 9 (U"' g+ At **AH" ((yt*)
Aa^ )** A 8e **A 17
+ a gAt a { *AH" [ a (y, n)
+xTI: ( II A }e [Y )I
a At a 1
+x7, (X'" ) A ]~ P* *p ( { .x ~ )
+1 a8 (Ux ] gAt + a H a n)
A Br p*** A I [**A 8a
t gAta H" n C (Y)
5A (e I, ] A c77 #***A (Y ")
+Xz (x,") + Ata 1 ( rw e'1,) (3.76)
Since none of the implicit terms contains any cross derivatives, the same proce
dure followed by Liu (1988) can be followed, breaking this equation in two, one for
the e derivative terms and the other for the r derivative terms:
A
2g (At)2 1
a f H"
# p***A
S 1
gAt8
Sa
; + a
#** 8) e1
+ X1 a (xI7
a *
H v r c")
P***A ["B (
A
2g (At)2
a H" a 1
+Xan (IXa) 1 g ***A "I (^T ")
+8 18^ 1 ^(^ )
 f***A [ (Ye A2) +
a H" } a
1+z (ux,") + A Hn
gArt 99 k *** 19 17 P***A [
a an #***A f"
(XA)] }
a (y,,,n)
(
[AY w XTwv)]
+q (3.78)
(3.79)
The spatial derivatives are approximated by centered differences, leading to the
following discrete equations:
+ 2
Sf.***, .A+ .
22
JA,,,A
y. i. i2
 q (3.77)
[ Ai,
2g (At)2
+ x"i+,ix.,i)
(y.,_j?. ,,, + x.,7f,.,,) YA,,
1Q
+ ra
Y17i+i jy,?i,j
+zE (x")
t, (A,_,iV ,_ + _, :,_, l,_x,,i
H ,j
8i***,, iA,
gAt ,3 /3 *J + 2.A._, Y .,+ +, U
*'**
y17 i+ yli,.i xi+i,i ++ x,
St'i 3 ( n2
ji 'y,7, tjaYj:
 ,+ fi1x ,i xi,ix _,l ai +,i n )
.+x l, i ~'j niA..+ i 1 1j
H2 iU 2
S n1 i+
+i ,iA + i ,,,i1
z x
2i e+ i
i ,_ 1i4 ,,,i / _s/,1 ,,
.i+ 1+,i ,31 },ii ,
+ 2_
[2g (At)2 1 A. i+ i + x +i
+ x ,l ,. r1 + 1 .2 2
* ,3 2 J
H7"
S i_,_ + x ,i
? A *.* A1.1 + x ,71,4 ai+ 2 .+
>t~^2 1,j I
1 'y+_li, *_,1
,+ ,3+
Hn.i, (14
SH". 12
231~ + 2 ~2,+
32
gAt i ** "+
_'_ 2 2
,, + ,i+ I i+ )
H.
A iJ I j+ + ,
sf.** A ,,+i , '
,i ii+ .j, + x,,+_; .+ +i+ ,+j+
,, ,,_,,, n I + A.. , i
t,3 t32
H"n
2+ V6.i
+ + x xEx,+i
+z.,_ t ,, a .,y_ z ,,iXiil)
/1+2
. (yI '.l. X~i i, + q,1 (3.81)
R*,* 1 in'ji in i
These equations can be written, in matricial form, as
H.
MAr = f NTq (3.82)
NAg = 0 (3.83)
where
M = 0 M2 (3.84)
N = [ N2 (3.85)
The partial matrices M1 and M2 are tridiagonal matrices that contain the left
hand side of equations 3.77 and 3.78 respectively, while, in this case, and since
A'I and A2 are computed in the same positions, N1 is an identity matrix and
N2 = N1.
The right hand side vector f is defined as
/=[f] (3.86)
where fi and f2 contain the right hand side terms (except for q) in equations 3.77
and 3.78.
Solving for A in equation 3.82 and substituting in equation 3.83, a new
equation on q is derived,
(NM1NT) q = NMlf (3.87)
This is the equation that is going to be solved using the iterative procedure
described in Hauguel (1979). This consists of choosing the optimum direction and
distance that the solution must move so that the error, NA reaches an absolute
34
minimum. If m denotes the iteration number, the solution after each iteration is
given by
qm+1 m qpm+pWm
(3.88)
The new variables in the right hand side are the distance the solution must be
moved, p and the direction, W. Hauguel (1979) defines a functional J, given by
J(q) = qT (NMNT ) q qT (NMif)
J (q = 2
(3.89)
which is the base for what he calls the descent method (choice of the optimal
distance for an arbitrary direction) and the conjugate gradient method (selection of
the optimal direction).
The optimal distance p is the one for which J (qm+') J (qm) has a zero
derivative. By substituting equation 3.88 in equation 3.89, this is found to be
S (Wm)T NAm
(Wm)T NM1NTWm
(3.90)
Similarly, an optimal direction W can be found by writing
Wm = gm M m1wm1
(3.91)
35
and finding the value of pm1 for which the difference J (qm+l) J (qm) has a zero
derivative.
This is called the gradient step method because it yields a solution for /L,
rnM1 (ASm)T NTNASm
(Am"1)T NTNAmI
(3.92)
for which Wm and Wm1 are conjugate relative to (NM'NT), since substituting
equation 3.92 into equation 3.91 leads to
WmlTNM1NTW" = 0
(3.93)
To make the computation of p and I simpler, a new vector V can be defined as
Vm = M1NTWm
(3.94)
and p is now
(Wm)T NASm
pm= )T
(W,)T NV,,
(3.95)
Multiplying both sides of equation 3.94 by M, a new equation can be derived
that looks similar to equation 3.82:
MVTm = NTWm (3.96)
The full iterative procedure can now be outlined:
1. Use the value of q at the last time step as qO.
2. Solve MAo = f NTqO.
3. Compute the error g = o A.
4. Define the initial direction WO = gO.
5. Solve MVm NTWm.
6. Compute pm = )
7. Advance q, using qm+l = qm + pmWm.
8. Solve MA^m+1 = f NTqm+l.
9. Compute the new error gm+l = A.m+1 A_2m+l and stop if it falls under a
specified control parameter.
10. Compute /zm = .AA+2.
11. Define the new direction Wm+1 = g+l + pmWm.
12. Increase m and go back to 5.
The two solutions, ASi and A2 are averaged point by point, to get the final At
solution. The new values of are given by
,n+1 = + A (3.97)
The new surface elevation values are the inserted into the discrete form of equa
tions 3.74 and 3.75 to get the new flux values.
+' '+ ,i gH+,,At [ +1 n+1
+X ,i (ii+, ,, ,n,i ini++ ),i ( ,
+j 1 yi, i+,j (,i + jji
,ini
in At
Li.+I_*n+1 + x vi+ i+/^.i+Ij
PAt,.* x j+ r.+ (3.98)
r *** + ,j sigHsAt
+1,+ .o**A + +n,+ +
a;'7 ,i,,+a )^, l (Ys,,++i^+ i,)
,_,+ (1 i_,^.j + X,,+  ,,+,++,,+ W+ +
At
+ p?**.* l/  *,,,+ r,,,j+ (3.99)
+ 2
3.8 The Boundary Conditions
The present model works only for closed basins, with noslip and noflow through
boundary condition. Both flux components are taken to be 0 at the boundary, and
the surface displacement is taken to be the same as at the nearest cell, or, which is
the same, the surface displacement gradient is taken to be 0 at the boundary.
38
3.9 Numerical Stability
The complex character of the NavierStokes equations and the non linearity of
some of its terms makes it very difficult to formally analyze the numerical stability of
the schemes used. As a guideline for the largest admissible time step, however, and
since the advection and diffusion terms are solved explicitly, the stability criteria
valid for linear equations can be used:
At < A (3.100)
< (Ax)2
At < (3.101)
4AH
where the minimum linear dimension of the grid cells and the maximum velocity in
the basin should be used.
Since the U" term in the UE Coriolis equation is also treated explicitly, the
Coriolis stability condition,
At < 1 (3.102)
can also be used as an indication. This condition is, in general, much less stringent
than 3.100, and it is unlikely to become a limitation on the possible time step.
Finally, the stability condition associated with the propagation step is
At < AX (3.103)
Since the crossderivative terms in the propagation step are treated explicitly, it
was not clear if this condition would effectively impose a limit to the possible time
step or not. However, several test cases were run with time steps up to five times
higher than this condition would allow, without any stability problems.
Equation 3.100 seems therefore to be, in the usual cases, the best indication of
the stability condition for the whole model.
The fact that the time derivatives are approached by forward differences and
the spatial derivatives by centered differences leads to firstorder accuracy in time
and secondorder accuracy in space.
CHAPTER 4
APPLICATIONS
After the development of a numerical model, it is necessary to validate it by
comparing its results with analytical solutions or with results obtained with well
known models whose accuracy has been, in due time, proved.
In the present case, a number of tests were run with the present model and with
the twodimensional version of the threedimensional model CH3D, which has been
well tested, as presented in, Sheng (1987), Sheng et al. (1988) and Sheng (1989).
4.1 Comparison with Analytical Solution
For the simple case of constant wind over a rectangular basin with flat bottom,
if only the pressure and wind terms are retained, the momentum equation in the
direction of the wind is reduced to
gH = r7 (4.1)
dx
which, when integrated, yields
o = x ( 0) + C (4.2)
gH
where the index 0 denotes a known point and C is the constant of integration. To
ensure mass conservation, and since the surface elevation gradient is uniform, the
41
constant of integration is 0 if ox is taken to be the center of the basin. Equation
4.2 then becomes
rW (x xo) (4.3)
gH
A skewed grid (figure 4.1) was designed for a square basin with 50km x 50km
with uniform depth of 3.0m. The model was run until it reached steady state, for a
uniform wind blowing from the west, corresponding to a wind stress of ldyne/cm2
and the results of the simulation were interpolated for a line parallel to the xaxis
through the center of the basin. Those results are shown in figure 4.2, together with
the analytical results.
4.2 Square Basin with Constant Slope
The same skewed grid as before (figure 4.2) was used, with a uniform bottom
slope in the southnorth direction. Depth is 7.5m at the north end and 2.5m at
the south end. The Coriolis factor f was taken as 0.O0009sec1, the eddy viscosity
coefficient AH as l0000cm2/sec and, in modeling the bottom friction, Manning's
n was taken as 0.040. The simulation was run for a constant and uniform wind
blowing from the east, with a wind stress r,, = Idyne/cm2, for 2.5 days, with a
timestep of 10min.
The CPU time needed for the computation compared favorably with CH3D
(12 min. and 44 sec. against 20 min. and 25 sec.) and the velocity and surface
displacement results at the last time step are shown in figures 4.4 and 4.6 for the
present model and in figures 4.5 and 4.7 for CH3D. Results at every timestep were
recorded at four different points (A through D in figure 4.1) and those results are
shown in figures 4.8 and 4.9 for the present model and in figures 4.10 and 4.11 for
CH3D.
L
Figure 4.1: Skewed grid for a square basin with uniform depth
MINO
3.0 M
 ANRLTTICRL SOLUTION
o NUMERICAL SOLUTION
2.5 7.5 12.5 17.5
22.5
X (KM)
Figure 4.2: Comparison with analytical solution
10.0 
0.0
5.0 @
1 .
10.0
27.5 32.5 37.5 12.5 47.5
I
Figure 4.3: Skewed grid for a square basin with uniform bottom slope
7.5 M
f N
MINO
4
, r .
ip/ ^ '
I. / 7 "
\ \ I
r r
 / / 1
I I
~, I
S / l 
.. .. / {
 / /
r s ,
I' >
Figure 4.4: Velocity results with constant slope present model
10.000 CM/SEC

# .
/ / 
It,
r t
I

II
/I

I,
I
I'
I
Figure 4.5: Velocity results with constant slope CH3D
10.000 CM/SEC
4. C r
 I I i I I I I I I
I I I4 I I I Ic I I I
Figure 4.6: Surface elevation with constant slope present model
III
Figure 4.7: Surface elevation with constant slope CH3D
I
Surface Elevation at Point A
0 12 24 36
Time ( Hours 1
48 60
OepthFveraged Current at Point A
0 12 24 36 48 60
Time ( Hours )
DepthAveroged Current at Point R
12 24 36
Time ( Hours
Surface Elevation at Point 8
\i
12 24 36
Time ( Hours
0
E
0
.3
o
u 5
10
10
48 60
0 12 24 36
Time ( Hours 1
48 60
48 60
DepthAveraged Current at Point B
J
0 12 24 36
Time ( Hours )
48 60
Figure 4.8: Time series results at points A and B present model

10
0
, r i i m i I I m I *
Surface Elevation at Point C
0 12 24 36 48
Time ( Hours )
DeothFveroaed Current at Point C
0 12 24 36 48
0 12 24 36 48
Time ( Hours )
DepthAveroged Current at Point C
12 24 36
Time ( Hours I
48 60
Surface Elevation at Point 0
10
0
a
,, ,m .. ... ..
10
0 12 24 36 48 60
Time ( Hours )
Oeptnhveraged Current at Point 0
10
I 5
0.
o
0
. 5
10
0 12 24 36 48 80
Time ( Hours )
DepthFveroged Current at Point 0
0 12 24 36
Time ( Hours I
48 60
Figure 4.9: Time series results at points C and D present model
Io
10
0
o
0 5
o
0
.
o 5
10
v
I
Surface Elevation at Point A
12 24 36 48
Time ( Hours ]
OeptnhFveroged Current at Point A
0 12 24 36 48
0 12 24 36 48 I
Time ( Hours
OepthAveraged Current at Point A
12 24 36
Tine ( Hours I
48 60
Surface Elevation at Point B
0 12 24 36 48 6C
Time ( Hours 1
DeothAveroged Current at Point B
10U
,
0
S5
to
I .
0 12 24 36 48 E
Time Hours )
DepthAveraged Current at Point B
0 12 24 36
Time ( HorTs
48 60
Figure 4.10: Time series results at points A and B CH3D
L
0
10
10
o
2 5
E
0
.o
o
.; 5
mn
*<  
. . . T _
A
Surface Elevation at Point C
0 12 24 36 48
Time ( Hours )
DepthAveraged Current at Point C
0 2 4 36 4
0 12 24 36 48
Time ( Hours )
DepthAveroged Current at Point C
0 12 24 36
Time ( Hours I
48 60
1u
o
)0
0
t
10
1
Surface Elevation at Point 0
0 12 24 36 48 60
Time ( ...Hours
DepthAveraged Current at Point 0
0 12 24 36 48 60
Time ( Hours )
DepthAveraged Current at Point 0
Oepl'hfveraged Current at Point 0
0 12 24 36
Time ( Hours )
48 60
Figure 4.11: Time series results at points C and D CH3D
10
10
U
U
0
o
 o
S5
10
r
53
The surface displacement results are practically the same for both models, while
the velocities computed by the present model are slightly higher than those com
puted by CH3D. This is consistent with all the other test cases run, and will be
discussed in the last chapter.
4.3 Square Basin with VShaped Bottom
The same skewed grid (figure 4.12) was used in this problem, with a Vshaped
bottom. The deepest part of the basin runs eastwest at the center of the basin.
Depth is 2.5m at the north and south ends and 5.0m at the center. The Coriolis fac
tor f was taken as 0.00009sec1, the eddy viscosity coefficient AH as l000cm2/sec
and, in modeling the bottom friction, Manning's n was taken as 0.04. The simula
tion was run for a constant and uniform wind blowing from the east, with a wind
stress Twz = ldyne/cm2, for 2.5 days, with a timestep of 10min.
The velocity and surface displacement results at the last time step are shown in
figures 4.13 and 4.15 for the present model and in figures 4.14 and 4.16 for CH3D.
Results at every timestep were recorded at the same four points and those results
are shown in figures 4.17 and 4.18 for the present model and in figures 4.19 and 4.20
for CH3D.
4.4 Lake Okeechobee with Constant Wind
A curvilinear grid (figure 4.21) was designed for Lake Okeechobee, in South
Florida. Lake Okeechobee is approximately circular, with an average diameter of
more than 50km and depth not exceeding 7m (figure 4.22). At the west end,
there is a vast area where depth is less than Im. The depth values used in the
computations were first digitized from a chart, interpolated for the grid nodes and
an offset was added to all depth values based on available field data, therefore
representing accurately the water level of the lake in the Spring of 1989.
Once again, the Coriolis factor f was taken as 0.00009sec1, the eddy viscosity
coefficient AH as l0000cm2/sec and, in modeling the bottom friction, Manning's n
Figure 4.12: Skewed grid for a square basin with Vshaped bottom
2.5 M
WIN N
4
7
t
* L
 c . . 
I w j
4 
. 4. .. r
r
I
/
10.000 CM/SEC
Figure 4.13: Velocity results with Vshaped bottom present model
56
tt
i 4 c bt
.1 ,/ ,
i I . I t
.  
S \

7 . 
I j
t .
10.000 CH/SEC
Figure 4.14: Velocity results with Vshaped bottom CH3D
57
Figure 4.15: Surface elevation with Vshaped bottom present model
Figure 4.15: Surface elevation with Vshaped bottom present model
58
Figure 4.16: Surface elevation with Vshaped bottom CH3D
Surface Elevation at Point A
A*1
0 12 24 36
Time ( Hours }
0 12 24 36
Time ( Hours )
48 60
48 60
DepthAveroged Current at Point A
12 24 36
Time ( Hours )
48 60
10
0
0 12 24 36
Time ( Hours
48 60
DepthAveroged Current at Point B
0 12 24 36 48 60
Time ( Hours
DepthAveroged Current at Point B
,,,,
0 12 24 36
Time { Hours I
48 60
Figure 4.17: Time series results at points A and B present model
0
J
4)
10
10 L
0
0
"
Surface Elevation at Point C
12 24 36
Time ( Hours )
48 60
0 12 24 36
Time ( Hours )
24 36
Time ( Hours I
Depthfveraged Current at Point C
. _ . 
DepthAveraged Current at Point 0
0
oi
5
iU
10
0 12 24 36 48 60
Time ( Hours I
OepthAveraged Current at Point 0
12 24 36 48 f
Time ( Hours )
12 24 36
Time ( Hours )
Figure 4.18: Time series results at points C and D present model
IV
10
0
48 60
10
o
S5
0
0
"0
S 5
10
5
10
0
48 60
. . i J I i I m
Surface Elevation at Point 0
A
5
10
0
0 12 24 36 48
Time I Hours )
Oepthfveraged Current at Point A
12 24 36 48
Time ( Hours )
OepthAveraged Current at Point A
0 12 24 36
Time ( Hours
48 60
10
 5
U
0
In
Surface Elevation at Point 8
0 12 24 36 48 60
Time ( Hours I
Depthveraged Current at Point B
A 1 24 38
0 12 24 36 48 6(
Time ( Hours )
DepthAveraged Current at Point 8
) . ,
0 12 24 36
Time ( Hours I
48 60
Figure 4.19: Time series results at points A and B CH3D
10
5
0
5
10
t U
$i
Surface Elevotion at Point C
10
0 12 24 36 48 6(
Time ( Hours i
DepthFveraged Current at Point C
Q 5
0
o
10
0 12 24 36
Time ( Hours
0 12 24 36
Time ( Hours
0 12 24 36
Tie ( Hours )
48 60
48 60
10
0
48 60
DepthAveroged Current at Point 0
3 12 24 36 48 60
Time ( Hours )
OepthAveroged Current at Point 0
12 24 36
Time ( Hours )
48 60
Figure 4.20: Time series results at points C and D CH3D
Surface Elevation at Point 0
MINO
IE
Figure 4.21: Curvilinear grid for Lake Okeechobee
Figure 4.22: Bottom contours for Lake Okeechobee
f/ j 
i "
T .\
/, I. ,
I ~~ I 1 
4 \ \ 4 . . *
.2
Figure 4.23: Velocity results for Lake Okeechobee present model
was taken as 0.04. The simulation was run for a constant and uniform wind blowing
from the east, with a wind stress ri = Idyne/cm2, for 2.5 days, with a timestep
of 10min.
The velocity and surface displacement results at the last time step are shown in
figures 4.23 and 4.25 for the present model and in figures 4.24 and 4.26 for CH3D.
Results at every timestep were recorded at six different points (A through E in
figure 4.21) and those results are shown in figures 4.27 through 4.29 for the present
model and in figures 4.30 through 4.32 for CH3D.
rl).000 tH/fE',
~ ~, 2.
I
 
tN 
4. 
Figure 4.24: Velocity results for Lake Okeechobee CH3D
10.000 Cn/SEC
Figure 4.25: Surface elevation for Lake Okeechobee present model
I I I I I I I I
Fig I4 : S e e n I I 
Figure 4.26: Surface elevation for Lake Okeechobee CH3D
Surface Elevation at Point A
0 12 24 36 48 60
Time ( Hours )
OepthRveroged Current at Point R
 :
12 24 36 48 60
Time ( Hours )
DepthFveroged Current at Point R
10
10
5
a
0
10
0I
0
10
10
o
o 5
0
s 5
Surface Elevation at Point B
0 12 24 36 48 60
10
n 5
u
0
t 5
10
48 60
Time ( Hours )
DepthAveroged Current at Point B
0 12 24 36 48 &
Time ( Hours )
Oepthlveroged Current at Point 8
0 12 24 36
Time ( Hous )
48 60
Figure 4.27: Time series results at points A and B present model
12 24 36
Time ( Hcurs )
. .
.... iIiiiiiimlllI
I
Surface Elevation at Point C
0 12 24 36 48 60
Time ( Hours )
DepthRveraged Current at Point C

0 12 24 36 48 64
Time ( Hours )
DepthAveroged Current at Point C
F
0 12 24 36
Time ( Hours J
48 60
Surfoce Elevation at Point 0
0 12 24 36 48 6
Time ( Hours )
OepthRveroged Current at Point 0
0 12 24 36 48 6C
Time Hours
OepthAveraged Current at Point 0
0 12 24 36
Time ( Hours I
48 60
Figure 4.28: Time series results at points C and D present model
10
a
10
10
; s
5
10
10
Surface Elevation at Point E
0 12 24 36 48
0 12 24 36 48
Time ( Hours )
DepthFveraged Current at Point E
0 12 24 36 48 6(
Time ( Hours )
OepthFveroged Current at Point E
0 12 24 36 48 6C
Time ( Hours I
60 0
0
10
10
5
0
to
10
0
U)
10
0
10 
0
12 24 36 48
Time ( Hours )
.hAveroqed Current at Point F
0 12 24 36 48
Time ( Hours )
OepthAveroged Current at Point F
12 24 36
Time ( Hours )
48 60
Figure 4.29: Time series results at points E and F present model
0
5
10
0
=
10
Surface Elevation at Point F
J
~
Surface Elevation at Point A
0 122364
0 12 24 36 48
Time ( Hours )
Oepthyveraged Current at Point R
0 12 24 36
Time ( Hours I
s6
48 60
10
C
10
5
0
5
10
0
Surface Elevation at Point 8
L
12 24 36 48
Time ( Hours I
DepthAveroged Current at Point 8
12 24 35 48
12 24 36 48 f
Time ( Hours 1
DepthAveroged Current at Point B
0 12 24 36
Time ( Hours I
48 60
Figure 4.30: Time series results at points A and B CH3D
3 12 24 38 48 6C
Time ( Hours )
DepthFveraged Current at Point A
W _, , __ .
1
Surface Elevation at Point C
1224 3648
12 24 36 48
Time ( Hours )
DeathRveroed Current at Point C
12 24 36 48
Time ( Hours )
OepthFveroged Current at Point C
V
0
10 
10
5
E
10
U
0
o
!o
i 5
0
s
10
48 60
Surface Elevation at Point 0
0 12 24 36 48 6(
Time ( Hours }
Depthfveraged Current at Point 0
0 12 24 36 48 6
Tine ( Hours
DepthAveroged Current at Point 0
0 12 24 36
Time ( Hours I
48 60
Figure 4.31: Time series results at points C and D CH3D
in
12 24 36
TiAe I Hours }
Surface Elevation at Point E
12 24 36 48
Time ( Hours
DeDthFveraqed Current at Point E
0
10
0
10
5
0
5
10
0
Surface Elevation at Point F
0 12 24 36 48 6(
Time I Hours )
Deothveraoed Current at Point F
10
5
0
5
10
0
0 12 24 36 48 6(
Time ( Hours I
12 24 36 48
Time Hours )
DepthAveroged Current at Point F
12 24 36
Time ( Hours )
48 60
Figure 4.32: Time series results at points E and F CH3D
12 24 36 48
Time ( Hours )
DepthAveroged Current at Point E
" l i m I .
75
4.5 Lake Okeechobee with Sinusoidal Wind
To test the longterm numerical stability of the model, the same problem as
before was run, with a uniform wind varying with time. The wind stress was taken
to follow a sine wave, of amplitude ldyne/cm2 and period 12hrs.
The simulation was run for 10 days, and the timeseries of the results at the
same six points as before are presented in figures 4.33 through 4.35.
10
0
0
to
Surface Elevation at Point A
O 0) N w0 0 w CO (M w
Time ( Hours
DepthAveroged Current at Point A
10 .. . .
S . . .
N ,, o N q 0)o o
Time ( Hours
DepthAveroged Current at Point A
O O ~T 0 ime OHours ( O
Tm o N N
Tlz'e ( Hours I
(WMyNWMVIWV
Iu
O u a0 N U) 0
N (N
Time ( Hours
Oepthveroged Current at Point B
u
5
i f
0
u
0
. S
10
10
0
to
. 5
u
O
4 5
10
o r 0) N (0 0 ^ 0) N^
N w r 0C) N a (W 0) 
Nu
Time ( Hours
ODethAveroged Current at Point B
o0 0 ) N (0 0 0) NO (D 0
T r' M 0 r' (0 0 
Tipme ( Hours I
Figure 4.33: Longterm simulation results points A and B
IVIwVVlV~fvwVVW
I V V A VV V A A
V VV V V\MIV
VrVVWWW\MlAMIWW
W\MVWMMMMMM
Surface Elevation at Point 8
Surface Elevation at Point C
0 VVVVVVVWMV VVEW
0
10 
IM
Time ( Hours I
OepthAveroged Current at Point C
to
Time ( Hours )
DepthAlveraged Current at Point C
10
10
 10   i . .. .
o a CD N Co o CD N D 0
Time ( Hours )
Surface Elevation at Point 0
S .
10
Time ( Hours
1 Depthveraged Current at Point 0
10
M 5
0
1 5
10 ,
Time ( Hours )
OepthAveraged Current at Point 0
10 . . .
N e 0S Hours 0
Time [ Hours I
Figure 4.34: Longterm simulation results points C and D
Surface Elevation at Point E
10
c
0 OD CO ) N C 0 0 4. C N 0
Time ( Hours 1
4. C) N (C
N NTi 0e
Tiwe
o0 H N 0 0
( Hours I
OepthAveraged Current at Point E
MMM
o C) N tO O 4 C) N (0 o
N N
Time [ Hours
OepthFveroged Current at Point E
1WPMAwhMW
o 4 C0 N C0 0 4. CO N (D
Tine I ours ]
10
t
0
j 5
10
N 4P4 qC C) N 4( (0 C) 
Tine ( Hours )
OeothAlveroged Current at Point F
10 . . .
I SMIWVWVVWW
I
oTi ( Hour 0
Time (Hours3
Figure 4.35: Longterm simulation results points E and F
m A tA I
iPV VV VV\Af VWVVWV
10
0 5
U
0
3
o
o 5
10
Surface Elevation at Point F
I
CHAPTER 5
CONCLUSIONS
The work leading to this thesis was started with the idea that a number of
techniques that have been developed for finite difference and finite element models
could be successfully applied together with a finite volume approach, therefore com
bining progress made in different directions in one single model. The time limits
involved cut the objectives of this work to a demonstration of the possibilities of the
method, instead of a full exploration of those possibilities. Therefore, this chapter
is a preview of possible future work more than an analysis of past work.
The advantages that derive from starting the development from the integral
form of the conservation equations instead of the differential form of the same equa
tions have been the object of Vinokur (1986) and, therefore, will not be further
analyzed here. Being naturally conservative, the resulting equations look consider
ably simpler than the corresponding finite difference equations, which seems to lead
to a considerable economy in terms of computation time. Recent work suggests that
a careful combination of cartesian and contravariant velocity components in some
terms in the equations might make them simpler. That is one line where future
work could bring some developments.
The fractional step method Yanenko (1971), allowing different terms in an equa
tion to be solved using different numerical techniques is a very powerful tool, used
before in finite difference models, as in Liu (1988) and in finite element models, as
in Baptista (1986). It opens the possibility of avoiding strict stability conditions, by
solving implicitly the terms associated with those conditions, allowing, at the same
time, other terms in the equations to be solved using simpler techniques, less time
80
consuming and, therefore, cheaper. It is also possible to use different time steps for
different terms in the equations, leading to considerable economy in terms of CPU
time.
Using a conjugate gradient method to solve the propagation step involved, as
was seen before, the linearization of the cross derivative terms. It is possible that
better ways to deal with the problem can be developed, leading to greater accuracy
without losing computational speed.
The more immediate need is, however, the extension of the model to allow for
open boundaries, needed to model any estuary or coastal area.
Finally, as was noted in the previous chapter, the velocity results obtained now
are consistently higher than those computed by CH3D. This is more evident for
smaller depths, as in the westernmost area in Lake Okeechobee. Further research
on why this happens is needed and, once again, the possibility to model open
boundaries can help, due to the large number of theoretical solutions available for
tidalforced circulation. It should also be noted that, although extreme care was
put into writing the computer code, and extensive debugging followed, it should be
kept in mind that, as in all new computer codes, it is still possible that errors might
remain in the code.
To conclude, the objective of showing the possibilities offered by the joint usage
of some of the latest developments in different numerical methods was met. The
time limitations implicit in a master's program did not allow for full exploration of
some of the features of the model and, therefore, considerable improvements should
still be possible. The fact that these first results look promising suggests that fur
ther development of the model is advisable.
BIBLIOGRAPHY
Baptista, A. E. M., "Accurate Numerical Modeling of AdvectionDominated
Transport of Passive Scalars A Contribution," Laborat6rio Nacional de
Engenharia Civil, Lisboa, 1986
Benque, J. P., J. A. Cunge, J. Feuillet, A. Hauguel and F. M. Holly,
"New Method for Tidal Current Computation," Journal of the Waterway,
Port, Coastal and Ocean Division, 108, 1982, pp. 396417
Hauguel, A., "Quelques Mots sur la Methode d'Eclatement d'Operateur avec
Coordination," Report HE042/79.39, Electricite de France, Chatou, 1979
Leendertse, J. J., "Aspects of a Computational Model for Long Water Wave
Propagation," Rept. RH5299RR, Rand Corp., 1967
Liu, Y., "TwoDimensional Finite Difference Moving Boundary Numerical
Model," Master's Thesis, University of Florida, 1988
Rosenfeld, M., D. Kwak and M. Vinokur, "A Solution Method for the Unsteady
and Incompressible NavierStokes Equations in Generalized Coordinate
Systems," AIAA 26th Aerospace Sciences Meeting, 1988
Sheng, Y. P., "On Modeling ThreeDimensional Estuarine and Marine
Hydrodynamics," in ThreeDimensional Models of Marine and Estuarine
Hydrodynamics, (J. C. G. Nihoul and B. M. Jamart, eds.), Elsevier, Am
sterdam, Oxford, New York, Tokyo, 1987, pp. 3554
Sheng, Y. P., "Currents and Contaminant Dispersion in the Nearshore Region
and Modification by a Jetpost," in Journal of Great Lakes Research, 2,
1976, pp. 402414
Sheng, Y. P., "Numerical Modeling of Coastal and Estuarine Processes Using
BoundaryFitted Grids," in River Sedimentation, vol. III, ed. Wang, Shen,
Ding, 1986, pp. 14161442
Sheng, Y. P., "Modeling WindInduced Mixing and Transport in Estuaries
and Lakes," in Estuarine Water Quality Management, ed. W. Michaelis,
SpringerVerlag, Berlin, Heidelberg, New York, Tokyo (in the press)
Sheng, Y. P. and H. L. Butler, "Modeling Coastal Currents and Sediment Trans
port," in Proc. from the 18th Coastal Engineering Conference, American
Society of Civil Engineers, 1982, pp. 11271148
82
Sheng, Y. P., T. S. Wu and P. F. Wang, Coastal and Estuarine Hydro
dynamic Modeling in Curvilinear Grids," in Proc. from the 21st Coastal
Engineering Conference, American Society of Civil Engineers, 1988, pp.
26552665
Spaulding, M. L., "A Vertically Averaged Circulation Model Using Boundary
Fitted Coordinates," in Journal of Physical Oceanography, 14, 1984, pp.
973982
Thompson, J. F., Z. U. A. Warsi and C. Wayne Mastin, "Numerical Grid Gen
eration Foundations and Applications," Elsevier, 1985
Vinokur, M., "An Analysis of FiniteDifference and FiniteVolume Formulations
of Conservation Laws," Contract Report 177416, NASA, Moffett Field,
1986
Yanenko, N. N., "The Method of Fractional Steps," SpringerVerlag, New York,
Heidelberg, Berlin, 1971
