UFL/COEL88/008
A TWODIMENSIONAL FINITEDIFFERENCE
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.
LiHwa 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 TWODIMENSIONAL 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 Nonlinear 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 xdirection at point
A . . . . . . . . .. .. 77
5.13 Transient variation of bottom friction in the ydirection 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 TWODIMENSIONAL FINITEDIFFERENCE 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 NavierStokes 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 threedimensional
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 threedimensional 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 multidimensional 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 sitedependent. 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 twodimensional finitedifference 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 twodimensional 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 TWODIMENSIONAL TIDAL CIRCULATION
MODEL
In this Chapter, a twodimensional implicit tidal circulation model is devel
oped using the method of fractional steps. The verticallyintegrated momentum
and continuity equations, which govern the twodimensional tidegenerated 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 threedimensional NavierStokes 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 righthanded 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 unitwidth discharges in the x and y
directions, respectively, u(z, y, t) and v(z, y, t) are the verticallyaveraged 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 ydirections, respectively, and rb, and rb, are the bottom
shear stresses in the x and ydirections. 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=ZZf
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(KU)
 + 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:
xadvection
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
yadvection
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 xadvection. This scheme is implicit
and unconditionally stable.
For the diffusion step, an alternating direction implicit method is used which
leads to
xsweep
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
ysweep
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 unitwidth discharge during the
computation of the diffusion step. Yanenko (1971) showed that the ADI scheme
was absolutely consistent and unconditionally stable for the twodimensibnal 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)(ghz)" 
'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
hAZ)} = 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 firstorder 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 + hAZ 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( )"V2/3A (2.46)
C2 h3
V"+1 = N{V+2/3 + VAZ 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 1a au Ma aU. ,, ,(h ) Ma(UU2+v)"+2/3
f = )" ( )n2/ +Ma 9X + C (2.52)
gAt Xs gAt az xz aZ xz
1a 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(ijl) + PC(i+,i) + Pz(i+l,y1))/4
(2.56)
xI
I
Figure 2.2: Schematic of finite difference mesh with variable rectangular cells
J2 t
Jl
I1 I
1*1
 ZGRID
 U
 A ( 1+
11 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,(i1) AX,(i)AY,(y)
4AX,(4 )AY,(y) P(i1,y+l) + 4AX(_..)AY() Pu(ly) (2.57)
Pz(i+/2,y) = 2AY (i) +,)+ 2A )Pu(ij) (2.58)
2AX,(i1) 2AXu(iI)
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(ii,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(ii,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
xadvection
n+1/6 n+1/6 n+1/6
i U + u"n i  = 0 (2.65)
At AX(i1) + AX,(i)
n+1/6 n n+1/6 n+1/6
vi, ,i + u Vi+1,) vii4 = 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(y1) + AYu(y)
where the subscripts i and j correspond to the (U,V) grid. After xadvection
and yadvection, 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:
xsweep
u  Un+1/
nAt
, v+1/3
,y i,j
K U lUi y U
K U +1.3i '1, ._* U1* I
AX',() AX( AX(OI)
Ay,(_i) AyX () AY.(I)
K s nis +/ v
iI i
AX,(i) X,(i)AX,(I)
K Vsi, + 1 iV iy id
AY,([1) AYU(y) Ayu,(1)
= 0
(2.71)
.(2.72)
ysweep
AX(i) AK,(i)
U:..j U:.j
14 Ui*Ij
r 1)
yadvection
(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,j1
Ay) AYu() AY.1)
= 0
K V*. V*. V* V.*,
AX.)W AX,() AX,(i1)
S "+2/s V+2/3 Vn2/3 T+2/3
S['ij+1 i,3 i i. ij1
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(i1,y) Z 
AiX[2y l + Xj A x, + A Z1(i+/2j) X
gAtAXu(il) 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,l1
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,, Ui1/2,jy a l+1/2,jy l/, U
gAt AXu(i1) gAt AXu(i1)
(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,31/2 Ay.(,)
a
"+ Yay [(QV rY),1 +1/2 (7 ,Y)!i1/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 Zgrid and the
discharges at the Zgrid 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, il/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 ( noslip condition ) or satisfy the zeroslope 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 twodimensional 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 twodimensional 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 onedimensional 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 xdirection, 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
S50
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
S50
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
d5
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 Nonlinear Theoretical Solution
When nonlinear terms are included in twodimensional shallow water equations,
it is impossible to obtain an analytical solution. In onedimensional 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 onedimensional 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 onedimensional 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

UI
.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 twodimensional 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[ky +ik(l z) iwt]
w
flk
+ AR exp[(b y) ik(l z) iwt]
+ {C,[cosKi,(by)+  sinKi,(by)]
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(by)
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 xdirections, 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
V50
.'>
s
h
S50
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 11 1 1 1 1 1
66 72 78 81 90 91
TIME (HR)
37
the flow velocities u in the xdirection at three points in the middle axis of the basin.
Figure 3.10 shows the flow velocities v in the ydirection 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 twodimensional shallow water equations can be
written as
+ gh a+ =0 (3.48)
at 9x p
0+ gha + = 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 + c2 +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{F2 [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
S100 I I I I fI
60 66 72 78 84 90 96
TIME (HR)
100
X=30KM. T=15KM
S50
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 ydirection, 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 verticallyintegrated shallow water problems is the ManningChezy 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 stairstep 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 broadcrested 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 912, 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 sitedependent.
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 ManningChezy 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 finiteelement 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
elementsplitting 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 stairsteps. 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 NavierStokes 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 noslip or freeslip 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 ZdZ (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 ManningChezy 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 twodimensional 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 NavierStokes 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 ManningChezy 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
i1/2 = i1 hi1 + (1 (4.8)
AZi+1/2 = AZ.i + (1 7i)AZi+i
AZi1/2 = fii AZiI + (1 'Y.i)AZ,
and
hj+1/2 = yj hi + (1 j)hj+l
h*1/2 = yl hj_ + (1 yiI)hj (4.9).
AZ+1/2 = j AZ, + (1 ,)AZ+1 (4)
Aj1/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 ManningChezy 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 xdirection, but in Eq. (4.19) they refer to
flow in the ydirection. 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+2ct
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 finitedifference 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 finitedifference
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 xdirection 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 ydirections,
p, is the density of air, W, and Wy are the wind speeds in the x and ydirections,
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 ManningChezy friction formulation.
As discussed previously, the computer programming becomes very complicated
when using the weir formulation to simulate the twodimensional 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 104 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 winddriven 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 ydirections 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
IO
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
S100
UJ
200
0 5 10 15 20 25
TIME (HR)
200
S POINT B
UJ
L 100
>
S100
LU
200 'iiI
0 5 10 15 20 25
TIME (MR)
200
POINT C
U
S100
> 0
s
S100
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 xdirection at point A
5 10 15 20 2
TIME (HR)
Figure 5.13: Transient variation of bottom friction in the ydirection at point A
CHAPTER 6
CONCLUSIONS
6.1 Conclusions
The objective of this study is to develop a twodimensional 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 verticallyintegrated
NavierStokes 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 sitedependent. It is also found that the
computer program for simulating the motion of the boundary is very complicated
for twodimensional 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 twodimensional 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 winddriven 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 winddriven 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 nonuniform 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 "Aconjugate" or "Aorthogonal" 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 + *.. + CNlPNi (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) + *
+ cNI(APN1,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, PN1) (.6)
h = Po + P + + PN1 (A.6)
(Apo,Po) (Api,pi) (ApN1,PN1)
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 Aconjugate 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 Aconjugate 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 GramSchmidt Orthogonalization procedure to obtain the
Aconjugate direction vectors {pi} They are expressed as
Pi+1 = ri+1 + ,ipp (A.11)
in which
r, = I r+12 (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 (BTA1B) 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 Af (BT A1 B) qm+l
= B (A f A1 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 A1B 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 doublesweep 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 doublesweep 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 doublesweep 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 xdirection;
KJO: Maximum number of grid points in the ydirection;
KX1(J), KX2(J), KY1(I), KY2(I): Grid numbers for the boundary's location;
CN: Manning coefficient;
F: Coriolis parameter;
G: Gravitational acceleration;
