|
UFL/COEL-89/023
A TWO-DIMENSIONAL HYDRODYNAMIC MODEL
USING A FINITE-VOLUME APPROACH
by
Joaquim Jose Areias Capitio
Thesis
1989
A TWO-DIMENSIONAL HYDRODYNAMIC MODEL USING A
FINITE-VOLUME 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. Pei-Fang 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 strong-willed 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 one-half 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 V-Shaped 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 V-shaped bottom . 54
4.13 Velocity results with V-shaped bottom present model . 55
4.14 Velocity results with V-shaped bottom CH3D . ... 56
4.15 Surface elevation with V-shaped bottom present model .. 57
4.16 Surface elevation with V-shaped 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
Long-term simulation results -
Long-term simulation results -
Long-term 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 TWO-DIMENSIONAL HYDRODYNAMIC MODEL USING A
FINITE-VOLUME APPROACH
By
JOAQUIM JOSE AREAS CAPITAL
December 1989
Chairman: Dr. Y. Peter Sheng
Major Department: Coastal and Oceanographic Engineering
A two-dimensional model of wind-driven circulation in a closed basin was developed
using a finite-volume 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 two-dimensional
version of the three-dimensional model CH3D. These included a square basin with
constant slope and with a V-shaped bottom and Lake Okeechobee, in South Florida.
To evaluate the long-term numerical stability of the model, a ten-day 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, well-publicized 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 Navier-Stokes equations are accepted as a good representation
of the physics involved, for instantaneous, three-dimensional 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 two-dimensional, 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 two-dimensional
vertically averaged Navier-Stokes 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 Navier-Stokes equations. The finite-element
methods have traditionally been considered to have a better ability to deal with
complex geometries, since they can use non-rectangular grid cells. The need for
rectangular grid cells was, therefore, the main limitation of finite-difference 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 boundary-fitted 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 non-conserva-
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 finite-volume 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 Navier-Stokes equations
into four separate equations, each one to be solved at a different intermediate step.
In the advection step, only the non-linear 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 two-dimensional, vertically-integrated flow in shal-
low water can be obtained from the three-dimensional Navier-Stokes 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 two-dimensional 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 + A-H-X2 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 vertically-integrated 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 cross-section 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 vertically-integrated
velocity vector U = Uci+ Vcf:
Ue= ytUc-xlVc
U,7 = -yc + XVc
(3.10)
(3.11)
Following the procedure used by Rosenfeld et al. (1988) for their three-dimensio-
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 vertically-integrated 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 non-linear Navier-Stokes 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 vertically-integrated velocity vector through a computational cell. Writing the
vertically-integrated 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 e-sweep equation and a tr-sweep equation to be solved
sequentially. If UeS and U'0 denote the intermediate flux components after the
e-sweep,
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 e-sweep 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.i-I
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 1-1 j 2
(3.39)
and the equations for the tr-sweep 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 .+,
d-A+ 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 four-point 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 vertically-integra-
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 vertically-integrated 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 E-yU, + 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
+ xl----U
+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
+ B-Y7 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 j-i^- ^
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 i-i,_ 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+-
t-pj .--21, '1 I XC XC
~i-Ij+ ,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 _id-j A.+, ,.+
A 2*'+,
(2
x2
Aij- 2 -
A 5 *
(3.57)
Ai,i
25
Ai Uirp (3.59)At
8t-2 1? I. Y7i+,,i
2^ t u (3.Y0)
+ Ai.'_ 24 U-n 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 four-point 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 non-linear, 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,") + At-a 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 ,3-1 }-,i-i -,
+ 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 ,,i-Xii-l)
/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,
(NM-1NT) q = NM-lf (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 (NM-if)
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 NM-1NTWm
(3.90)
Similarly, an optimal direction W can be found by writing
Wm = gm M m-1wm-1
(3.91)
35
and finding the value of pm-1 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,
rnM-1 (ASm)T NTNASm
(Am"1)T NTNAm-I
(3.92)
for which Wm and Wm-1 are conjugate relative to (NM-'NT), since substituting
equation 3.92 into equation 3.91 leads to
Wm-lTNM-1NTW" = 0
(3.93)
To make the computation of p and I simpler, a new vector V can be defined as
Vm = M-1NTWm
(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 = .A----A+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 no-slip and no-flow 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 Navier-Stokes 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 cross-derivative 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 first-order accuracy in time
and second-order 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 two-dimensional version of the three-dimensional 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 x-axis
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 south-north direction. Depth is 7.5m at the north end and 2.5m at
the south end. The Coriolis factor f was taken as 0.O0009sec-1, 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
time-step 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 time-step 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
Oepth-Fveraged Current at Point A
0 12 24 36 48 60
Time ( Hours )
Depth-Averoged 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
Depth-Averaged 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 )
Deoth-Fveroaed Current at Point C
0 12 24 36 48
0 12 24 36 48
Time ( Hours )
Depth-Averoged 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 )
Oeptnh-veraged Current at Point 0
10
I 5
0.
o
0
. -5
-10
0 12 24 36 48 80
Time ( Hours )
Depth-Fveroged 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 ]
Oeptnh-Fveroged Current at Point A
0 12 24 36 48
0 12 24 36 48 I
Time ( Hours
Oepth-Averaged 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
Deoth-Averoged Current at Point B
10U
,
0
S-5
-to
I .
0 12 24 36 48 E
Time Hours )
Depth-Averaged 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 )
Depth-Averaged Current at Point C
0 2 4 36 4
0 12 24 36 48
Time ( Hours )
Depth-Averoged 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
Depth-Averaged Current at Point 0
0 12 24 36 48 60
Time ( Hours )
Depth-Averaged Current at Point 0
Oepl'h-fveraged 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
S-5
-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 V-Shaped Bottom
The same skewed grid (figure 4.12) was used in this problem, with a V-shaped
bottom. The deepest part of the basin runs east-west 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.00009sec-1, 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 time-step 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 time-step 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.00009sec-1, 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 V-shaped 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 V-shaped 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 V-shaped bottom CH3D
57
Figure 4.15: Surface elevation with V-shaped bottom present model
Figure 4.15: Surface elevation with V-shaped bottom present model
58
Figure 4.16: Surface elevation with V-shaped 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
Depth-Averoged Current at Point A
12 24 36
Time ( Hours )
48 60
10
0
0 12 24 36
Time ( Hours
48 60
Depth-Averoged Current at Point B
0 12 24 36 48 60
Time ( Hours
Depth-Averoged 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
Depth-fveraged Current at Point C
. _--- .-- -
Depth-Averaged Current at Point 0
0
-oi
-5
-iU
10
0 12 24 36 48 60
Time ( Hours I
Oepth-Averaged 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 )
Oepth-fveraged Current at Point A
12 24 36 48
Time ( Hours )
Oepth-Averaged 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
Depth-veraged Current at Point B
A 1 24 38
0 12 24 36 48 6(
Time ( Hours )
Depth-Averaged 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
Depth-Fveraged 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
Depth-Averoged Current at Point 0
3 12 24 36 48 60
Time ( Hours )
Oepth-Averoged 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 time-step
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 time-step 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 )
Oepth-Rveroged Current at Point R
- :
12 24 36 48 60
Time ( Hours )
Depth-Fveroged 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 )
Depth-Averoged Current at Point B
0 12 24 36 48 &
Time ( Hours )
Oepthl-veroged 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 )
Depth-Rveraged Current at Point C
-
0 12 24 36 48 64
Time ( Hours )
Depth-Averoged 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 )
Oepth-Rveroged Current at Point 0
0 12 24 36 48 6C
Time Hours
Oepth-Averaged 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 )
Depth-Fveraged Current at Point E
0 12 24 36 48 6(
Time ( Hours )
Oepth-Fveroged 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 )
.h-Averoqed Current at Point F
0 12 24 36 48
Time ( Hours )
Oepth-Averoged 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 12-2-36-4
0 12 24 36 48
Time ( Hours )
Oepth-yveraged 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
Depth-Averoged Current at Point 8
12 24 35 48
12 24 36 48 f
Time ( Hours 1
Depth-Averoged 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 )
Depth-Fveraged Current at Point A
W _, , __ .
1
Surface Elevation at Point C
12--24 -36-48
12 24 36 48
Time ( Hours )
Death-Rveroed Current at Point C
12 24 36 48
Time ( Hours )
Oepth-Fveroged 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 }
Depth-fveraged Current at Point 0
0 12 24 36 48 6
Tine ( Hours
Depth-Averoged 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
DeDth-Fveraqed 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 )
Deoth-veraoed 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 )
Depth-Averoged 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 )
Depth-Averoged Current at Point E
" l i m I .
75
4.5 Lake Okeechobee with Sinusoidal Wind
To test the long-term 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 time-series 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
Depth-Averoged Current at Point A
10 .. . .
S- . . .
N ,, o N q 0)o o-
Time ( Hours
Depth-Averoged 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
Oepth-veroged 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
ODeth-Averoged 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: Long-term 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
Oepth-Averoged Current at Point C
to
Time ( Hours )
Depth-Alveraged 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 Depth-veraged Current at Point 0
10
M 5
0
1 -5
-10 ,
Time ( Hours )
Oepth-Averaged Current at Point 0
10 . . .
N e 0S Hours 0
Time [ Hours I
Figure 4.34: Long-term 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
Oepth-Averaged Current at Point E
MMM
o C) N tO O 4 C) N (0 o
N N
Time [ Hours
Oepth-Fveroged 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 )
Oeoth-Alveroged Current at Point F
10 . . .
I SMIWVWVVWW
-I
oTi ( Hour 0
Time (Hours3
Figure 4.35: Long-term 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
tidal-forced 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 Advection-Dominated
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. 396-417
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. RH-5299-RR, Rand Corp., 1967
Liu, Y., "Two-Dimensional 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 Navier-Stokes Equations in Generalized Coordinate
Systems," AIAA 26th Aerospace Sciences Meeting, 1988
Sheng, Y. P., "On Modeling Three-Dimensional Estuarine and Marine
Hydrodynamics," in Three-Dimensional Models of Marine and Estuarine
Hydrodynamics, (J. C. G. Nihoul and B. M. Jamart, eds.), Elsevier, Am-
sterdam, Oxford, New York, Tokyo, 1987, pp. 35-54
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. 402-414
Sheng, Y. P., "Numerical Modeling of Coastal and Estuarine Processes Using
Boundary-Fitted Grids," in River Sedimentation, vol. III, ed. Wang, Shen,
Ding, 1986, pp. 1416-1442
Sheng, Y. P., "Modeling Wind-Induced Mixing and Transport in Estuaries
and Lakes," in Estuarine Water Quality Management, ed. W. Michaelis,
Springer-Verlag, 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. 1127-1148
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.
2655-2665
Spaulding, M. L., "A Vertically Averaged Circulation Model Using Boundary-
Fitted Coordinates," in Journal of Physical Oceanography, 14, 1984, pp.
973-982
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 Finite-Difference and Finite-Volume Formulations
of Conservation Laws," Contract Report 177416, NASA, Moffett Field,
1986
Yanenko, N. N., "The Method of Fractional Steps," Springer-Verlag, New York,
Heidelberg, Berlin, 1971
|
RSS
HIDE
