A two-dimensional finite-difference model for moving boundary hydrodynamic problems

UFL/COEL-88/008

A TWO-DIMENSIONAL FINITE-DIFFERENCE
MODEL FOR MOVING BOUNDARY
HYDRODYNAMIC PROBLEMS

by

Yuming Liu

Thesis

1988

ACKNOWLEDGEMENTS

I would like to express my sincere appreciation and gratitude to my advisor

Dr. Y. Peter Sheng, Professor of Coastal and Oceanographic Engineering, for his

continuous guidance and encouragement throughout this study. I would also like to

extend my thanks and appreciation to my thesis committee members, Dr. James

T. Kirby, Associate Professor of Coastal and Oceanographic Engineering, and Dr.

D. Max Sheppard, Professor of Coastal and Oceanographic Engineering, for their

patience in reviewing this paper.

I would like to thank Dr. Yixin Yan, Dr. Tienshun Wu, Dr. Weichi Wang, Dr.

Li-Hwa Lin, and Steven Peene for their help and suggestions in this study.

Finally, I would like to express my deepest appreciation and gratitude to my

financial sponsor, K. C. Wong Education Foundation Ltd. of Hong Kong, for offering

me the opportunity to study at the University of Florida.

TABLE OF CONTENTS

ACKNOWLEDGEMENTS ...................

LIST OF FIGURES ..... ..................

ABSTRACT ...........................

CHAPTERS

1 INTRODUCTION ......................

2 DEVELOPMENT OF A TWO-DIMENSIONAL TIDAL

M ODEL ..........................

2.1 Governing Equations .................

2.2 Numerical Scheme ...................

2.3 Grid System .................. .....

2.4 Finite Difference Approximation ...........

2.5 Boundary and Initial Conditions ...........

. . . ii

. . . vii

CIRCULATION

2.6 Consistency and stability

3 COMPARISONS BETWEEN THEORETICAL AND NUMERICAL SOLU-

T IO N S . . . . . . . . . .

3.1 Comparison With a Linear Theoretical Solution ............

3.2 Comparison With a Non-linear Theoretical Solution .........

3.3 Comparison to a Theoretical Solution with Coriolis Effect ......

3.4 Comparison to a Theoretical Solution with Friction Effect . .

4 NUMERICAL SIMULATION OF FLOW OVER TIDAL FLATS ....

4.1 Properties of a Moving boundary ....................

4.2 Past Study . . . . . . . .

4.3 Numerical Treatment of a Moving Boundary . . . .

4.3.1 Method to Treat a Moving Boundary with the Weir Formu-

lation . . . . . . . .

4.3.2 Method to Treat a Moving Boundary With a Thin Water

Layer . . . . . . . .

4.4 Theoretical Solution of Wave Propagation on a Sloping Beach .

4.5 Comparison of Theoretical Solution with Numerical Solution .

5 APPLICATION TO LAKE OKEECHOBEE . . . .

6 CONCLUSIONS ................................

6.1 Conclusions .................. .............

6.2 Future Study ...............................

APPENDICES

A APPLICATION OF THE CONJUGATE GRADIENT METHOD TO THE

PROPAGATION STEP ..........

A.1 Conjugate Direction Method .

A.2 Conjugate Gradient Method . .

A.3 Application to the Propagation Step

B DERIVATION OF Z2 AND U2 .

C PROGRAM LISTING ..........

C.1 Flow Chart ..............

C.1.1 Main Routine: T2D .....

C.1.2 Subroutine: PROP ......

C.2 Program listing ............

C.2.1 Description of Parameters .

C.2.2 Program listing ........

BIBLIOGRAPHY ...............

BIOGRAPHICAL SKETCH .........

. . . . 81

. . . . 81

. . . . 82

. . . . 83

. . . . 88

. . . . 91

. . . . 91

. . . . 91

. . . . 92

. . . . 92

. . . . 92

. . . . 94

. . . . 125

. . . . 127

LIST OF FIGURES

2.1 Definition sketch for tidal equations . . . . 5

2.2. Schematic of finite difference mesh with variable rectangular cells 13

2.3 Computational grid definition . . . .... 13

3.1 Schematic diagram of a rectangular basin . . .... 23

3.2 Computational grid ........................ .. 23

3.3 Comparison between theoretical and numerical solutions for wa-
ter surface elevation ......................... 24

3.4 Comparison between theoretical and numerical solutions for ve-
locity ................... .......... .. 25

3.5 Comparison of theoretical solution to numerical results with dif-
ferent time steps ........................... 26

3.6 Comparison between theoretical and numerical solutions for wa-
ter surface elevation with nonlinear effects . . .... 29

3.7 Comparison between theoretical and numerical solutions for ve-
locity with nonlinear effects . . . .. ... 30

3.8 Comparison between theoretical and numerical solutions for wa-
ter surface elevation with coriolis effect . . ... 34

3.9 Comparison between theoretical and numerical solutions for ve-
locity U with coriolis effect ..................... 35

3.10 Comparison between theoretical and numerical solutions for ve-
locity V with coriolis effect ............ ......... 36

3.11 Comparison between theoretical and numerical solutions for wa-
ter surface elevation with friction effect . . ... 39

3.12 Comparison between theoretical and numerical solutions for ve-
locity U with friction effect ..................... 40

4.1 Definition sketch for wave propagation on a sloping beach 55

4.2 Computational grid ......................... 59

4.3 Comparison between wave profiles as predicted by theory and
the numerical model ( t7 = 7*w'2/1g, x = z*w2/4g ) ...... .62

4.4 Comparison between wave profiles as predicted by theory and
the numerical model ( 7 = '*w2/02g, x = z*w2/g ) . ... 63

4.5 Comparison between wave profiles near moving boundary as
predicted by theory and the numerical model ( T = rl*w2/42g,
z = x*w2/4g ) ............................ 64

4.6 Comparison between theoretical and numerical solutions of wa-
ter elevation ( 7 = ,7*w2/42g,t = ut* ) . . . 65

4.7 Comparison between theoretical and numerical solutions of ve-
locity ( u = u*w/4g,t = wt* ) .................... 66

5.1 The sketch of Lake Okeechobee . . . ... 68

5.2 Computational grid ......................... 68

5.3 Wind driven circulation with moving boundary . ... 70

5.4 Wind driven circulation with fixed boundary . . ... 70

5.5 Variation of wind speed with time . . . ..... 72

5.6 Locations of the moving boundary at four different time .. 72

5.7 Extra mass introduced by considering the moving boundary .73

5.8 Transient variation of dry area . . . .... 73

5.9 Comparisons of water surface elevation with moving boundary
and fixed boundary .......................... 74

5.10 Comparisons of velocity U with moving boundary and fixed
boundary ... .... .... .... ... ... ... ... 75

5.11 Comparisons of velocity V with moving boundary and fixed
boundary ................... .......... 76

5.12 Transient variation of bottom friction in the x-direction at point
A . . . . . . . . .. .. 77

5.13 Transient variation of bottom friction in the y-direction at point
A .... ..................... ......... .. 77

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 FINITE-DIFFERENCE MODEL FOR MOVING
BOUNDARY HYDRODYNAMIC PROBLEMS

By

YUMING LIU

December 1988

Chairman: Dr. Y. Peter Sheng
Major Department: Coastal and Oceanographic Engineering

STo predict the hydrodynamics of lakes, estuaries and shallow seas, a two di-

mensional numerical model is developed using the method of fractional steps. The

governing equations, i.e., the vertically integrated Navier-Stokes equations of fluid

motion, are solved through three steps: advection, diffusion and propagation. The

characteristics method is used to solve the advection, the alternating direction im-

plicit method is applied to compute the diffusion, and the conjugate gradient it-

erative method is employed to calculate the propagation. Two ways to simulate

the moving boundary problem are studied. The first method is based on the weir

formulation. The second method is based on the assumption that a thin water layer

exists over the entire dry region at all times. A number of analytical solutions are

used to validate the model. The model is also applied to simulate the wind driven

circulation in Lake Okeechobee, Florida.

CHAPTER 1
INTRODUCTION

In the past two decades, numerical modeling has been widely applied to the
study of the hydrodynamics of lakes, estuaries, coastal regions, etc. Many numerical

models ( Reid and Bodine, 1968, Leendertse, 1970, Yeh and Chou, 1979, etc.) have

been developed from the shallow water equations, i.e. vertically integrated Navier-

Stokes equations of fluid motion, using the finite difference technique.

Many numerical schemes have been proposed to solve the shallow water equa-
tions in the development of numerical models. They all have advantages and dis-

advantages. The explicit scheme is computationally simple, but the time step must

be sufficiently small such that the Courant number is less than 1 (Smith, 1969) in

order to attain numerical stability. The implicit scheme does not have this limita-

tion, but it requires solving the matrix equations. For. two- and three-dimensional
flows, it is very difficult to overcome the computational difficulties resulting from

the sheer size of the matrices. The alternating direction implicit method, or ADI

method, avoids solving complex matrix equations, but can obtain accurate solutions

only for Courant numbers less than 5 (Weare, 1979). Furthermore the ADI method

is not applicable to three-dimensional problems (Yanenko, 1971). The method of

fractional steps developed by Yanenko (1971), on the other hand is known to be

an effective method for solving complicated multi-dimensional problems in several

variables. In this scheme, the computation from one time level to the next is divided

into a series of intermediate steps. For each intermediate step, the computational

procedure is relatively simple, an exact solution can be obtained in some cases and

the time step can be quite large. However, this scheme still has the disadvantage

2
that its consistency has not completely been justified theoretically. Nevertheless, it

has been used to solve the shallow water equations (Benque et al., 1982).

In the development of numerical hydrodynamic models, the treatment of the

shoreline boundary is very important. Most existing numerical hydrodynamic mod-

els were developed based on the assumption of a fixed boundary with a vertical wall

located at the shoreline defined by the mean water depth. However, the shoreline

boundary can actually move with time in such problems as wave runup on a beach

and coastal flooding into a dry coastal region due to tides or storm surges. Some

researchers have tried to simulate this problem numerically. Reid and Bodine (1968)

considered the motion of the shoreline according to the water elevation and used

the empirical weir formulation to compute the velocity in which water flows into or

out of the dry land region. The disadvantage of this method is that the empirical

coefficients in the weir formulation are very site-dependent. Yeh and Chou (1979)

treated the shoreline boundary as a discrete moving boundary, but the impulsive

motion of the boundary could introduce serious numerical problems. Lynch and

Gray (1980) simulated the shoreline boundary as a continuous moving boundary
using continuous grid deformation. This method, however, can not be applied to

problems with complex topography.

In this study, a two-dimensional finite-difference numerical hydrodynamic model
is developed from the shallow water equations by using the method of fractional

steps. Two ways to simulate the moving boundary problem are studied. For some

simplified flow situations, analytical solutions are compared with numerical solu-

tions to verify the consistency of the fractional step method and the accuracy of

the numerical model. The numerical model is used to investigate the wind driven

circulation in Lake Okeechobee, Florida.

In Chapter 2, a careful derivation of the numerical model from the shallow

water equations will be presented. The computational mesh consists of rectangular

cells with variable sizes. The governing equations are solved through three steps:

advection, diffusion and propagation. The characteristics method is used in the

advection step, the ADI method is applied to the diffusion step, and the conjugate

gradient method is applied to the propagation step.

In Chapter 3, linear and nonlinear analytical solutions to the one- dimensional
long wave propagation are developed and used to verify the accuracy of the numer-

ical model. Theoretical solutions of tidal responses inside a rectangular basin with

Coriolis effect and bottom friction obtained by Rahman (1981) will be used to test

the ability of the model to simulate two-dimensional problems.

In Chapter 4, two ways to simulate the moving boundary problem are discussed

in detail. One way is similar to that used by Reid and Bodine (1968). Another way

is to include the dry land region into the computational domain by assuming a thin

water layer on the dry land region at all times. The theoretical solutions for wave

propagation on a linearly sloping beach developed by Carrier and Greenspan (1958)
are then presented and compared to the numerical solutions to verify the second

method for the moving boundary problem.

In Chapter 5, the wind driven circulation in Lake Okeechobee, Florida, including

the effects of the Coriolis force and a moving boundary, is investigated numerically.

In Chapter 6, conclusions are drawn and suggestions are made towards further
studies on numerical simulation of the moving boundary problem.

CHAPTER 2
DEVELOPMENT OF A TWO-DIMENSIONAL TIDAL CIRCULATION

MODEL

In this Chapter, a two-dimensional implicit tidal circulation model is devel-

oped using the method of fractional steps. The vertically-integrated momentum

and continuity equations, which govern the two-dimensional tide-generated cur-

rents, are solved through three steps which include an advection step, a diffusion

step and a propagation step. Momentum advection is solved using the method of
characteristics. An alternating direction implicit (ADI) method is applied to calcu-

late momentum diffusion. The wave propagation is calculated using the conjugate

gradient method.

2.1 Governing Equations

The mathematical equations describing tidal flow in shallow water can be ob-

tained by vertically integrating the three-dimensional Navier-Stokes equations of

fluid motion. Generally, it is assumed that the density of water over the depth is

constant and the vertical pressure variation is hydrostatic, thus leading to the fol-

lowing continuity and momentum equations in a right-handed Cartisian coordinate
system shown in the Fig. 2.1

az aU av

au auu avU z rbz r,, a au a a u
+ + +gh fV + [(K x ) + (K ) = 0 (2.2)
av auv avv az p bz a avy
V uV vV Z rV r (K )+V (K )V= (2.3)
+ + + gh + U+ ) = 0 (2.3)
xt 'z y ay P ax ax ay ay

z

FREE SURFACE

Z (X, .T) TV

u.U

ZB( X.T)

x

Figure 2.1: Definition sketch for tidal equations

where t is time, x and y are the spatial coordinates, Z(x,y,t) is the free surface

elevation, U(x,y,t) and V(x,y,t) are the unit-width discharges in the x- and y-

directions, respectively, u(z, y, t) and v(z, y, t) are the vertically-averaged velocities

corresponding to U(x,y,t) and V(x,y,t), h(z,y,t) is the water depth, f is the

Coriolis acceleration parameter, g is the gravitational acceleration, p is the density

of water, K is the horizontal turbulent diffusion coefficient, r,, and r,, are the wind

shear stresses in the x- and y-directions, respectively, and rb, and rb, are the bottom

shear stresses in the x- and y-directions. rb2 and rby can be expressed as

pgU U +V2
rTb = pg 2 (2.4)
C2h2

and

pgVVU2 + V2
Tbv = (2.5)
C2h2

where C is the Chezy coefficient which is

R1/6
C = 8.21-- cm l2/sec. (2.6)
n

or

R1/6
C = 1.49-- ft/2/sec. (2.7)
n

where R is the hydraulic radius and n is the Manning coefficient.
Obviously we have
U=uh

h=Z-Zf
where Zb is the bottom bed elevation of an estuary which is only a function of x
and y. Using these relationships and Eq. (2.1), the nonlinear terms in Eqs. (2.2)
and (2.3) can be expressed as

auU 9vU au au az
-+ = h(u + v-) u- (2.9)
xz dy xz 9y t"

and

auV avV v 9v az
S+ = h(ua + v-) v-- (2.10)
9z By 9x 9y .9t
Substituting them into Eqs. (2.2) and (2.3), we obtain

aU au au aZ aZ rb r7, 8(KK,) a(K-U)
- + h(u + v ) "u + gh- OV + + ] = 0
St Tz Ty at ax p az ay
(2.11)
av av av z hZ u + rbv r) a(K~ =
w h(u4 ^- -v- .g -+ [a(K)a + ay
+h(u +v-)-v-+gh + U + Tel=
9t TX Ty 9t 9y p zx ay
(2.12)
Using the method of fractional steps, Eqs. (2.1), (2.11) and (2.12) can be solved
through three steps which are called advection step, diffusion step and propagation
step (Benque et al., 1982). In order to represent working equations at each step,
we use the subscript n to denote the model variables at time nAt and n+1 for the
model variables at time (n + 1)At.
The working equations for the advection step are as follows:

u n+1/3 u" au au
t- +-- + v-=0 (2.13)
At 9z 9y

v"+l/3 _v av av
+ u-~ + v = 0 (2.14)
At xz ay
U"+1/S = un+1/3h" (2.15)

Vn+l/S = v+l/3 h (2.16)

where the subscript n+1/3 is a symbolic used to denote the result at time (n+ 1)At

due to the advection step alone.

The working equations for the computation of the diffusion step are

U"+2/s Un+1/3 a au a aU
n V [-(K + -(K = 0 (2.17)
At Lz 8x Qy +y
V"+2/S Vn+1/s a 8 v a av
At + n U [-(K ) + -(K )] = 0 (2.18)
At 2x Ox ay ay
where the subscript n + 2/3 denotes the result after the diffusion step.

The working equations for the propagation step are as follows:

Zn+1 Z" au av
t + + = 0 (2.19)
At zx Ty

Un"+ U"+2/3 aZ Z Tbz o,,
t 7 + gh-- + -- = 0 (2.20)
At at rz p
V"+1 V"+2/3 az + Z Iby r7,y =0 (2.21)
v + gh = 0 (2.21)
At at ay p
The terms ug- and vM- come from the nonlinear terms. Compared with other

terms in Eqs. (2.20) and (2.21), their values are small and could be neglected in

the propagation step.

2.2 Numerical Scheme

Just like the splitting of a single time step into three steps as shown above, three

steps could be further split into two directional steps. For example, the advection

step can be treated with the following two steps:

x-advection
Un+1/6 Un un+1/6
S + u" = 0 (2.22)
At ax

Un+1/6 Vn n-"+l1/6
+ ,U =0 (2.23)
At Ox
y-advection
U"+1/3 un+1+1/ 6OUn+1/3
At + V U+/6 y = 0 (2.24)
Aht By
V'+1/3 Vn+1/6 0n+1/3
+ vn+1/6 a = 0 (2.25)
lt aBy
where n + 1/6 represents the state after the x-advection. This scheme is implicit
and unconditionally stable.

For the diffusion step, an alternating direction implicit method is used which

leads to

x-sweep

U* Un+'1/3 a au a aBU* BUn+
S (K- () = V"+1/S (2.26)
'At ax (z ay (y

V* V"+1/3 9 av* a aV"+1/3
--- (K ) (K ) = 0 (2.27)
'At 8x 8z y ay

y-sweep

U' u a u a auv2" 3
--(K (K ) = 0 (2.28)

V+2/3 a a* a aVn+2/3
At -(K -) -(K ) = -n U* (2.29)
At Bx x Y 9y By

where U* and V* are the intermediate values of unit-width discharge during the
computation of the diffusion step. Yanenko (1971) showed that the ADI scheme
was absolutely consistent and unconditionally stable for the two-dimensibnal heat

conduction equations, which become the working equations for the diffusion step if
the source terms are added.

The propagation step is the most important of the three steps and its work-

ing equations are more complicated than the other two steps. By introducing a

coefficient a (0 < a < 1) for the spatial derivatives, the working equations for the
propagation step can be written as

Un+1 Un+2/3
At

Zn aZ U"+2/3 Zn+1 Zn
+ ac(gh )"n+ + (1 c)(gh ) -
az az h" At
+ a(,,( r")"+' + (1 a)( ( ")" = 0 (2.30)
P P

+ ci(gh-)"n +
ay

(Z V"+2/3 Zn+1 Zn
(1- a)(gh-)" hn(
By h" at

(2.31)

+ a(Ty ")+1 + (1 )( ")" = 0
P P

aU
+ (1i- )(-)X)
8z

+ V l
+ a( ) +1
ay

av)
+ (1- a)C( )" = 0 (2.32)
By

It is clear that the scheme is fully implicit when a is equal to 1 and fully explicit
when a is equal to 0. From Eqs. ( 2.30) and ( 2.31), we can obtain the formulae
for U"+1 and Vn"+ as follows

aZ aZ Un"
- At{a(gh()9 + (1- a)(gh")n- -
- At{a( )n+1 + (1 a)(Tab re)"}
P P

+2/3 Zn+1 Z"
-( At )}
(2.At
(2.33)

- At{a(gh-)"n+ +
Bay

(- a)(gh-z)" -
'y

Vn+2/3 Zn+1 Zn
h" At

- At{a(rs, r.,).+ + (1l- )(TbV -,,)n}
p p

(2.34)

Treating the bottom friction and surface friction explicitly for clarity, and substi-
tuting U"+1 and V"+1 into equation ( 2.32) yield

a hn AZ aZ"
2{(h" + AZ-) +
Xz x az
a 8 U"+2/3 a V
+ {A ( zAZ) +
gA t TX hn TY

a t U 5 +/s
(At -)

a(h"
ay ay

az"
+ AZ )}
ay

n+2/3
h-AZ)} = f +
h"

(2.35)

a (hZ) + a r, r,,
+ aa (h ) + a s)" (2.36)
x -zx g9 Tx p

Vn+1 Vn+2/3
At

Z"+' Z" U
At z

Un+l = Un+2/3

Vn+l = V"+2/3

AZ
g(At)2

where

S(1- a) a8U
1it a ( ) a
gAt 8%

10
(1 av)l a +2/3 a <9Z" a a ( r, -
(1 )(_ )V + af (h )++ ( )" (2.37)
gAt ay -gAt y ay ( y ga- y 37)

and

AZ = Z"'+ Z" (2.38)

It is very complicated to solve for AZ directly from Eq. (2.35) because of the
existence of second derivative terms with respect to x and y. Some researchers
(Benque et al, 1982) proposed solving this problem by splitting the Eq. (2.35) into
the following two equations

a -(h" + AZ, ) + AZ) = f, q (2.39)
2g(At)2 ax ax ax gAt az h"

TAZ2 2a aaz az A a vn+2/5
Sa (h" + AZ2 ) + ( Z2) = f2 + q (2.40)
2g(At)2 ay ay ay gAt ay h A

in which q(z, y) is the coordination term. If AZi = AZ2, the solution of Eq. (2.35),
AZ, will be exactly the same as the solution of Eq. (2.39) or (2.40), AZI or AZ2.
The details of solving Eqs. (2.39) and (2.40) for AZi, AZ2 and q are presented in
appendix A.
If the bottom friction terms in Eqs. ( 2.30) and ( 2.31) are treated implicitly,
the first-order series expansion in terms of velocity (U, V) and water depth h can
be used to linearize them as

(-)n+l = (rbs)n+2/3 + a T(bn+2/3Ah + a )rb)n+2/3A (2.41)
p p ah p aU p

y)n+l = ()n+2/3 + (T)n+2/3A + (v)n+2/3V (2.42)
p p dh p av p

where
Ah = h"+ h" = Z
AU = Un+1 Un+2/3 (2.43)
AV = V"+1 V"+2/3

Substituting Eqs. (2.4) and (2.5) into the above equations, we get
n+1 = )n+2/3 _g 2UvU +V +2/3
P P C2 h3
+ g ( VU2+ )n+2/3(U"l Un+2/3) (2.44)
C2 h2 + hvU2 + V2

Tbn+1 ( Vn+2/3 g (2V/U +V2 n 3
(-) = ( 3
P P +V2 h
+ (V+ )"n+2/3(Vn"+'- V"+2/3) (2.45)
C2 h2 h~U + V2)V4
Plugging them into Eqs. ( 2.30) and ( 2.31), expressions for U"+* and V"+1 can be
extracted as
U*n+I/S Z aZ
_h" az az
U"+1 = M{U2+/13 + h-----AZ aAt(gh -)"+ (1 a)At~gh )"

g UVWTV1 g 2U2 + V +
At UVU 2 )n+2/3} + aAt( 2 + U
C2 h2 C2 h2U2 + V2
g U/U2 + V )l
+ 2aAt-( )"V-2/3A (2.46)
C2 h3

V"+1 = N{V+2/3 + V-AZ aAt(gh a )"f+ (1 a)At(gh )"
AZ- ay ay
V- At (V ) U+2/3} +a }+at ( U + 2V 2/3V+2/
C2 h2 C2 h/ U2 +'V2
g 2VU + V2
+ 2aAt-(Vt )n+2/ (2.47)
CZ h3
where

M =1/,1 + aAt 2 ) 3] (2.48)
C=2 ,h2VU+V2 (2V4
and
N + agAt U2 + 2V2 )n+2/s (2.49)
N = 1/[1 + + ] (2.49)+
C2 hC 2 U2 vTV2
Substituting U""+ and V"+1 into Eq. ( 2.32), an equation for the single unknown
AZ(x, y) will be obtained. It can be solved by using a splitting scheme. The two
split equations are
AZ1 a2(h"na ) (AZaz) Maa(U n(+2 AZi-)
SMa{ + q (2.0)
2g(At)2 dx ax gAt 8x

12
a, a(-)) Na a(Vn3 AZ)
2(A)2 ,Na ( + + a N --- f2 + q (2.51)
2g(At)2 ay ay gAt ay
where

f 1-a au Ma aU. ,, ,(h ) Ma(UU2+v)"+2/3
f = )" ( )n2/ +Ma -9X + C (-2.52)
gAt Xs gAt az xz aZ xz

1-a aV Ma +V 2(h2 )" Ma(V h )n*+2/3
f, = -) ) -(y + M a y + ,(2.53)
gAt -y gAt 9y +y C2 ay

2.3 Grid System

A mesh of rectangular cells with variable sizes is established for the development
of the finite difference approximation. This is shown in Fig. 2.2 in which AX,
expresses the size for ith column and AYI represents the size for jth row.
The grid system is shoWn in Fig. 2.3. In this grid, each (U, V) grid point is
located at the center of rectangles made by four Z grid points, but each Z grid
point is not at the center of rectangles marked by (U, V) points. The cell sizes in
the (U, V) grid system and in the Z grid system have the following relationship

AX,(i) A ( AX+ (2.54)
2

AY,(,) = AY-() + A (2.55)
2
where AX,(j) and AY,(j) present the sizes of the cell (i,j) in the (U, V) grid system,
AX,(,) and AY(ij) express the sizes of the cell (i,j) in the Z grid system.
To illustrate the relationship between the numerical values in the Z and (U, V)
systems, a property denoted by P is introduced. Its value at the (U, V) grid points
is denoted by Pu,(,j) and at the Z grid points by Pz(ij). P,(i+i/2,j) is used to denote
the value at the middle point between the Z grid points (i,j) and (i + 1,j), and

Pz(i,j+1/2) represents the value at the middle point between the Z grid points (i,j)
and (i,j + 1). Given these, we have

Pu,(i) = [P(ij) + Pz(ij-l) + PC(i+,i) + Pz(i+l,y-1))/4

(2.56)

xI

I

Figure 2.2: Schematic of finite difference mesh with variable rectangular cells

J2 t

J-l

I-1 I

1*1

----- Z-GRID

-------- ----U---

-- A ( 1+

1-1 I I*1 1*2
1Xznnc

-J*2

"'T
AYU U)
.jI

Figure 2.3: Computational grid definition

AX,(_i_)_Ay--) AX, (,-)AY, )
P(iy) 4AXU,(-..)AY+,() PU(iy+l) + P4AX((i,j)

+ AX,()AY,(i-1) AX,(i)AY,(y)
4AX,(4- )AY,(y) P(i-1,y+l) + 4AX(_..-)AY() Pu(-ly) (2.57)

Pz(i+/2,y) = 2AY (i) +,)+ 2A )Pu(ij) (2.58)

2AX,(i-1) 2AXu(i-I)

P.(i+l/2,j) and Ps(i,j+l/2) can also be evaluated by

Pz(i+1/2,i) = [Pz(i,,) + Pz(i+lj)]/2 (2.60)

P.(i,,+1/2) = [Pz(i,) + P,(ij+1)]/2 (2.61)

For a spatially uniform grid system, Eqs. (2.57), (2.58) and (2.59) become

P.(ij) = [Pu(i+l)) + P,(i.) + Pu(i-.ij+1) + Pu(i-i,j)]/4 (2.62)

Pz(i+i/2j) = [Pu(i,,+1) + Pu(i,)]/2 (2.63)

P,(i,y+1/2) = [Pu(i,j+i) + Pt(i-i,i+,)]/2 (2.64)

2.4 Finite Difference Approximation

Suppose that the state of the model is known at time nAt. Then the state of
the model at time (n + 1)At can be obtained by successively solving the advection,
diffusion and propagation steps. From Eqs. (2.22) to (2.29), it is seen that the water
surface elevation Z does not appear in the advection and diffusion steps. Therefore
these steps can be executed solely on the (U, V) grid.

The finite difference equations for the advection step are
x-advection
n+1/6 n+1/6 n+1/6
i U + u"n i -- = 0 (2.65)
At AX(i-1) + AX,(i)

n+1/6 n n+1/6 n+1/6
vi, -,i + u Vi+1,) vi-i4 = 0
At + AX(-1 + AX ()

n+1/3 n+1/6 n+1/3 n+1/3
Uy u~., n+1/6 ui,+l i, -1
At* Ay.(y_+ + Ay.(y)

n+1/3 n+1/6
Ati,

n+1/3 n+1/3
+ 7+1/6 i,+1 i,,-1 = 0
'I' AY(y-1) + AYu(y)

where the subscripts i and j correspond to the (U,V) grid. After x-advection

and y-advection, we get the modified velocities u"n+1/ and v+'1/3. The discharges

corresponding to u"+1/3 and v"f+1/ can be calculated by

(2.69)

(2.70)

id = ui hI.

iV+l/s "+1/sl h
J j "i "

in which h,". is water depth on the (U, V) grid at time nAt.

The finite difference equations for the diffusion step are as follows:

x-sweep

u -- Un+1/
nAt

, v-+1/3
,y i,j

K U l-Ui y U-
K U +1.3i '1, ._* U1*- I
AX',() AX( AX(O-I)

Ay,(_i) AyX () AY.(I)

K s nis -+/ v
iI i

AX,(i) X,(i)AX,(I)
K Vs-i, + 1 iV iy id-

AY,(-[1) AYU(y) Ayu,(-1)
= 0

(2.71)

.(2.72)

y-sweep

AX-(i) AK,(i)

U:..j U:.-j-
14 Ui*Ij
r 1)

y-advection

(2.66)

(2.67)

(2.68)

Uj+2/3 -
2At

n+2/3

'At
2

16

+K U 2- U+2/3 U -+2/ Un+2/3
K [',3+1 '_ 1 i,j-1
Ay-) AYu() AY.-1)
= 0

K V*. V*. V* V.*,
AX.)W AX,() AX,(i-1)
S "+2/s V+2/3 V-n2/3 T-+2/3
S['ij+1 i,3 i i. ij-1
AY.(_-1) AY(Y) Ay(-_)
= -aU
iivl

where (i,j) also corresponds to the (U,V) grid.

After the diffusion step, the discharges U"+2/3 and V"+2/3 are known on the

(U, V) grid. The propagation step is performed in terms of these modified discharges
and consists of two computations. One is the computation of Z"+1 which is executed

on the Z grid. The other is the evaluation of Un+1 and V"+' which is performed on

the (U, V) grid. The finite difference equations for the computation of Z"+x can be

expressed as

a2 AZl(i+l,) AZ(,.) .
h AZ+(ii) AZ(i-1,y) Z -
AiX-[2y l + Xj- A x,- + A Z1(i-+/2j) --X

gAtAXu(i-l) hi A(/2,j) _-_,_(i )
= f- q (2.

a- A AZ2(ij+l) AZ2(ij) Z"j+1 Zi
AY ) [hi,+12 AY*(i) + AZ(i,1+1/2)
h AZ2,(i,) AZ2(iI -1,l-1

Un+2/3 Vn+2/3

+ gAtAY h, ^h "-AZl(i*+l/2) iT-- (2
Siu() j+1/2 ij -1/2A
= -+ q (2.
f2 2+q (2.-

75)

76)

in which

.1 rr-)+2/3 Un+2/3
(1 a) U+1l/2,, Ui-1/2,jy a l+1/2,jy -l/, U
gAt AXu(i-1) gAt AXu(i-1)

(2.73)

(2.74)

A2((it)
2g(At)2

AZ2(i.,)
2g(At)2

a. 7n _7n nn
S Zin
+ n hi+ j i hj d X 1j -l,
AX ,) hi+1/2,j A ,(i) / AX

+ pgAX (u,)[(- r)+1/2, (rT b r~ ),-1/2,jl

S- "n+2/3 _Vn+2/3
f (1 a) V,.i+1/2 .J_-1/2 a J1+/_ -. /2
sgAt AYU(,) gAt AYu(,)
S Z" Z" Z" Z"
+ aY [hl ij+1 J_ n *,J '.-1
+ AY( ,i+1/2 AY(,) 1,3-1/2 Ay.(-,)
a
"+ Yay [(QV rY),1 -+1/2 (7 ,Y)!i-1/2]
pg j J Id(

hill/2 = (hil,, + hi,y)/2
hijl/2 = (hijl + hij)/.2
AZis/2,i = (AZi~t,i + AZ~j,)/2
AZi.1/2 = (iAZi,j + AZ.,)/2

Note that the subscripts i, j, i 1/2 and j 1/2 refer to the Z-grid and the
discharges at the Z-grid can be calculated by solving Eqs. (2.57), (2.58) and

(2.59) in terms of the discharges on the (U, V) grid.

The finite difference forms for U"+* and V"n+ are

Un+l
'3

Zn+1 7n+1
= Ui+2 At_ gh1 Z+l1/2, i-l/2,
'ij 't ,AX (i)

Z" Zn
- At(1 a)ghn Z+1/2,j i/2,j
S AX.(i)

Un+2/3
+ .---- (
(z)3

- Ata At( a)
o p,,,I (-- ),,),

fn+l 7n+1
= V 2/3 Atagh +1/2 i,-1/2

Anii+1/2 Yij+/2
n n v n+2/3
- At(1 a)gh, i A"+ 2l) + +n+ (Z 1 Z

- ta At(1 a)
(si, (,)(rb, r,), by
P P

in which the indices i, j, i 1/2 and j 1/2 refer to the (U, V) grid.

(2.77)

and

(2.78)

(2.79)

v,'+
iJ

(2.80)

(2.81)

18
2.5 Boundary and Initial Conditions

Two types of boundary conditions are normally encountered in the numerical

simulation of tidal circulation. One is the open boundary which is an artificial

termination of the computational system. The other is the shoreline boundary. At

the open boundary, the water surface elevation and/or the current velocities are

specified. In this model, the velocity gradients in the direction normal to the open

boundary is set to zero, and the water surface elevation is specified along the open

boundary.

The shoreline boundary can be treated as either a moving boundary or a fixed

boundary. The moving boundary problem will be discussed in detail in Chapter

4. Along a fixed boundary, the normal velocity is zero, but the tangential velocity

may either be set to zero ( no-slip condition ) or satisfy the zero-slope in the normal

direction ( slip condition).

In this model, boundary conditions must be imposed at each of the fractional

steps. Since the water surface elevation does not appear in the advection and

diffusion steps, only velocity boundary conditions are imposed for these steps. For

the propagation step, however, boundary conditions in terms of surface elevation

and velocities must be specified along the open and closed ( fixed or moving )

boundaries.

The specification of the initial conditions requires the knowledge of the free

surface position and the flow field at t = 0. Usually this is impossible and the

model is started with the still water condition that Z =constant and U = V = 0 at

t =0.

2.6 Consistency and stability

As seen previously, the numerical scheme for each of the three fractional steps

is relatively simple. However, we need to show that the combination of these three

steps is consistent with the original governing equations. This is hard to prove

19
theoretically because of the existence of. the two-dimensional nonlinear terms and

diffusion terms. This is a drawback in applying the method of fractional steps to the

problems of fluid dynamics in several variables. It will be seen in the next chapter,

however, that this drawback does not appear to affect the accuracy of the numerical

model.

Since each fractional step is unconditionally stable, the complete scheme is un-

conditional stable. However, for numerical accuracy, it is desirable to keep the time

step sufficiently small such that the Courant number based on the velocities does

not become exceedingly larger than 1.

CHAPTER 3
COMPARISONS BETWEEN THEORETICAL AND NUMERICAL SOLUTIONS

In this Chapter, results of the two-dimensional tidal circulation model will be

compared to four different theoretical solutions. The first comparison will primarily

investigate the model's capability to simulate the propagation step. The second

comparison will validate the model's ability to compute the nonlinear effects. The

other two comparisons will focus on the simulation of the Coriolis and friction forces.

From these comparisons, we can then assess the model's accuracy and consistency.

3.1 Comparison With a Linear Theoretical Solution

Neglecting the advective terms, diffusion terms, bottom friction and wind sur-

face stress, the one-dimensional shallow water equations are

au aZ
+ gh = 0 (3.1)
t 8z

BZ aU
a+ u= (3.2)
at az

where U = uh represents the vertically integrated velocity in the x-direction, u is

the vertically averaged velocity, h is the mean water depth and Z is water surface

elevation. Assuming a closed boundary at z = I and an open boundary at z = 0,

the boundary conditions and initial conditions associated with Eqs. (3.1) and (3.2)

are:

Boundary conditions

U(, t) =0 (3.3)

Z(0,t) = Zo + asinwt (3.4)

Initial conditions

Z(,0) =Zo (3.5)

U(x,0) =0 (3.6)

where I is the length of an estuary, Zo is the still water level elevation, and a

and w are the amplitude and frequency of the forcing tide at the open boundary,

respectively.

Combining Eq. (3.1) with Eq. (3.2), we get

c Z (3.7)

where c is the wave speed. which can be expressed as c = VgT. Let the solution of

Eq. (3.7) take the form of

acosk(I z) +
z, t) cos kl sin wt + E A, sin k, sin wt + Zo (3.8)
n=O
cos kl .=o

which apparently satisfies the boundary condition (3.4) and initial condition (3.5).

Substituting Eq. (3.8) into Eq. (3.7), we can obtain the following dispersion rela-

tionship

w2 = k2 c2 (3.9)

Skn c2 (3.10)

where k and k, represent the wave numbers. From Eq. (3.1), it is seen that

boundary condition (3.3) leads to

BZ
S l,= = (3.11)

In order to satisfy this boundary condition, it follows

00
E Akn, cos kl sin wt = 0 (3.12)
n=0
Therefore

cos knl = 0 (3.13)

So
(2n + 1)7r
k = 2 (3.14)

where n is an integer. By using initial condition (3.6), we can determine the coeffi-

cient A, as

-2aw If cos k(1 x) s
An = sin k,/ dx
1,n Jo cos kl
-2aw
(3.15)
lc(kk k2)

Consequently, we obtain

acosk( x) 0[ -2aw
Z(z,t) = sint + (k k) sin kx sinwt] + Zo (3.16)

Substituting the above equation into Eq. (3.1), we get

Sac sink(l x) -2aw
U(, kt) = cost + E(k cos k,z cos wt (3.17)
cos knl n k 2 kS)

To compare the numerical solution with a linear theoretical solution, a rectangu-

lar basin shown in Fig. 3.1 with a constant water depth of 10 meters is considered.

Assume that a periodic tide with an amplitude of 50 centimeters and a period of

12.4 hrs is forced along the mouth of the basin. The still water level elevation Zo is

set to be 0. The numerical solutions can be obtained by discretizing the rectangular

basin with 15 x 30 grid points as shown in Fig. 3.2 and the use of a time step of 30

minutes.

It is noted that in the numerical computation, extra diffusion is introduced due

to the use of the implicit numerical scheme. Thus the numerical solution should

correspond to the first mode of the theoretical solution, i.e., the first terms in the

Eqs. (3.16) and (3.17). The comparison of them is shown in Figs. 3.3 and 3.4.

Figure 3.3 shows the water surface elevation near the mouth (x = 10km), at the

middle point of the estuary (x = 30km), and near the closed boundary (x = 60km).

Figure 3.4 presents the velocity near the mouth (x = 10km), at the middle point

3

0

CLOSED BOUNDRRT
K ////////////////////////
V=0 O

U

Y=O I

/0 CLOSED BOUNDRR/ 60KM///
0 CLOSED BOUNDRRT 60KH

Figure 3.1: Schematic diagram of a rectangular basin

AX= AT= 2KM

Figure 3.2: Computational grid

-bL

- THEORETICAL SOLUTION
S NUMERICAL SOLUTION

c
NJ
z 50
z o

c -50
Z-

-100

S100
-

50
z
60

ID

-J
> 0

S-50
UJ
t-

-100
60

66 72 78 84 90
TIME (HR)

66 72 78 84 90 96
TIME (HR)

Figure 3.3: Comparison between theoretical and numerical solutions for water sur-
face elevation

66 72 78 84 90 96
TIME (HR)

SIUU

N2
z 50
0
I-

6 0
LU
_J
UJ
c -5 -
I-
O

60

THEORETICAL SOLUTION
NUMERICAL SOLUTION
100
X=1OKM. Y=15KM

so

& &

u -50
-LJ

-100
60 66 72 78 84 90 96
TIME (HR)

100
X=30KM. T=15KM
so
50 -

N
0

S-50

-100
60 66 72 78 84 90 96
TIME (HR)

100
X=5OKM. T=15KM

so
so -

U -50
-lo

-100
60 66 72 78 84 90 96
TIME (HR)

Figure 3.4: Comparison between theoretical and numerical solutions for velocity

I

----- THEORETICAL SOLUTION
------- NUMERICAL SOLUTION WITH TIME STEP: 30MIN.
---------- NUMERICAL SOLUTION WITH TIME STEP: 15MIN.

100

Sso
50
z
0
Cr
> 0

UJ

I-
c -50'

z

-100
60

78
TIME (HR)

Figure 3.5: Comparison
time steps

of theoretical solution to numerical results with different

66 72 78 84 90
TIME (HR)

50
U

w
- 0

d-5
ua -50

X=30KM. T=I5KM

SI _____ 1 I

-100-
60

I I f I III I I I I...

27
(z = 30km), and near the closed boundary (z = 50km). Both Fig. 3.3 and Fig. 3.4
show that there is a reasonably good agreement between theoretical and numerical
solutions, even though a slight phase shift between the theoretical and numerical
results exists.

To investigate the origin of the phase shift, different time steps for the numerical
computation are used and the numerical results for water surface elevation and
velocity are presented in Fig. 3.5. From these results, it is seen that the phase
shift between theoretical and numerical solutions tends to diminish as the time step
decreases.

3.2 Comparison With a Non-linear Theoretical Solution

When nonlinear terms are included in two-dimensional shallow water equations,
it is impossible to obtain an analytical solution. In one-dimensional nonlinear tidal
motion, however, we can use harmonic analysis to develop a theoretical solution.
Obtaining exact theoretical solutions, however, is still difficult because the high
order terms are difficult to solve for. In this study, only the zeroth and first order

harmonic solutions are developed.

Without the consideration of Coriolis and bottom friction forces, the governing
equations for one-dimensional nonlinear tidal motion can be written as

aU aU az
+ u 5 + gh = O (3.18)

az aU
-t + = 0 (3.19)

where h is assumed to be constant. The bounary conditions are

Z(0,t) = a sinwt (3.20)

U(, t) =0 (3.21)

Based on the idea of harmonic analysis, the theoretical solution of Eqs. (3.18) and

(3.19) can be expressed as

Z = Z1 + Z2 (3.22)

U = UI + U2 (3.23)
Ux + U2
u -- + (3.24)

where Z1 and U1 are the theoretical solutions of one-dimensional linear shallow
water equations which have been obtained in the preceding section:

Sacosk(I z)
Z1 = cos kl( sin wt (3.25)
cos ki
Sac sin k( z)
Ucos = cos wt (3.26)
cos ki
Subsituting Z, U and u into Eqs. (3.18) and (3.19) and neglecting terms higher
than the first order ones, we can get the governing equations for Z2 and U2:

aU aZ2 a2c2k
+ gh O 2 k {sin[2k(1 x) + 2wt]
at ax 8h cos2 ki
+ sin[2k(l x) 2wt] + 2sin2k(l- x)} (3.27)

az2+ au 0 (3.28)
at ax
The boundary conditions for Z2 and U2 can be obtained from Eqs. (3.20), (3.21),

(3.22) and (3.23) as follows

U2(, t) =0 (3.29)

Z2(0,t) =0 (3.30)

Subject to these boundary conditions, the solutions of Eqs. (3.27) and (3.28) are

a2k I
Z2= 8h 2k l[sin2k(- x) + in sin2k(l + x)
8h cos2 ki cos4ki
s -- tan 2kl cos 2k(l x)] cos 2wt (3.31)
cos4kl

a_ 1
U2 c= 2 Icos 2k(l X) + sin 2k(l z)
8h cos2 ki 2k
I 1
cos4l cos 2k( + ) + 4 tan 2kl sin 2k(1 x)] sin 2wt (3.32)
cos4ki cos4ki

THEORETICAL SOLUTION
NUMERICAL SOLUTION

X=10KM, T=15KM

I I I I I I I I I

100

50

0

-50

-100
6

66 72 78 84 90
TIME (HR)

UU
X=60KM. T=15KM

0 -

0f0 I I I I I I I

78
TIME (HR)

Figure 3.6: Comparison between theoretical and numerical solutions for water sur-
face elevation with nonlinear effects

66 72 78 84 90
TIME (HR)

0

- IUU
E
Nj
z 50
0
I-
> 0
LJ
-I

- -50

z
h-

-100
60

z

CC
Nj

o -5

-
LU

-1

60

-- THEORETICAL SOLUTION
NUMERICRL SOLUTION

100

u
U 50

>-

so
,-
U

-100
60

78
TIME (HR)

66 72. 78 84 90
TIME (MR)

Figure 3.7: Comparison between theoretical and
with nonlinear effects

numerical solutions for velocity

66 72 78 84 90
TIME (HR)

D00

Uj

so
5 0

-50
U.'

-100
60

u
--

U-I
.U

U -50
w
I-
-so

-100
60

X=30KM. T=15KM

I I I I 1 I I 1 1 I
N<^ ^r^ ^

- I -- -- 1 -- -- I --

31
Details on the derivation of Z2 and U2 is presented in appendix B.
The same basin and forcing conditions of the preceding section will now be used
to compare theoretical solutions with numerical results obtained with a time step
of 30 minutes. The comparisons between the numerical results and the theoretical
solutions are shown in Figs. 3.6 and 3.7. Figure 3.6 presents water surface elevations
near the mouth, at the middle point, and near the closed boundary of the rectangular
basin. Figure 3.7 shows velocities near the mouth, at the middle point, and near
the closed boundary. From Figs. 3.6 and 3.7, it can be seen that the numerical
results are reasonably similar to the theoretical solutions.

3.3 Comparison to a Theoretical Solution with Coriolis Effect

Neglecting convection, diffusion, and bottom friction, the two-dimensional linear
shallow water equations can be written as

BU BZ
+au V =0 (3.33)
at xc
av 2az
+ c- + fU = 0 (3.34)
at ay
az av av (
+ --+ = (3.35)

where fl is the Coriolis parameter, and c = V/gi in which h is assumed to be
constant. Referring to the system shown in Fig. 3.1, the boundary conditions
associated with Eqs. (3.33), (3.34) and (3.35) are

Z(0,y,t) = Zm (3.36)

U(,y, t)= 0 (3.37)

V(x,0,t) 0 (3.38)

V(, b, t)=0 (3.39)

where Z, is the forced water surface elevation at the mouth of the basin, I and b
are the length and width of the basin, respectively.

I

32
By using the concepts of Kelvin waves and spectrum of Poincare waves, Rahman

(1982) proposed the theoretical solutions of Eqs. (3.33), (3.34) and (3.35) as follows:

Z = Aexp[k-y +ik(l- z) iwt]
w
flk
+ AR exp[-(b y) ik(l z) iwt]

+ -{C,[cosKi,(b-y)+ -- sinKi,(b-y)]
n=l wKln
exp[-iK2n,( z) iwt]} (3.40)

S Akc2 [lk
U -- exp[-y + ik(l x) iwt]
w w
ARkc2 ,k
+ -- exp[- (b y) ik(l z) iwt]
w w
+ -{Cn[K2cosKi(b y) + s (b
n=l KInC
ezp[-iK2n,( z) iwt]} (3.41)

ic2w 00 n2K2
V = -{CKn(1 + ) sin Kin(b-y)
W2 n2 n=1 w02Kl
n=1 in
exp[-iK2,(l z) iwt]} (3.42)

where A is the Kelvin wave amplitude at x = I and y = 0, R is the reflection

coefficient of Kelvin waves, k and w are wave number and frequency of Kelvin

waves, respectively, Kin and K2, are wave numbers of the nth Poincare mode with

respect to the y- and x-directions, respectively, and Cn is the amplitude of the nth

Poincare mode. The dispersion relationship for Kelvin waves is

w2 = k c2 (3.43)

For Poincare waves, we have

w2 22
K1. + K = 2 (3.44)

From the boundary condition of V (x,0,t) = 0, it is easy to obtain

n7r
Ki, = b (3.45)

The unknown coefficients R and C, are obtained when the boundary condition

U(l, y, t) = 0 is satisfied. Thus we have

-Ak exp(-y) + ARk exp[-(b y)]
w w
oo w
+ C,[K2 cos K,(b y) + sin K,(b y)] = 0 (3.46)
n=l K1nc

The simple way to solve for the unknown coefficients R and C, is through the matrix

method. Let us truncate the sum in Eq. (3.46) at the Nth term; the number of

unknown coefficients is then N + 1 (R and C1, ..., C,). If Eq. (3.46) is made to

hold at N +1 points on (1, y), we then have N+1 nonhomogeneous linear equations

for the same number of unknowns. Solutions are readily obtained by inversing the

matrix of the coefficients.

The rectangular basin defined previously is used to compare the tidal responses

inside the basin calculated from the theory and the numerical model. The Coriolis

parameter i is chosen to be 10". The period of the forcing tide at the mouth of

the basin is 12 hrs and the amplitude of the Kelvin wave is 50 cm. Given this basic

information, the theoretical solutions for the tidal response inside the basin can be

obtained by choosing the real parts of Eqs. (3.40), (3.41) and (3.42).

For the numerical computation, the time step is chosen to be 30 minutes and

the grid system shown in Fig. 3.2 is employed. The amplitude of the forcing tide

at the mouth of the basin is

flk flk
Z = Re{Aexp(-y + ikl) + ARexp[-(b y) ikl]

+ Z C,[cos Kin(b y) + "K sin K,,(b y)]exp(-iK2zl)} (3.47)
n=l wKin

which is obtained from Eq. (3.40) by setting t = 0. The comparison between the

numerical results and the theoretical solutions is shown in Figs. 3.8, 3.9 and 3.10.

Figure 3.8 shows the variation of water surface elevation with time near the

mouth, at the middle point, and at the closed end of the basin. Figure 3.9 presents

- THEORETICAL SOLUTION
S NUMERICAL SOLUTION

100

50

50

-50

-100
60

UU
X=30KM. T=15KM

0 -

0

0

ri l I II I I I I

66 72 78 84 90
TIME (HR)

66 72 78 84 90 96
TIME (HR)

Figure 3.8: Comparison between theoretical and numerical solutions for water sur-
face elevation with coriolis effect

66 72 78 84 90
TIME (HR)

- ,
S 1
NJ
z 5

-5
_1
0--

Cr
c -c

-1I"

60
60

iUU

S-
z 50
0
o o

I-
LU
_j1
LU
c -50
LU
c-.
-100
60

----- THEORETICAL SOLUTION
NUMERICRL SOLUTION

X=10KM. T=15KM

I I I I I

78
TIME (HR)

66 72 78 84 90
TIME (HR)

78
TIME (HR)

Figure 3.9: Comparison between theoretical and numerical solutions for velocity U
with coriolis effect

u

V-50
.'>

-s

h-
S-50
LU

-100
6

0

C-

100)

so
>-
-J

-100
60

50

0

-50

-100
60

X=50KM. T=15KM

1 11111

I

- THEORETICAL SOLUTION
NUMERICAL SOLUTION

66 72 78 84 90
TIME (MR)

25

u0
U

u -25

-so
E

50

> 0
I-

>-

" -25

-50

Figure 3.10: Comparison between theoretical and numerical solutions for velocity
V with coriolis effect

THEORETICAL SOLUTION
NUMERICAL SOLUTION

X=60KM. T=15KM

-- ,1 1-1 1 1 1 1 1
66 72 78 81 90 91
TIME (HR)

37
the flow velocities u in the x-direction at three points in the middle axis of the basin.
Figure 3.10 shows the flow velocities v in the y-direction at the middle point and at
the closed boundary. Both the theoretical and the numerical solutions at x = 10km

are very small, thus are not shown in the figures. These figures illustrate that the

model results agree quite well with the theoretical solutions.

3.4 Comparison to a Theoretical Solution with Friction Effect

With the bottom frictions, the two-dimensional shallow water equations can be
written as

+ gh a+ =0 (3.48)
at 9x p

0+ gh-a + = 0 (3.49)
Qt 9y p

az au av
BZ QU QV
-- + + = 0 (3.50)
9t ax 9y

where p is the density of water. Assume that the water depth h is constant and the
bottom frictions can be calculated by linear friction formulae

rt = pFU (3.51)

rby = pFV (3.52)

where

g 9v, + V,
FC2h2 (3.53)

Given these, Eqs. (3.48), (3.49) and (3.50) can be simplified as

QU QZ
a- + c-2 +FU=0 (3.54)

0V aZ
C- + FV = 0 (3.55)
at ay

38
az aU av
-t + + = 0 (3.56)
9t ax ay

Subject to the boundary conditions defined by Eqs. (3.36), (3.37), (3.38) and (3.39),

Rahman (1982) proposed the theoretical solutions of Eqs. (3.54), (3.55) and (3.56)

as follows

Z = Re {Ao(coskx + tanklsinkx)exp(-iwt)
00
+ E [A cos Ki(b y)
n=l
(cos K2,x + tan K2,l sin K2nx)]ezp(-iwt)} (3.57)

-C2
U = Re{ Aok(-sinkx+ tanklcoskx)exp(-iwt)
F iuw
-c2 oo
+ F= Z [AnK2n cos Kin(b y)
F --
Sn=1
(- sin K2,x + tan K2nl cos K2,n)] exp(-iwt)} (3.58)

C2 oo
V = Re{F-2 [AnKx sin Ki(b y)
n=1
(cos K2x, + tan Kz2l sin K2nz)]exp(-iwt)} (3.59)

where b and I are the width and length of the basin, respectively, and

nwr
Kin = n (3.60)

2 F
S= (1 + i-) (3.61)

K = + ) ) (3.62)
c w b
The coefficients Ao,..., An can be obtained using the condition

Z(, y, t) = Z,,exp(-wt) (3.63)

at x = 0. They are given as

Ao = Z,,dy (3.64)

39

-- THEORETICAL SOLUTION
NUMERICAL SOLUTION
100
SX=10KM, T=15KM

S50

-J

z
% o

S-100 I I I I f-I---
60 66 72 78 84 90 96
TIME (HR)

100
X=30KM. T=15KM

S-50
I-

00

60 66 72 78 84 90 96
TIME (HR)

S100
S. X=60KM. T=15KM

S50
I,-,

c -50

z
-100
60 66 72 78 84 90 96
TIME (HR)

Figure 3.11: Comparison between theoretical and numerical solutions for water
surface elevation with friction effect

----- THEORETICAL SOLUTION
NUMERICAL SOLUTION
00
X=10KM. T=15KM

0

0

nn I I I I I

66 72 78 I8 90 96
TIME (HR)

66 72 78 8 90
TIME (HR)

78
TIME (HR)

Figure 3.12: Comparison between theoretical and numerical solutions for velocity
U with friction effect

30

S-

0 -S
-s
u5
-14

>-

o
-J

(

0

0

0

-100'
60

100

0
U' SO
50-

0
= 0

I-
U -50
UJ
>

X=50KM. T=15KM

S1 1
-

-

-100
60

L-

A, = Z, cos K,,(b y)dy (3.65)

To compare the theory with the numerical results, the same basin defined pre-

viously is used. The amplitude of the forcing tide along the mouth of the basin

is assumed to be constant (50 cm), which leads to no flow in the y-direction, i.e.

V = 0, and zero value for the coefficients A1, ..., A,. The period of the forcing

tide is 12 hrs. The Manning coefficient is 0.02 and the maximum velocity is 50

cm/sec.. Thus the theoretical solutions for the tidal responses inside the basin can

be easily obtained from Eqs. (3.57) and (3.58). With a time step of 30 minutes,

the numerical results can be obtained by the use of the grid system shown in Fig.

3.2. The comparisons between the theory and the model results are shown in Figs.

3.11 and 3.12.

Figure 3.11 shows the variation of the water surface elevations with time at
three points inside the basin. And Fig. 3.12 presents the comparison for the veloc-

ities. From both figures, it is seen that the numerical results are very close to the

theoretical ones.

CHAPTER 4
NUMERICAL SIMULATION OF FLOW OVER TIDAL FLATS

In the computation of flow over tidal flats, the difficulty is to simulate prop-

erly the shoreline boundary which moves with time. .The work associated with it

includes the determination of the instantaneous location of the shoreline and the

implementation of boundary conditions. In this chapter, two ways to treat this

problem will be discussed in detail, also the wave propagation on a sloping beach

will be theoretically and numerically studied to investigate the ability of these two

methods to simulate the moving boundary problems.

4.1 Properties of a Moving boundary

The moving boundary has the basic property that it can move with time. Its

movement is mainly controlled by the topography of the coastal region, the tidal am-

plitude, the water level associated with a storm surge, etc. Based on the Lagrangian

description of fluid motion, the movement of a boundary can be mathematically ex-

pressed as

?(t) = go + 6b dt (4.1)

where X(t) is the location of the boundary at the time t, Xo is the initial location,

and 6'b is the velocity of the boundary motion. The basic boundary conditions on

the moving boundary are

h = 0 (4.2)

v = Vi, (4.3)

where h is total water depth and 6T is the velocity of a fluid particle on the moving

boundary.

43
Although the mathematical expressions for the movement of the boundary are

simple, it is difficult to couple numerically the boundary movement to the main

numerical model. The first problem which arises is the difficulty in simulating the

continuous boundary motion with conventional finite difference techniques. Without

a generalized finite difference formulation for irregular grids, it is impossible to

consider the continuous boundary motion.

The second problem is that it is difficult to compute the velocity of the moving

boundary, which is governed by the topography and the water surface gradient. Usu-

ally the topography is given, but the water surface gradient on the moving boundary

is unknown and is related to the velocity of the boundary motion implicitly.

The third problem arises from the computation of the bottom friction in the

vicinity of a moving boundary. The most common bottom friction formulation

in vertically-integrated shallow water problems is the Manning-Chezy formulation.

This formulation provides finite bottom friction throughout the interior of the com-

putational domain. However, it breaks down near the moving boundary since the

computed bottom friction approaches infinity as the water depth tends to zero. Us-

ing the conventional finite difference model, the computed strong bottom friction

causes problems both in the flooding and drying. During the flooding, the strong

bottom friction leads to very sharp surface slope which may lead to very large veloc-

ity and cause numerical instability. During the drying, the strong bottom friction

leads to very slow water motion locally.

4.2 Past Study

Several numerical models for simulating the moving boundary problem have

been proposed in the past. Most of them are developed using the finite difference

technique and finite element technique with fixed grids or deforming grids. The

technique that makes use of the fixed grids usually treats the moving boundary by

turning cells on or off at the boundary based on the mass conservation. Although

44
this technique is simpler to implement than the technique that makes use of de-

forming grids, it possesses other problems. The impulsive filling of a cell with fluid

often leads to numerical problems unless treated very carefully.

Reid and Bodine (1968) investigated the transient storm surges in Galveston

Bay, Texas. Omitting nonlinear advection and Coriolis terms, they developed a

finite difference numerical model with the inclusion of the tidal flat. A uniform

Cartesian mesh and staggered grid system were used. The elevation of the sea bed

or land was represented by a constant value at each grid point within a square

grid. Hence, the actual topography was approximated in a stair-step fashion. The

movement of the boundary was controlled by the water elevation. If the water

elevation in a flooded square was less than the base elevation of an adjacent dry

square, then a zero normal flow boundary condition was applied along their com-

mon boundary. However, if the water elevation in a flooded square was greater

than that of an. adjacent dry square, then the water was permitted to flow into the

dry square. The flow rate between two squares was determined using an empirical

equation for flow over a broad-crested barrier. The overtopping of a barrier could

be treated also. However, the model could treat only barriers aligned along the grid

mesh division. Flow across the barrier was permitted when the water height on one

side exceeded the barrier height. If the water height exceeded the barrier height on

both sides, then the flow rate was determined using an empirical equation for flow

over a submerged weir. The empirical coefficients in the model were determined

by iteration, comparing the model with tidal data and data from hurricane Carla

(September 9-12, 1963). The gross features of the inundation were predicted rea-

sonably well. However, it should be noted that the empirical coefficients used could

be very site-dependent.

Yeh and Yeh (1976) developed a nonlinear nondispersive moving boundary

model for simulating storm surge using an ADI technique. The shoreline in this

45
numerical model moved as the flow inundated low lying land. However, details

of the treatment of the boundary were not given. It appears that the shoreline

advanced or retreated in discrete increments of grid cells. They reported good

agreement with field data.

Yeh and Chou (1979) developed a nonlinear nondispersive finite difference surge

model using an explicit technique with reference to a fixed grid system. The bound-

ary between dry land and the water was simulated as a discrete moving boundary,

i.e. the boundary moves in discrete jumps. It advances or retreats according to the

rise or recession of the surge level. During the rising surge, a new grid was added to

the computations if the surge elevation of any of its neighbors was above the base

elevation of that grid point. During the receding surge, a grid point was taken out

of the computations if its total water depth was decreased to a preset value. If any

of its neighboring points still has a surge elevation above its bottom, this grid point

was kept in the computations. They compared the model to the field results and

also with a similar numerical model, which used a fictitious vertical wall instead of

a sloping shoreline. They reported that their model showed much better agreement

with field data than the fixed boundary model. The fixed boundary model pre-

dicted up to 30 percent higher surge levels than their moving boundary model. The

discrepancy was greater for higher surge values. They explained the discrepancy as
being due to the storage effect of the inland region where water can accumulate but

which is not part of the computational domain of the fixed boundary model.

Hirt and Nichols (1981) proposed a method of treating the free boundary which

was similar to the marker particle method. They called this method the volume of

fluid (VOF) technique. According to this technique, they defined a step function,

F, whose value is unity at any point occupied by fluid and zero otherwise. The

average value of F in a cell would then represent the fractional volume of the cell

occupied by fluid. In particular, a unity value of F would correspond to a cell full

46
of fluid, while a zero value would indicate that the cell contains no fluid. Cells with

F values between zero and one must then contain a free surface. The location of
J
the free boundary in a boundary cell is determined in terms of the F value and

the normal direction of the boundary. The normal direction to the boundary is

thought to be the direction in which the value of F changes most rapidly. The F

field is governed by the equation which states that F moves with the fluid at any

time. This equation is solved by the donator- acceptor iteration. With this method

to simulate the free boundary, they developed a finite difference free surface model

with a variable rectangular mesh using an iteration scheme. The model was applied

to the broken dam, undular bore and breaking bore problems, etc. They reported

good agreements with the experimental results.

Some French researchers ( Benque et al., 1982 ) developed a finite difference

numerical model with the inclusion of the tidal flat using the method of fractional

steps. They solved shallow water equations through three steps, which are advec-

tion, diffusion and propagation steps. A different numerical scheme was applied at

each computational step. The treatment of the boundary motion was considered

at the propagation step. The dry land region was assumed to be covered by a thin

water layer. The flow in this region was governed by the bottom friction. During

the computation of the propagation, the shallow water equations were first applied

to the whole computational domain including the region of the thin water layer.

Then the solution in the vicinity of the moving boundary needed to be adjusted

to satisfy the resistance flow. The Manning-Chezy friction formulation was used

to compute the bottom friction. They reported that the moving boundary model

developed in this way violated the continuity equation slightly. Good agreement of

the numerical results with the measured data were presented based on the model

application to the Bay of Saint Brieuc and the River Canche Estuary, France.

Lynch and Gray (1980) outlined a general technique whereby a moving boundary

47
can be treated by finite element Eulerian models. The finite element basis functions

are chosen to be functions of time so that the element boundaries track the moving

shoreline. They showed how this motion generates extra terms which, if treated

properly, reduce the problem to one that can be treated by standard finite element

procedures. They showed how to apply the method to treat the propagation and

runup of long waves. They looked at two simple problems involving the runup of

waves on plane beaches. They showed that estimating the runup by extrapolating

the wave height at a vertical wall could introduce significant errors. They also

pointed out that treating deforming elements is computationally more expensive

than fixed ones, and recognized that potential problem can arise if the mesh becomes

geometrically too distorted.

Gopalakrishnan and Tung (1983) described a finite-element nonlinear long wave

runup model valid only for one horizontal dimension. The model contained terms

that accounted for vertical accelerations. The moving shoreline was handled by

allowing the shoreline element to deform so that the beach node always tracked the

shoreline. If the shoreline element became too stretched, it split into two elements.

The element containing the shoreline node continued to deform but the other new

element created by splitting stayed fixed. They showed some plots that detailed

the runup process, but they did not present any results about the rundown process.

The technique outlined by the authors seems applicable to tsunami runup, but

they did not present a thorough or convincing argument to show that their model
could be used reliably for such studies. It should be noted that the technique in

their work cannot be extended easily to include two horizontal dimensions since the

element-splitting procedure would be very complex.

4.3 Numerical Treatment of a Moving Boundary

In this section, two ways to simulate the moving boundary problem will be

studied in detail. The first method was proposed by Reid and Bodine, and the

48
second method was proposed by Benque et al. (1982).

4.3.1 Method to Treat a Moving Boundary with the Weir Formulation

Reid and Bodine (1968) treated the continuous moving boundary as a discrete

moving boundary. In their scheme, a discrete Cartesian grid system is used and the

actual topography in the vicinity of a moving boundary is approximated by two-

dimensional stair-steps. Thus the elevation of the sea bed or land can be regarded

as uniform over each grid square.

The boundary motion is controlled by the water level. During the flooding, a

grid point is added to the computational system if its water depth is greater than

a preset value h, which is the minimum water depth for the effective application of

the Navier-Stokes equations. The value of hi also depends on the topography and

the grid system and is usually taken to be 10 to 20 cm. During the.drying, a grid

point is taken out of the computation if its water depth is decreased to a preset

value, h2, which is usually slightly different from hi.

In the computation of flooding, the condition of no normal mass flux is applied

at the moving boundary when the water elevation is less than that of the adjacent

dry land, i.e., the normal component of flow, Q,, at the juncture of a flooded cell

and a dry cell is taken as zero, while the tangential component of flow may satisfy

the no-slip or free-slip condition. However, if the water elevation is greater than

that of the dry land, flow will be allowed to flood into the dry cell until the water

depth in this cell reaches hi. The mass flux per unit width of the dry cell, Q,, can

be calculated by the weir formulation ( Reid and Bodine, 1968) as:

Q. =CoDD gI Zd-Z (4.4)

where Zd and Z, are the water surface elevation in the donator and acceptor cells,

respectively; Co is an empirical dimensionless coefficient which was suggested to be

49
less than 0.5 by Reid and Bodine (1968); and

D = Zd Zb (4.5)

where Zb is the bottom elevation of the acceptor cell. During the period of At, the

increment of water level in the acceptor cell, AZ,, is

AZ, = (4.6)
A,
where A, is the area of the acceptor cell and B is the width of the acceptor cell.

The decrement of water level at the donator cell, AZd, is

Q,, AtB
AZd = A (4.7)
Ad
where Ad is the area of the donator cell. During the flooding computation of each

time step, each donator cell must be examined to see if too much water is flooded

from the donator cell to the acceptor cell. If the water level in the donator cell is

less than that in the acceptor cell, the mass flux computed by Eq. (4.4) must be

adjusted to make the water level in the donator cell the same as that in the acceptor

cell.

When the water depth becomes small during the drying, an artificial water

depth has to be considered so that the Manning-Chezy friction formulation can still

be used to compute the bottom friction. The value of this artificial water depth is

determined empirically. Usually we can set it to be 20 cm.

The development of a two-dimensional moving boundary model in this manner

is very cumbersome since the location of the new boundary must be determined at

each time step. But if the grids in the vicinity of the moving boundary and the time

step are kept small, we can obtain reasonable numerical solutions for the moving

boundary.

4.3.2 Method to Treat a Moving Boundary With a Thin Water Layer

Based on an assumption that the water flow in the region of shallow water

is dominated by the bottom friction, Benque et al. (1982) proposed a way to

50
treat the moving boundary by controlling the water flow discharges through the

adjustment of the water depth. In this scheme, a grid system is first established

in the computational domain which includes the entire tidal flat region. To ensure

numerical accuracy, grid sizes in the tidal flat region should be smaller than those

in the main computational domain.

A thin water layer is assumed to exist at all times over the dry land region so

that all grids in the computational domain are always wet. Thus, no grid needs to

be taken out of the computational domain during the computation and we do not

need to consider the motion of the boundary. The shoreline boundary can thus be

treated simply as a fixed boundary, so the Navier-Stokes equations of fluid motion

can be applied to the entire computational domain so long as the bottom friction

is adequately resolved for small water depths. From the Manning-Chezy friction

formulation, we see that the water depth plays an important role in the computation
of the bottom friction, so it is possible to "control" the bottom friction by adjusting

the water depth.

As seen in Chapter 2, the free surface elevation Z does not appear during the

advection and diffusion steps. Since the motion of the boundary is mainly controlled

by the water elevation, it is only necessary to consider the propagation step. From

the governing equations, Eqs. (2.75) and (2.76), of the propagation step, it is seen

that the water depth, h, and the increment of water surface elevation during one

time step, AZ, at the points of i 1/2 and j 1/2 must be calculated. They

are governed by the water depth and the surface elevation at their upstream and

downstream grid points, and can be calculated by the following formulae:

hi+1/2 = Ti hi + (1 yi)hi+l
i-1/2 = i-1- hi-1 + (1 (4.8)
AZi+1/2 = AZ.i + (1 7i)AZi+i
AZi-1/2 = -fi-i AZi-I + (1 'Y.-i)AZ,

and

hj+1/2 = yj hi + (1 -j)hj+l
h*-1/2 = y-l hj-_ + (1 yi-I)hj (4.9).
AZ+1/2 = -j AZ, + (1 ,)AZ+1 (4)
Aj-1/2 = 'Y,-i AZ,_1 + (1 '-Y-)AZ, )
where is a weighting coefficient. Normally the shape of the water surface between
two successive grid points can be assumed to be a straight line, so that 7 can be set

to 0.5. However, in the regions where the water flow is dominated by the bottom
resistance, -7 cannot be equal to 0.5 since the shape of the water surface is greatly
curved. Therefore we need to look for a formulation to calculate it.
Neglecting the effect of the interior force, we can obtain the governing equations
for the flow dominated by the bottom resistance as follows

49Z rb.
gh + =0 (4.10)
aZ p

ghf +- = 0 (4.11)
By P
Using the Manning-Chezy friction formulation,

pgUVU2 + V2(4
rb = C2 h (4.12)

pgVVU2 + V2
v = C2 h2 (4.13)

we can rewrite Eqs. (4.10) and (4.11) as

9Z UVU2 +V2
-x + C2h3 = 0 (4.14)
Z VVU2 + V2
a- C2h3 = 0 (4.15)
For flow controlled by the bottom resistance, it can be assumed that the discharge
should always increase when the downstream level decreases, i.e,

aU
< z 0 (4.16)

av
a < 0 (4.17)
8Zd, -

52
Where the subscript ds refers to the downstream. Rewriting Eqs. (4.14) and (4.15)

in the discrete form, we have

Zd, Zu, UVU2 +V2
A + = 0 (4.18)
Az~ C2 hi+1/2

Za, Z., V VU + V2
+ = 0 (4.19)
Ay, C2 hC+112
where the subscript us refers to the upstream. In Eq. (4.18), the upstream and

downstream correspond to flow in the x-direction, but in Eq. (4.19) they refer to

flow in the y-direction. Introducing a coefficient f, the water depth h in Eqs. (4.18)

and (4.19) can be expressed as

hi+1/2 = P hu. + (1 /)hd, (4.20)

hj+l/2 = P hu, + (1 ?)hd, (4.21)

where

hu, = Z,, ZB,, (4.22)

hd, = Zda ZBd, (4.23)

in which ZB,, and ZBd, are the bottom elevations upstream and downstream of

the point i + 1/2 or j + 1/2, respectively. Differentiating Eq. (4.18) with respect to

Zd, and bearing in mind that U is a function of Z,, and Zd,, and V is constant, we

obtain

aU U2 C2h ah
(U + V2 + = ) = [-h + 3(Z,, Zd,) ] (4.24)
az, (U + U= +V) = z' a--,

Applying the condition -- < 0, we get
Oh
h + 3(Z,, Zd.) -Z, < 0 (4.25)
ozd,

From Eq. (4.20), we know
ah
1 (4.26)
8 at,

53
Substituting h- and h into inequality (4.25) yields

[fh,, + (1 P)hd} ] + 3(1 0)(Z,, Zd,) < 0 (4.27)

So we obtain
S 3(Z,, Za,) hd, .
P> (4.28)
3(Z,, Zd,) + h.u hda
Differentiating Eq. (4.19) with respect to Zda, treating U as a constant and following
the same procedure as above, we can get the same inequality for f as in (4.28).
The weighting coefficient y in Eq. (4.8) for the propagation step must be re-
placed by either f or 1 in the shallow water region. Since the direction of flow
during flooding is different from that during drying, the concepts of upstream and
downstream are changed. Thus during flooding, = )3 and during drying y = 1- P.
If 3 computed from Eq. (4.28) is less than 0, it can be set equal to 0. If it is greater
than 1, it is set equal to 1. In general, 3 is 0.5 in the computational domain except
in the transitional region between the thin water layer and the deep water region.
For the flow controlled by the bottom friction, we can obtain the maximum
discharge, Umx and V,,,, from Eqs. (4.14) and (4.15)

aZ
Uma, = sign(u) C2h3 | 1 (4.29)
aZ
Vma, = sign(v) C'h3 -z (4.30)
ay
in which
aZ aZ
sign(u) = -z/ | I (4.31)
8- ax
dZ BZ
sign(v) = -z/ I -a (4.32)
-y ay
After the propagation step, the new state of the model is known. Corresponding to
the water surface gradient at this state, we can calculate the maximum discharge
from the above two equations. Also we have discharges U and V at this state which
are computed from the momentum and continuity equations. If U or V is greater

54
than Uma, br Vma, U or V must be replaced by Uma, or Vma, before the computation
proceeds to the next time step.
It should be noted that using this way to develop a moving boundary model

leads to relatively simple computer program since we do not need to simulate the

motion of the boundary. But the continuity equation is slightly violated by always
maintaining a thin water layer on the dry land region and replacing the discharges
U and V by Uma and Vm,. The thickness of the water layer on the dry land region

is determined by the requirement that the computational cell cannot become dry

during one time step.

4.4 Theoretical Solution of Wave Propagation on a Sloping Beach

In this section, the theoretical solution for the wave propagation on a linearly
sloping beach obtained by Carrier and Greenspan (1958) will be presented briefly.

Referring to the system shown in Fig. 4.1, the one- dimensional nonlinear shallow
water equations can be written as

t9 + a[(1* + h')u'] = 0 (4.33)

au* *.u* 8n*
+ + g = 0 (4.34)

where the asterisks denote dimensional quantities, 1t is the displacement of water
surface above the mean water level, h is the still water depth which varies linearly
with z, u is the velocity in the x direction. Let L be the characteristic horizontal
length scale of the wave motion. Then we can define a time scale, T, and velocity

scale, uo, as follows

T = /L/g (4.35)

uo = V L (4.36)

where 4 is the beach angle. Let us choose the following nondimensionalization

z* t* r 7
L T eL

Figure 4.1: Definition sketch for wave propagation on a sloping beach

h* U*
h = x = (4.37)
4L Uo

and define
c2 h = h + = x + r (4.38)
L
With these definitions, Eqs. (4.33) and (4.34) become

rt + [(77 + x)u], = 0 (4.39)

ut + u uu + ri = 0 (4.40)

In terms of u and c, Eqs. (4.39) and (4.40) become

2ct + 2u c, + c = 0 (4.41)

ut + uu + 2c c = 1 (4.42)

Through the elegant series of transformations, Carrier and Greenspan (1958) were

able to transform this problem, with two coupled nonlinear equations, into a new

problem with only one linear equation. A brief derivation will be presented in the

following.

56
If Eqs. (4.41) and (4.42) are added and subtracted, we obtain

d dz
(u 2c -t) = 0 along d= u C
at dt

Let us define the characteristic variables and ( by

S=u+2c-t

S= -u + 2c + t

Hence, Eq. (4.43) becomes

= const

S= const

Now let us consider x and t to be

S= const we get
dxz ax at
dt a ea
dx x at
-n,
dti w gT
From above two equations, we get

along

along

functions

if

if

e = te (u + c)

X = tr (u c)

From Eqs. (4.44) and (4.45), we can obtain

S+ c = (3 e)/4 + t

u c = ( 3)/4 + t

Substituting u + c and u c into Eqs. (4.50) and

relationship between (x, t) and (,, )) as follows:

C = te(3 ()/4 + (t2/2)e

dx
=u + c (4.46)
dt

d = u c (4.47)
of and Then for = on or
of and e. Then for S = cost or

S= const

S= const

(4.48)

(4.49)

(4.50)

(4.51)

(4.52)

(4.53)

(4.51) yields the transform

(4.54)

(4.43)

(4.44)

(4.45)

57
x, = t,(- 3()/4 + (t2/2), (4.55)

Eliminating x from Eqs. (4.54) and (4.55), we get

2( + ))t+ + 3(t, + t) = 0 (4.56)

This is a linear partial differential equation for t((, (). It is convenient to introduce

new variables a and A by

A = (- = 2(t u) (4.57)

a = + = 4c (4.58)

Then Eq. (4.56) becomes

txx = t~, + 3t (4.59)
a
Since from Eq. (4.57), t = A/2 + u, u must also satisfy Eq. (4.59)

3
Uxn = Uao + -uo (4.60)

If we introduce a "potential" p(a, A) so that

S= (4.61)

then Eq. (4.60) reduces to
1
Px = Poo + -Po (4.62)

This is a single linear partial differential equation. The boundary condition at the

shoreline for Eq. (4.62) is

a= 0 (4.63)

which corresponds to the condition c = 0, i.e., the total water depth must be

identically zero at the shoreline for all time.

In terms of the variables a, A and the potential p(a, A), Carrier and Greenspan

(1958) proposed the following expressions for t, x, r7, u and c

t A + +A (4.64)
2 2 a

1 2 1 1 2o 1 o2
S= -u + c (+- = )2 + -= + (4.65)
2 4 2 a 4 16
22 1 a2
77 = c = = --PA (4.66)
16 4 16
= (4.67)

c = 4 (4.68)

Although Eq. (4.62) is certainly much simpler to solve than the two original coupled

nonlinear equations (4.39) and (4.40), it is difficult to obtain rl or u as explicit

functions of z and t. If p(a,A) is given, then Eqs. (4.64)-(4.68) will give t, x, rl

and u all parametrically in terms of the variables a and A. In general, it is very

difficult to eliminate a and A to obtain direct functional relationships for l7 and u

in terms of x and t.

4.5 Comparison of Theoretical Solution with Numerical Solution

A simple solution of Eq. (4.62) pointed out by Carrier and Greenspan (1958) is

p(a, ) = -8AoJo( ) sin (4.69)

where Ao is an arbitrary amplitude parameter and Jo is the Bessel function of

the first kind of order zero. This potential corresponds to a standing wave solution

resulting from the perfect reflection from the shore of a wave of unit frequency. With

p(a, A) given by Eq. (4.69), Eqs. (4.64) to (4.68) will implicitly give the solution of
this standing wave. To evaluate r(z,,t) and u(x,t) for arbitrary x and t, Eqs. (4.64)

to (4.68) must be solved numerically. For specific values of z and t, a and A can

be obtained from Eqs. (4.64) and (4.65) by using Newton's Method, so that rl(x,t)

and u(x,t) can easily be obtained from Eqs. (4.66) and (4.67), respectively.

To test the ability of the finite-difference model to simulate the moving boundary

problem, a numerical simulation of a long wave will be performed in a rectangular

basin with linearly varying water depth. All quantities used in the finite-difference

Figure 4.2: Computational grid

model are dimensional since the model is developed based on the dimensional gov-
erning equations, but the final solution obtained by the numerical model will be
converted to dimensionless form in order to be compared with the theoretical solu-
tion. In the rest of this section, all dimensional quantities will be presented with
units and dimensionless quantities will be presented without units.

The length of the rectangular basin is 60 km as measured from the mean water
level and the width of the basin is 10 km. The slope of the bottom, 0, is 1 : 2500.
The still water depth at the offshore boundary is 24 meters. The period of the long
wave is set equal to 1 hour. Thus the frequency, w', is 0.00174sec.-'. We may now

define the characteristic horizontal length scale by

L = = 129500 cm (4.70)

Therefore we have the velocity scale

uo = VgL = 224 cm/sec. (4.71)

and the time scale

T = v =L/g = sec.

(4.72)

60
The computational domain is shown in the Fig. 4.2. Here we see that the grid

density in the x-direction in the vicinity of the moving boundary is higher than in

the offshore region. From the grid line of 50 km to the offshore boundary, the grid

space increment is held at 1000 meters. Starting at this grid line, the grid space

increment decreases 10 percent successively for about 6 km until the grid space

increment of 200 meters is obtained. In the region of the moving boundary, the

grid space increment is fixed at 200 meters. The grid density in the y- direction is

unchanged, and the space increment is 2000 meters.

The moving boundary is simulated using the way of assuming a thin water layer

to cover the dry region all the time. The thickness of this water layer is set to be

5 cm. To decrease the violation of mass conservation resulting from this layer, the

friction with the Manning coefficient of 0.03 is considered in the region of thin water

layer. However, there is no friction considered in the main computational domain

since we did not consider the effect of friction in the derivation of the theoretical

solution.

A time step of At = 30 seconds is chosen. At t = 0, the fluid is quiescent. For

t > 0, the wave amplitude at the offshore boundary, r?*, can be given by

r7*(t) = t(t) 4L cm (4.73)

where r (t) is the dimensionless value of the wave amplitude which can be obtained

from the theoretical solution. After about three periods, we can obtain the state

numerical solution of the standing wave which will be compared with the theoretical

solution.

Figures 4.3 and 4.4 shows the comparisons between numerical and theoretical

wave profiles over the entire length of the basin at 7 time instant during half of a

period. Fig. 4.5 presents the amplification of these comparisons near the moving

boundary region. From Figs. 4.3, 4.4 and 4.5, it is seen that the agreement between

the numerical solution and the theoretical solution is good.

61
Figure 4.6 presents the comparisons of water surface displacement near the

moving boundary and the offshore boundary. Figure 4.7 presents the comparisons

of velocities at the same points as in Fig. 4.6. Both Figs. 4.6 and 4.7 show that the

agreement between the numerical solution and the theoretical solution is good.

62
-- THEORETICAL SOLUTION
A NUMERICAL SOLUTION

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

5

5

20 X

0 10 20 X 30 o0 50

Figure 4.3: Comparison between wave profiles as predicted by theory and the nu-
merical model ( t7 = rT77*w1/g, z = z'*w2/g )

5 0 10 20 X 30 l0 50

TIME: T1/6

------ -------- --- -------- -------- -----------------

5 0 10 20 X 30 40 50

TIME: TT/3

---__-I---_------ -__----I_--- --------------
IMT/

TIME: 0

- ------- ---------- ---- ------------------ -- -- - -

THEORETICAL SOLUTION
& NUMERICAL SOLUTION

TIME: 21T/3

- ------- ---------- ----------- ------------- --

1.5

1.0
0.5

0.0
-0.5
-1.0

-1.5
-i

1.5
1.0
0.5

0.0

-0.5
-1.0

-1.5

1.5

1.0
0.5
0.0

-0.5

-1.0

-1.5

5 0 10 20 30 40 50
X

0 10 20 30 40 5

Figure 4.4: Comparison between wave profiles as predicted by theory and the nu-
merical model ( 7 = 7'w2/xg, z = x*w2/Ig )

5 0 10 20 30 40 50
X

TIME: 57T/6

\

TIME: TT

5

-
-
-

-
-

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

-2 0 2 X 4

1.5

1.0.5

0.0 -

-0.5

-1.0 -

-1.5
-2

THE
A NUM

-_ \ TIME: 0

-
---1- w-----------------

-2 0 2 X 4 6

,,, TIME: TT/3

I I I

-2 0 2 X 4 6

TIME: 2Tr/3
\

1

---- -----,w --------------

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

64

ORETICRL
ERICAL SO
1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

TIME: TT

0 2 6----

Figure 4.5: Comparison between wave profiles near moving boundary as predicted
by theory and the numerical model ( 77 = q*wl2/2g, z = z*'w2/g )

SOLUTION
LUTION

TIME: Ir/6

...- -_------------------

-2 0 2 X 4 6

TIME: IT/2

I I I

-2 0 2 X 4 6

TIME: 5T/6

-
-\

-2 0 2 X 4 6

-.5w

T 2TT 3T 4WT

TIME

T 2T 3Tn 47T

TIME

X=7.5

- ----- THEORETICAL SOLUTION
NUMERICAL SOLUTION
I f I
T 2W 3W IMW
TIME

Figure 4.6: Comparison between theoretical and numerical solutions of water ele-
vation ( 7 = r*w'2/k2g,t = wt* )

1.0

0.5

0.0

-0.5

-1.0

1.0

0.5

-0.5

-1.0

1.0

0.5

0.0

-0.5

-1.0

X=15.

Tr 21T 3TT

TIME

)T 2TT 3T 4n1

TIME

X=7.5

STHEORETICRL SOLUTION
A NUMERICAL SOLUTION

TT 2Tr 3T 41
TIME

Figure 4.7: Comparison between theoretical and numerical solutions of velocity (
u = u*w/lg,t = wt* )

-0.2

-0.4

U. 4
X=15.
0.2

0.0

-0.2

-0.4

0.2

-0.2

-0.4

CHAPTER 5
APPLICATION TO LAKE OKEECHOBEE

In this chapter, the moving boundary numerical model developed in the previous

chapter will be applied to simulate wind driven circulation in Lake Okeechobee,
Florida.

As shown in Fig. 5.1, the western part of Lake Okeechobee is the grass region

where the water depth is about 30 to 100 centimeters. During a seiche induced by

a hurricane, this region may be inundated or dried due to the temporal variation
of water level. The water depth outside the grass area is about 4 meters on the

average.

A 31 x 50 grid is constructed as shown in Fig. 5.2. The north to south grid

spacing is uniform, but the east to west grid spacing is variable with smaller spacing

in the grass area. The north to south grid spacing is 1120 meters. The maximum
east to west grid spacing is 2240 meters and the minimum is 1120 meters.

Water depth values at the grid points are obtained by adding 1.2 meters to the

numbers given in the 1987 contour map of Lake Okeechobee which is for low water

conditions.

The wind shear stress acting on the lake surface can be calculated by the fol-

lowing formulae

Tz=PaCda W+W- W (5.1)

r=p. = Cd W/ +W2 W, (5.2)

where r, and ry denote respectively the wind shear stresses in the x- and y-directions,

p, is the density of air, W, and Wy are the wind speeds in the x- and y-directions,

respectively, and Cda is the wind shear stress coefficient. The value of Cda can be

Figure 5.1: The sketch of Lake Okeechobee

Figure 5.2: Computational grid

I A

S I I T I I T

69
obtained from the following formulation proposed by Garrett (1977)

Cd, = 0.001 x (0.75 + 0.00067 x W_ + W2) (5.3)

where the unit of W, and W, is cm/sec.. The maximum value of Cda is 0.003. The

bottom friction can be calculated from the Manning-Chezy friction formulation.

As discussed previously, the computer programming becomes very complicated

when using the weir formulation to simulate the two-dimensional moving boundary

problem. Thus the method proposed by Reid and Bodine (1968) is not applied to

study the wind driven circulation of Lake Okeechobee. The method to simulate the

moving boundary problem with a thin water layer on the dry land area is employed

here. To illustrate the significance of the moving boundary to the actual problem,

the simulations will be conducted with and without the moving boundary. For the

moving boundary case, the grass area is included in the computational domain and

the thickness of the thin water layer on the dry area is assumed to be 10 cm. For

the fixed boundary case, the grass area is taken out of the computational domain

and a vertical wall is assumed to be located at the edge of the grass region. Thus

the boundary between the grass area and open water can be considered as a closed

boundary.

A numerical time step of 3 minutes is used here. The Coriolis parameter, f, in

Lake Okeechobee is 0.73 x 10-4 sec.-1 (Schmalz,1986).

A spatially uniform wind from the east to the west is applied over the lake

surface for 9 hours. The wind speed is 20 m/sec.. After the wind stops, a seiche

in the lake will be produced. When the wind is applied over the lake, the Manning

coefficient is chosen to be 0.02 ( Schmalz, 1986). But after the wind stops, the

Manning coefficient is taken as 0.005 in order that the seiche can last for a longer

time for studying the shoreline motion.

The steady state of wind-driven circulation in the lake is reached after the wind

is applied for 7 hours. Figures 5.3 and 5.4 show the circulations for the moving

MIND
150 CM/SEC.

Fi/gue .,< d e c i .-- .
'1 t

Figure 5.3: Wind driven circulation with moving boundary

r n HiINDa
150 CM/SEC, ).\ -. ,\

.: ., ,,, \ \

<'..- / ; '

\ '.-' ',, 5 /

71
boundary case and the fixed boundary case respectively after the wind has been

applied for 8 hours. It is seen that the circulation for the two cases are quite

similar, although the velocities in the moving boundary case appear to be larger

than those in the fixed boundary case. This can be explained by the fact that

the area of the lake in the moving boundary case is greater than that in the fixed

boundary case.

Figure 5.5 shows the variation of the wind speed with time. Figure 5.6 shows

the motion of the western shoreline during the first cycle of seiche oscillation in the

lake after the wind has stopped.

As mentioned previously, mass conservation is violated by assuming that a thin

water layer exists on the dry area at all times. Figure 5.7 shows the extent of the

mass conservation violation during the entire computational period. It is found that

an extra mass of 0.2 percent of the total water mass in the lake is induced due to

the consideration of the shoreline motion. Figure 5.8 presents the variation of the

total dry area in the lake with the time.

Figure 5.9 shows the comparisons of the transient variations of the water surface

elevation at points A, B and C between the moving boundary case and the fixed

boundary case. The locations of points A, B and C in the lake are shown in the

Fig. 5.1. Figures 5.10 and 5.11 present the comparisons for the velocities U and V,

respectively, between the moving boundary case and the fixed boundary case. From

Figs. 5.9, 5.10 and 5.11, it is seen that there is a phase lag between the results for

the moving boundary case and the fixed boundary case.

Figures 5.12 and 5.13 present the transient variation of the bottom friction in

the x- and y-directions at point A, respectively.

20 ----- -----------------------
0.

LO
=n 20.
z:

0 .
(0

0 -20.
z

-40. --- -- --- -- --- -- 1 -- -- --- -
0 5 10 15 20 25
TIME (HR)

Figure 5.5: Variation of wind speed with time

Figure 5.6: Locations of the moving boundary at four different time

0.5

O. -

0.3

V
0.2

0.1

0.0
0 5 10 15 20 25
TIME (HR)

Figure 5.7: Extra mass introduced by considering the moving boundary

5.0

4.0 -

3.0
x100

2.0

1.0

0.0
0 5 10 15 20 25
TIME (HR)

Figure 5.8: Transient variation of dry area

-- SOLUTION WITH MOVING BOUNDRRT
-------- SOLUTION WITH FIXED BOUNDARY
200
POINT A
150

100

50 -

-50

-100
0 5 10 15 20 2!
TIME (HR)

5 10 15 20 25
TIME (HR)

5 10 15 20 25
TIME (MR)

Figure 5.9: Comparisons of water surface elevation with moving boundary and fixed
boundary

U

r
z 50
0

0
.J
UJ
e -50

-100
O

- UU

N

> 0
UJ
-I
UJ
a- -50
UJ
I--O
cc
-100
0

75

--- SOLUTION WITH MOVING BOUNDRRT
---------. SOLUTION WITH FIXED BOUNDRART
200
POINT R
UJ
100

0 --------------------------------

U -100
-j
-200
0 5 10 15 20 25
TIME [HR)

200
POINT B
U
100
. -..

L -100
o
-J

-200
0 5 10 15 20 25
TIME (MR)

POINT C

S100

0
100
0 o -/- -- -
I-

L -100
0
-J

-200
0 5 10 15 20 25
TIME (MR)

Figure 5.10: Comparisons of velocity U with moving boundary and fixed boundary

----- SOLUTION WITH MOVING BOUNDARY
--------- SOLUTION WITH FIXED BOUNDRRT
200
POINT R

C 100

S-100

UJ
-200
0 5 10 15 20 25
TIME (HR)

200
S POINT B
UJ
L 100

>-

S-100-

LU
-200 'iiI
0 5 10 15 20 25
TIME (MR)

200
-POINT C
U
S100-

> 0
s-

S-100
-J
LU
-200
0 5 10 15 20 25
TIME (HR)

Figure 5.11: Comparisons of velocity V with moving boundary and fixed boundary

10

0

-10

0 5 10 15 20 25
TIME (HR)

Figure 5.12: Transient variation of bottom friction in the x-direction at point A

5 10 15 20 2
TIME (HR)

Figure 5.13: Transient variation of bottom friction in the y-direction at point A

CHAPTER 6
CONCLUSIONS

6.1 Conclusions

The objective of this study is to develop a two-dimensional moving boundary

numerical model which can predict the hydrodynamics in lakes, estuaries and shal-

low seas with the effect of a moving boundary The study consists of derivation,

verification, and application of the model.

The finite difference technique is used to develop the model in terms of a non-

uniform rectangular grid. The governing equations, which are vertically-integrated

Navier-Stokes equations of water motion, are solved using the method of fractional

steps. The transition from one stage of the computation to the next is divided into

three steps which are advection, diffusion, and propagation. Different numerical

schemes are employed for each computational step. The method of characteristics

is used for the advection, the ADI method is applied to the diffusion, and the

conjugate gradient iterative method is used for the propagation. The numerical

schemes used in the model are implicit so that there is no limitation for the choice

of the numerical time step.

Two methods for simulating the moving boundary problem are discussed in

this study. The first, which was proposed by Reid and Bodine (1968), is examined

briefly. It is found that it is difficult to determine the values of empirical coefficients

in the weir formulation since they are very site-dependent. It is also found that the

computer program for simulating the motion of the boundary is very complicated

for two-dimensional problems. The second, which was proposed by Benque et al.

(1982), is studied in detail. The advantage of this method is that the computer

78

program for simulating two-dimensional moving boundary problems is very simple.

The disadvantage is that the mass conservation is violated slightly.

For the verification of the model, four analytical solutions are used to compare

with the numerical solutions. From these comparisons, it can be concluded that the

consistency and the accuracy of the model are acceptable. It is also found that the

method of fractional steps is a powerful means of solving the complicated problems

in several variables although its consistency has not been completely proved. In

order to validate the model's ability to simulate the motion of the boundary, wave

propagation on a linearly sloping beach is studied theoretically and numerically. It
is found that the moving boundary model developed using this method can simu-

late moving boundary problems reasonably well although the mass conservation is

violated slightly.

The model is applied to the wind-driven circulation in Lake Okeechobee, Florida.

Comparison is made between the results obtained from the moving boundary model

and the fixed boundary model. From this comparison, it is seen that the boundary

motion can not be neglected when studying wind-driven circulation in Lake Okee-

chobee. It is found that the violation of mass conservation in the moving boundary

model can be negligible.

6.2 Future Study

In the moving boundary region, the water depth is usually very small and the

water motion is controlled by the bottom friction. In this case, the Manning-

Chezy formulation cannot give a correct estimation of the bottom friction since it
breaks down when the water depth approaches zero. Therefore, basic hydrodynamic

research is needed in this area. For example, Yih (1963) investigated the stability

of liquid flow down an inclined plane and Melkonian (1987) studied nonlinear waves

in thin films.

When using a uniform or non-uniform rectangular grid to approximate the com-

80
plex geometry of water bodies, a large number of grid points is generally required.

As a result, the computational cost is increased greatly. To avoid this, it is useful

to modify this model to work in a curvilinear grid.

In addition, we can use the method of fractional steps to develop a three-

dimensional numerical model.

APPENDIX A
APPLICATION OF THE CONJUGATE GRADIENT METHOD TO THE

PROPAGATION STEP

In this appendix, the algorithm of the conjugate gradient method for the so-

lution of simultaneous linear algebraic equations will be presented briefly and the

procedure of the application of this method to the propagation step will be described

in detail.

A.1 Conjugate Direction Method

In order to understand the algorithm of the conjugate gradient method, the

conjugate direction method needs to be reviewed briefly since the conjugate gradient

method is a special case of the conjugate direction method. It is assumed that there

is a solution vector h existed to a system of simultaneous linear algebraic equations,

Ax=k (A.1)

in which A is an N x N matrice of coefficients, x is an N x 1 vector of unknowns,
and k is an N x 1 vector of constants. In the following description, the symbol (p, q)

presents the inner product of the vectors p and q, or the value of pTq.

Let us suppose that a set of N "A-conjugate" or "A-orthogonal" vectors {pi},

i = 0,1,2,***,N 1, is available to us. This means that the inner product

(APi,pj) = 0 where i # j. If A is positive definite, then (Ap,,pi) > 0. In this
case, since the vectors {p,} are necessarily linearly independent and span the N-

dimensional space, the solution vector h can be written as

h = coPo + cll + C2P2 + *.. + CN-lPN-i (A.2)

82
where {ci} are coefficients. If we can determine {c,}, the solution h can be quickly

calculated. Since

(k, P) = (Ah, p) (A.3)

(Ah, p,) = co(Apo,Pi) + ci(Api,pi) + + ci(Api,p i) + *

+ cN-I(APN-1,Pi) = c,(Ap,,pi) (A.4)

we get
S(k,p) (A.5)
Ci = (Ap,p) (A.5)
(Api, pi)
Consequently

(k, o) (k, pi) (k, PN-1) (.6)
h = Po + -P + + PN-1 (A.6)
(Apo,Po) (Api,pi) (ApN-1,PN-1)
This computation of h defined by equation ( A.6) also can be described as the

following iterative scheme

x1 = 0opo (A.7)

zX+l = zi + Pip (A.8)

where

P = (k, p) (A.9)
(Api,p,)

and Po, Pi,*" ,PN- is the given set of N A-conjugate vectors. The iteration can
be terminated when XM = h where M < N.

A.2 Conjugate Gradient Method

In conjugate gradient method, a particular set of A-conjugate vectors, {pi}, is

developed and a solution to equation (A.1) formed in terms of these. Introducing

residual vectors,

ri+l = k A xi+i

(.A.10)

83
Beckman (1958) used Gram-Schmidt Orthogonalization procedure to obtain the
A-conjugate direction vectors {pi} They are expressed as

Pi+1 = ri+1 + ,ipp (A.11)

in which

r, = I r+1-2 (A.12)
I ri 12

where I r+i I and I ri represent the length of the vectors ri,+ and r,. The fun-
damental conjugate gradient iterative procedure leading to a solution of equation
(A.1) can be defined by following formulae:

Po = ro = k A xo (A.13)

(p(Pi, r) (A.14)
S (p,, Api)

X+1 = Xi + a;ip (A.15)

ri+l = k A zi+l (A.16)

A = -, 12 (A.17)
I ri I

Pi+i = ri+1 + ?iP, (A.18)

After M iterations, with M < N, XM will be equal to the solution h if all conputa-
tions are done with no loss of accuracy.
A.3 Application to the Propagation Step

The key to the prapagation step presented in Chapter 2 is to solve for AZi (z, y),
AZ2(z,y) and q(z,y) which satisfy

Z_ 8(h" ') a(AZIaz ) I a (Ell AZ)
Az cc 2{ ax 9a )+ Z = f- q (A.19)
2g(At)2 ax + gAt dx

84
AZ2_ a(hnaZ2) a(AZ2W) O a( 2AZ2) -f+q
_C,2 a2 + + h 2-- + q (A.20)
2g(At)2 ay ay gt ay

AZi(z, y) = AZ2(, y) (A.21)

and subject to the open boundary and fixed boundary or moving boundary condi-
tions. In order to present the way to solve above equations clearly, we need to write
the above equations in the matrix form as

[A,] [AZ1] = [I(] [B], [q] (A.22)

[Ay] [AZ2] = [f2] [By] [q] (A.23)

[By] [AZ1] + [B.] [AZ2] = 0 (A.24)

in which [A,] and [A,] are symetrical coefficient matrix of Eqs. (A.19) and (A.20),
respectively, [AZ1] and [AZ] are vectors of unknowns, [fil and [f2] are vectors of
known terms at the right sides of Eqs. (A.19) and (A.20), respectively, [q] is the
vector of coordinate terms, and [B.] is identity matrice.
It should be noted that the arrangement of elements in [AZI] is different from
that in [AZ2]. Equation (A.21) is for the same grid point. Therefore, a matrice
needs to be defined to make the same arrangement of elements in [AZ1] and [AZ2].
[By] is defined to do it. Actually for large grid system, it is difficult to obtain the
explicit form of [By]. Fortunately it will be seen in later description that we do not
really need to know [By].
Eqs. ( A.22), (A.23) and (A.24) can be rewritten as

S]= B [q] (A.25)

B 0 (A.26)
Byz Az IZ =

Denoting A = A, A = AZ, B = B, and [q] = q,
Eqs. (A.25) and (A.26) become

A Z = f B q (A.27)

BT AZ = 0 (A.28)

where BT is transposed matrice of B. From Eq. (A.27), we can get

AZ = A-' (f Bq) (A.29)

Substitution of AZ into Eq. (A.28) yields

(BT A- B) q =BT A-' f (A.30)

where (BTA-1B) is a symetric matrice.
The conjugate gradient method is used to solve Eq. (A.30) for q. The iteration
on q consists of looking for qm+I by given qm where m means the mth iteration.
The residual vector is

r"+1 = BT A-f (BT A-1 B) qm+l

= B (A- f A-1 B qm+l)

= BT'AZm+1

= AZ+' _AZ+1 (A.31)

Coefficient 6m can be calculated by

I rm+l 2
Srm 12

= ZE[Az+) AZm ] z2/ E[AZ ) AZM,,)2 (A.32)
I,J ij
The direction vector is

pm+1 = rm+1 + Om Pm

(A.33)

The coefficient a, is

(pm, rm)
S (pm,BT A-1B pm)

[pm]T [BT AZm]
[pm]T [BT A-' B pm A
Defining A-'Bpm = Hm, therefore we have

A Hm = Bpm (A.35)

Comparing Eq. (A.35) with Eq. (A.27), it is seen that H" can be obtained by
using the double-sweep method to solve

Hx, (h"nl{) a+(H,)} a 8(tnr" H,)
2g ) a { + e} + g" = B, pm (A.36)
2g(At)2 ax a gAt Ox

H a(h"i) (H, + a 8(Vn+2/3 r2)
2 -( ao(%{ + ) } +a = By pm (A.37)
2g(At)2 ay ay gAt ay

where [ ] = H. Consequently

am = [pn]T (BT Zm
[pm]T [BT Hm)
I,J I,J
= .[A2Zi,,)- ArZ, i)] P,/)/ E[Hmir ) H,\ij] pmj (A.38)
ij Ijd
The iteration procedure for q can be sumarized as:
(1) Let the final value of q at the previous time step be the initial approximation
to the solution of q at the new time step.
(2) Apply the double-sweep method to solve Eqs. (A.19) and (A.20) for AZ
and AZ0.
(3) Calculate the residual vector by ro = AZ? AZo.
(4) Let po = ro
(5) Use the double-sweep method to solve Eqs. (A.36) and (A.37) for H? and

87
(6) Compute the coefficient ao using Eq. (A.38).

(7) Advance q by q' = qO + ao pO

(8) Determine successively

rm+1 = AZ1+ AZm+1 (A.39)

I,J I,J
S= -[AZit) AZJ)]2/ /[AZi,) Az, (A.40)
i,i i'j

pm+l = rm+l + m pm (A.41)

I,J I,J
am+1 = AZ ) AZPmn)+ / E[Hi'+) H+l P+)1 (A.42)
l -j) ij) 2(i,)J (ij)
dj *,i

qm+2 = qm+1 + am+1 pm+1 (A.43)

(9) Repeat Step 8 with m +1 replacing m and continue until m = N- 1 or until

the residual vector becomes sufficiently small, whichever condition may be satisfied

first.

(10) Let AZ = (AZ1 + AZ2)/2

APPENDIX B
DERIVATION OF Z2 AND U2

The governing equations are

+ = = o (B.1)
at ax

+ gh a 2 k2 {sin[2k(l z) + 2wt]
at + a 8h cos2 kl
+ sin[2k(l z) 2wt] + 2 sin 2k(l z)} (B.2)

and the boundary conditions are

U2(, t) =0 (B.3)

Z2(0, t) =0 (B.4)

Eliminating U2 from Eqs. (B.1) and (B.2), we obtain

a2Z2 2Z2 a'c2 k2
2 4cos k{cos[2k(i x) + 2wt]
at z 4h cos2 ki
+ cos(2k(I z) 2wt] + 2k cos 2k(l x)} (B.5)

Since we want to obtain the solutions which vary with the time, the third term in
the right hand side of Eq. (B.5) can be neglected. Let the general solution of Eq.
(B.5) take the form of

Z2 = F{sin[2k(l z) + 2wt] + sin[2k(l x) 2wt]}

+ G{cos[2k(l x) + 2wt] + cos[2k(l z) 2wt]} (B.6)

89
where F and G are fountion of z only. Substitution of Z2 into Eq. (B.5) yields

[-c'F" c24kc']{sin[2k(1 x) + 2wt] + sin[2k(l x) 2wt]}

+ [c'4kF' c'G"]{cos[2k(l x) + 2wt] + cos[2k(l z) 2wt]}
a2 c2k2
=4hc2 k {cos[2k(l z) + 2wt] + cos[2k(l ) 2wt]} (B.7)

Obviousely we get
c2F" c24kG' = 0 (B.8)
a2c2k2
c24kF' C2G" = c (B.9)
4h cos2 kl
Integrating Eq. (B.9) with respect to z and set the integrating coefficient as zero,

we have
a2k2
4kF G' = (B.10)
4h cos2 kl
Eliminating G' from Eqs. (B.8) and (B.10), we obtain

F" a2k2
+ 4kF = x (B.11)
4k 4hcos2kl

A general solution for F of Eq. (B.11) is

F, = C1 cos 4kz (B.12)

where C1 is a constant to be determined from the boundary conditions. A particular

solution of Eq. (B.11) is

ac2k
Fv = k (B.13)
16h cos2 kl B.13)

Therefore

F = F,+F,
a2 k
= C1 cos 4kx + 6 k (B.14)
16h cos2 kl

Substituting F into Eq. (B.8), we can obtain

G = C1 sin 4kx + C2 (B.15)

90
where C2 is a integration constant. Plugging F and G into Eq. (B.6) and simplifying

it, we have

a2k
Z = [CI cos 4kx + zs2 k] sin 2k(1 x) cos 2wt
8h cos2 kl
+ [C1 sin 4k + C2] cos 2k( z) cos2wt (B.16)

From the boundary condition (B.3), we obtain

a2kl
[CI cos4kl + s kl ](-2k) + C14k cos 4kl = 0 (B.17)

Thus

a2ki
Ca (B.18)
S= 8h cos2 k cos 4kl (

Similarly from boundary condition (B.4), we can get

a2kl
C2 = a _tan 2ki (B.19)
h cos2 kl cos 4k .

Substituting C1 and C2 into Eq. (B.16) and simplifying it, we consequently obtain

a2k I
Z2 = [zsin2k(l z) + sin2k(l + z)
8h cos2 ki cos 4ki
Stan 2kl cos 2k( z)] cos 2wt (B.20)
cos 4kl

Substituting Z2 into Eq. (B.1) and integrating with respect to z, we have

a2w 1
U2 = s [cos 2k(1 ) + sin2k(I z)
8hcos2kl 2k
1 1
cos 2k(1 + z) + tan 2k sin 2k(l- z)]sin 2wt (B.21)
cos 4kl cos 4k

APPENDIX C
PROGRAM LISTING

C.1 Flow Chart

C.1.1 Main Routine: T2D

C.1.2 Subroutine: PROP

C.2 Program listing

C.2.1 Description of Parameters

KIO: Maximum number of grid points in the x-direction;

KJO: Maximum number of grid points in the y-direction;

KX1(J), KX2(J), KY1(I), KY2(I): Grid numbers for the boundary's location;

CN: Manning coefficient;

F: Coriolis parameter;

G: Gravitational acceleration;

MILLISECOND	CLASS.METHOD	MESSAGE
0	sobekcm_page_globals.constructor
0	sobekcm_page_globals.constructor	Application State validated or built
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.constructor	Navigation Object created from URI query string
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.display_item	Retrieving item or group information
0	sobekcm_page_globals.get_entire_collection_hierarchy	Retrieving hierarchy information
0	sobekcm_assistant.get_entire_collection_hierarchy
0	cached_data_manager.retrieve_item_aggregation
0	cached_data_manager.retrieve_item_aggregation	Found item aggregation on local cache
0	item_aggregation_builder.get_item_aggregation	Found 'all' item aggregation in cache
0	system.web.ui.page.page_load (ufdc.page_load)
0	sobekcm_page_globals.constructor.on_page_load
0	html_echo_mainwriter.add_style_references	Adding style references to HTML
0	html_echo_mainwriter.add_text_to_page	Reading the text from the file and echoing back to the output stream
120	html_echo_mainwriter.add_text_to_page	Finished reading and writing the file