• TABLE OF CONTENTS
HIDE
 Front Cover
 Title Page
 Acknowledgement
 Table of Contents
 List of Figures
 Abstract
 Introduction
 Brief review
 Numerical model
 Applications
 Conclusions
 Bibliography






Group Title: UFL/COEL (University of Florida. Coastal and Oceanographic Engineering Laboratory) ; 89/023
Title: A two-dimensional hydrodynamic model using a finite-volume approach
CITATION PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00076122/00001
 Material Information
Title: A two-dimensional hydrodynamic model using a finite-volume approach
Series Title: UFLCOEL
Physical Description: viii, 83 leaves : ill. ; 28 cm.
Language: English
Creator: Capitäao, Joaquim Josâe Areias, 1959-
University of Florida -- Coastal and Oceanographic Engineering Laboratory
Publication Date: 1989
 Subjects
Subject: Lakes -- Circulation -- Mathematical models   ( lcsh )
Lakes -- Circulation -- Measurement   ( lcsh )
Hydrodynamics -- Mathematical models   ( lcsh )
Coastal and Oceanographic Engineering thesis M.S
Coastal and Oceanographic Engineering -- Dissertations, Academic -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
 Notes
Thesis: Thesis (M.S.)--University of Florida, 1989.
Bibliography: Includes bibliographical references (leaves 81-82).
Statement of Responsibility: by Joaquim Josâe Areias Capitäao.
General Note: Typescript.
General Note: Vita.
Funding: This publication is being made available as part of the report series written by the faculty, staff, and students of the Coastal and Oceanographic Program of the Department of Civil and Coastal Engineering.
 Record Information
Bibliographic ID: UF00076122
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved, Board of Trustees of the University of Florida
Resource Identifier: oclc - 22254287

Table of Contents
    Front Cover
        Front Cover
    Title Page
        Title Page
    Acknowledgement
        Acknowledgement 1
        Acknowledgement 2
    Table of Contents
        Table of Contents 1
        Table of Contents 2
    List of Figures
        List of Figures 1
        List of Figures 2
    Abstract
        Abstract
    Introduction
        Page 1
        Page 2
        Page 3
    Brief review
        Page 4
        Page 5
        Page 6
    Numerical model
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
    Applications
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
    Conclusions
        Page 79
        Page 80
    Bibliography
        Page 81
        Page 82
Full Text




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




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs