Citation |

- Permanent Link:
- http://ufdc.ufl.edu/UF00076135/00001
## Material Information- Title:
- A two-dimensional finite-difference model for moving boundary hydrodynamic problems
- Series Title:
- UFLCOEL
- Creator:
- Liu, Yuming, 1965- (
*Dissertant*) University of Florida -- Coastal and Oceanographic Engineering Laboratory Sheng, Y. Peter (*Thesis advisor*) Kirby, James T. (*Reviewer*) Sheppard, D. Max (*Reviewer*) - Place of Publication:
- Gainesville, Fla.
- Publisher:
- University of Florida
- Publication Date:
- 1988
- Copyright Date:
- 1988
- Language:
- English
- Physical Description:
- vii, 128 leaves : ill. ; 28 cm.
## Subjects- Subjects / Keywords:
- Advection ( jstor )
Boundary conditions ( jstor ) Friction ( jstor ) Modeling ( jstor ) Shallow water ( jstor ) Shorelines ( jstor ) Subroutines ( jstor ) Surface water ( jstor ) Velocity ( jstor ) Water depth ( jstor ) Coastal and Oceanographic Engineering -- Dissertations, Academic -- UF Coastal and Oceanographic Engineering thesis M.S Hydrodynamics -- Mathematical models ( lcsh ) Navier-Stokes equations -- Numerical solutions ( lcsh ) Numerical analysis ( lcsh ) - Genre:
- bibliography ( marcgt )
non-fiction ( marcgt ) theses ( marcgt ) - Spatial Coverage:
- United States -- Florida -- Lake Okeechobee
## Notes- Abstract:
- To predict the hydrodynamics of lakes, estuaries and shallow seas, a two dimensional numerical model is developed using the method of fractional steps. The governing equations, i.e., the vertically integrated Navier-Stokes equations of fluid motion, are solved through three steps: advection, diffusion, and propagation. The characteristics method is used to solve the advection, the alternating direction implicit method is applied to compute the diffusion, and the conjugate gradient iterative 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.
- Thesis:
- Thesis (M.S.)--University of Florida, 1988.
- Bibliography:
- Includes bibliographical references.
- General Note:
- Typescript.
- General Note:
- Vita.
- Funding:
- This publication is being made available as part of the report series written by the faculty, staff, and students of the Coastal and Oceanographic Program of the Department of Civil and Coastal Engineering.
- Statement of Responsibility:
- by Yuming Liu.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- All rights reserved, Board of Trustees of the University of Florida
- Resource Identifier:
- 20117334 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

UFL/COEL-88/008 A TWO-DIMENSIONAL FINITE-DIFFERENCE MODEL FOR MOVING BOUNDARY HYDRODYNAMIC PROBLEMS by Yuming Liu Thesis 1988 ACKNOWLEDGEMENTS I would like to express my sincere appreciation and gratitude to my advisor Dr. Y. Peter Sheng, Professor of Coastal and Oceanographic Engineering, for his continuous guidance and encourgement throughout this study. I would also like to extend my thanks and appreciation to my thesis committee members, Dr. James T. Kirby, Associate Professor of Coastal and Oceanographic Engineering, and Dr. D. Max Sheppard, Professor of Coastal and Oceanographic Engineering, for their patience in reviewing this paper. I would like to thank Dr. Yixin Yan, Dr. Tienshun Wu, Dr. Weichi Wang, Dr. Li-Hwa Lin, and Steven Peene for their help and suggestions in this study. Finally, I would like to express my deepest appreciation and gratitude to my financial sponsor, K. C. Wong Education Foundation Ltd. of Hong Kong, for offering me the opportunity to study at the University of Florida. TABLE OF CONTENTS ACKNOWLEDGEMENTS ................... LIST OF FIGURES ........................ ABSTRACT ........................... CHAPTERS 1 INTRODUCTION ...................... 2 DEVELOPMENT OF A TWO-DIMENSIONAL TIDAL M ODEL ........................... 2.1 Governing Equations .................. 2.2 Numerical Scheme ................... 2.3 Grid System ....................... 2.4 Finite Difference Approximation ........... 2.5 Boundary and Initial Conditions ........... v . . vii CIRCULATION 2.6 Consistency and stability 3 COMPARISONS BETWEEN THEORETICAL AND NUMERICAL SOLUT IO N S . . . . . 3.1 Comparison With a Linear Theoretical Solution ............ 3.2 Comparison With a Non-linear Theoretical Solution ......... 3.3 Comparison to a Theoretical Solution with Coriolis Effect ...... 3.4 Comparison to a Theoretical Solution with Friction Effect 4 NUMERICAL SIMULATION OF FLOW OVER TIDAL FLATS . 4.1 Properties of a Moving boundary . . . 4.2 Past Study . 4.3 Numerical Treatment of a Moving Boundary . . 4.3.1 Method to Treat a Moving Boundary with the Weir Formulation . 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 . . CA 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 water surface elevation. .. .. .. ... ... ... ... ... ....24 3.4 Comparison between theoretical and numerical solutions for velocity .. .. .. ... ... ... ... ... ... ... ... ....25 3.5 Comparison of theoretical solution to numerical results with different time steps .. .. .. .. .... ... ... ... ... .....26 3.6 Comparison between theoretical and numerical solutions for water surface elevation with nonlinear effects. .. .. .. .. ... ..29 3.7 Comparison between theoretical and numerical solutions for velocity with nonlinear effects .. .. .. .. ... ... ... ... ..30 3.8 Comparison between theoretical and numerical solutions for water surface elevation with coriolis effect. .. .. ... ... .....34 3.9 Comparison between theoretical and numerical solutions for velocity U with coriolis effect .. .. .. ... ... ... ... ....35 3.10 Comparison between theoretical and numerical solutions for velocity V with coriolis effect .. .. .. ... ... ... ... ....36 3.11 Comparison between theoretical and numerical solutions for water surface elevation with friction effect. .. .. ... ... .....39 3.12 Comparison between theoretical and numerical solutions for velocity U with friction effect .. .. .. ... ... ... ... .....40 4.1 DefiniLion sketch for wave propagation on a sloping beach 55 I 4.2 Computational grid ...... ......................... 59 4.3 Comparison between wave profiles as predicted by theory and the numerical model ( j7 = r7*w'10'g, x = z*wI2/g ) ....... .62 4.4 Comparison between wave profiles as predicted by theory and the numerical model ( ,7 = 7*7w2/02g, x X*w2/og ) ....... .63 4.5 Comparison between wave profiles near moving boundary as predicted by theory and the numerical model ( 77 = 77*w2/2g, X = x*w2g ) ............................... 64 4.6 Comparison between theoretical and numerical solutions of water elevation ( 7 = ,.*w2/02gt = wt- ) ............... 65 4.7 Comparison between theoretical and numerical solutions of velocity ( U = u*w/4g, t = wt* ) .................... 66 5.1 The sketch of Lake Okeechobee ..... .................. 68 5.2 Computational grid ...... ......................... 68 5.3 Wind driven circulation with moving boundary ............ 70 5.4 Wind driven circulation with fixed boundary .............. 70 5.5 Variation of wind speed with time ..................... 72 5.6 Locations of the moving boundary at four different time .. 72 5.7 Extra mass introduced by considering the moving boundary 73 5.8 Transient variation of dry area ..... .................. 73 5.9 Comparisons of water surface elevation with moving boundary and fixed boundary ................................ 74 5.10 Comparisons of velocity U with moving boundary and fixed boundary ....... ............................... 75 5.11 Comparisons of velocity V with moving boundary and fixed boundary ....... ............................... 76 5.12 Transient variation of bottom friction in the x-direction at point A ........ .................................... 77 5.13 Transient variation of bottom friction in the y-direction at point A ........ .................................... 77 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science A TWO-DIMENSIONAL FINITE-DIFFERENCE MODEL FOR MOVING BOUNDARY HYDRODYNAMIC PROBLEMS By YUMING LIU December 1988 Chairman: Dr. Y. Peter Sheng Major Department: Coastal and Oceanographic Engineering .To predict the hydrodynamics of lakes, estuaries and shallow seas, a two 'dimensional numerical model is developed using the method of fractional steps. The governing equations, i.e., the vertically integrated Navier-Stokes; equations of fluid motion, are solved through three steps: advection, diffusion and propagation. The characteristics method is used to solve the advection, the alternating direction implicit method is applied to compute the diffusion, and the conjugate gradient iterative 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 NavierStokes equations of fluid motion, using the finite difference technique. Many numerical schemes have been proposed to solve the shallow water equations in the development of numerical models. They all have advantages and disadvantages. 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 limitation, but it requires solving the matrix equations. For two- and three-dimensional flows, it is very difficult to overcome the computational difficulties resulting from the sheer size of the matrices. -The alternating direction implicit method, or ADI method, avoids solving complex matrix equations, but can obtain accurate solutions only for Courant numbers less than 5 (Weare, 1979). Furthermore the ADI method is not applicable to three-dimensional problems (Yanenko, 1971). The method of fractional steps developed by Yanenko (1971), on the other hand is known to be an effective method for solving complicated multi-dimensional problems in several variables. In this scheme, the computation from one time level to the next is divided into a series of intermediate steps. For each intermediate step, the computational procedure is relatively simple, an exact solution can be obtained in some cases and the time step can be quite large.. However, this scheme still has the disadvantage 2 that its co sistency 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 models were developed based on the assumption of a fixed boundary with a vertical wall located at the shoreline defined by the mean water depth. However, the shoreline boundary can actually move with time in such problems as wave runup on a beach and coastal flooding into a dry coastal region due to tides or storm surges. Some researchers have tried to simulate this problem numerically. Reid and Bodine (1968) considered the motion of the shoreline according to the water elevation and used the empirical weir formulation to compute the velocity in which water flows into or out of the dry land region. The disadvantage of this method is that the empirical coefficients in the weir formulation are very site-dependent. Yeh and Chou (1979) treated the shoreline boundary as a discrete moving boundary, but the impulsive motion of the boundary could introduce serious numerical problems. Lynch and Gray (1980) simulated the shoreline boundary as a continuous moving boundary using continuous grid deformation. This method, however, can not be applied to problems with complex topography. In this study, a two-dimensional finite-difference numerical hydrodynamic model is developed from the shallow water equations by using the method of fractional steps. Two ways to simulate the moving boundary problem are studied. For some simplified flow situations, analytical solutions are compared with numerical solutions 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 numerical model. Theoretical solutions of tidal responses inside a rectangular basin with Coriolis effect and bottom friction obtained by Rahman (1981) will be used to test the ability of the model to simulate two-dimensional problems. In Chapter 4, two ways to simulate the moving boundary problem are discussed in detail. One way is similar to that used by Reid and Bodine (1968). Another way is to include the dry land region into the computational domain by assuming a thin water layer on the dry land region at all times. The theoretical solutions for wave propagation on a linearly sloping beach developed by Carrier and Greenspan (1958) are then presented and compared to the numerical solutions to verify the second method for the moving boundary problem. In Chapter 5, the wind driven circulation in Lake Okeechobee, Florida, including the effects of the Coriolis force and a moving boundary, is investigated numerically. In Chapter 6, conclusions are drawn and suggestions are made towards further studies on numerical simulation of the moving boundary problem. CHAPTER 2 DEVELOPMENT OF A TWO-DIMENSIONAL TIDAL CIRCULATION MODEL In this Chapter, a two-dimensional implicit tidal circulation model is developed using the method of fractional steps. The vertically-integrated momentum and continuity equations, which govern the two-dimensional tide-generated currents, 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 calculate 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 obtained by vertically integrating the three-dimensional Navier-Stokes equations of fluid motion. Generally, it is assumed that the density of water over the depth is constant and the vertical pressure variation is hydrostatic, thus leading to the following continuity and momentum equations in a right-handed Cartisian coordinate system shown in the Fig. 2.1 az au av -t- + + 5Y = 0 (2.1) au auu avu az bV 'r., a au a au - +.-- +-5y+ghax-flV + P lax(K-ax + (K-- y)] =0 (2.2) av auv avv az Trb TrOY a 8V a av -y g+ ---+ Uh- +v[1UK+ -+-(- )]=0 (2.3) P [1 (K x + y(K--) 23 x a-aT a x y a Z " FREE SURFACE Z (X, Y. T V u.U ZB (X.T) X Figure 2.1: Definition sketch for tidal equations where t is time, x and y are the spatial coordinates, Z(x,y,t) is the free surface elevation, U(x,y,t) and V(x,y,t) are the unit-width discharges in the x- and ydirections, respectively, u(x, y, t) and v(z, y, t) are the vertically-averaged velocities corrresponding to U(z,y,t) and V(x,y,t), h(z,y,t) is the water depth, f is the Coriolis acceleration parameter, g is the gravitational acceleration, p is the density of water, K is the horizontal turbulent diffusion coefficient, r,, and r,, are the wind shear stresses in the x- and y-directions, respectively, and rb, and rb, are the bottom shear stresses in the x- and y-directions. rb and rb, can be expressed as pgU'Uz + V2 rb = pgUI2 (2.4) C2h2 and pgVVU2 + V2 (2.5) Tby = c2h(2.5) 'C2h2 where C is the Chezy coefficient which is R1/6 C = 8.21- cm 12/sec. (2.6) n or R1/6 C = 1.49- f tl2/slec. (2.7) n where R is the hydraulic radius and n is the Manning coefficient. Obviously we have U=uh V = v h (2.8) h=Z-Z where Z 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 avU au au az -+ = hu- + v-) u-- (2.9) ax ay = z By +t and auV avV v av az + a = h(ua + v-) v-5- (2.10) Substituting them into Eqs. (2.2) and (2.3), we obtain au Bu au Uaz z rb r0, 8(K) 8(K--) -+ h(u + v ) u- + g v+ [ + ] = 0 t Tz y at ax p az ay (2.11) aV av av az + Z flu r+,v r (K)) ) = 0 + h(u- +v-)-v- +g + Y -[ ax ay = 8t 8z By( + t -B, +g y +r+ p 8z By u,-T, (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: un+1/3 u" au au +ut + +~-=0 (2.13) v"la- v 8v 8v S + u~- + va- = 0 (2.14) At T8z Bay U -+1/ = un+1/3h" (2.15) V"+I/ = v"+I/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/3 Un+1/3 8 au a aU - n V [y-(K-) + -(K--)] = 0 (2.17) At z8a ay y V"+2/+ Vu+1/s a 8 V a av At + 0 U [- (K -) + -y (K y )] = 0 (2.18) At 2: 8 ay' aBy 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 At + + = 0 (2.19) At az y Un+1 U"**/3 aZ z Z az r,, At -7 + gh- + = 0 (2.20) At at 8z p Vra+1 V"**2/3 az az 9 r,(221 At at ay P The terms u- and vM come from the nonlinear terms. Compared with other terms in Eqs. (2.20) and (2.21), their values are small and could be neglected in the propagation step. 2.2 Numerical Scheme Just like the splitting of a single time step into three steps as shown above, three steps could be further split into two directional steps. For example, the advection step can be treated with the following two steps: x-advection Un+1/6 Un un+1/6 S + u" = 0 (2.22) At axz Un+1/6 Vn -n+1/6 t + Un = 0 (2.23) At Ox y-advection Un+1/3 Un+1+1/OUn+1/3 At + Un+V/6 y = 0 (2.24) At By Vn+1/3 n+1/6 Bvn+1/3 At + vn+1/6 a = 0 (2.25) At By where n + 1/6 represents the state after the x-advection. This scheme is implicit and unconditionally stable. For the diffusion step, an alternating direction implicit method is used which leads to x-sweep U* Un+/s a a BU* a BU(n+13 S- (K ay ) = V+/ (2.26) V* V+1/s 9 8V* a aV"**/3 (K ) (K) = 0 (2.27) 'At y y y-sweep Un+2/3 U* u U a Un+2/3 1 --(K -)- -(K l ) = 0 (2.28) At az az By By V _/- V* a av* a aV+2/3 'At / V (K ) ) (K ) = -U (2.29) At B 8x B y By where U* and V* are the intermediate values of unit-width discharge during the computation of the diffusion step. Yanenko (1971) showed that the ADI scheme was absolutely consistent and unconditionally stable for the two-dimensibnal heat conduction equations, which become the working equations for the diffusion step if the source terms are added. The propagation step is the most important of the three steps and its working 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' Un+2/3 At 8Z nB Zn U"**/3 Zn+1 zn + a(gh )+" + (1- a)(gh )- un ,3 ( aaz h At + a(rs. r,")"+, + (1 a) ( ")" = 0 (2.30) P P + a(gh-)" + ay Z(gh az V"+2/3 ( zn+ Zn (1- a)g ) ( At (2.31) + a( Tby T-. )-+ +(1 )(-b V)n = 0 p p aU + (1- a)( ) ay + Vn~ aV + (1 a).( )" = 0 (2.32) + 5y 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 Un+I and Vn+l as follows aZ aZ Un - At{a(gh -)"+ + (1- a)(gh.)" - At{( r r8 )n+1 + (1 a)(b- _7~esz)"} P p +2/3 zn+l zn ,"( zt )} n At (2.33) - At{a(gh- z)n+1 + B y (1 a)(gh-z)" B y Vn+2/3 Zn+1 g Zn hn At - At{a(Trsv r,).+l +(1_ a)( r,)"} p p (2.34) Treating the bottom friction and surface friction explicitly for clarity, and substituting Un+1 and Vn+1 into equation (2.32) yield 2a (na~z +Zazn+ 8 8 AZ BZ." a2{ X(h" + AZ ) + a 8 U"+*Is 8 V + { a ( U--2/3Az) + a( gA t Tz hY a BU -+,7/s -- ( a ) gd-t 8z X 8 ,AZ a(h n ay ay az" + AZ -)} ay n+2/3 n AZ)} = f + f2 (2.35) 8 8Z+)" aa 8rb-r, + a a(hn az a ( -) (2.36) Tz z g-) 8 x p Vn+l Vn+2/3 At Z'+1 Z" aU ,+ At a(x) Un+1 = Un+2/ Vn+1 = V +2/ AZ g(At)2 where _ (1 a) (aU t u ) gAt 8% 10 (1 a) aV n a Z" a a rs, ( r, (1 ) a( )" + a-(+ )Z + ( )" (2.37) 2 gAt ay -At ay) y (hBy ay(.37 gy and AZ = Za+1 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'Zx 8= O AZI BZ" a 8 U"+*/s AZ1 2 a 8 ,AZ1 az" a a vU+2Is - a -(ha + AZ -) + a( AZI) = f1 q (2.39) 2g(At)2 ax 8: ax gAta8: h" TAZ 2 a 2 a a n2/ A a (ha + AZ- ) + ( AZ2) = f2 + q (2.40) 2g(At)2 ay ay ay gAt y h in which q(z, y) is the coordination term. If AZ = 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 AZ, AZ2 and q are presented in appendix A. If the bottom friction terms in Eqs. ( 2.30) and ( 2.31) are treated implicitly, the first-order series expansion in terms of velocity (U, V) and water depth h can be used to linearize them as ()n+1 = (rb)n+2/3 +a(bz)n+2/3Ah + a ()rb)n+2/3AU (2.41) p p ah p aUp (y)n+1 = ()n+2/3 +a (TbV)n+2/3Ah + -v) n+2/3AV (2.42) p p BA p aV p where Ah = h"+ -h"= AZ) AU = Un+1 Un+2/3 (2.43) AV = Vn+1 Vn+2/3 Substituting Eqs. (2.4) and (2.5) into the above equations, we get n+1 = n g ),+2/s g 2UvU + V2) n+2/3z g VU2 + V2 U )"+2/3(U"+l Un+2/3) (2.44) C2( h2 + h U2 +V2 n+1 .. n+2/3 _g 2vU V2n2VU + 23AZ (-)V- = C3 P ~ +V2 V + V-- + V2 2)"+/3(V+'- V+2/3) (2.45) C2 h2 h2. /U + V2)( Plugging them into Eqs. (2.30) and ( 2.31), expressions for U"+*' and V'+1 can be extracted as U"+1 = M{U"n+13 + h-----AZ aAt(gh- )n+ (1 a)At(gh )n g UV TV1g r2U2+ V 2 +lr+l - At-( UVU2 + V)+2/3} + aAt-( 2U + V n+U C2 h C2 h2V/U2 + V2 ..g U U U2 + V)72 lA + 2aAt- ( UV) V (2.46) C2 h3 (2.46) Vn+1 = N{V"+2/3 + Vn2/3 AZ aAt(gh +1 (1 a)At(gh z)n AZ- ayg-- -(1 C2 h2 C2 htdU2 +V2 g VVUs + V2 +e2at( )n+2/3 C'h3 where M = 1/1 + ag 2U V 2+ 2 ] (2.48) C2 ,h2V/U2 + V" and N 1/ a + gAt U2 + 2V2 )n+2/s] (2.49) N = 1/[1 + ( )"*U-- (2.49 C2 h2VU2 +V2 Substituting Un+1 and Vn+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 AZa 8(hnaAZ) a(AZz") Ma a(U+2/ AZ ) SMa{ aZ + -- f q (2.50) 2g(At)2 ax ax gAt 8x 12 2g(At)2 Na h { + az N+ =--- ) f2 + q (2.51) 2g(At)2 ay ay gAt ay where f= 1-a au n Ma au + 8(hM)" Ma(UU2+v2)n+2/3 f g 5 ) "n+- ( )" + Ma a= + C2 (2.52) gAt 8X gAt 8z azx28 1-c a V Ma V +2 (h )" Maa(V )n*+2/3 12 = (- ( ) -(y) + Me O + C h (2.53) gAt y) gAt 9Y +y C2 y 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 Ayl 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) AX(i) + AX(i+l) (2.54) 2 A Xu(t) = 2(2.54) AY,(,) = AY 2(u_.) + AY,(.) (2.55) 2 where AX,(j) and AY,(,) present the sizes of the cell (i,j) in the (U, V) grid system, AX,(i) and AY(i) 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,(i,) and at the Z grid points by Pz(ij). P,(i+i/2,y) is used to denote the value at the middle point between the Z grid points (i,j) and (i + 1,j), and Pz(ij.i+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,) = [Pz(id) + Pz(i,.-1) + P(i+',i) + Ps(i+l,y-1))/4 (2.56) xI Figure 2.2: Schematic of finite difference mesh with variable rectangular cells J*2 J-I I-1 I 1Xu I I I I - A----- Z-GRID -------- ---. (U.V) -GRID --+ ----.--_____ (it 1 L, A. 1-I I II l*2 1Xz c1 -J*2 "T AU U)J .jI Figure 2.3: Computational grid definition P(iy) 4AXU(i1)AY,(j) PU(iy+) + 4AX()AYz() P(i,) + AX,()AY,(+-1) AX.(i)AY(y) 4AX,(-)AY,(y),.,(i-1,i+1) + 4AXe(4-I)AYu)Pu(4-Ly) (2.57) AY.(y-1)AY (y) Pz(i+1/2,y) = AYz(i)p(i+,)+ 2A )Pu(ij) (2.58) and Pz(i,j+1/2) AXz(i-l1) P(i,+l) + AX.(i) P(i-1,jy+i) (2.59) 2AXu(i-,) 2AXu(i-I) Pa(i+l/2,j) and Pz(i,j+l/2) can also be evaluated by Pz(i+1/2,.i) = [Pz(,,ij) + Pz(i+li)]/2 (2.60) P.(i,j+1/2) = [P(,ij) + P.(,j+1)]/2 (2.61) For a spatially uniform grid system, Eqs. (2.57), (2.58) and (2.59) become PF(id) = [Pu(ii+l) + P(ii) + P u(i-.ij+1) + Pu(i-ij)]/4 (2.62) Pz(i+i/2 i) = [Pu(i,+1)+ P.(id)]/2 (2.63) PF(ii+1/2) = [Pu(i,j+I) + Ptu(i-,ji+,)]/2 (2.64) 2.4 Finite Difference Approximation Suppose that the state of the model is known at time nAt. Then the state of the model at time (n + 1)At can be obtained by successively solving the .advection, diffusion and propagation steps. From Eqs. (2.22) to (2.29), it is seen that the water surface elevation Z does not appear in the advection and diffusion steps. Therefore these steps can be executed solely on the (U, V) grid. The finite difference equations for the advection step are x-advection n+1/6 n+1/6 n+1/6 '7~ s, + Un i l -U -~ S U + U i+,j 0 (2.65) At "'AXe(i-1) + AX,(i) n+1/6 n n+1/6 n+1/6 vij ,i + u + Vi+l,i vi-1i 0 At AXU(-1) + AX o) n+1/3 -n+1/6 n+1/3 n+1/3 Uy ui, + n+1/6 ui,y+ iy-1 A t 0y(~ ) + A .y n+1/3 n+1/6 Vi. Vi, At n+1/3 n+1/3 + n+1/6 i,+ -- i,i-1 = 0 +II AY(y-1) + AY.(y) where the subscripts i and j correspond to the (U, V) grid. After x-advection and y-advection, we get the modified velocities un+1/s and vn+1/3. The discharges corresponding to u"+1/3 and V"+l/s can be calculated by (2.69) (2.70) Un+1/3 n+1/3 n." id = ui I. +1/3 +1/3 Vi +I/s 10I/3h in which h". is water depth on the (U, V) grid at time nAt. The finite difference equations for the diffusion step are as follows: x-sweep U Un+1/3 1At U~j III *. n+1/3 At K U! -U!. U: K U *+ 1. U-*1, U U*- AX,(W AXo AXO(i-1) K U,"3+/ UV+13 U.+"ls U+1/s K [ i+1 i i, ij-1 Ay,(I) Ay() X AY.) K Vn+1/s i +13g~i V., IT K V V* V _i 3Xzw AX,,(i) AX,,(i1)] V-+l/S n.+1/s F-+1/3 V.-+1/3 ,Ij+1 ,; s,y d-1 AY(y-1) AYU(y) AYu(-1) = 0 (2.71) .(2.72) y-sweep AX-(i) AxK,) U:..j U:.-j14 i*-I) y-advection (2.66) (2.67) (2.68) U+2/3 , 2At n+2/3 i t Vii 'At 16 K U+2/ U '+2/3 U"+2/3 U.+2/3 KA [1+1 s 1f ,j-1 Y- AYu() AYC,.y_1) = 0 K V. V. Vj,* V-1, AX.([ AX.,W AX,(i-1) K2-_ Vn+2/ V-+2/3 V-+2/s K .i+ i i" i"-1 AY.=-1) AYU(*) Ayg_ ) = Uij where (i, j) also corresponds to the (U,V) grid. After the diffusion step, the discharges U"+2/3 and V"n+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 Zn+x can be expressed as a2 Ail(1+l_ ) S h AZi(i++) AZ(ii) A+1 AX,(i-~) hi+1/2j AX- + AZ(i+1/2j) AX h" AZI(i,.i) AZx(i-I,y) y. ."1, -/2 A X (i ) --A Z x(i- / i A X ,(p:" ] U +2/3 rn+2/3 gdtX,(-)[h~,t AZ(i+1/2.i) AZ1(i-1/2,i)] ' ~x(-1 i+1/2,y h Fll2,y = fi-q (2. a2 n AZ2(ij+l) AZ2(ij) Z"+1 Ay) [hi,+1/2 AY(i) + AZ2(i,+1/2) AYYUCY) h l AZ2(ij,) AZ2(i,-1) y AZ2(ij-1/2) ,-1 J-1/2 AYz(-1) (id-1/2) AYz(-1) t n+2/3 Vn+2/3 + cc Yi) AZ2(i,f+1/2) -2(i-1/2) =() j+1/2 ij-1/2 = 2+ q (2. 75) 76) in which . (1n+2/3 U)+2/3 (1 a) Us+1/2,y Ui-1/2,y a +1 /2,y i-1/ , gAt AXu(i-1) gAt AXu(i-1) (2.73) (2.74) AZI(ii) 2g(At)* A Z2(,i) 2g(At)2 "l Z Zin Z +n + Zi i id 1 J XX( hi+1/2,j aXi) a + pgAX(-1) [(r: rS)+1/2, (rb r,,)-1/2,( 1n + 2 / 3 v n + 2 /3 f (1 a) Vi'3+1/2 V.-1/2 a J+1/ yi.Y-1/2 gat AYU(j) gAt AYu(j) a Z" Z" Z. Z" " Y [h ij+1 s n i.jI-1I Y(A ,3+1/2 Y() 1,i-1/2 ( a y ay [(QrV r,), +1/2 (r5 1,)-1/2] p ~ j J Id,=y hi/2; = (hil,, + hi,y)/2 hi/2 = (hjdt + hi,2)/.2 AZit/2,i = (AZit,i + AZ~j,)/2 AZi.1/2 = (AZi,j + AZti)/2 Note that the subscripts i, j, i 1/2 and j 1/2 refer to the Z-grid and the discharges at the Z-grid can be calculated by solving Eqs. (2.57), (2.58) and (2.59) in terms of the discharges on the (U, V) grid. The finite difference forms for U"+1 and Vn+l are U n+l '3 Zn+1 7n+1 = Ui, At!gh.i1 Z+(1/2 1/2, '. 2.3 A . Z" Zn - At(1 a)ghn y +1/2,; i-1/2,j i AX.() Un+2/3 + (Z'+ hl (,y n )3 - Ata( -T )'- A t(1 a) (b ( o -r,,, I (r. p),. P $13 p i fn+l 7n+1 SV"+2/3 Ataght' "i+1Z/2 i- 1/2 A, n ( 1/ ~) n n Vn+2/3 - At(1 a)gh-, in,++ 1 + 2/ 1 z ) W~ Ay,(yI) j', (Y - Ata At(1 a) (rs r )"f (rbs r,,Y p p in which the indices i, j, i 1/2 and 5 + 1/2 refer to the (U, V) grid. (2.77) and (2.78) (2.79) '.+1 (2.80) (2.81) 18 2.5 Boundary and Initial Conditions Two types of boundary conditions are normally encountered in the numerical simulation of tidal circulation. One is the open boundary which is an artificial termination of the computational system. The other is the shoreline boundary. At the open boundary, the water surface elevation and/or the current velocities are specified. In this model, the velocity gradients in the direction normal to the open boundary is set to zero, and the water surface elevation is specified along the open boundary. The shoreline boundary can be treated as either a moving boundary or a fixed boundary. The moving boundary problem will be discussed in detail in Chapter 4. Along a fixed boundary, the normal velocity is zero, but the tangential velocity may either be set to zero ( no-slio condition ) or satisfy the zero-slope in the normal direction ( slip condition). In this model, boundary conditions must be imposed- at each of the fractional steps., Since the water surface elevation does not appear in the advection and diffusion steps, only velocity boundary conditions are imposed for these steps. For the propagation step, however, boundary conditions in terms of surface elevation and velocities must be specified along the open and closed ( fixed or moving boundaries. The specification of the initial conditions requires the knowledge of the free surface position and the flow field at t = 0. Usually this is impossible and the model is started with the still water condition that Z =constant and U = V = 0 at t = 0. 2.6 Consistency and stability As seen previously, the numerical scheme for each of the three fractional steps is relatively simple. However, we need to show that the combination of these three steps is consistent with the original governing equations. This is hard to prove theoretically because of the existence of. the two-dimensional nonlinear terms and diffusion terms. This is a drawback in applying the method of fractional steps to the problems of fluid dynamics in several variables. It will be seen in the next chapter, however, that this drawback does not appear to affect the accuracy of the numerical model. Since each fractional step is unconditionally stable, the complete scheme is unconditional stable. However, for numerical accuracy, it'is desirable to keep the time step sufficiently small such that the Courant number based on the velocities does not become exceedingly larger than 1. CHAPTER 3 COMPARISONS BETWEEN THEORETICAL AND NUMERICAL SOLUTIONS In this Chapter, results of the two-dimensional tidal circulation model will be compared to four different theoretical solutions. The first comparison will primarily investigate the model's capability to simulate the propagation step. The second comparison will validate the model's ability to compute the nonlinear effects. The other two comparisons will focus on the simulation of the Coriolis and friction forces. From these comparisons, we can then assess the model's accuracy and consistency. 3.1 Comparison With a Linear Theoretical Solution Neglecting the advective terms, diffusion terms, bottom friction and wind surface stress, the one-dimensional shallow water equations are au az + gh- = 0 (3.1) az au (3.2) where U = uh represents the vertically integrated velocity in the x-direction, u is the vertically averaged velocity, h is the mean water depth and Z is water surface elevation. Assuming a closed boundary at z = I and an open boundary at x = 0, the boundary conditions and initial conditions associated with Eqs. (3.1) and (3.2) are: Boundary conditions U(, t) =0 (3.3) Z(O,t) = Zo +asinwt (3.4) Initial conditions Z(z,O) Zo (3.5) U(XO) =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 a2z= C (3.7) where c is the wave speed. which can be expressed as c = Vg-. Let the solution of Eq. (3.7) take the form of Z~~t =acosk(I ) +0 (coskl sinwt + E A,,sink,,xsinwnt + Zo (3.8) n0O 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 relationship w2 = k2 c2 (3.9) n kn 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 az 0 (3.11) In order to satisfy this boundary condition, it follows E A kcos kl sin wt = 0 (3.12) n=O Therefore cos knl = 0 (3.13) So k = (2n+ 1)r (3.14) where n is an integer. By using initial condition (3.6), we can determine the coefficient An as -2aw f I cos k( x) sink lw, Jo cos ki - 2aw ,(k )(3.15) lc(k2 k2) Consequently, we obtain acosk(Q-x) s0[ -2aw Zcoskl s + c(k k2) sin +2x sinw t] + Z0 (3.16) Substituting the above equation into Eq. (3.1), we get U~x,t) = ac sink(l x) 00 -2aw Cos cost ==o(kn k2) cos kX cos wnt] (3.17) To compare the numerical solution with a linear theoretical solution, a rectangular 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 i CLOSED BOUNORRY KM ///////////////////////////// V=0Or U V=O I 0 CLOSED BOUNDRRT 60KM Figure 3.1: Schematic diagram of a rectangular basin &X= &T= 2KM Figure 3.2: Computational grid -bL - THEORETICAL 50LUTION A NUMERICAL SOLUTION 100 NJ z 00 > 0o LU = -50 -100L 60 -100 NJ z s ID. > 0 =-SO LU -100L60 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 surface elevation 66 72 78 84 90 96 TIME (HR) - I.UU 50 a 0 -50 60 LU -j LU -100 60 THEORETICRL SOLUTION NUMERICAL SOLUTION 100 X=IOKM. Y=15KM w50 *** so & .& L -50 -LJ -100 I I I 60 66 72 78 84 90 96 TIME (HR) 100 - X=3OKM. T=15KM o 50 N. 0 S-50 -LJ -100 I I I I I I 60 66 72 78 84 90 96 TIME (HR) 100 X=5OKM. T=15KM W so U -50 Uj 00 -lo - 1 0I I I I I I 60 66 72 78 84 90 96 TIME (HR) Figure 3.4: Comparison between theoretical and numerical solutions for velocity - THEORETICAL SOLUTION ------- NUMERICAL SOLUTION WITH TIME STEP: 30MIN. ----------- NUMERICAL SOLUTION WITH TIME STEP: 15MIN. 66 72 78 84 90 TIME (HR) 78 TIME (HR) Figure 3.5: Comparison time steps of theoretical solution to numerical results with different 2 o so z 0 - 50 ..J u-1 Ir -SO -100 L 60 S50 ,l 0 m -5 -S X=3OKM. T=I5KM I _____ ~ ~ I -100 60 I I f III I I I I. 27 (x = 30km), and near the closed boundary (x = 50km). Both Fig. 3.3 and Fig. 3.4 show that there is a reasonably good agreement between theoretical and numerical solutions, even though a slight phase shift between the theoretical and numerical results exists. To investigate the origin of the phase shift, different time steps for the numerical computation are used and the numerical results for water surface elevation and velocity are presented in Fig. 3.5. From these results, it is seen that the phase shift between theoretical and numerical solutions tends to diminish as the time step decreases. 3.2 Comparison With a Non-linear Theoretical Solution When nonlinear terms are included in two-dimensional shallow water equations, it is impossible to obtain an analytical solution. In one-dimensional nonlinear tidal motion, however, we can use harmonic analysis to develop a theoretical solution. Obtaining exact theoretical solutions, however, is still difficult because the high order terms are difficult to solve for. In this study, only the zeroth and first order harmonic solutions are developed. Without the consideration of Coriolis and bottom friction forces, the governing equations for one-dimensional nonlinear tidal motion can be written as au au az az au where h is assumed to be constant. The bounary conditions are Z(0, t) = a sin wt (3.20) U(l'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 + U"2 u I + U2 (3.24) h where Z1 and U1 are the theoretical solutions of one-dimensional linear shallow water equations which have been obtained in the preceding section: =a cos k(1- z) Z acosk( x) sin wt (3.25) cos ki U ac sin k(I z) cosk 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 Z2 a2c2k + gh O 8hcos2 k {sin[2k(1 x) + 2wt] at 8x 8& cos2 ki + sin[2k(l x) 2wt] + 2sin2k(l- x)} (3.27) az2+ 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 (l, 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[xsin2k(I- x) + sin2k(l + x) 8h cos2ki cos 4ki 1 - tan 2kl cos 2k(I x)] cos 2wt (3.31) cos4kl ab_____ 1 U2 = 2 [ xcos 2k(l X) + sin 2k(1 x) 8h cos2 ki 2k l 1 o 4Il cos 2k(1 + x) + o I tan 2ki sin 2k(1 x)] sin 2wt (3.32) c: os4ki c os 4ki - 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 (MR) UU X=60KM. T=15KM 0 f0 I I I I I I I 78 TIME (HR) Figure 3.6: Comparison between theoretical and numerical solutions for water surface elevation with nonlinear effects 66 72 78 84 90 TIME (HAR) 0 R_ luu- LUU. z 50 0 0- > 0 LU -J LU - -50 LU z -100 60 N .I o z 1j u-I LU P-1 60 THEORETICAL SOLUTION NUMERICRL SOLUTION 100 50 0 , U -so -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) 0DI Li .U50 0 a-50 os ,-I U.' -100 60 U 50 1 0 .-50 I-, -100 60 X=30KM, T=15KM I I I I 1 I I I 31 Details on the derivation of Z2 and U2 is presented in appendix B. The same basin and forcing conditions of the preceding section will now be used to compare theoretical solutions with numerical results obtained with a time step of 30 minutes. The comparisons between the numerical results and the theoretical solutions are shown in Figs. 3.6 and 3.7. Figure 3.6 presents water surface elevations near the mouth, at the middle point, and near the closed boundary of the rectangular basin. Figure 3.7 shows velocities near the mouth, at the middle point, and near the closed boundary. From Figs. 3.6 and 3.7, it can be seen that the numerical results are reasonably similar to the theoretical solutions. 3.3 Comparison to a Theoretical Solution with Coriolis Effect Neglecting convection, diffusion, and bottom friction, the two-dimensional linear shallow water equations can be written as au 2az- nv = 0 (3.33) at+ c -av 2az t+c-y flU =0 (3.34) az au av (.5 S+ -+ = (3.35) where fl is the Coriolis parameter, and c = V'- 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(O,y,t) Z. (3.36) U(I,y,t) =0 (3.37) V(,otO 0 (3.38) V(x,b, = o (3.39) where Z, is the forced water surface elevation at the mouth of the basin, 1 and b are the length and width of the basin, respectively. 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[-y +ik(1- z) iwt] Glk + AR exp[l (b y) ik(1 z) iwt] + -{C.[cosKin(b-y)+ K'sinKi,(b-y)] n=1 wKl exp[-iK2,(l z) iwt]} (3.40) _ Akc2 [k U = exp[ -y + ik(I z) iwt] ARkc2 + ezp[-(b y) ik(I z) iwt] w w C 2 o l~w snK.b ) + E {C.[K cosKin.(b- y) + KC2sinKn(by)] n=1 in exp[-iK2n(I x) iwt]} (3.41) ic2w 00 l2K2 V = E- {C.Kn(1 + w 2, sin KIn(b y) W2 n2 n= 2Kl n=1 n exp[-iK2n(1 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, Kl, and K2, are wave numbers of the nth Poincare mode with respect to the y- and x-directions, respectively, and Cn is the amplitude of the nth Poincare mode. The dispersion relationship for Kelvin waves is w = k2c2 (3.43) For Poincare waves, we have K2 2 K1 + K2. = (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(- ky) + ARC exp[-k(b -y) W W + C,[ K2cosK,,(b y) + K-- sin K,,(b y)] =0 (3.46) n=1 K1c The simple way to solve for the unknown coefficients R and Cn 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 fi 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 Z = Re{Aexp(-y + ikl) + AR exp[-(b y) iki] 1flK2, + C,,[cosKi,(b y) + "K22 sin Kl,,(b y)]exp(-iK2,l)} (3.47) n=1wK 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 NUMERICAL SOLUTION 100 50 0 -50 -100 60 UU .X=30KM. r=15KM 0 0 00 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 surface elevation with coriolis effect 66 72 78 84 90 TIME (HR) - I S1 z -5 LU . 1U cc -5c LU. -I. 60 IUU z 50 o > 0 LU -j LU S-50 LU -.r -100 60 THEORETICRL SOLUTION NUMERICAL SOLUTION X=10KM. T=15KM I I I 78 TIME (HR) 66 72 78 84 90 TIME (HR) 78 TIME (HMR) Figure 3.9: Comparison between theoretical and numerical solutions for velocity U with coriolis effect U.J #- -50 -J -100 6 0 so 0 u, 50 -50 C-) U.' O -100 60o 50o 0 -50 -100 60 X=50KM. T=15KM I l l l I I I THEORETICAL SOLUTION NUMERICAL SOLUTION 66 72 78 84 90 TIME (MHR) 25 U_) ). 0 Iul -25 -50 -so E 50 U 25 UI L) U > 0 0 -25 -50 Figure 3.10: Comparison between theoretical and numerical solutions for velocity V with coriolis effect THEORETICAL SOLUTION NUMERICAL SOLUTION X=6OKM. T=I5KM I I I I I I II I I I 66 72 78 84 90 91 TIME (HR) 37 the flow velocities u in the x-direction at three points in the middle axis of the basin. Figure 3.10 shows the flow velocities v in the y-direction at the middle point and at the closed boundary. Both the theoretical and the numerical solutions at x = 10km are very small, thus are not shown in the figures. These figures illustrate that the model results agree quite well with the theoretical solutions. 3.4 Comparison to a Theoretical Solution with Friction Effect With the. bottom frictions, the two-dimensional shallow water equations can be written as BU BZ n,= + gh a + = 0 (3.48) at 8x p 8 V B Z 4~ + gh-= + = 0 (3.49) az au av Bt y p BZ BU BV - + 5 + 7= 0 (3.50) 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 r = pFU (3.51) rbv = pFV (3.52) where g /2. + v2 = C2h2 (3.53) Given these, Eqs. (3.48), (3.49) and (3.50) can be simplified as BU B2Z a + C +FU = 0 (3.54) BV ax av 2az - + c- + FV = 0 (3.55) at ay az aU av 5- + + = 0 (3.56) 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) o00 + E [A,, cos Kin,,(b y) n= 1 (cos K2nx + tan K2,l sin K2,z)]exp(-iwt)} (3.57) -C2 U = Re{_-. Aok(-sinkx+ tanklcoskz)exp(-iwt) -C 2 oo + F- E [AK2, cos Kin(b y) rUn=1 (- sin K2,x + tan K2nl cos K2,z)]exp(-iwt)} (3.58) C2 o V = Re{F. F, [AKix sin Kin(b y) F n=1 (cos K2,x + tan K2n1 sin K2nz)]exp(-iwt)} (3.59) where b and I are the width and length of the basin, respectively, and w~r KIn = b (3.60) =w' F K = (1 + i-) (3.61) ;n = 2-(1+i) () (3.62) The coefficients A0,..., An can be obtained using the condition Z(0,y,t) = Z,, exp(-iwt) (3.63) at x = 0. They are given as Ao = Z,,dy (3.64) 39 THEORETICAL SOLUTION NUMERICAL SOLUTION - 100 X=10KM, T=15KM 2 50 so -J ,-50 LU I- I Cz -100 I I I I 60 66 72 78 84 90 96 TIME (HR) 100 x=30KM. T=15KM > 0 - 100 I I I I I I I I I I 60 66 72 78 84 90 96 TIME (HR) 100 X=60KM. T=15KM " r z 50 -I o AU c -50 -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 THEORETICRL SOLUTION NUMERICAL SOLUTION 00 X=10KM. T=15KM 0 0 0n I I I I I 66 72 78 BLI8 90 96 TIME (HAR) 66 72 78 8L 90 TIME (HAR) 78 TIME (HR) Figure 3.12: Comparison between theoretical and numerical solutions for velocity U with friction effect 5 LIJ 0 -5 ) -5 u. -ii 60 50 o -so 0 -J ..J -100 60 in. 50 0 ) LAJ S-50 o 0 .-I u-I X=50KM. Y=15KM I I I I -100 60 A jb Zm cos Ki b- y) dy (3.65) To compare the theory with the numerical results, the same basin defined previously is used. The amplitude of the forcing tide along the mouth of the basin is assumed to be constant (50 cm), which leads to no flow in the y-direction, i.e. V = 0, and zero value for the coefficients A,, ..., 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 velocities. 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 properly 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 boundaryThe 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 amplitude, 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 expressed as '?(t = o + o 6 dt(4.1) where X(t) is the location of the boundary at the time t, Yo~ is the initial location, and 6b is the velocity of the boundary motion. The basic boundary conditions on the moving boundary are h=O0 (4.2) 6 = 6b (4.3) where h is total water depth and tT 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. Usually the topography is given, but the water surface gradient on the moving boundary is unknown and is related to the velocity of the boundary motion implicitly. The third problem arises from the computation of the bottom friction in the vicinity of a moving boundary. The most common bottom friction formulation in vertically-integrated shallow water problems is the Manning-Chezy formulation. This formulation provides finite bottom friction throughout the interior of the computational domain. However, it breaks down near the moving boundary since the computed bottom friction approaches infinity as the water depth tends to zero. Using the conventional finite difference model, the computed strong bottom friction causes problems both in the flooding and drying. During the flooding, the str ng bottom friction leads to very sharp surface slope which may lead to very large velocity 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 deforming grids, it possesses other problems. The impulsive filling of a cell with fluid often leads to numerical problems unless treated very carefully. Reid and Bodine (1968) investigated the transient storm surges in Galveston Bay, Texas. Omitting nonlinear advection and Coriolis terms, they developed a finite difference numerical model with the inclusion of the tidal flat. A uniform Cartesian mesh and staggered grid system were used. The elevation of the sea bed or land was represented by a constant value at each grid point within a square grid. Hence, the actual topography was approximated in a stair-step fashion. The movement of the boundary was controlled by the water elevation. If the water elevation in a flooded square was less than the base elevation of an adjacent dry square, then a zero normal flow boundary condition was applied along their common boundary. However, if the water elevation in a flooded square was greater than that of an. adjacent dry square, then the water was permitted to flow into the dry square. The flow rate between two squares was determined using an empirical equation for flow over a broad-crested barrier. The overtopping of a barrier could be treated also. However, the model could treat only barriers aligned along the grid mesh division. Flow across the barrier was permitted when the water height on one side exceeded the barrier height. If the water height exceeded the barrier height on both sides, then the flow rate was determined using an empirical equation for flow over a submerged weir. The empirical coefficients in the model were determined by iteration, comparing the model with tidal data and data from hurricane Carla (September 9-12, 1963). The gross features of the inundation were predicted reasonably well. However, it should be noted that the empirical coefficients used could be very site-dependent. Yeh and Yeh (1976) developed a nonlinear nondispersive moving boundary model for simulating storm surge using an ADI technique. The shoreline in this 45 numerical model moved as the flow inundated low lying land. However, details of the treatment of the boundary were not given. It appears that the shoreline advanced or retreated in discrete increments of grid cells. They reported good agreement with field data. Yeh and Chou (1979) developed a nonlinear nondispersive finite difference surge model using an explicit technique with reference to a fixed grid system. The boundaxy 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 predicted 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 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 w-ith the experimental results. Some French researchers ( Benque et al., 1982 ) developed a finite difference numerical model with the inclusion of the tidal fiat using the method of fractional steps. They solved shallow water equations through three steps, which are advection, diffusion and propagation steps. A different numerical scheme was applied at each computational step. The treatment of the boundary motion was considered at the propagation step. The dry land region was assumed to be covered by a thin water layer. The flow in this region was governed by the bottom friction. During the computation of the propagation, the shallow water equations were first applied to the whole computational domain including the region of the thin water layer. Then the solution in the vicinity of the moving boundary needed to be adjusted to satisfy the resistance flow. The Manning-Chezy friction formulation was used to compute the bottom friction. They reported that the moving boundary model developed in this way violated the continuity equation slightly. Good agreement of the numerical results with the measured data were presented based on the model application to the Bay of Saint Brieuc and the River Canche Estuary, France. Lynch and Gray (1980) outlined a general technique whereby a moving boundary 47 can be treated by finite element Eulerian models. The finite element basis functions axe chosen to be functions of time so that the element boundaries track the moving shoreline. They showed how this motion generates extra terms which, if treated properly, reduce the problem to one that can be treated by standard finite element procedures. They showed how to apply the method to treat the propagation and runup of long waves. They looked at two simple problems involving the runup of waves on plane beaches. They showed that estimating the runup by extrapolating the wave height at a vertical wall could introduce significant errors. They also pointed out that treating deforming elements is computationally more expensive than fixed ones, and recognized that potential problem can arise if the mesh becomes geometrically too distorted. Gopalakrishnan and Tung (1983) described a finite-element nonlinear long wave runup model valid only for one horizontal dimension. The model contained terms that accounted for. vertical accelerations. The moving shoreline was handled by allowing the shoreline element to deform so that the beach node always tracked the shoreline. If the shoreline element became too stretched, it split into two elements. The element containing the shoreline node continued to deform but the other new element created by splitting stayed fixed. They showed some plots that detailed the runup process, but they did-not present any results about the rundown process. The technique outlined by the authors seems applicable to tsunami runup, but they did not present a thorough or convincing argument to show that their model could be used reliably for such studies. It should be noted that the technique in their work cannot be extended easily to include two horizontal dimensions since the element-splitting procedure would be very complex. 4.3 Numerical Treatment of a Moving Boundary In this section, two ways to simulate the moving boundary problem will be studied in detail. The first method was proposed by Reid and Bodine, and the 48 second method was proposed by Benque et al. (1982). 4.3.1 Method to Treat a Moving Boundary with the Weir Formulation Reid and Bodine (1968) treated the continuous moving boundary as a discrete moving boundary. In their scheme, a discrete Cartesian grid system is used and the actual topography in the vicinity of a moving boundary is approximated by twodimensional stair-steps. Thus the elevation of the sea bed or land can be regarded as uniform over each grid square. The boundary motion is controlled by the water level. During the flooding, a grid point is added to the computational system if its water depth is greater than a preset value h, which is the minimum water depth for the effective application of the Navier-Stokes equations. The value of h, 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 hl. In the computation of flooding, the condition of no normal mass flux is applied at the moving boundary when the water elevation is less than that of the adjacent dry land, i.e., the normal component of flow, Q,,, at the juncture of a flooded cell and a dry cell is taken as zero, while the tangential component of flow may satisfy the no-slip or free-slip condition. However, if the water elevation is greater than that of the dry land, flow will be allowed to flood into the dry cell until the water depth in this cell reaches hl. The mass flux per unit width of the dry cell, Q,,, can be calculated by the weir formulation ( Reid and Bodine, 1968) as: Q. =CoD g I Zd-ZC, (4.4) where Zd and Z, are the water surface elevation in the donator and acceptor cells, respectively; C0 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' A (4.6) 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 AZd=Q, Atd (4.7) where Ad is the area of the donator cell. During the flooding computation of each time step, each donator cell must be examined to see if too much water is flooded from the donator cell to the acceptor cell. If the water level in the donator cell is less than that in the acceptor cell, the mass flux computed by Eq. (4.4) must be adjusted to make the water level in the donator cell the same as that in the acceptor cell. When the water depth becomes small during the drying, an artificial water depth has to be considered so that the Manning-Chezy friction formulation can still be used to compute the bottom friction. The value of this artificial water depth is determined empirically. Usually we can set it to be 20 cm. The development of a two-dimensional moving boundary model in this manner is very cumbersome since the location of the new boundary must be determined at each time step. But if the grids in the vicinity of the moving boundary and the time step are kept small, we can obtain reasonable numerical solutions for the moving boundary. 4.3.2 Method to Treat a Moving Boundary With a Thin Water Layer Based on an assumption that the water flow in the region of shallow water is dominated by the bottom friction, Benque et al. (1982) proposed a way to 50 treat the inoving boundary by controlling the water flow discharges through the adjustment of the water depth. In this scheme, a grid system is first established in the computational domain which includes the entire tidal flat region. To ensure numerical accuracy, grid sizes in the tidal flat region should be smaller than those in the main computational domain. A thin water layer is assumed to exist at all times over the dry land region so that all grids in the computational domain are always wet. Thus, no grid needs to be taken out of the computational domain during the computation and we do not need to consider the motion of the boundary. The shoreline boundary can thus be treated simply as a fixed boundary, so the Navier-Stokes equations of fluid motion can be applied to the entire computational domain so long as the bottom friction is adequately resolved for small water depths. From the Manning-Chezy friction formulation, we see that the water depth plays an important role in the computation of the bottom friction, so it is possible to "control" the bottom friction by adjusting the water depth. As seen in Chapter 2, the free surface elevation Z does not appear during the advection and diffusion steps. Since the motion of the boundary is mainly controlled by the water elevation, it is only necessary to consider the propagation step. From the governing equations, Eqs. (2.75) and (2.76), of the propagation step, it is seen that the water depth, h, and the increment of water surface elevation during one time step, AZ, at the points of i 1/2 and j 1/2 must be calculated. They are governed by the water depth and the surface elevation at their upstream and downstream grid points, and can be calculated by the following formulae: hi+1/2 = Ii hi + (1 yi)hi+l y1 = Z1 h,_ + (1 -y -)Z i (4.8) Azi1/2 =i-1 Azi1 + (1 -Yi,.)AZi Ai12= 'y,-1 AZ.... + (1 Y-) i and h+1/2= yji hi + (1 y.)h,+ hy-112 = yj-l hy-1 + (1 qy_)hj (4.9). AZ+1/2 = -I AZ, + (1 Y,)AZ+1( A -12= -i AZj..1 + (1- Y-tl)AZj where -y 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 y can be set to 0.5. However, in the regions where the water flow is dominated by the bottom resistance, -y 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-g + = 0 (4.10) 8:p BZ n ghY +- = 0 (4.11) Using the Manning-Chezy friction formulation, pgUVU2V (4.12)+V2 Tb = C2 h2 (4.12) pgVv/U2 + V2 rb = C2 h2 (4.13) we can rewrite Eqs. (4.10) and (4.11) as BZ UVU2 +V2 -x- + C2h3 0 (4.14) aZ VV'U2 +V2 - + 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, BU -a, < 0 (4.16) aZda 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, U-VU2 +V2 +A = 0 (4.18) A + C2 hi+1/2 Zd, Z,+ V VU2 + V2 + = 0 (4.19) Ayi C2 h+112 where the subscript us refers to the upstream. In Eq. (4.18), the upstream and downstream correspond to flow in the x-direction, but in Eq. (4.19) they refer to flow in the y-direction. Introducing a coefficient f, the water depth h in Eqs. (4.18) and (4.19) can be expressed as hi+1/2 = 8 ha. + (1 0)ha, (4.20) hi+1/2 = 1 hu, + (1 9)hd, (4.21) where hu, = Zu, ZBu, (4.22) hd, = Zda, ZBd. (4.23) in which ZBu, 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 Zds and bearing in mind that U is a function of Z,, and Zd,, and V is constant, we obtain B U U2 C2h2 ah (U + V2 +) = [-h + 3(Z,, Zd,) ] (4.24) B a, (V' U= + V) =Azx---T, a~du Applying the condition a < 0, we get Oh - h + 3(Z,. Zd,) -Z, <0 (4.25) From Eq. (4.20), we know 8h - 1 (4.26) 8 34,d 53 Substituting 3h and h into inequality (4.25) yields -[Oh,, + (1 P)hd.] + 3(1 0)(Z Zd,) < 0 (4.27) So we obtain P _> 3(Zu. Zd,) h, (4.28) 3(Z. Zds) + h,. hda. Differentiating Eq. (4.19) with respect to Zd,, 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 replaced by either f or 1 ,3 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, -y = )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, U, and Vm,,, from Eqs. (4.14) and (4.15) U,..,= sign(u) C2h3 fa 1 (4.29) ax V a = sign(v) C2h3 a- I (4.30) in which sign(u) = az/ -z (4.31) 8z- ax sign(v) = -z / az- (4.32) 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 Umaz br Vm.., U or V must be replaced by Umz or Vmaz 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 Umz 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 49-- + ax [(17 + h')u'] = 0 (4.33) agu* (u" ag/ + U 7 + On- = 0 (4.34) where the asterisks denote dimensional quantities, 17 is the displacement of water surface above the mean water level, h is the still water depth which varies linearly with x, 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'Cg (4.35) U0 = 0 (4.36) where 0 is the beach angle. Let us choose the following nondimensionalization x* t"r* X- t* 77 X L. T eL Figure 4.1: Definition sketch for wave propagation on a sloping beach h x u = (4.37) L 210 and define C2 = h* + l = h + x + n (4.38) L With these definitions, Eqs. (4.33) and (4.34) become 17, + [(77 + x)uJ = 0 (4.39) ut + u u. + 77z = 0 (4.40) In terms of u and c, Eqs. (4.39) and (4.40) become 2ct + 2u c, + c u 0 (4.41) ut + uu, + 2cc, = 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 dx (u =2c -t) = 0 along dt= u Let us define the characteristic variables ( and ( by = u + 2c t = -u + 2c + t along along functions Hence, Eq. (4.43) becomes = const = const Now let us consider x and t to be = const we get dx ax at dt 8a e8( dz ax at From above two equations, we get From above two equations, we get e = te (u + c) z = tr (u c) From Eqs. (4.44) and (4.45), we can obtain u + c = (3 e)/4 + t u c = (c 3e)/4 + t Substituting u + c and u- c into Eqs. (4.50) and relationship between (x, t) and (g, () as follows: z = te(3( ()/4 + (t/2)e dx -t = u + c (4.46) dx d- = u c (4.47) of and Then for = const or of ( and e. Then for costnt or = const = 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 xf = t 3 )/4 + (t2/2) (4.55) Eliminating x from Eqs. (4.54) and (4.55), we get 2( + )tr + 3(tt + te) = 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) ar + = 4c (4.58) Then Eq. (4.56) becomes txx = t"" + 3-t" (4.59) a Since from Eq. (4.57), t = A/2 + u, u must also satisfy Eq. (4.59) 3 uAA = uao + -u (4.60) If we introduce a "potential" (a, A) so that (4.61) U then Eq. (4.60) reduces to 1 AX= Poo + P (4.62) This is a single linear partial differential equation. The boundary condition at the shoreline for Eq. (4.62) is or=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 A A + (4.64) 2 2 C 1 2 +2 + 1 o2 1 a2 -u +c +-A =-(__)+ + (4.65) 2 4 2 a 4 16 2 a2 1 a2 77 = C X, --PA = (4.66) 77C 16 4 16 (.6 S-- (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 r7 or u as explicit functions of x and t. If p(a,A) is given, then Eqs. (4.64)-(4.68) will give t, x, t7 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 T7 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,A) = -8AoJ(2) sin ?2 (4.69) where A0 is an arbitrary amplitude parameter and J0 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 t7(X,t) and u(x,t) for arbitrary x and t, Eqs. (4.64) to (4.68) must be solved numerically. For specific values of x and t, a and A can be obtained from Eqs. (4.64) and (4.65) by using Newton's Method, so that 7(x, t) and u(x,t) can easily be obtained from Eqs. (4.66) and (4.67), respectively. To test the ability of the finite-difference model to simulate the moving boundary problem, a numerical simulation of a long wave will be performed in a rectangular basin with linearly varying water depth. All quantities used in the finite-difference Figure 4.2: Computational grid model axe dimensional since the model is developed based on the dimensional governing equations, but the final solution obtained by the numerical model will be converted to dimensionless form in order to be compared with the theoretical solution. 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 -g = 129500 cm (4.70) Therefore we have the velocity scale uo = VigL = 224 cm/sec. (4.71) and the time scale T = v-L/g = sec. (4.72) 60 The computational domain is shown in the Fig. 4.2. Here we see that the grid density in the x-direction in the vicinity of the moving boundary is higher than in the offshore region. From the grid line of 50 km to the offshore boundary, the grid space increment is held at 1000 meters. Starting at this grid line, the grid space increment decreases 10 percent successively for about 6 km until the grid space increment of 200 meters is obtained. In the region of the moving boundary, the grid space increment is fixed at 200 meters. The grid density in the y- direction is unchanged, and the space increment is 2000 meters. The moving boundary is simulated using the way of assuming a thin water layer to cover the dry region all the time. The thickness of this water layer is set to be 5 cm. To decrease the violation of mass conservation resulting from this layer, the friction with the Manning coefficient of 0.03 is considered in the region of thin water layer. However, there is no friction considered in the main computational domain since we did not consider the effect of friction in the derivation of the theoretical solution. A time step of At 30 seconds is chosen. At t = 0, the fluid is quiescent. For t > 0, the wave amplitude at the offshore boundary, i?*, can be given by 77*(t) = 77(t) OL Cm (4.73) where 77(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 6 NUMERICAL SOLUTION 1.5 1.0 70"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 0o 50 Figure 4.3: Comparison between wave profiles as predicted by theory and the numerical model ( rt = r;*wx/4'g, z = z*w'2/g ) 5 0 10 20 X 30 0 so50 SI I I 5 0 10 20 x 30 4 0 50 TIME: If/3 I I 1 I I - TIME: T1/3 -. . . - 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 numerical model ( 77 = x'w'/'g, = X*w/g ) 5 0 10 20 30 40 50 X TIME: 57T/6 \ TIME: T 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 0.0 -0.5 -1.0 -1.5 -2 THE A NUM - TIME: 0 - I - ---1- w---------------2 0 2 X 4 6 ,,,, I ME: TT/3 I I I -2 0 2 X 4 6 - TIME: 2H/3 - \ - 1 -,TIME: 2Tr/3 ---- -----,w -------------1 I 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 64 ORETICRL ERICRL SO 1.5 1.0 0.5 0.0 -0.5 -1.0 -1.5 TIME: TT ------- ------------------0 2 X 4 6 Figure 4.5: Comparison between wave profiles near moving boundary as predicted by theory and the numerical model ( r = q*w'/lg, z = z*w'2/g ) SOLUTION LUTION - TIME: I/6 - -------------------2 0 2 X 4 6 ,, TIME: TT/2 I I I -2 0 2 X 4 6 TIME: 5n/6 I I! -2 0 2 X 4 6 -1.5 T 2WT 3WT 4WTr TIME Tr 2T 3r 4WT TIME X=7.5 ..... = 5 - ------THEORETICAL SOLUTION NUMERICAL SOLUTION I f I S2 TW MW TIME Figure 4.6: Comparison between theoretical and numerical solutions of water elevation ( 7 = r*w'2/kg,t = wt* ) 1.0 0.5 0.0 -0.5 -1.0 1.0 0.5 '7 0.0 -0.5 -1.0 1.0 0.5 '7 0.0 -0.5 -1.0 X=15. Tr 21T 3Tr TIME )T 2T 3T 41P TIME X=7.5 -THEORETICAL SOLUTION NUMERICAL SOLUTION I I I Tr 2 Tr 3 41 TIME Figure 4.7: Comparison between theoretical and numerical solutions of velocity ( u = u*w/g,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 following formulae Tzd= P W+WYW (5.1) r.= Cd. W+W WY (5.2) where r, and ry denote respectively the wind shear stresses in the x- and y-directions, p, is the density of air, W, and Wy are the wind speeds in the x- and y-directions, respectively, and Cd is the wind shear stress coefficient. The value of Cd0 can be Figure 5.1: The sketch of Lake Okeechobee Figure 5.2: Computational grid . . . I A I I T F I T 69 obtained from the following formulation proposed by Garrett (1977) Cd,, = 0.001 X (0.75 + 0.00067 x VW+ W,2) (5.3) where the unit of W. and W. is cm/sec.. The maximum value Of Cd. is 0.003. The bottom friction can be calculated from the Manning-Chezy friction formulation. As discussed previously, the computer programming becomes very complicated when using the weir formulation to simulate the two-dimensional moving boundary problem. Thus the method proposed by Reid and Bodine (1968) is not applied to study the wind driven circulation of Lake Okeechobee. The method to simulate the moving boundary problem with a thin water layer on the dry land area is employed here. To illustrate the significance of the moving boundary to the actual problem, the simulations will be conducted with and without the moving boundary. For the moving boundary case, the grass area is included in the computational domain and the thickness of the thin water layer on the dry area is assumed to be 10 cm. For the fixed boundary case, the grass area is taken out of the computational domain and a vertical wall is assumed to be located at the edge of the grass region. Thus the boundary between the grass area and open water can be considered as a closed boundary. A numerical time step of 3 minutes is used here. The Coriolis parameter, fl, in Lake Okeechobee is 0.73 x 10-1 sec.-' (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 in/sec.. After the wind stops, a seiche in the lake will be produced. When the wind is applied over the lake, the Manning coefficient is chosen to be 0.02 ( Schmalz, 1986). But after the wind stops, the Manning coefficient is taken as 0.005 in order that the seiche can last for a longer time for studying the shoreline motion. The steady state of wind-driven circulation in the lake is reached after the wind is applied for 7 hours. Figures 5.3 and 5.4 show the circulations for the moving WIND 150 CM/SEC. //,,if,,-,;<.--- . , .-1 ( # ., , Figure 5.3: Wind driven circulation with moving boundary Fge IND 150 CM/SEC. F i g u r e 5 .4 : W i n d d ri v e n c i r u f x b Figure 5.4: Wind driven circulation with fixed boundary 71 boundary case and the fixed boundary case respectively after the wind has been applied for 8 hours. It is seen that the circulation for the two cases are quite similar, although the velocities in the moving boundary case appear to be larger than those in the fixed boundary case. This can be explained by the fact that the area of the lake in the moving boundary case is greater than that in the fixed boundary case. Figure 5.5 shows the variation of the wind speed with time. Figure 5.6 shows the motion of the western shoreline during the first cycle of seiche oscillation in the lake after the wind has stopped. As mentioned previously, mass conservation is violated by assuming that a thin water layer exists on the dry area at all times. Figure 5.7 shows the extent of the mass conservation violation during the entire computational period. It is found that an extra mass of 0.2 percent of the total water mass in the lake is induced due to the consideration of the shoreline motion. Figure 5.8 presents the variation of the total dry area in the lake with the time. Figure 5.9 shows the comparisons of the transient variations of the water surface elevation at points A, B and C between the moving boundary case and the fixed boundary case. The locations of points A, B and C in the lake are shown in the Fig. 5.1. Figures 5.10 and 5.11 present the comparisons for the velocities U and V, respectively, between the moving boundary case and the fixed boundary case. From Figs. 5.9, 5.10 and 5.11, it is seen that there is a phase lag between the results for the moving boundary case and the fixed boundary case. Figures 5.12 and 5.13 present the transient variation of the bottom friction in the x- and y-directions at point A, respectively. 0. LO =f -20. 0 u-i 0. 0 -20.i 2 2 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 0.4 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 POINT A 2DO 150 10050 -50 -100 0 5 10 15 20 2! TIME (HR) 5 10 15 20 25 TIME (MR) 5 10 15 20 25 TIME (MHR) Figure 5.9: Comparisons of water surface elevation with moving boundary and fixed boundary oo r -50 0-1 S 0 -J UJ =', -SO LU I. -100 0 LUU z so > 0 LU -J LU c- -50 -50 100 0 75 SOLUTION WITH MOVING BOUNDRRT ---------- SOLUTION WITH FIXED BOUNDRAT 200 ~POINT R Li 100 0---------------------------------------------------------------U -100 0 -j -00 0 5 10 15 20 25 TIME [MR) 200 -- POINT I U 100 Li200 0 -100 U -200 0 5 10 15 20 25 TIME [MR) 200PONC L; LU 200 0i -100 22 TIME (HR) Figure 5.10: Comparisons of velocity U with moving boundary and fixed boundary - SOLUTION WITH MOVING BOUNDARY .--------- SOLUTION WITH FIXED BOUNDRRY 200 POINT A La 100 0 ----- ------------------------------- ---- ----- ------------- -------------U -100 -j ,=J -200' 0 5 10 15 20 25 TIME (HR) 200 - POINT B L 100 >0 -&J O-10 LU -200i 0 5 10 15 20 25 TIME (MR) 200 - POINT C 100 > 0 - 200 0 5 10 15 20 25 TIME (MR) Figure 5.11: Comparisons of velocity V with moving boundary and fixed boundary La-00 -200 T 0 5 1 15 20 2S TIME (HR) Figure 5.11: Comparisons of velocity V with moving boundary and fixed boundary 10 0 -10 0 5 10 15 20 25 TIME (HR) Figure 5.12: Transient variation of bottom friction in the x-direction at point A 5 10 15 20 2 TIME (HR) Figure 5.13: Transient variation of bottom friction in the y-direction at point A CHAPTER 6 CONCLUSIONS 6.1 Conclusions The objective of this study is to develop a two-dimensional moving boundary numerical model which can predict the hydrodynamics in lakes, estuaries and shallow 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 nonuniform rectangular grid. The governing equations, which are vertically-integrated Navier-Stokes equations of water motion, are solved using the method of fractional steps. The transition from one stage of the computation to the next is divided into three steps which are advection, diffusion, and propagation. Different numerical schemes are employed for each computational step. The method of characteristics is used for the advection, the ADI method is applied to the diffusion, and the conjugate gradient iterative method is used for the propagation. The numerical schemes used in the model are implicit so that there is no limitation for the choice of the numerical time step. Two methods for simulating the moving boundary problem are discussed in this study. The first, which was proposed by Reid and Bodine (1968), is examined briefly. It is found that it is difficult to determine the values of empirical coefficients in the weir formulation since they are very site-dependent. It is also found that the computer program for simulating the motion of the boundary is very complicated for two-dimensional problems. The second, which was proposed by Benque et al. (1982), is studied in detail. The advantage of this method is that the computer program for simulating two-dimensional moving boundary problems is very simple. The disadvantage is that the mass conservation is violated slightly. For the verification of the model, four analytical solutions are used to compare with the numerical solutions. From these comparisons, it can be concluded that the consistency and the accuracy of the model are acceptable. It is also found that the method of fractional steps is a powerful means of solving the complicated problems in several variables although its consistency has not been completely proved. In order to validate the model's ability to simulate the motion of the boundary, wave propagation on a linearly sloping beach is studied theoretically and numerically. It is found that the moving boundary model developed using this method can simulate moving boundary problems reasonably well although the mass conservation is violated slightly. The model is applied to the wind-driven circulation in Lake Okeechobee, Florida. Comparison is made between the results obtained from the moving boundary model and the fixed boundary model. From this comparison, it is seen that the boundary motion can not be neglected when studying wind-driven circulation in Lake Okeechobee. 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 ManningChezy formulation cannot give a correct estimation of the bottom friction since it breaks down when the water depth approaches zero. Therefore, basic hydrodynamic research is needed in this area. For example, Yih (1963) investigated the stability of liquid flow down an inclined plane and Melkonian (1987) studied nonlinear waves in thin films. When using a uniform or non-uniform rectangular grid to approximate the com- 80 plex geometry of water bodies, a large number of grid points is generally required. As a result, the computational cost is increased greatly. To avoid this, it is useful to modify this model to work in a curvilinear grid. In addition, we can use the method of fractional steps to develop a threedimensional 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 solution 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 simultaneouse linear algebraic equations, Ax=k (A.1) in which A is an N x N matrice of coefficients, x is an N x 1 vector of unknowns, and k is an N x 1 vector of constants. In the following description, the symbol (p, q) presents the inner product of the vectors p and q, or the value of pTq. Let us suppose that a set of N "A-conjugate" or "A-orthogonal" vectors {pi}, i = 0, 1,2,... N 1, is available to us. This means that the inner product (Api,pj) = 0 where i # j. If A is positive definite, then (Api,pi) > 0. In this case, since the vectors {pi} are necessarily linearly independent and span the Ndimensional space, the solution vector h can be written as h = cOP + C1pi + C2P2 + "". + CN-lPN-i (A.2) 82 where {ci} are coefficients. If we can determine {ci}, the solution h can be quickly calculated. Since (k,p,) = (Ah, p) (A.3) (Ah,p,) = co(Apo,p,) + ci(Api, pi) + + ci (Api,p i) + . + CN-1(APN-1,Pi) = c,(Ap,,pi) (A.4) we get S(k, p) (A.5) S(Ap, pi) Consequently (k, po) + (k,pi) (k, PN-1) PN1 (A.6) (Apo,Po) (Api,pi) (ApN-1,PN-1) This computation of h defined by equation (A.6) also can be described as the following iterative scheme x1 = opo (A.7) X+ = zi + Pip; (A.8) where S= (k, p,) (A.9) (Api,p,) and Po, PI,* ,PN-1 is the given set of N A-conjugate vectors. The iteration can be terminated when zM = h where M < N. A.2 Conjugate Gradient Method In conjugate gradient method, a particular set of A-conjugate vectors, {pi}, is developed and a solution to equation (A.1) formed in terms of these. Introducing residual vectors, ri+1 = k A xi+ (.A.10) 83 Beckman (1958) used Gram-Schmidt Orthogonalization procedure to obtain the A-conjugate direction vectors {pi} They are expressed as Pi+1 = ri+1 + fip; (A.11) in which I 1 2 = ri+12 (A.12) where I ri+ I and I ri represent the length of the vectors ri+1 and ri. The fundamental conjugate gradient iterative procedure leading to a solution of equation (A.1) can be defined by following formulae: Po = ro = k A zo (A.13) (p, r') (A.14) (p,, Api) xi+ = Xi + aipg (A.15) ri+l = k A z~+1 (A.16) S= ,+, 12 (A.17) | ri 1I Pi+, = ri+1 + iP, (A.18) After M iterations, with M < N, ZM will be equal to the solution h if all conputations 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 AZ1 (z, y), AZ2(z,y) and q(z,y) which satisfy Aa (h"M&Z) 8(dZIa9z") C a (El! A Zx) Azi cc 2{. az aU + Z f q (A.19) 2g(At)2 8x + xz gAt 8x 84 e,, ) aVn+2/3 r - {(h) } a (AZ ) v--- )- f2 + q (A.20) 2g(At)2 ay ay gAt ay AZi(X, y) = AZ2(z, y) (A.21) and subject to the open boundary and fixed boundary or moving boundary conditions. In order to present the way to solve above equations clearly, we need to write the above equations in the matrix form as [A,] [AZ] = [f] [B,] [q] (A.22) [Au] [AZ] = [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, [AZ] 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 [AZI] 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 A 0 B]Z[Az] B2 0 B z[q] (A.25) 0 Ay A-Z2 f2 By B ~ 0Z (A.26) Bz Az Z2 = Denoting 0 A = A, AZ2 = AZ, B B, and [q] = q, Eqs. (A.25) and (A.26) become AAZ = 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 B q) (A.29) Substitution of AZ into Eq. (A.28) yields ( BT A-' B) q = BT A-' f (A.30) where (BTA-1B) is a symetric matrice. The conjugate gradient method is used to solve Eq. (A.30) for q. The iteration on q consists of looking for qm+I by given qm' where m means the mth iteration. The residual vector is r m+1 = BT A-1 f (BT A-1 B) qm+l = BT (A- f A-1 B qm+l) = BT AZm+l 1 = AZm+'- +AZ 1 (A.31) Coefficient 6m can be calculated by m I r- m+1 2 I rm 2 I,J I,J Th+[Az i Az 1 AM 2 (A.32) I., (A (+,j -- 2 (i) ij ii T ,3 s The direction vector is pm+1 = rm+1 + OmP (A.33) The coefficient a, is (pm, rm) S (pm,BT A-1 Bpm) [pm]T [BT AZm( [pm] T [BT A-' B pm (A.34) Defining A-'Bpm = Hm, therefore we have A HM = B pm (A.35) Comparing Eq. (A.35) with Eq. (A.27), it is seen that Hm can be obtained by using the double-sweep method to solve H, a(h"E) +8(H,) a 8("""" H,) a'{ ae + + g h B, pm (A.36) 2g(At)2 ax x gAt Ox H2 2(h" 8) a (Vn8(+2/ H2) H { a2) + a(H-) +a "- = By pm (A.37) 2g(At)2 ay ay gAt y where [ ] = H. Consequently am = [pn]T (BT AZmJ S= [pm]r [BT Hmi I,J I,J = .[21Z(,)- AZ2(i )]p,/j) E[Hmi,.) H2,g] pmj (A.38) The iteration procedure for q can be sumarized as: (1) Let the final value of q at the previouse time step be the initial approximation to the solution of q at the new time step. (2) Apply the double-sweep method to solve Eqs. (A.19) and (A.20) for AZ and AZ2. (3) Calculate the residual vector by ro = AZ AZ. (4) Let po = ro (5) Use the double-sweep method to solve Eqs. (A.36) and (A.37) for H and 87 (6) Compute the coefficient ao using Eq. (A.38). (7) Advance q by q' = qo + ao po (8) Determine successively r~m+1 = Am+ Azm+1 (A.39) I,J I,J M-- ~ m 2 (A.40) An = '-[AZI'(+,) AZ2(dj)]2/ "[AZ(,i) (A.40) i,i i' pm+1 = rm+1 + m Pm (A.41) I,J I,J am+1 = m-')[A A +m )l/-[Hi7 ) H,+lI P+)1 (A.42) 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 az2 au2 + 2 = 0 (B.1) at ax au+ ghZ ao {sin[2k(l ) + 2wt] at 8 8= hcos2 ki + sin[2k(l z) 2wt] + 2sin2k(l z)} (B.2) and the boundary conditions are U2 (1,t) =0 (B.3) Z2(0, t) = 0 (B.4) Eliminating U2 from Eqs. (B.1) and (B.2), we obtain a2Z2 2a2 a2 c k2 t2 4hcos k{cos[2k(i x) + 2wt] Bt 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(1 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" c'4kG']{sin[2k(1 x) + 2wt] + sin[2k(1 x) 2wt]} + [c'4kF'- c'G"]{cos[2k(l x) + 2wt] + cos[2k(l z) 2wt]} a2 c2 k2 S4hcos2 kl{cos[2k(1 z) + 2wt] + cos[2k(l x) 2wt]} (B.7) Obviousely we get - c2F" c24kG' = 0 (B.8) a'c'k2 c24kF' ckG" = c k(B.9) 4h cos2 ki Integrating Eq. (B.9) with respect to z and set the integrating coefficient as zero, we have a'k2 4kF G' = kx (B.10) 4h cos2 kl Eliminating G' from Eqs. (B.8) and (B.10), we obtain F" a2 k2 - + 4kF = h x (B.11) 4k 4hcos2ki A general solution for F of Eq. (B.11) is F, = C1cos 4kz (B.12) where C1 is a constant to be determined from the boundary conditions. A particular solution of Eq. (B.11) is a2k F = 16h cos (B.13) 16-h cos2 kl Therefore F = F, + F, a2 k = C1 cos 4kx + 6hcos kl x (B.14) 16h cos2 k1 Substituting F into Eq. (B.8), we can obtain G = C1 sin 4kz + 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 + cos2 kx]sin 2k(1 x) cos 2wt 8h cos2 kl + [Ci sin4kx + C2] cos2k(l z) cos2wt (B.16) From the boundary condition (B.3), we obtain [CI os kI +jh-a 2 kI_1-7 [CI cos4k + cos2 kl ](-2k) + C 4kcos4kl = 0 (B.17) Thus a ki C = k (B.18) C = 8h cos2 ki cos 4kl (8) Similarly from boundary condition (B.4), we can get a2kl C2 = k cos 4k tan 2ki (B.19) 8h cos2 kl cos 4ki. Substituting C1 and C2 into Eq. (B.16) and simplifying it, we consequently obtain ask I Z2 = k[zsin2k(l z) + sin2k(l + z) 8h cos2 ki cos 4ki 1 Stan 2kl cos 2k(1 z)] cos 2wt (B.20) cos 4kl Substituting Z2 into Eq. (B.1) and integrating with respect to z, we have a'w 1 U2 = I[zcos2k(1 z) + -sin2k( z) 8hcos2kl 2k 1 1 cos 2k(1 + z) + tan2kl sin2k(l z)]sin2wt (B.21) cos 4kl cos 4ki APPENDIX C PROGRAM LISTING C.1 Flow Chart C.1.1 Main Routine: T2D C.1.2 Subroutine: PROP C.2 Program listing C.2.1 Description of Parameters KIO: Maximum number of grid points in the x-direction; KJO: Maximum number of grid points in the y-direction; KX1(J), KX2(J), KY1(I), KY2(I): Grid numbers for the boundary's location; CN: Manning coefficient; F: Coriolis parameter; G: Gravitational acceleration; |

Full Text |

xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EXQKWDCLU_UPDD1S INGEST_TIME 2017-07-14T23:39:03Z PACKAGE UF00076135_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES PAGE 2 UFL/COEL-88/008 A TWO-DIMENSIONAL FINITE-DIFFERENCE MODEL FOR MOVING BOUNDARY HYDRODYNAMIC PROBLEMS by Yuming Liu Thesis 1988 PAGE 3 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 encourgement throughout this study. I would also like to extend my thanks and appreciation to my thesis committee members, Dr. James T. Kirby, Associate Professor of Coastal and Oceanographic Engineering, and Dr. D. Max Sheppard, Professor of Coastal and Oceanographic Engineering, for their patience in reviewing this paper. I would like to thank Dr. Yixin Yan, Dr. Tienshun Wu, Dr. Weichi Wang, Dr. Li-Hwa Lin, and Steven Peene for their help and suggestions in this study. Finally, I would like to express my deepest appreciation and gratitude to my financial sponsor, K. C. Wong Education Foundation Ltd. of Hong Kong, for offering me the opportunity to study at the University of Florida. ii PAGE 4 TABLE OF CONTENTS ACKNOWLEDGEMENTS ............................ ii LIST OF FIGURES ................................ v ABSTRACT .................................... vii CHAPTERS 1 INTRODUCTION ................... ............ 1 2 DEVELOPMENT OF A TWO-DIMENSIONAL TIDAL CIRCULATION M ODEL ... .. .. ... .... .... .... .... ... ... .. ... 4 2.1 Governing Equations ................... ........ 4 2.2 Numerical Scheme ................... ......... 7 2.3 Grid System ................................ 12 2.4 Finite Difference Approximation ................. .... 14 2.5 Boundary and Initial Conditions ................ .... 18 2.6 Consistency and stability ........................ 18 3 COMPARISONS BETWEEN THEORETICAL AND NUMERICAL SOLUT IO N S .....................................20 3.1 Comparison With a Linear Theoretical Solution ............ 20 3.2 Comparison With a Non-linear Theoretical Solution ......... 27 3.3 Comparison to a Theoretical Solution with Coriolis Effect ...... 31 3.4 Comparison to a Theoretical Solution with Friction Effect ...... 37 4 NUMERICAL SIMULATION OF FLOW OVER TIDAL FLATS .... 42 4.1 Properties of a Moving boundary ........................ 42 4.2 Past Study .... ................... ......... .43 iii PAGE 5 4.3 Numerical Treatment of a Moving Boundary ............. 47 4.3.1 Method to Treat a Moving Boundary with the Weir Formulation ...............................48 4.3.2 Method to Treat a Moving Boundary With a Thin Water Layer ...............................49 4.4 Theoretical Solution of Wave Propagation on a Sloping Beach ....54 4.5 Comparison of Theoretical Solution with Numerical Solution ....58 5 APPLICATION TO LAKE OKEECHOBEE ................ 67 6 CONCLUSIONS ................................ .78 6.1 Conclusions ................................ 78 6.2 Future Study ............................... 79 APPENDICES A APPLICATION OF THE CONJUGATE GRADIENT METHOD TO THE PROPAGATION STEP ............................ .. .81 A.1 Conjugate Direction Method ...................... 81 A.2 Conjugate Gradient Method ....................... ..82 A.3 Application to the Propagation Step .................. 83 B DERIVATION OF Z2 AND U2 .. .................. 88 C PROGRAM LISTING ............................ 91 C.1 Flow Chart ...................... .......... .. 91 C.1.1 Main Routine: T2D ....................... .. 91 C.1.2 Subroutine: PROP ........................ 92 C.2 Program listing .............................. 92 C.2.1 Description of Parameters .................... 92 C.2.2 Program listing .......................... ..94 BIBLIOGRAPHY ................ .................. ..125 BIOGRAPHICAL SKETCH ........................... ..127 iv PAGE 6 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 water surface elevation ......................... 24 3.4 Comparison between theoretical and numerical solutions for velocity ................. ... ........... .25 3.5 Comparison of theoretical solution to numerical results with different time steps ........................... 26 3.6 Comparison between theoretical and numerical solutions for water surface elevation with nonlinear effects ............. 29 3.7 Comparison between theoretical and numerical solutions for velocity with nonlinear effects ................ ... .30 3.8 Comparison between theoretical and numerical solutions for water surface elevation with coriolis effect ............. 34 3.9 Comparison between theoretical and numerical solutions for velocity U with coriolis effect ..................... 35 3.10 Comparison between theoretical and numerical solutions for velocity V with coriolis effect ............ ......... 36 3.11 Comparison between theoretical and numerical solutions for water surface elevation with friction effect ............. 39 3.12 Comparison between theoretical and numerical solutions for velocity U with friction effect ..................... 40 4.1 Definition sketch for wave propagation on a sloping beach ...55 v PAGE 7 4.2 Computational grid ......................... 59 4.3 Comparison between wave profiles as predicted by theory and the numerical model ( t7 = 7*w'2/1g, z = z*w2/4g ) ...... ..62 4.4 Comparison between wave profiles as predicted by theory and the numerical model ( 77 = '*w2/02g, x = z*w2/og ) .......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 water elevation ( 7 = r7*w2/42g,t = ut* ) ..............65 4.7 Comparison between theoretical and numerical solutions of velocity ( u = u*w/4g,t = wt* ) .................... 66 5.1 The sketch of Lake Okeechobee .................. 68 5.2 Computational grid ......................... 68 5.3 Wind driven circulation with moving boundary ......... 70 5.4 Wind driven circulation with fixed boundary ........... 70 5.5 Variation of wind speed with time ................. 72 5.6 Locations of the moving boundary at four different time .... 72 5.7 Extra mass introduced by considering the moving boundary ..73 5.8 Transient variation of dry area ................. 73 5.9 Comparisons of water surface elevation with moving boundary and fixed boundary .......................... 74 5.10 Comparisons of velocity U with moving boundary and fixed boundary ..... .... .... .... ... .... .... ... 75 5.11 Comparisons of velocity V with moving boundary and fixed boundary ................... ........... 76 5.12 Transient variation of bottom friction in the x-direction at point A .................................. .. 77 5.13 Transient variation of bottom friction in the y-direction at point A .... ..................... ......... .. 77 vi PAGE 8 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science A TWO-DIMENSIONAL FINITE-DIFFERENCE MODEL FOR MOVING BOUNDARY HYDRODYNAMIC PROBLEMS By YUMING LIU December 1988 Chairman: Dr. Y. Peter Sheng Major Department: Coastal and Oceanographic Engineering STo predict the hydrodynamics of lakes, estuaries and shallow seas, a two dimensional numerical model is developed using the method of fractional steps. The governing equations, i.e., the vertically integrated Navier-Stokes equations of fluid motion, are solved through three steps: advection, diffusion and propagation. The characteristics method is used to solve the advection, the alternating direction implicit method is applied to compute the diffusion, and the conjugate gradient iterative 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. vii PAGE 9 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 NavierStokes equations of fluid motion, using the finite difference technique. Many numerical schemes have been proposed to solve the shallow water equations in the development of numerical models. They all have advantages and disadvantages. 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 limitation, but it requires solving the matrix equations. For. twoand three-dimensional flows, it is very difficult to overcome the computational difficulties resulting from the sheer size of the matrices. -The alternating direction implicit method, or ADI method, avoids solving complex matrix equations, but can obtain accurate solutions only for Courant numbers less than 5 (Weare, 1979). Furthermore the ADI method is not applicable to three-dimensional problems (Yanenko, 1971). The method of fractional steps developed by Yanenko (1971), on the other hand is known to be an effective method for solving complicated multi-dimensional problems in several variables. In this scheme, the computation from one time level to the next is divided into a series of intermediate steps. For each intermediate step, the computational procedure is relatively simple, an exact solution can be obtained in some cases and the time step can be quite large. However, this scheme still has the disadvantage 1 PAGE 10 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 models were developed based on the assumption of a fixed boundary with a vertical wall located at the shoreline defined by the mean water depth. However, the shoreline boundary can actually move with time in such problems as wave runup on a beach and coastal flooding into a dry coastal region due to tides or storm surges. Some researchers have tried to simulate this problem numerically. Reid and Bodine (1968) considered the motion of the shoreline according to the water elevation and used the empirical weir formulation to compute the velocity in which water flows into or out of the dry land region. The disadvantage of this method is that the empirical coefficients in the weir formulation are very site-dependent. Yeh and Chou (1979) treated the shoreline boundary as a discrete moving boundary, but the impulsive motion of the boundary could introduce serious numerical problems. Lynch and Gray (1980) simulated the shoreline boundary as a continuous moving boundary using continuous grid deformation. This method, however, can not be applied to problems with complex topography. In this study, a two-dimensional finite-difference numerical hydrodynamic model is developed from the shallow water equations by using the method of fractional steps. Two ways to simulate the moving boundary problem are studied. For some simplified flow situations, analytical solutions are compared with numerical solutions 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 PAGE 11 3 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 onedimensional long wave propagation are developed and used to verify the accuracy of the numerical model. Theoretical solutions of tidal responses inside a rectangular basin with Coriolis effect and bottom friction obtained by Rahman (1981) will be used to test the ability of the model to simulate two-dimensional problems. In Chapter 4, two ways to simulate the moving boundary problem are discussed in detail. One way is similar to that used by Reid and Bodine (1968). Another way is to include the dry land region into the computational domain by assuming a thin water layer on the dry land region at all times. The theoretical solutions for wave propagation on a linearly sloping beach developed by Carrier and Greenspan (1958) are then presented and compared to the numerical solutions to verify the second method for the moving boundary problem. In Chapter 5, the wind driven circulation in Lake Okeechobee, Florida, including the effects of the Coriolis force and a moving boundary, is investigated numerically. In Chapter 6, conclusions are drawn and suggestions are made towards further studies on numerical simulation of the moving boundary problem. PAGE 12 CHAPTER 2 DEVELOPMENT OF A TWO-DIMENSIONAL TIDAL CIRCULATION MODEL In this Chapter, a two-dimensional implicit tidal circulation model is developed using the method of fractional steps. The vertically-integrated momentum and continuity equations, which govern the two-dimensional tide-generated currents, 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 calculate 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 obtained by vertically integrating the three-dimensional Navier-Stokes equations of fluid motion. Generally, it is assumed that the density of water over the depth is constant and the vertical pressure variation is hydrostatic, thus leading to the following continuity and momentum equations in a right-handed Cartisian coordinate system shown in the Fig. 2.1 az aU av -+ 5+ -Y = 0 (2.1) tt 8x yy au auu avU 8Z rbz -r,, a au a au + + + gh fV + [-(K ) +((K ) = 0 (2.2) at Bz By az p dz az ay By av auv avv aZ rTb -r 8 av a av S+ + gh+ U + -[ (K )+ (K ) = (2.3) at 9z 9y ay p zx azx y ay 4 PAGE 13 5 Z FREE SURFRCE Z (X, .T) TV u.U ZB( X.T) x Figure 2.1: Definition sketch for tidal equations where t is time, x and y are the spatial coordinates, Z(x,y,t) is the free surface elevation, U(x,y,t) and V(x,y,t) are the unit-width discharges in the xand ydirections, respectively, u(z, y, t) and v(z, y, t) are the vertically-averaged velocities corrresponding 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 xand y-directions, respectively, and rb, and rb, are the bottom shear stresses in the xand y-directions. rb2 and rby can be expressed as pgU U +V2 rTb = pg 2 (2.4) C2h2 and pgVVU2 + V2 b = 2h (2.5) where C is the Chezy coefficient which is R1/6 C = 8.21-cm l2/sec. (2.6) n PAGE 14 6 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 V = v h (2.8) h=Z-Zf I 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 + 5 = h(u+ v-) -u(2.9) xz dy axz y 9t and auV avV v 9v az S+ -= h(u+ v-) -v(2.10) Qz Qy x By .Qt Substituting them into Eqs. (2.2) and (2.3), we obtain aU au au aZ aZ rb -fr, 8(Ke,) a(K-) +h(, + -)h -"' + gh -+ ] = 0 ot Tz Ty at ax p az ay (2.11) av av av +z V aZ Z rbV -r h (Kn) a(Kj+) + h(u-+v-) -v--+g a + h [ a aU ]= Qt T T9y Qt 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 +u + v-=0 (2.13) At 9z ay PAGE 15 7 v"+l/3 -v u 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+1 -U"+2/3 aZ aZ 7bb -r,, At -+ gh + = 0 (2.20) V"+1 -V"+2/3 aZ Z by -r (2.21) -v + gh +-= 0 (2.21) At at ay p The terms ugand vMcome from the nonlinear terms. Compared with other terms in Eqs. (2.20) and (2.21), their values are small and could be neglected in the propagation step. 2.2 Numerical Scheme Just like the splitting of a single time step into three steps as shown above, three steps could be further split into two directional steps. For example, the advection step can be treated with the following two steps: x-advection Un+1/6 Un un+1/6 S + u" = 0 (2.22) At axz PAGE 16 8 Un+1/6 Vn n-"+l1/6 -+ ,U -=0 (2.23) At Ox y-advection U"+1/3 un+1+1/Oaun+1/3 At + Un+1V/6 = 0 (2.24) at By v -+1/3 yn+1/6 un+l1/3 + vn+1/6 a = 0 (2.25) At By where n + 1/6 represents the state after the x-advection. This scheme is implicit and unconditionally stable. For the diffusion step, an alternating direction implicit method is used which leads to x-sweep U* -Un+'1/3 a au a aBU* BUn+ ---y(K ay ) = V+1/ (2.26) rAt ax oz ay By V* -V"+1/ a avV* a aV"+1/3 (K) (K ) = 0 (2.27) 'At 8x 8z y ay y-sweep U' -u a u a auv2" 3 --(K --(K ) = 0 (2.28) V+2/3 a a* a aVn+2/3 At --(K -) -(K ) = -n U* (2.29) At Bx x Y 9y By where U* and V* are the intermediate values of unit-width discharge during the computation of the diffusion step. Yanenko (1971) showed that the ADI scheme was absolutely consistent and unconditionally stable for the two-dimensibnal heat conduction equations, which become the working equations for the diffusion step if the source terms are added. The propagation step is the most important of the three steps and its working equations are more complicated than the other two steps. By introducing a PAGE 17 9 coefficient a (0 < a < 1) for the spatial derivatives, the working equations for the propagation step can be written as Un+l -U"n+2/3 aZ az U"-+2/3 n+l + aC(gh )"+' + (1 -a)(gh )" t z+ z -h" At + a(Tb_ -")"n+ + (1a)(ult ")" = 0 (2.30) P P Vn+1 -V"+2/3 .) aZ V"n+/3 Zn+1 -Z + a(gh5-)"Y + (1 -a)(gh )" -( At At -y By h" (t + a(Ty ")+1 + (1 -a)( ")" = 0 (2.31) P P Z"+* -Z" U au av avU V V S + a+(-)" + (1a)(( )" + a( )+ (1 -a)( ) =0 (2.32) At az 8z ay 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 U"*2/1 Zn+" Zn U"+' = U"+/ -At{a(gh.)"-' + (1a)(gh-)" -(" )) az xz h" At -At{a(rs" -,=) n+l + (1 -a)('bz -r*z)"} (2.33) P P V-t+'f=V"+ at -h -n+1+ aZ V"/3 Z+'-1 Z 1 +2/3 -At {(gh+ (1 -a)(gh-)" -( )} By Ty h" -At -At{a(Tbv -r,)n+l + (1 -a)(r'b -"")"} (2.34) P P Treating the bottom friction and surface friction explicitly for clarity, and substituting U"+1 and V"+1 into equation ( 2.32) yield AZ a a AZ aZ" a QAZ az" gt 2 a (h" + AZ--) + + (h" + AZ-y)} g(At)2 Xz 8 az 9y ay ay Sa U"+'2/3 a V"+2/3 + {Z) + ( A)} = f + f (2.35) where (1 -a) aU -a aU +v2/S a, ( aZ" a a r,-r,, Sfit ) )gn+2 + a (h) + ( ) (2.36) gAt XZ git 5z TX z gTx p PAGE 18 10 (1av) L aa +2/3 a az aZ" a ra -r (1 -)(_ )aV + -(h2 )++ ( -_ )" (2.37) gAt ay -gAt y ay ( ayh gy 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 AZ1 2 a .AZ1 azn a a vU+23 a -(h + AZI ) + -AZI) = fi -q (2.39) 2g(At)2 ax ax ax gAt az h" Az22 a az2 az a a Vn+2/5 -a (h" + AZ2 ) + ( Z2) = f2 + q (2.40) 2g(At)2 ay ay ay gAt ay h A in which q(z, y) is the coordination term. If AZi = AZ2, the solution of Eq. (2.35), AZ, will be exactly the same as the solution of Eq. (2.39) or (2.40), AZI or AZ2. The details of solving Eqs. (2.39) and (2.40) for AZi, AZ2 and q are presented in appendix A. If the bottom friction terms in Eqs. ( 2.30) and ( 2.31) are treated implicitly, the first-order series expansion in terms of velocity (U, V) and water depth h can be used to linearize them as (-)n+l = (rLb)n+2/3 +, a ) n+2/3Ah + a rb)n+2/3A (2.41) p p ah p au p )n+1 = ()n+2/3 (TV)n+2/3Ah Tb)n+2/3AV (2.42) p p dh p av p where Ah = h"+ -h" = Z AU = Un+1 -UU 2/3 (2.43) AV = V"+1-V"+2/3 PAGE 19 Substituting Eqs. (2.4) and (2.5) into the above equations, we get / n+1 = ,),n+2/3 9 2UvU + V2) +2/3 p p C2 h3 + gvU2 T+ V -)n+2/3(U"+2 -U+/3) (2.44) C2 h2 + h2VU2 + V2) n(y+l ( +2/3_ 9 VVU +V2 n+2/3vZ P P +V2 h + --+ V )"+2/(V"+' -v"+2/3) (2.45) C2 h2 h. U + V2 Plugging them into Eqs. ( 2.30) and ( 2.31), expressions for U"+* and V"+1 can be extracted as U*n+2/S a az U"+ = M{U"+23 + ---Az -aAt(gh )"+Z -(1 -a)At(gh )" h" ax ax g UVWTV1 g 2U2 + V + -t ( UVU2 )n+2/3} + aAt( + U C2 h2 C2 h2U2 + V2 + 2aAtt ( U)2/3'A (2.46) V"+ = N{V"+2/3 + V--3 AZ -aAt(gh + -1 )t(gh ()" -At( /U )+ V2/)" } + aACt ( U T + 2V +2/3V+2 C2 h2 C2 h/ U2 -+V2 g 2VU + V2 + 2aAtY -( V)n+2/ (2.47) CZ h3 where M =1/,1 + agA( 2)3] (2.48) CM2 h2V/U2 + VV2 and N = 1/[ + + 2V )L+2/s] (2.49) C2 h2zU2 + V2 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"Z) (AZ1az") M+ a a(Un-+2/ AZ,) -Ma2{ as + a )+ h" = f1 -q (2.50) 2g(At)2ax ax gAt 8x PAGE 20 12 At, (2) 2a(AM ) Na a(Vn+213 A Z) 2(A Na a{ + + = Nay "I----n f2 + q (2.51) 2g(At)2 ay ay gAt ay where f 1-a au Ma aU. (hZ)" Ma(U~U+v)"+2/3 f, = --( )" + Ma + C)+2 h (2. 52) 1c av M +a 2V (hMC,)" a(V' )n+2/3 f, = -i (y -( + Ma Oy ) ( +M) (2.53) f gAt 5y 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 A -AX(i) + AX(i+1) (2.54) AXu(t) = 2-(2.54) AY,(,) = AY-() + A (2.55) 2 where AXt(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 P,(i,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,i+I/2) represents the value at the middle point between the Z grid points (i,j) and (i,j + 1). Given these, we have Pu(,,) = [P,(ij) + P,(i,-l + P.,(+ij) + P,(ilj-i)l/4 (2.56) PAGE 21 13 J III I Figure 2.2: Schematic of finite difference mesh with variable rectangular cells l(I I II I 141 -A--AZ-GRID ------------(U.V) -GRID ..---.........0..-.......V)-GBID J. J---------a,, Zn (i-.-(......I J* I J-l _------------I-I I 1*1 I*2 Figure 2.3: Computational grid definition Figure 2.3: Computational grid definition PAGE 22 14 AX,(_i_)_Ay--) AX, (,-)AY, ) P(i) -4AX,(..)AY,()PU +4AX ,(-)AY,(,y) AX,()AY,(i-1) AX.,()AY(,) 4AX(i-)AY,,(y).(i-1,+1) + 4AX.(iI)AY(j) Pu(i-L,j) (2.57) Pz(i+1/2,y) = AYz(i)p(i+,)+ 2A )Pu(ij) (2.58) 2AK"(i) 2A*'u(j) and Pz(i,j+1/2) P(i.X+l) + AX.(i) Pu(i-,y+i) (2.59) 2AXu(,i-) 2AXu(i-I) P.(i+l/2,j) and Ps(i,j+l/2) can also be evaluated by Pz(i+1/2,i) = [Pz(i,,) + Pz(i+lj)]/2 (2.60) P.(i,,+1/2) = [Pz(i,) + P(i.j+1)]/2 (2.61) For a spatially uniform grid system, Eqs. (2.57), (2.58) and (2.59) become P.(ij) = [Pu(i.+l) + P,(i.) + Pu(i-ij+1) + Pu(i-i,j)]/4 (2.62) Pz(+i+/2j) = [Pu(i,j+) + P.u(i,)]/2 (2.63) Pz(i,j+1/2) = [Pu(i,j+i) + Pu(i-,,i+)]/2 (2.64) 2.4 Finite Difference Approximation Suppose that the state of the model is known at time nAt. Then the state of the model at time (n + 1)At can be obtained by successively solving the advection, diffusion and propagation steps. From Eqs. (2.22) to (2.29), it is seen that the water surface elevation Z does not appear in the advection and diffusion steps. Therefore these steps can be executed solely on the (U, V) grid. The finite difference equations for the advection step are x-advection n+1/6 n+1/6 n+1/6 S ,+ Un. il -i -,i + u" i+,,j -, -0 (2.65) At AX(i-1) + AX,(i) PAGE 23 15 n+1/6 n n+1/6 n+1/6 V',J -" + ty Vi+li -v -/'p i + u" i+ i-1 -= 0 (2.66) At AXu(t-1) + AXu(i) y-advection n+1/3 n+1/6 n+1/3 n+1/3 Ui -t. 1 n+1/6U+l UUif27 A +1/6 ,+1 i, -0 (2.67) At Ay ) + AY(,) n+1/S n+1/6 n+1/3 n+1/3 vi -i. + ,7+1/6 i+1 -i,-1 = 0 (2.68) At l AY,(-) + AY,(y) where the subscripts i and j correspond to the (U, V) grid. After x-advection and y-advection, we get the modified velocities u"+1/3 and v"+1/3. The discharges corresponding to u"+1/3 and v"f+1/ can be calculated by U /S u+/, h", (2.69) id = j So, V+1'/3 "+1/3 (2.70) in which h,. is water depth on the (U, V) grid at time nAt. The finite difference equations for the diffusion step are as follows: x-sweep Ui -/S_ K U+ -U!, U!, -U.-.] at AX,(..) AX(,) AX(i-1) K UT+'13/ -U"n+1/ Uf+/1 Urn+1/s K ij+l rij i,j jI1 Ay,(_){ Ay() AY.(-IY)() = V+l/s (2.71) -1/3 K *+l -,i V --1,. ,At AX,(i) AX() AX,(_) K V"+t/s v"+/3s V.+1/3 7"+1/ K ,+l i,j 'i,iy -id= 0 .(2.72) y-sweep U"+2/3 -K U/,-Un,-1, Jt Ax,(,) Ax,(,) a .2-) .) PAGE 24 16 Uin+2/s U-n2/3 Uin+2/S Ui:+2/S K u" +lU 3 Uy iU 1 Ay(__-i) Ay+1 ) Ay 4 _-1) = 0 (2.73) v+M/3 -,_ K V, -V V ,V..-, at AX,() AX,() AX.(i-) V. -"+2/3 V+2/3 Vn+2/3 TV"+2/S ['Kj+1 i ._ i.j ij-1 AY.(_-1) Ay,(Y) Ay-(y-1) SUi.* (2.74) where (i,j) also corresponds to the (U,V) grid. After the diffusion step, the discharges U"+2/3 and Vn+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"+' which is executed on the Z grid. The other is the evaluation of Un+1 and V"+I which is performed on the (U, V) grid. The finite difference equations for the computation of Z"+' can be expressed as Z(i_) a2 h, AZ,(i+l,) -AZ(i,) Z iZ. 2g(At)2 AX(,-_) [hi-+11/2 AX.+ AZ1(i+1/2J) AX h!" AZi(i,) -AZl(i-,) ZP"-i, i-1/2,j--y,--Al-/J-r --S-1/2,1 AX,(i-i) A-AZ(i-/2)AX Un+2/S rrn+2/S + gAtAX,[) i+/2-AZ /2,)] tan (i-+1/2,j h( ) -1/2 1(i-1 = iq (2.75) AZ2(i) a hn AZ2(ij+l) -AZ2(ij) Zj+1 -Z 2g(At)2 AY,,) 1/2 AY,() aY,2(i Z.n _Z h lAZ2(i,y) -AZ2(i,y-1i) Z1 'I-Zj_ a n+2/3 n+2/3 + gAtAY u'AZ2(i,j+1/2) -AZ2(i,-1/2) At+u(j) /,j+1/2 i,y-i/2 = f2+q (2.76) in which (1a) U. U rr+2/3 Un+2/3 (1a) Uin+x/2,y -Un-/2,y a Ul+1 /2,jy ---l/ gAt AXu(-i) gAt AXu(i-1) PAGE 25 17 Sh Zin+1, i -_ zd nz1,, + n Z VZ _i 7 'j AXl) [ hi+1/2,j A, / AX, + pgAX(-)[(rT -rS,) +1/2, -(rbz -r),-1/2j] (2.77) S"n+2/3 _rn+2/3 f (1 -a) Vi+1/2 -J-1/2 J+1/ -i-1/2 gAt AYU(j) gAt AY(j,) Z" Z" Z" -Z" S-Y [h" +-ij+1 J n h ij-1 Y+A (j) ',i1/ AY(,) I.-1/2 y ,) "+ y [(obV -rY),y'+1/2 -(rb -,',).'-1/2] (2.78) and hil/2j = (hil,, + h,j)/2 hij/2 = (hij + h,j)/2 (2.7) AZi,i/2 = (AZi,,it + AZ)/2 Note that the subscripts i, j, i 1/2 and j 1/2 refer to the Z-grid and the discharges at the Z-grid can be calculated by solving Eqs. (2.57), (2.58) and (2.59) in terms of the discharges on the (U, V) grid. The finite difference forms for U"+* and V"+1 are Zn+1 n4+1 UP+i = U. +2/A .ag.+ +/2, -i-1/2, ti 3d -tuy.,j A3X (' Z" -Z" U"12/3 t(1 -)gh +/2j i-1/2,j + i.n hill aX.(i) + h.". ,j ',j) At-a) + At(l --Za) -at -,,7, I ( -,.2),n (2.80) P P Z-n+ Zn+1 1= + -Atagh T+12 "j-12 zn z_ V+2/3 -At(1 -a)gh", in,+/2 +1/2 + (Zn+1 -At -At(1 -) (2.81)a Swih t iis i, j, i -(2.81) rr P P in which the indices i, j, i i 1/2 and 5 i 1/2 refer to the (U, V) grid. PAGE 26 18 2.5 Boundary and Initial Conditions Two types of boundary conditions are normally encountered in the numerical simulation of tidal circulation. One is the open boundary which is an artificial termination of the computational system. The other is the shoreline boundary. At the open boundary, the water surface elevation and/or the current velocities are specified. In this model, the velocity gradients in the direction normal to the open boundary is set to zero, and the water surface elevation is specified along the open boundary. The shoreline boundary can be treated as either a moving boundary or a fixed boundary. The moving boundary problem will be discussed in detail in Chapter 4. Along a fixed boundary, the normal velocity is zero, but the tangential velocity may either be set to zero ( no-slip condition ) or satisfy the zero-slope in the normal direction ( slip condition). In this model, boundary conditions must be imposed at each of the fractional steps. Since the water surface elevation does not appear in the advection and diffusion steps, only velocity boundary conditions are imposed for these steps. For the propagation step, however, boundary conditions in terms of surface elevation and velocities must be specified along the open and closed ( fixed or moving ) boundaries. The specification of the initial conditions requires the knowledge of the free surface position and the flow field at t = 0. Usually this is impossible and the model is started with the still water condition that Z =constant and U = V = 0 at t =0. 2.6 Consistency and stability As seen previously, the numerical scheme for each of the three fractional steps is relatively simple. However, we need to show that the combination of these three steps is consistent with the original governing equations. This is hard to prove PAGE 27 theoretically because of the existence of. the two-dimensional nonlinear terms and diffusion terms. This is a drawback in applying the method of fractional steps to the problems of fluid dynamics in several variables. It will be seen in the next chapter, however, that this drawback does not appear to affect the accuracy of the numerical model. Since each fractional step is unconditionally stable, the complete scheme is unconditional 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. PAGE 28 CHAPTER 3 COMPARISONS BETWEEN THEORETICAL AND NUMERICAL SOLUTIONS In this Chapter, results of the two-dimensional tidal circulation model will be compared to four different theoretical solutions. The first comparison will primarily investigate the model's capability to simulate the propagation step. The second comparison will validate the model's ability to compute the nonlinear effects. The other two comparisons will focus on the simulation of the Coriolis and friction forces. From these comparisons, we can then assess the model's accuracy and consistency. 3.1 Comparison With a Linear Theoretical Solution Neglecting the advective terms, diffusion terms, bottom friction and wind surface stress, the one-dimensional shallow water equations are au aZ + gh = 0 (3.1) at ax aZ aU -+ --= 0 (3.2) where U = uh represents the vertically integrated velocity in the x-direction, u is the vertically averaged velocity, h is the mean water depth and Z is water surface elevation. Assuming a closed boundary at z = I and an open boundary at z = 0, the boundary conditions and initial conditions associated with Eqs. (3.1) and (3.2) are: Boundary conditions U(, t) =0 (3.3) Z(0,t) = Zo + asinwt (3.4) 20 PAGE 29 21 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 zc (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) = os ) sin wt + E A, sin k, sin wat + Zo (3.8) cos kl .= 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 relationship w2 = k2 c2 (3.9) S= k 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 SAkn cos kl sin wt = 0 (3.12) n=0 Therefore cos kl = 0 (3.13) PAGE 30 22 So (2n + l)7r k_ = 2(3.14) where n is an integer. By using initial condition (3.6), we can determine the coefficient An as -2aw If cos k(1 -x) s An T= -k sin kxdx 1,n Jo cos kl -2aw (3.15) lc(kk -k2) Consequently, we obtain acosk( -x) 0[ -2aw Z(z, t) = so inw t + E ck -k s2in k sinwt] + Zo (3.16) Substituting the above equation into Eq. (3.1), we get = ac sin k(l -) + -2aw U(o, t) = cos wt + (k ) cos c,,z cos w,t] (3.17) cos kl ,= (k2 -k2) To compare the numerical solution with a linear theoretical solution, a rectangular 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 (z = 10km), at the middle point of the estuary (z = 30km), and near the closed boundary (x = 60km). Figure 3.4 presents the velocity near the mouth (x = 10km), at the middle point PAGE 31 23 CLOSED BOUNDRRT 30K //////////////////////////// V=O Y=0 / M I z o a 0 D II 0 CLOSED BOUNDRRY 60KM Figure 3.1: Schematic diagram of a rectangular basin z Lu Lu 0 CLOSED BOUNDRRT 60KM Figure 3.1: Schematic diagram of a rectangular basin AX= AT= 2KM Figure 3.2: Computational grid PAGE 32 24 -THEORETICRL SOLUTION NUMERICAL SOLUTION 100 S X=10KM. T=15KM N z50 LU -J CC --1 0 0 I I I I I I I I I I 60 66 72 78 84 90 96 TIME (HR) -100 5 -X=30KM. T=15KM -50 cc -100 60 66 72 78 84 90 96 TIME (HR) 100 S -X=60KM, T=15KM -50 60 66 72 78 84 90 96 TIME (HR) = -00 LJ 60 66 72 78 8u 90 96 TIME (HR) 0 Figure 3.3: Comparison between theoretical and numerical solutions for water surface elevation -_J -100 60 66 72 78 64 90 96 TIME (HR) Figure 3.3: Comparison between theoretical and numerical solutions for water surface elevation PAGE 33 25 THEORETICRL SOLUTION NUMERICRL SOLUTION 100 X=1OKM. Y=15KM so & & u -50 -LJ -100 60 66 72 78 84 90 96 TIME (HR) 100 X=30KM. T=15KM so 50 N 0 S-50 -100 60 66 72 78 84 90 96 TIME (HR) 100 X=5OKM. T=15KM so so U -50 -100 60 66 72 78 84 90 96 TIME (HR) Figure 3.4: Comparison between theoretical and numerical solutions for velocity PAGE 34 26 -THEORETICAL SOLUTION ------NUMERICRL SOLUTION WITH TIME STEP: 30MIN. ---------NUMERICAL SOLUTION WITH TIME STEP: 15MIN. 100 X=30KM, T=15KM .100 --------1-1---I--1-I---1-z 60 66 72 78 84 90 96 TIME IHR) S50 U -100 S-50 60 66 72 78 84 90 96 TIME (HR) Figure 3.5: Comparison of theoretical solution to numerical results with different time steps U* -j -100 I I I I 60 66 72 78 8(4 90 96 Figure 3.5: Comparison of theoretical solution to numerical results with different time steps PAGE 35 27 (z = 30km), and near the closed boundary (z = 50km). Both Fig. 3.3 and Fig. 3.4 show that there is a reasonably good agreement between theoretical and numerical solutions, even though a slight phase shift between the theoretical and numerical results exists. To investigate the origin of the phase shift, different time steps for the numerical computation are used and the numerical results for water surface elevation and velocity are presented in Fig. 3.5. From these results, it is seen that the phase shift between theoretical and numerical solutions tends to diminish as the time step decreases. 3.2 Comparison With a Non-linear Theoretical Solution When nonlinear terms are included in two-dimensional shallow water equations, it is impossible to obtain an analytical solution. In one-dimensional nonlinear tidal motion, however, we can use harmonic analysis to develop a theoretical solution. Obtaining exact theoretical solutions, however, is still difficult because the high order terms are difficult to solve for. In this study, only the zeroth and first order harmonic solutions are developed. Without the consideration of Coriolis and bottom friction forces, the governing equations for one-dimensional nonlinear tidal motion can be written as aU aU az + u 5 + gh = O (3.18) az aU -t + -= 0 (3.19) where h is assumed to be constant. The bounary conditions are Z(0,t) = a sinwt (3.20) U(, t) =0 (3.21) Based on the idea of harmonic analysis, the theoretical solution of Eqs. (3.18) and PAGE 36 28 (3.19) can be expressed as Z = Z1 + Z2 (3.22) U = U + U2 (3.23) Ux + U2 u -+ (3.24) where Z1 and U1 are the theoretical solutions of one-dimensional linear shallow water equations which have been obtained in the preceding section: Sacosk(I -z) os = k) sin wt (3.25) cos kli U ac sin k(I -z) cos = 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 + g h cos k{sin[2k(1 -x) + 2wt] + sin[2k(l -x) -2wt] + 2sin2k(lx)} (3.27) a + -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 kl[xsin2k(x) + in sin2k(I + x) 8h cos2 ki cos4ki s -tan 2kl cos 2k(l -x)] cos 2wt (3.31) cos4kl a2w 1 U2 I[z cos 2k(l -X) + sin 2k( -z) 8h cos2 ki 2k I I 4l cos 2k(l + z) + tan 2kI sin 2k(l -x)] sin 2wt (3.32) cos4kc cos4kli PAGE 37 29 -THEORETICRL SOLUTION NUMERICAL SOLUTION -100 SX=10KM, T=15KM z 50 0 LU tLU -SO 60 66 72 78 84 90 96 TIME (HR) -100 u X=30KM. T=15KM so 0 66 72 78 90 96 TIME (MR) 100 60 66 72 78 84 90 96 TIME (HR) X=60KM. T=15KM NJ 50 M -50 -100 I I I I I I I 60 66 72 78 84 90 96 TIME (HR) Figure 3.6: Comparison between theoretical and numerical solutions for water surface elevation with nonlinear effects PAGE 38 30 -THEORETICRL SOLUTION SNUMERICRL SOLUTION 100 X=10KM, T=15KM u, 50 Mu -s50 -100 ,-J -100 I -I I I I I I I I I 60 66 72 78 84 90 96 TIME (HR) 100 X=30KM. T=15KM 0 -50 -J UJ -100 I I I I I 60 66 72 78 BLI 90 96 TIME (HR) 100 X=50KM. T=15KM u so 00 -50 -100 60 66 72. 78 84 90 96 TIME (HR) Figure 3.7: Comparison between theoretical and numerical solutions for velocity with nonlinear effects PAGE 39 31 Details on the derivation of Z2 and U2 is presented in appendix B. The same basin and forcing conditions of the preceding section will now be used to compare theoretical solutions with numerical results obtained with a time step of 30 minutes. The comparisons between the numerical results and the theoretical solutions are shown in Figs. 3.6 and 3.7. Figure 3.6 presents water surface elevations near the mouth, at the middle point, and near the closed boundary of the rectangular basin. Figure 3.7 shows velocities near the mouth, at the middle point, and near the closed boundary. From Figs. 3.6 and 3.7, it can be seen that the numerical results are reasonably similar to the theoretical solutions. 3.3 Comparison to a Theoretical Solution with Coriolis Effect Neglecting convection, diffusion, and bottom friction, the two-dimensional linear shallow water equations can be written as BU BZ au+2az n =0 (3.33) av 2az + -+ flU = 0 (3.34) at ay az av av -+ -+ --= 0 (3.35) where fl is the Coriolis parameter, and c = V/gWi 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. PAGE 40 32 By using the concepts of Kelvin waves and spectrum of Poincare waves, Rahman (1982) proposed the theoretical solutions of Eqs. (3.33), (3.34) and (3.35) as follows: Z = Aexp[k-y +ik(lz) -iwt] w flk + AR exp[-(b -y) -ik(l -z) -iwt] + {C,[cosKi,(b -y)+ -n 'sinKi(b -y)] n=l wKln exp[-iK2n,( -z) -iwt]} (3.40) S Akc2 [lk U --exp[-y + ik(l -x) -iwt] w w ARkc2 ,k + -exp[(b -y) -ik(l -z) -iwt] w w + -E{Cn[K2 cosKi -y)(by) + k n=l KInC ezp[-iK2n,( -z) -iwt]} (3.41) ic2w2l 00 n2K2 V = -{CnKn(1 + ,) sin Kin(b-y) W2 -n2 n=1 w02Kl n=1 in exp[-iK2,(l -z) -iwt]} (3.42) where A is the Kelvin wave amplitude at x = I and y = 0, R is the reflection coefficient of Kelvin waves, k and w are wave number and frequency of Kelvin waves, respectively, Kin and K2, are wave numbers of the nth Poincare mode with respect to the yand x-directions, respectively, and Cn is the amplitude of the nth Poincare mode. The dispersion relationship for Kelvin waves is w = k2c2 (3.43) For Poincare waves, we have K2 22 K1. + K = 2 (3.44) From the boundary condition of V (x,0,t) = 0, it is easy to obtain n7r Ki. = -(3.45) PAGE 41 33 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) + s-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] + C,,[cosKin(b -y) + 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 PAGE 42 34 ----THEORETICAL SOLUTION SNUMERICAL SOLUTION -100 y -X=1OKM. T=15KM N 50 a-50 icU S-100 60 66 72 78 84 90 96 TIME (HR) 100 -X=30KM. T=15KM TIME (HR)o -100|-------------------------50 -J = -50 -100 I I I I I 60 66 72 78 84 90 96 TIME (HR) 100 5 X=60KM. T=15KM 0 -j c -50 -100 I I I I I 60 66 72 78 8U4 90 96 TIME (HR) Figure 3.8: Comparison between theoretical and numerical solutions for water surface elevation with coriolis effect PAGE 43 35 ----THEORETICRL SOLUTION NUMERICRL SOLUTION 100 X=10KM. T=15KM so ,-S 0 U-J -100 't '' I 60 66 72 78 84 90 96 TIME (HR) 100 -X=30KM. T=15KM u so C-j U. n O U -50 -J -100 I I I I I 60 66 72 78 84 90 96 TIME (HR) 100 -X=50KM. T=15KM u 50 so -s 0 S-50 -J -100 60 66 72 78 84 90 96 TIME (HR) Figure 3.9: Comparison between theoretical and numerical solutions for velocity U with coriolis effect PAGE 44 36 -THEORETICRL SOLUTION SNUMERICAL SOLUTION 50 X=30KM. T=15KM L 25 UJ 0 S-25 >-25 -50I I I I I I I I I 60 66 72 78 84 90 96 TIME (HR) THEORETICAL SOLUTION NUMERICAL SOLUTION 50 X=60KM. T=15KM J 25 U L) Iw -25 -50 I I I I I 60 66 72 78 811 90 96 TIME (HR) Figure 3.10: Comparison between theoretical and numerical solutions for velocity V with coriolis effect PAGE 45 37 the flow velocities u in the x-direction at three points in the middle axis of the basin. Figure 3.10 shows the flow velocities v in the y-direction at the middle point and at the closed boundary. Both the theoretical and the numerical solutions at x = 10km are very small, thus are not shown in the figures. These figures illustrate that the model results agree quite well with the theoretical solutions. 3.4 Comparison to a Theoretical Solution with Friction Effect With the bottom frictions, the two-dimensional shallow water equations can be written as + gh a+ =0 (3.48) at 9x p 0+ gh-a + -= 0 (3.49) Qt 9y p az au av BZ QU QV -+ + = 0 (3.50) 9t ax 9y where p is the density of water. Assume that the water depth h is constant and the bottom frictions can be calculated by linear friction formulae rt = pFU (3.51) rby = pFV (3.52) where /U2 + V2 F g U a (3.53) C2h2 (3.53) Given these, Eqs. (3.48), (3.49) and (3.50) can be simplified as QU QZ -+ -x + FU = 0 (3.54) V aZ( + c+ FV = 0 (3.55) at ay PAGE 46 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 K(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 -iu -c2 oo + F= Z [AnK2, cos Kin(b -y) Sn=1 (sin K2,x + tan K2nl cos K2,n)] exp(-iwt)} (3.58) C2 oo V = Re{FF[AnKx sin Ki(b-y) F -iwn=l (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 = (1 + i-) (3.61) c w b k 2 = F (36r 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) PAGE 47 39 THEORETICAL SOLUTION NUMERICAL SOLUTION -100 SX=10KM, T=15KM -50 -J IU S-100 I I I I -I I -I I 60 66 72 78 84 90 96 TIME (HR) SX==30KM. T=15KM 50 0 100 TIME (HR) -100 z 10o Â£ -,*Â•. .,., X =60 KM, T = 1S KM ... I> 0 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 -00 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 PAGE 48 40 -THEORETICRL SOLUTION NUMERICAL SOLUTION 100 -X=10KM. T=15KM U SO Sso 9 -50 -100 60 66 72 78 84 90 96 TIME (HR) 100 .X=30KM. T=15KM S50 UJ o so 0 > _io ----Y-'-'-'--'-' -U -50 -100 I I I I I 60 66 72 78 81L 90 96 TIME (HR) 100 X=50KM. T=ISKM 0 ui S50 10 0 -50 -J --I00 I I I l I i l I 60 66 72 78 84i 90 96 TIME (HR) Figure 3.12: Comparison between theoretical and numerical solutions for velocity U with friction effect PAGE 49 41 A, = b Z, cos Ki,(b -y)dy (3.65) To compare the theory with the numerical results, the same basin defined previously is used. The amplitude of the forcing tide along the mouth of the basin is assumed to be constant (50 cm), which leads to no flow in the y-direction, i.e. V = 0, and zero value for the coefficients A1, ..., A,. The period of the forcing tide is 12 hrs. The Manning coefficient is 0.02 and the maximum velocity is 50 cm/sec.. Thus the theoretical solutions for the tidal responses inside the basin can be easily obtained from Eqs. (3.57) and (3.58). With a time step of 30 minutes, the numerical results can be obtained by the use of the grid system shown in Fig. 3.2. The comparisons between the theory and the model results are shown in Figs. 3.11 and 3.12. Figure 3.11 shows the variation of the water surface elevations with time at three points inside the basin. And Fig. 3.12 presents the comparison for the velocities. From both figures, it is seen that the numerical results are very close to the theoretical ones. PAGE 50 CHAPTER 4 NUMERICAL SIMULATION OF FLOW OVER TIDAL FLATS In the computation of flow over tidal flats, the difficulty is to simulate properly 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 amplitude, 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 expressed 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. 42 PAGE 51 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. Usually the topography is given, but the water surface gradient on the moving boundary is unknown and is related to the velocity of the boundary motion implicitly. The third problem arises from the computation of the bottom friction in the vicinity of a moving boundary. The most common bottom friction formulation in vertically-integrated shallow water problems is the Manning-Chezy formulation. This formulation provides finite bottom friction throughout the interior of the computational domain. However, it breaks down near the moving boundary since the computed bottom friction approaches infinity as the water depth tends to zero. Using 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 velocity 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 PAGE 52 44 this technique is simpler to implement than the technique that makes use of deforming grids, it possesses other problems. The impulsive filling of a cell with fluid often leads to numerical problems unless treated very carefully. Reid and Bodine (1968) investigated the transient storm surges in Galveston Bay, Texas. Omitting nonlinear advection and Coriolis terms, they developed a finite difference numerical model with the inclusion of the tidal flat. A uniform Cartesian mesh and staggered grid system were used. The elevation of the sea bed or land was represented by a constant value at each grid point within a square grid. Hence, the actual topography was approximated in a stair-step fashion. The movement of the boundary was controlled by the water elevation. If the water elevation in a flooded square was less than the base elevation of an adjacent dry square, then a zero normal flow boundary condition was applied along their common boundary. However, if the water elevation in a flooded square was greater than that of an. adjacent dry square, then the water was permitted to flow into the dry square. The flow rate between two squares was determined using an empirical equation for flow over a broad-crested barrier. The overtopping of a barrier could be treated also. However, the model could treat only barriers aligned along the grid mesh division. Flow across the barrier was permitted when the water height on one side exceeded the barrier height. If the water height exceeded the barrier height on both sides, then the flow rate was determined using an empirical equation for flow over a submerged weir. The empirical coefficients in the model were determined by iteration, comparing the model with tidal data and data from hurricane Carla (September 9-12, 1963). The gross features of the inundation were predicted reasonably well. However, it should be noted that the empirical coefficients used could be very site-dependent. Yeh and Yeh (1976) developed a nonlinear nondispersive moving boundary model for simulating storm surge using an ADI technique. The shoreline in this PAGE 53 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 boundary 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 predicted 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 PAGE 54 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 donatoracceptor 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 advection, diffusion and propagation steps. A different numerical scheme was applied at each computational step. The treatment of the boundary motion was considered at the propagation step. The dry land region was assumed to be covered by a thin water layer. The flow in this region was governed by the bottom friction. During the computation of the propagation, the shallow water equations were first applied to the whole computational domain including the region of the thin water layer. Then the solution in the vicinity of the moving boundary needed to be adjusted to satisfy the resistance flow. The Manning-Chezy friction formulation was used to compute the bottom friction. They reported that the moving boundary model developed in this way violated the continuity equation slightly. Good agreement of the numerical results with the measured data were presented based on the model application to the Bay of Saint Brieuc and the River Canche Estuary, France. Lynch and Gray (1980) outlined a general technique whereby a moving boundary PAGE 55 47 can be treated by finite element Eulerian models. The finite element basis functions are chosen to be functions of time so that the element boundaries track the moving shoreline. They showed how this motion generates extra terms which, if treated properly, reduce the problem to one that can be treated by standard finite element procedures. They showed how to apply the method to treat the propagation and runup of long waves. They looked at two simple problems involving the runup of waves on plane beaches. They showed that estimating the runup by extrapolating the wave height at a vertical wall could introduce significant errors. They also pointed out that treating deforming elements is computationally more expensive than fixed ones, and recognized that potential problem can arise if the mesh becomes geometrically too distorted. Gopalakrishnan and Tung (1983) described a finite-element nonlinear long wave runup model valid only for one horizontal dimension. The model contained terms that accounted for vertical accelerations. The moving shoreline was handled by allowing the shoreline element to deform so that the beach node always tracked the shoreline. If the shoreline element became too stretched, it split into two elements. The element containing the shoreline node continued to deform but the other new element created by splitting stayed fixed. They showed some plots that detailed the runup process, but they did not present any results about the rundown process. The technique outlined by the authors seems applicable to tsunami runup, but they did not present a thorough or convincing argument to show that their model could be used reliably for such studies. It should be noted that the technique in their work cannot be extended easily to include two horizontal dimensions since the element-splitting procedure would be very complex. 4.3 Numerical Treatment of a Moving Boundary In this section, two ways to simulate the moving boundary problem will be studied in detail. The first method was proposed by Reid and Bodine, and the PAGE 56 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 twodimensional stair-steps. Thus the elevation of the sea bed or land can be regarded as uniform over each grid square. The boundary motion is controlled by the water level. During the flooding, a grid point is added to the computational system if its water depth is greater than a preset value h, which is the minimum water depth for the effective application of the Navier-Stokes equations. The value of hi also depends on the topography and the grid system and is usually taken to be 10 to 20 cm. During the.drying, a grid point is taken out of the computation if its water depth is decreased to a preset value, h2, which is usually slightly different from hi. In the computation of flooding, the condition of no normal mass flux is applied at the moving boundary when the water elevation is less than that of the adjacent dry land, i.e., the normal component of flow, Q,, at the juncture of a flooded cell and a dry cell is taken as zero, while the tangential component of flow may satisfy the no-slip or free-slip condition. However, if the water elevation is greater than that of the dry land, flow will be allowed to flood into the dry cell until the water depth in this cell reaches hi. The mass flux per unit width of the dry cell, Q,, can be calculated by the weir formulation ( Reid and Bodine, 1968) as: Q. =CoDD gI Zd-Z (4.4) where Zd and Z, are the water surface elevation in the donator and acceptor cells, respectively; Co is an empirical dimensionless coefficient which was suggested to be PAGE 57 49 less than 0.5 by Reid and Bodine (1968); and D = Zd -Zb (4.5) where Zb is the bottom elevation of the acceptor cell. During the period of At, the increment of water level in the acceptor cell, AZ,, is AZ, = (4.6) A, where A, is the area of the acceptor cell and B is the width of the acceptor cell. The decrement of water level at the donator cell, AZd, is Q,, AtB AZd = A (4.7) Ad where Ad is the area of the donator cell. During the flooding computation of each time step, each donator cell must be examined to see if too much water is flooded from the donator cell to the acceptor cell. If the water level in the donator cell is less than that in the acceptor cell, the mass flux computed by Eq. (4.4) must be adjusted to make the water level in the donator cell the same as that in the acceptor cell. When the water depth becomes small during the drying, an artificial water depth has to be considered so that the Manning-Chezy friction formulation can still be used to compute the bottom friction. The value of this artificial water depth is determined empirically. Usually we can set it to be 20 cm. The development of a two-dimensional moving boundary model in this manner is very cumbersome since the location of the new boundary must be determined at each time step. But if the grids in the vicinity of the moving boundary and the time step are kept small, we can obtain reasonable numerical solutions for the moving boundary. 4.3.2 Method to Treat a Moving Boundary With a Thin Water Layer Based on an assumption that the water flow in the region of shallow water is dominated by the bottom friction, Benque et al. (1982) proposed a way to PAGE 58 50 treat the moving boundary by controlling the water flow discharges through the adjustment of the water depth. In this scheme, a grid system is first established in the computational domain which includes the entire tidal flat region. To ensure numerical accuracy, grid sizes in the tidal flat region should be smaller than those in the main computational domain. A thin water layer is assumed to exist at all times over the dry land region so that all grids in the computational domain are always wet. Thus, no grid needs to be taken out of the computational domain during the computation and we do not need to consider the motion of the boundary. The shoreline boundary can thus be treated simply as a fixed boundary, so the Navier-Stokes equations of fluid motion can be applied to the entire computational domain so long as the bottom friction is adequately resolved for small water depths. From the Manning-Chezy friction formulation, we see that the water depth plays an important role in the computation of the bottom friction, so it is possible to "control" the bottom friction by adjusting the water depth. As seen in Chapter 2, the free surface elevation Z does not appear during the advection and diffusion steps. Since the motion of the boundary is mainly controlled by the water elevation, it is only necessary to consider the propagation step. From the governing equations, Eqs. (2.75) and (2.76), of the propagation step, it is seen that the water depth, h, and the increment of water surface elevation during one time step, AZ, at the points of i 1/2 and j 1/2 must be calculated. They are governed by the water depth and the surface elevation at their upstream and downstream grid points, and can be calculated by the following formulae: hi+1/2 = Ti hi + (1 -yi)hi+l i-1/2 = i-1hi-1 + (1 --(4.8) AZi+1/2 = AZ.i -+ (1 -7i)AZi+i AZi-1/2 = -fi-1 AZi-I + (1 -'Y.-i)AZ. PAGE 59 51 and hj+1/2 = yj h, + (1 --j)hj+ h*-1/2 = y-l hj-_ + (1 -y--I)hj 4.9). AZ,+1/2 = -j AZ, + (1 -Y)AZ+ (4) Aj-1/2 = 'Y,-i AZ,_1 + (1 -'-Y-)AZ, ) where -is a weighting coefficient. Normally the shape of the water surface between two successive grid points can be assumed to be a straight line, so that 7 can be set to 0.5. However, in the regions where the water flow is dominated by the bottom resistance, -7 cannot be equal to 0.5 since the shape of the water surface is greatly curved. Therefore we need to look for a formulation to calculate it. Neglecting the effect of the interior force, we can obtain the governing equations for the flow dominated by the bottom resistance as follows 49Z rb. gh + =0 (4.10) aZ p gh += 0 (4.11) By P Using the Manning-Chezy friction formulation, pgUVU2 + V2(4 rb = C2 h (4.12) pgVVU2 + V2 v = C2 h2 (4.13) we can rewrite Eqs. (4.10) and (4.11) as 9Z UVU2V + x + C2h = 0 (4.14) Z VVU2 + V2 y + C2h3 = (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 S< 0 (4.17) azd, PAGE 60 52 Where the subscript ds refers to the downstream. Rewriting Eqs. (4.14) and (4.15) in the discrete form, we have Zda -Zu, U/VU2 +V2 A + =0 (4.18) Az C2 hi+1/2 Za, -Z, V V'U + V Zd+ =V0 (4.19) Ay, C2 h,+112 where the subscript us refers to the upstream. In Eq. (4.18), the upstream and downstream correspond to flow in the x-direction, but in Eq. (4.19) they refer to flow in the y-direction. Introducing a coefficient f, the water depth h in Eqs. (4.18) and (4.19) can be expressed as hi+1/2 = P hu. + (1 -/)hd, (4.20) hj+l/2 = P hu, + (1 -?)hd, (4.21) where hu, = Z,, -ZB,, (4.22) hd, = Zda -ZBd, (4.23) in which ZB,, and ZBd, are the bottom elevations upstream and downstream of the point i + 1/2 or j + 1/2, respectively. Differentiating Eq. (4.18) with respect to Zd, and bearing in mind that U is a function of Z,, and Zd,, and V is constant, we obtain aU U2 C2h ah (U + V2= -[-h + 3(Z,, -Zd,) ] (4.24) az, (U + U= +V) = z' a--, Applying the condition -< 0, we get Oh -h + 3(Z,, -Zd.) --Z < 0 (4.25) zd,' From Eq. (4.20), we know ah -1 -(4.26) 8 4, PAGE 61 53 Substituting hand 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 replaced 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 = 1P. 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) 8Z Uma, = sign(u) C2h3 | 1 (4.29) aZ Vma, = sign(v) C'h3 -z (4.30) ay in which sign(u) = az/ | a (4.31) 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 PAGE 62 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 onedimensional nonlinear shallow water equations can be written as t + a.[('* + h*)u'] = 0 (4.33) au* *u* 8r*a + + = 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* *77 X=t -r L L T eL PAGE 63 55 xx MWL Figure 4.1: Definition sketch for wave propagation on a sloping beach h x u = -(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 us + ri = 0 (4.40) In terms of u and c, Eqs. (4.39) and (4.40) become 2ct + 2u c. + c u, = 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. PAGE 64 56 If Eqs. (4.41) and (4.42) are added and subtracted, we obtain d dz c(u 2c -t) = 0 along dt = u c (4.43) Let us define the characteristic variables and ( by S= u + 2c -t (4.44) S= -u + 2c + t (4.45) Hence, Eq. (4.43) becomes dx S = const along -= u + c (4.46) dx ( = const along -= u -c (4.47) dt Now let us consider x and t to be functions of and e. Then for S = const or S= const we get -= if a = const (4.48) dt a( e ( dx x at -Tt -/ if = const (4.49) From above two equations, we get xe = te (u + c) (4.50) xt = t (u -c) (4.51) From Eqs. (4.44) and (4.45), we can obtain u c = (3" -e)/4 + t (4.52) u -c = (3)/4 + t (4.53) Substituting u + c and u -c into Eqs. (4.50) and (4.51) yields the transform relationship between (x, t) and (, )) as follows: xe = te(3 -E)/4 + (t2/2)e (4.54) PAGE 65 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 n = (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=O (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 = + = (4.64) 2 2 a PAGE 66 58 1 1 1 o 1 o2 S= + += ( )2 + -p + -(4.65) 2 4 2a 4 16 2 2 1 a2 77 = c -z = -P= -(4.66) 16 4 16 = -(4.67) c (4.68) Although Eq. (4.62) is certainly much simpler to solve than the two original coupled nonlinear equations (4.39) and (4.40), it is difficult to obtain rl or u as explicit functions of z and t. If p(a,A) is given, then Eqs. (4.64)-(4.68) will give t, x, rl and u all parametrically in terms of the variables a and A. In general, it is very difficult to eliminate a and A to obtain direct functional relationships for l7 and u in terms of x and t. 4.5 Comparison of Theoretical Solution with Numerical Solution A simple solution of Eq. (4.62) pointed out by Carrier and Greenspan (1958) is p(a, ) = -8AoJo( ) sin (4.69) where Ao is an arbitrary amplitude parameter and Jo is the Bessel function of the first kind of order zero. This potential corresponds to a standing wave solution resulting from the perfect reflection from the shore of a wave of unit frequency. With p(a, A) given by Eq. (4.69), Eqs. (4.64) to (4.68) will implicitly give the solution of this standing wave. To evaluate r(z,,t) and u(x,t) for arbitrary x and t, Eqs. (4.64) to (4.68) must be solved numerically. For specific values of z and t, a and A can be obtained from Eqs. (4.64) and (4.65) by using Newton's Method, so that rl(x,t) and u(x,t) can easily be obtained from Eqs. (4.66) and (4.67), respectively. To test the ability of the finite-difference model to simulate the moving boundary problem, a numerical simulation of a long wave will be performed in a rectangular basin with linearly varying water depth. All quantities used in the finite-difference PAGE 67 59 Figure 4.2: Computational grid model are dimensional since the model is developed based on the dimensional governing equations, but the final solution obtained by the numerical model will be converted to dimensionless form in order to be compared with the theoretical solution. 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) (W*)2 Therefore we have the velocity scale uo = -gL = 224 cm/sec. (4.71) and the time scale T = vLlg = sec. (4.72) PAGE 68 60 The computational domain is shown in the Fig. 4.2. Here we see that the grid density in the x-direction in the vicinity of the moving boundary is higher than in the offshore region. From the grid line of 50 km to the offshore boundary, the grid space increment is held at 1000 meters. Starting at this grid line, the grid space increment decreases 10 percent successively for about 6 km until the grid space increment of 200 meters is obtained. In the region of the moving boundary, the grid space increment is fixed at 200 meters. The grid density in the ydirection 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. PAGE 69 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. PAGE 70 62 -THEORETICAL SOLUTION A NUMERICAL SOLUTION 1.5 1.0 -TIME: 0 -1.5 ------------------1-----------------------05 0.0 --------0.5 -1.0 -1.5 -5 0 10 20 X 30 40 50 1.5 1.0 TIME: T/6 0.5 -1.5 I --------------------5 0 10 20 X 30 40 50 1.5 1.0 TIME: TT/3 0.5 0.O -" "=--------'.,--............>. .. .. ..,---0.5 -1.0 -1.5 -5 0 10 20 X 30 10 50 1.5 1.0 -TIME: T7/2 1.0 0.5 0.0 -0.5 -1.0 -1.5 -5 0 10 20 X 30 40 50 Figure 4.3: Comparison between wave profiles as predicted by theory and the numerical model ( i7 = r7*wx/12g, z = z*w2/4g ) PAGE 71 63 -THEORETICAL SOLUTION & NUMERICAL SOLUTION 1.5 1.0 TIME: 2T/3 0.5 -1.5 -----------------i-----------0.0 ----------------------------------------0.5 -1.0 -1.5 -5 0 10 20 30 40 50 1.0 -\ TIME: 5T/6 0.5 0.0 -1.0 -----------------------------------------------------------------0.5 -1.0 -1.5 -5 0 10 20 30 40 50 X 1.5 1.0 \ TIME: TT -0.5 -1.0 -1.5 I I I I I -5 0 10 20 30 40 50 X Figure 4.4: Comparison between wave profiles as predicted by theory and the numerical model ( 7 = x'w2/42g, z = z*'2/4g ) PAGE 72 64 -THEORETICRL SOLUTION A NUMERICAL SOLUTION 1.5 1.5 1.0 -TIME: 0 1.0 -TIME: T/6 0.5 -0.5 -10.50 ------_ ..... -----1.5 0.-0 ---1-------------2 0 2 X I 6 -2 o 2 X I 6 1.5-, -----] 1.5 0.0 1 -\TIM TE: /3 1.0 \ TIME: TT/2 -1.0 -1.0 -1.5 I I I -1.5 I I 1.5 1.5 1.5 ,1.5s 1.0 -TIME: 2T/3 1.0 TIME: 5n/6 0. -0. 0.0 ~-rn3-------------------0.0 -I -0.5 --0.5 -1.0 --1.0 -1.5 I I -1.5 I -2 0 2 X 4 6 -2 0 2 X 6 1.5 1.5 1.01.0 -TIME: 2/3 1.0 -TIME: 5T/6 0. -0.5 0.0 0.O -----,w-----------0 .0 -------------------------0.5 --5 -1.0 -1.0 1.5 1.0 TIME: T 0.5 0.0 --------------------------0.5 -1.0 -2 0 2 X 4 6 Figure 4.5: Comparison between wave profiles near moving boundary as predicted by theory and the numerical model ( 77 = q*w2/2'g, x = x'w2/og ) PAGE 73 65 1.0 X=30. 0.5 r7 0.0 -0.5 -1.0 0 T 2T 3T 4T TIME 1.0 X=15. 0.5 0.0 -0.5 -1.0 I I 0 IT 2T 3TT 4wT TIME 1.0 X=7.5 0.50.0 -0.5 ---THEORETICL SOLUTION A NUMERICAL SOLUTION -1.0 I 1 I 0 w 2T 3T 4U TIME Figure 4.6: Comparison between theoretical and numerical solutions of water elevation ( 7 = r*w'2/2g,t = wt* ) PAGE 74 66 0.4 X=30. 0.2 U 0.0 -0.2 -0.4 0 T 21T 3Tr 4T TIME 0.4 X=15. 0.2 0.0 -0.2 -0.4 0 T 2T 3T 4TW TIME 0.4 X=7.5 0.2 0.0 -0.2 --THEORETICRL SOLUTION NUMERICRL SOLUTION -0.4 I I 0 T 2T 3T 4Tr TIME Figure 4.7: Comparison between theoretical and numerical solutions of velocity ( u = u*u/g,t = wt* ) PAGE 75 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 following formulae Tz=P.Cda W+WW (5.1) rY= .o Cd W/ +W W, (5.2) where r, and ry denote respectively the wind shear stresses in the xand y-directions, p, is the density of air, W, and Wy are the wind speeds in the xand y-directions, respectively, and Cda is the wind shear stress coefficient. The value of Cda can be 67 PAGE 76 68 S\ GRASS / S REGION I o / &C EAST Figure 5.1: The sketch of Lake Okeechobee Figure 5.2: Computational grid I I I J* PAGE 77 69 obtained from the following formulation proposed by Garrett (1977) Cd, = 0.001 x (0.75 + 0.00067 x W_ + W2) (5.3) where the unit of W, and W, is cm/sec.. The maximum value of Cda is 0.003. The bottom friction can be calculated from the Manning-Chezy friction formulation. As discussed previously, the computer programming becomes very complicated when using the weir formulation to simulate the two-dimensional moving boundary problem. Thus the method proposed by Reid and Bodine (1968) is not applied to study the wind driven circulation of Lake Okeechobee. The method to simulate the moving boundary problem with a thin water layer on the dry land area is employed here. To illustrate the significance of the moving boundary to the actual problem, the simulations will be conducted with and without the moving boundary. For the moving boundary case, the grass area is included in the computational domain and the thickness of the thin water layer on the dry area is assumed to be 10 cm. For the fixed boundary case, the grass area is taken out of the computational domain and a vertical wall is assumed to be located at the edge of the grass region. Thus the boundary between the grass area and open water can be considered as a closed boundary. A numerical time step of 3 minutes is used here. The Coriolis parameter, f, in Lake Okeechobee is 0.73 x 10-4 sec.-1 (Schmalz,1986). A spatially uniform wind from the east to the west is applied over the lake surface for 9 hours. The wind speed is 20 m/sec.. After the wind stops, a seiche in the lake will be produced. When the wind is applied over the lake, the Manning coefficient is chosen to be 0.02 ( Schmalz, 1986). But after the wind stops, the Manning coefficient is taken as 0.005 in order that the seiche can last for a longer time for studying the shoreline motion. The steady state of wind-driven circulation in the lake is reached after the wind is applied for 7 hours. Figures 5.3 and 5.4 show the circulations for the moving PAGE 78 70 S WIND 150 CM/SEC. //Figure 5.3: Wind driven circulation with moving boundary W.IND Figure 5.4: Wind driven circulation with fixed bound--ary ,, ,\ Figure 5.3: Wind driven circulation with moving boundary _-^ ^" ^ ^ < "< > *' -' PAGE 79 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 xand y-directions at point A, respectively. PAGE 80 72 40. ( 20. z: .J 0o 0 -20. z -40. Â• -1 --I --I I I I 0 5 10 15 20 25 TIME (HR) Figure 5.5: Variation of wind speed with time ---............ 05 MR --11.0 HR -----11.25 HR --Figure 5.6: Locations o1.5 the moving boundary at four different time Figure 5.6: Locations of the moving boundary at four different time PAGE 81 73 0.5 0. 0.3 V x&100 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 PAGE 82 74 ----SOLUTION WITH MOVING BOUNDRRY -------SOLUTION WITH FIXED BOUNDARY 200 u POINT A S150 o 100 ( 50 ... a: > -50 3j -0 ------------------------------------------------------100 0 5 10 15 20 25 TIME (HR) 100 yu -POINT B z 50 0 Cr > O e -50 -100 0 5 10 15 20 25 TIME (HR) 100 POINT C N 0 0 ILU .0 5 10 s15 20 25 TIME (MHR) Figure 5.9: Comparisons of water surface elevation with moving boundary and fixed boundary PAGE 83 75 -SOLUTION WITH MOVING BOUNDRRT .--------SOLUTION WITH FIXED BOUNDRART 200 -POINT R UJ w 100 U -100 -J -200 1 I -I -I 1 I 0 5 10 15 20 25 TIME [HR) 200 POINT B -100 0 u -100 o .J Lj -200 0 5 10 15 20 25 TIME (MR) 200 POINT C UJ S100 0 I-. U -100 -200 I I I I 0 5 10 15 20 25 TIME (MR) Figure 5.10: Comparisons of velocity U with moving boundary and fixed boundary PAGE 84 76 ----SOLUTION WITH MOVING BOUNDART --------SOLUTION WITH FIXED BOUNDRRT 200 POINT R C 100 -J 0 --------------------------------------------------------------L -100 =J -200 0 5 10 15 20 25 TIME (HR) 200 S -POINT B U >0 S-100-200 'iiI 0 5 10 15 20 25 TIME (MR) 200 -POINT C U S100> 0 sS-100 -J LU -200 I 0 5 10 15 20 25 TIME (HR) Figure 5.11: Comparisons of velocity V with moving boundary and fixed boundary PAGE 85 77 10 x x N. <-0 LOJ IL -10 -20 0 o -20 I-----I ------I-----I -I LO 10 E-20 I-,, 0 5 10 15 20 25 TIME (HR) Figure 5.13: Transient variation of bottom friction in the y-direction at point A 10[----------------------Figure 5.13: Transient variation of bottom friction in the y-direction at point A PAGE 86 CHAPTER 6 CONCLUSIONS 6.1 Conclusions The objective of this study is to develop a two-dimensional moving boundary numerical model which can predict the hydrodynamics in lakes, estuaries and shallow 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 nonuniform rectangular grid. The governing equations, which are vertically-integrated Navier-Stokes equations of water motion, are solved using the method of fractional steps. The transition from one stage of the computation to the next is divided into three steps which are advection, diffusion, and propagation. Different numerical schemes are employed for each computational step. The method of characteristics is used for the advection, the ADI method is applied to the diffusion, and the conjugate gradient iterative method is used for the propagation. The numerical schemes used in the model are implicit so that there is no limitation for the choice of the numerical time step. Two methods for simulating the moving boundary problem are discussed in this study. The first, which was proposed by Reid and Bodine (1968), is examined briefly. It is found that it is difficult to determine the values of empirical coefficients in the weir formulation since they are very site-dependent. It is also found that the computer program for simulating the motion of the boundary is very complicated for two-dimensional problems. The second, which was proposed by Benque et al. (1982), is studied in detail. The advantage of this method is that the computer 78 PAGE 87 79 program for simulating two-dimensional moving boundary problems is very simple. The disadvantage is that the mass conservation is violated slightly. For the verification of the model, four analytical solutions are used to compare with the numerical solutions. From these comparisons, it can be concluded that the consistency and the accuracy of the model are acceptable. It is also found that the method of fractional steps is a powerful means of solving the complicated problems in several variables although its consistency has not been completely proved. In order to validate the model's ability to simulate the motion of the boundary, wave propagation on a linearly sloping beach is studied theoretically and numerically. It is found that the moving boundary model developed using this method can simulate moving boundary problems reasonably well although the mass conservation is violated slightly. The model is applied to the wind-driven circulation in Lake Okeechobee, Florida. Comparison is made between the results obtained from the moving boundary model and the fixed boundary model. From this comparison, it is seen that the boundary motion can not be neglected when studying wind-driven circulation in Lake Okeechobee. 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 ManningChezy formulation cannot give a correct estimation of the bottom friction since it breaks down when the water depth approaches zero. Therefore, basic hydrodynamic research is needed in this area. For example, Yih (1963) investigated the stability of liquid flow down an inclined plane and Melkonian (1987) studied nonlinear waves in thin films. When using a uniform or non-uniform rectangular grid to approximate the com- PAGE 88 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 threedimensional numerical model. PAGE 89 APPENDIX A APPLICATION OF THE CONJUGATE GRADIENT METHOD TO THE PROPAGATION STEP In this appendix, the algorithm of the conjugate gradient method for the solution 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 simultaneouse linear algebraic equations, Ax=k (A.1) in which A is an N x N matrice of coefficients, x is an N x 1 vector of unknowns, and k is an N x 1 vector of constants. In the following description, the symbol (p, q) presents the inner product of the vectors p and q, or the value of pTq. Let us suppose that a set of N "A-conjugate" or "A-orthogonal" vectors {pi}, i = 0,1,2,***,N -1, is available to us. This means that the inner product (APi,pj) = 0 where i # j. If A is positive definite, then (Ap,,pi) > 0. In this case, since the vectors {p,} are necessarily linearly independent and span the Ndimensional space, the solution vector h can be written as h = coPo + cpl + C2P2 + *. + CN-lPN-1 (A.2) 81 PAGE 90 82 where {ci} are coefficients. If we can determine {c,}, the solution h can be quickly calculated. Since (k, P) = (Ah, p) (A.3) (Ah, p,) = co(Apo,Pi) + ci(Api,pi) + + ci(Api,p -i) + + cN-I(APN-1,Pi) = c,(Ap,,pi) (A.4) we get S(k,p) (A.5) Ci = (Ap,p)A.5) (Api, pi) Consequently (k, o) (k, pi) (k, PN-1) (.6) h = Po + -P + + PN-1 (A.6) (Apo,Po) (Api,pi) (ApN-1,PN-1) This computation of h defined by equation ( A.6) also can be described as the following iterative scheme x1 = 0opo (A.7) zX+l = zi + Pip (A.8) where S= (k, p) (A.9) (Api,p,) and Po, Pi,*" ,PNis the given set of N A-conjugate vectors. The iteration can be terminated when XM = h where M < N. A.2 Conjugate Gradient Method In conjugate gradient method, a particular set of A-conjugate vectors, {pi}, is developed and a solution to equation (A.1) formed in terms of these. Introducing residual vectors, r+l = k -A i+ (A.10) PAGE 91 83 Beckman (1958) used Gram-Schmidt Orthogonalization procedure to obtain the A-conjugate direction vectors {pi} .They are expressed as Pi+1 = ri+1 + ,ip3 (A.11) in which I 12 where I rI+1 | and I ri represent the length of the vectors ri,+ and r,. The fundamental 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= I, 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 conputations 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(x,y) and q(z,y) which satisfy a (h"_'&) a(AZIa9z) I a (Ell AZ) Az -cc { a 9 z + A( Z f -q (A.19) 2g(At)2ax + gAt dx PAGE 92 84 AZ2 a(hfna-z2) a(AZ2 ) ,O a(n/AZ3) f+ A 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 conditions. 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] = [A] -[B], [q] (A.22) [Ay] [AZ2] = [f2] -[By] [q] (A.23) [By] [AZ1] + [B.] [AZ2] = 0 (A.24) in which [A,] and [Au] are symetrical coefficient matrix of Eqs. (A.19) and (A.20), respectively, [AZ1] and [AZ] are vectors of unknowns, [fi] 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 = PAGE 93 85 Denoting[ = A, A2 = A, 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 AB) q =BT A-' f (A.30) where (BTA-1B) is a symetric matrice. The conjugate gradient method is used to solve Eq. (A.30) for q. The iteration on q consists of looking for qm+I by given qm where m means the mth iteration. The residual vector is r"+1 = BT A-lf -(BTA-B) qm+l = B (Af -A-1 B qm+l) = BT AZm+l = AZ"+ -AZm+1 (A.31) Coefficient 8m can be calculated by m I rm+l 2 I rm 12 I,J I,J = E[Azm+l) -AZ 2 [ 2 A ,,/E) -AZM)i2 (A.32) I., (.j) --(ij) 2 (ij) iij ij 1,3 sj The direction vector is pm+1 rm+1 +. .p" (A.33) PAGE 94 86 The coefficient a, is (pm, rm) S (pm,BT A-1B pm) [pm]T [BT AZm] [pm]T [BT A-' B pmj (A Defining A-'Bpm = Hm, therefore we have A Hm = Bpm (A.35) Comparing Eq. (A.35) with Eq. (A.27), it is seen that H" can be obtained by using the double-sweep method to solve Hx a(h" ) a(R, ) a 8(~U" Hi,) 2g aa {-a } + + ag(= B, pm (A.36) 2g(At)2 ax a gAt Ox H2 2 a(hn ) aH ) a-8("Vn+" H2) 2 -ao(n{ 2) + a( -} + = B, p (A.37) 2g(At)2 ay ay gAt ay where [ = H. Consequently [pm]r [BT AZm) am [pm [BT Hi] I,J I,J = .[AZmi,,)AZ2'i,)] P(,)/ E[Hi.i) -H \j,] p( (A.38) The iteration procedure for q can be sumarized as: (1) Let the final value of q at the previouse time step be the initial approximation to the solution of q at the new time step. (2) Apply the double-sweep method to solve Eqs. (A.19) and (A.20) for AZ and AZ0. (3) Calculate the residual vector by ro = AZ? -AZ2. (4) Let po = ro (5) Use the double-sweep method to solve Eqs. (A.36) and (A.37) for HO and H,. PAGE 95 87 (6) Compute the coefficient ao using Eq. (A.38). (7) Advance q by q' = qO + ao pO (8) Determine successively rm+1= AZ+ -Am+1 (A.39) I,J I,J = -[AZI ) -AZ )]2/ /"[AZ ,) -A 1]2 (A.40) i,i ij pm+l = rm+1 +Pm pm (A.41) I,J I,J am+1 = )[A]m) -Az )1/ E[Hi7 ) -H,a+ P)1+ (A.42) -( i) / .L ij) 2(ij)J Yij) j *,y qm+2 = qm+1 + m+1 m+1 (A.43) (9) Repeat Step 8 with m +1 replacing m and continue until m = N1 or until the residual vector becomes sufficiently small, whichever condition may be satisfied first. (10) Let AZ = (AZ1+ AZ2)/2 PAGE 96 APPENDIX B DERIVATION OF Z2 AND U2 The governing equations are az2 aU, + = 0 (B.1) at ax + gh a 2 k2 {sin[2k(l -z) + 2wt] at ax 8h cos2 ki + 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 a Z2 a'c2k2 a -2 C 4hcos {cos[2k(i -x) + 2wt] at2 z 4h cos2 ICI + 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) 88 PAGE 97 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" = k (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 a2k F,= 1h cosk (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) PAGE 98 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 a2kl a =ki (B.18) 8h cos2 ki cos 4kI Similarly from boundary condition (B.4), we can get a2kl C2 = c2 k tan 2ki (B.19) 8h cos2 klcos4ki 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(lz)]sin 2wt (B.21) cos 4kl cos 4k PAGE 99 APPENDIX C PROGRAM LISTING C.1 Flow Chart C.1.1 Main Routine: T2D Start ) T2DIN Input Parameters Initiate Variables Yes No Advection ADVT1 ADVT2 Yes -No Diffusion DIFN1 ---DIFN2 PROP Perform Propagation CUV Advance Velocities T2DOT No NT=Nmax Output Variables Yes Stop 91 PAGE 100 92 C.1.2 Subroutine: PROP Initiate Q(I,J) CIWSE Compute DZ1(I,J) DZ2(I,J) < TE< c Yes Evaluate DZ(I,J) No Return R(I,J) Compute P(IJ) COYM1 Compute YM(I,J) Update Q(I,J) CIWSE1 Compute DZ1(I,J) DZ2(I,J) C.2 Program listing C.2.1 Description of Parameters KIO: Maximum number of grid points in the x-direction; KJO: Maximum number of grid points in the y-direction; KX1(J), KX2(J), KY1(I), KY2(I): Grid numbers for the boundary's location; CN: Manning coefficient; F: Coriolis parameter; G: Gravitational acceleration; PAGE 101 93 NT: Number of time steps; HMIN: Minimum water depth; DELT: Time increment; EPSLON: Accuracy of the conjugate gradient iteration; NCG: Number of the conjugate gradient iterations; ALFA: Coefficient of implicitization; DC: Horizontal diffusion coefficient; ATA: Amplitude of the focing tide; WOMAGA: Frequency of the forcing tide; AX(I,J), BX(I,J), CX(I,J), DDX(I,J): Coefficients of the matrix equations for the propagation in the x-direction; AY(I,J), BY(I,J), CY(I,J), DDY(I,J): Coefficients of the matrix equations for the propagation in the y-direction; DELX(I): Space increment of (U, V) grid in the x-direction; DELY(J): Space increment of (U, V) grid in the y-direction; ZDELX(I): Space increment of Z grid in the x-direction; ZDELY(J): Space increment of Z grid in the y-direction; Z(I,J): Water surface elevation; ZB(I,J): Bed elevation; H(I,J): Water depth at Z grid; HUV(I,J): Water depth at (U, V) grid; ZP(I,J): Water surface elevation at the previous time; HP(I,J): Water depth of Z grid at the previouse time; HUVP(I,J): Water depth of (U, V) grid at the previous time; U(I,J): Unit-width discharge in the x-direction; V(I,J): Unit-width discharge in the y-direction; PU(I,J): Velocity corresponding to U(I,J); PAGE 102 94 PV(I,J): Velocity corresponding to V(I,J); UD(I,J): Unit-width discharge in the x-direction after the diffusion; VD(I,J): Unit-width discharge in the y-direction after the diffusion; UAD(I,J): Unit-width discharge in the x-direction after the advection; VAD(I,J): Unit-width discharge in the y-direction after the advection; W(I,J): Wind speed; DZ(I,J): Increment of water surface elevation during one time step; DZ1(I,J): Increment of water surface elevation due to the propagation in the x-direction; DZ2(I,J): Increment of water surface elevation due to the propagation in the y-direction; -GAMA(I,J): Weighting coefficient; P(I,J): Direction vector; R(I,J)i Residual vector; Q(I,J): Lagrange coordinates; HX(I,J): Water depth at the middle point of (I+1,J) and (I,J) on Z grid; HY(I,J): Water depth at the middle point of (I,J) and (I,J+1) on Z grid; C.2.2 Program listing C************** MAIN PROGRAM: T2D ************ C** THIS PROGRAM IS USED TO SIMULATE THE TWO-DIMENSIONAL C WATER MOTION IN ESTUARIES, LAKES AND SHALLOW SEAS. C IT IS DEVELOPED USING THE METHOD OF FRACTIONAL STEPS C INCLUDE 'DIM.FOR' C C--INPUT PARAMETERS AND INITIATE VARIABLES C CALL T2DIN C C--ADVACE VARIABLES C DO 1000 NT-1,600 C IF(NT.GT.200) THEN DO 10 I-i.KIO PAGE 103 95 DO 10 J-1,KJO W(IJ)-0. 10 CONTINUE ENDIF C C--DETERMINE THE POSITION OF THE MOVING BOUNDARY C CALL FMOVB C C--EVALUATE WATER DEPTH AT 'UV' GRID C HYDRAULIC RADIUS HX, AND HY C CALL UVDEP C C--PERFORM CALCULATION OF THE ADVETION C CALL ADVT1 I FOR CASE OF ADVECTION C CALL ADVT2 I FOR CASE OF NO ADVECTION C C--CALCULATE THE DIFFUSION C CALL DIFN1 FOR CASE OF DIFFUSION C CALL DIFN2 FOR CASE OF NO DIFFUSION C C--COMPUTE THE PROPAGATION C CALL PROP(NCG,EPSLON) C C--ADVANCE THE UNIT-WIDTH DISCHARGES C CALL CUV C C--OUTPUT THE RESULTS C CALL T2DOT C 1000 CONTINUE C STOP END C*********.******** SUBROUTINE T2DIN *************** C*** THIS SUBROUTINE IS DEVELOPED TO INPUT ALL NECESSARY C PARAMETERS AND DATA FOR RUNNING THE MAIN PROGRAM T2D C SUBROUTINE T2DIN C INCLUDE 'DIM.FOR' C OPEN(UNIT-1,NAME='KX.DAT' ,STATUS-'OLD') OPEN(UNIT-2,NAME-'KY.DAT',STATUS-'OLD') C C--EVALUATE CONSTANT PARAMETERS PAGE 104 96 C ALFA-1.0 G-980. I CM/SECOND**2 CN-0.02 CM**(1/6.) DC=0. I CM**2/SEC. F-0.73E-4 1/SEC. EPSLONI .E-4 HMIN=20. CM DELT-180. I SECOND ATA-50. CM WOMAGA=0.14E-3 1/SEC. NCG=20 KIO-32 KJO-51 C C--INPUT THE POSITION OF THE BOUNDARY KX1(J),KX2(J),KY1(I),KY2(I) C DO 2 J-1,KJO READ(1,3)KX1(J),KX2(J) 2 CONTINUE 3 FORMAT(6X, I2,5X, 2) CLOSE(1) C DO 5 I-1.KIO READ(2,6)KY( I) ,KY2(I) 5 CONTINUE 6 FORMAT(7X, I2,5X, 2) CLOSE(2) C C--INPUT GRID STRUCTURE CALL GRID C C--INITIATE VARIABLES CALL INIT C C--EVALUATE WEIGHTING COEFFICIENTS C DO 10 I-I,KIO DO 10 J-1,KJO IF(I.EQ.1) THEN GAMA(I, J)-0.5 ELSE GAMA(I, J)-DELX(I)/(DELX(I) +DELX(I-1)) ENDIF 10 CONTINUE C RETURN END C C*************SUBROUTINE GRID **"************** C** PERFORM CONSTRUCTION OF THE GRID C SUBROUTINE GRID PAGE 105 97 C INCLUDE 'DIM.FOR' C C--EVALUATE THE SCALE OBTAINED FROM THE MAP OF LAKE OKEECHOBEE. C C (1+7/16)IN=2 MILES C 1 MILE-1.609 KM= 1609M -160900 CM SCALE=2*160900./(1.+7/16.) CM/INCH C C--INPUT THE GRID SIZE IN THE X-DIRECTION C AT (U,V) GRID C DO 10 I-1,KIO-1 IF(I.LE.15) THEN DELX(I)-0.5*SCALE ELSE IF(I.LE.20) THEN DELX(I)-DELX(I-1)+0.1*SCALE ELSE DELX(I)1. *SCALE ENDIF ENDIF 10 CONTINUE DELX(KIO)-DELX(KIO-1) C C AT Z GRID C DO 20 I-2,KIO-1 ZDELX(I)-(DELX(I)+DELX(I-1))/2. 20 CONTINUE ZDELX(1)-DELX(I) ZDELX(KIO)-ZDELX(KIO-1) C C--INPUT GRID SIZE IN THE Y-DIRECTION C AT (U,V) GRID C DO 30 J-1, KJO-1 DELY(J)-0.5*SCALE 30 CONTINUE DELY(KJO)-DELY(KJO-1) C C AT Z GRID C DO 40 J-2,KJO-1 ZDELY(J)-(DELY(J)+DELY(J-1))/2. 40 CONTINUE ZDELY(I)-DELY(l) ZDELY(KJO)-ZDELY(KJO-1) C RETURN END C C***************** SUBROUTINE INIT ************* C** INPUT INITIAL VALUES OF VARIABLES PAGE 106 98 C SUBROUTINE INIT C INCLUDE 'DIM.FOR' DIMENSION WDP(100,100),WDP1(100,100) C OPEN (UNIT=1.NAME-'DEPTH1.DAT' STATUS-'OLD') OPEN (UNIT-2.NAME-'DEPTH2.DAT'.STATUS-'OLD') OPEN (UNIT-3, NAME-'DEPTHS.DAT',STATUS-'OLD') C C--READ WATER DEPTH AT Z GRID POINTS C DO 10 J-2,KJO READ(1,11) WDP(2,J),WDP(3,J),WDP(4,J),WDP(5.J),WDP(6,J). & WDP(7.J).WDP(8.J).WDP(9,J).WDP(10.J),WDP(11,J) READ(2,11) WDP(12,J) ,DP(13,J),WDP(14,J),WDP(15.J),WDP(16.J), & WDP(17,J),WDP(18,J),WDP(19,J),WDP(20,J). DP(21,J) READ(3,12) WDP(22,J),WDP(23, J),DP(24,J),WDP(25,J),WDP(26,J), & WDP(27,J).WDP(28,J),WDP(29,J),VDP(30,J).WDP(31,J), & WDP(32,J) 10 CONTINUE 11 FORMAT (4X,F4.1,2X,F4.1,2XF4.12X,F.12X,F4.1,2X,F4.1,2X,F4.1, & 2X.F4.1.2X,F4.1,2X,F4.1) 12 FORMAT (4XF4.1,2XF4.1.2XF4.1.2X.F4.1.2X,F4.1.2XF4.1,2X.F4.1, & 2X,F4.1,2X,F4.1,2X,F4.1,2X.F4.1) C CLOSE(I) CLOSE(2) CLOSE(3) C C--SMOOTH WATER DEPTH C DO 5 I-2,KIO DO 5 J-2.KJO WDP1(IJ)-(WDP(I-1,J)+WDP(I,J)+VDP(I+1.J)+WDP(I-1.J+1)+WDP(I,J+1) & +WDP(I+1,J+1)+WDP(I-1,J-1)+WDP(I,J-1)+WDP(I+1,J-1))/9. 5 CONTINUE C C--TRANSFORM THE UNIT C DO 15 J-2,KJO DO 15 I-2,KIO H(I,J)-WDP1(I.J)*30.48 CM 15 CONTINUE C DO 16 J-2,KJO H(1,J)-H(2,J) 16 CONTINUE C DO 17 I-I,KIO H(I.1)-H(I,2) 17 CONTINUE C DO 18 I-1.KIO PAGE 107 99 DO 18 J-1,KJO IF(H(I,J).LT.HMIN) THEN H(IJ)-HMIN ENDIF 18 CONTINUE C C--EVALUATE THE WATER SURFACE ELEVATION AT 'Z' GRID C DO 20 I-1,KIO DO 20 J-1,KJO Z(I,J)-2000. CM 20 CONTINUE C C--EVALUATE THE BED ELEVATION ZB(I,J) C DO 30 I-1,KIO DO 30 J-1,KJO ZB(I,J)-Z(IJ)-H(I,J) I CM 30 CONTINUE C C--EVALUATE WATER DEPTH AT (U,V) GRID C DO 40 J-1,KJO-1 DO 40 I-1,KIO-1 HUV(IJ)-(H(I.J)+H(I, J+1)+H(I+1, J)+H(I+1, J+1))/4. 40 CONTINUE DO 41 J-1,KJO HUV(KIO,J)-HUV(KIO-1,J) 41 CONTINUE DO 42 I-1,KIO HUV(IKJO)-HUV(I,KJO-1) 42 CONTINUE C C--INITIATE UNIT-WIDTH DISCHARGES AND VELOCITIES C DO 50 I-1,KIO DO 50 J-1,KJO U(I,J)-0. I CM**2/SEC. V(I,J)-0. UD(I,J)-0. VD(IJ)-0. UAD(I,J)-0. VAD(I,J)-0. PU(IJ)-0. CM/SEC. PV(I.J)-0. 50 CONTINUE C C--INITIATE THE INCREMENT OF WATER SURFACE ELEVATION C DO 60 I=1,KIO DO 60 J-1,KJO DZ(I,J)-0. I CM 60 CONTINUE C PAGE 108 100 C--INPUT THE WIND SPEED C DO 70 J-2,KJO IF(J.LT.25) THEN N1-KX1(J-1)+1 N2KX2 (J-1) ELSE Nl-KXI(J)+1 N2-KX2(J) ENDIF DO 70 I-N1,N2-1 W(I,J)-2000. CM/SEC. 70 CONTINUE C RETURN END C************* SUBROUTINE ADVTI *******l***** ******* C**THIS SUBROUTINE IS USED FOR SOLVING THE ADVECTION. USING THE METHOD C OF FRACTIONAL STEPS, THE TWO-DIMENSIONAL ADVECTION IS DIVIDED AS TWO C ONE-DIMENSIONAL ADVCTIONS WHICH ARE X-ADVECTION AND Y-ADVECTION. THE C IMPLICIT CHARACTERISTIC NUMERICAL SCHEME IS USED TO SOLVE EACH C ADVECTION. IT IS NOTED THAT THE FOLLOWING PROGRAM IS DEVELOPED C ACCORDING TO A RECTANGULAR BASIN. UX(I,J),VX(I,J): VELOCITY C COMPONENTS IN THE XAND Y-DIRECTIONS AFTER X-ADVECTION. UY(I,J), C VY(I,J): VELOCITIES COMPONENTS IN THE XAND Y-DIRECTIONS AFTER C Y-ADVECTION. C SUBROUTINE ADVT C INCLUDE 'DIM.FOR' DIMENSION A(100),B(100),C(100),D(100),X(100) DIMENSION UX(100,100),VX(100,100) UY(100,100),VY(100,100) C C--X-ADVECTION C COMPUTATION OF UX(100,100) C KI2-KIO-1 KJ2-KJO-1 DO 10 J-2,KJ2 DO 20 1=2,KI2 C A(I)-1. B(I)-PU(I,J)*DELT/(DELX(I)+DELX(I-1)) C(I)--PU(I,J)* DELT/(DELX(I)+DELX(I-1)) D(I)-PU(IJ) C IF(I.EQ.2) THEN A(I) A(I)+C(I) C(I)O0. ENDIF C IF (I.EQ.KI2) THEN PAGE 109 101 B(I)=0. ENDIF C 20 CONTINUE C KM=KI2-1 CALL CUGA(X,A,B,C,D,KM) DO 30 I=2,KI2 UX(I, J)-X(I) 30 CONTINUE 10 CONTINUE C DO 40 J-2,KJ2 UX(1,J)-UX(2,J) 40 CONTINUE C C COMPUTATION OF VX(100,100) C DO 50 J-2,KJ2 DO 60 I-2,KI2 A(I)-1. B(I)-PU(I, J)*DELT/(DELX(I)+DELX(I-1)) C(I)--PU(I,J)*DELT/(DELX(I)+DELX(I-1)) D(I)-PV(I.J) IF(I.EQ.2) THEN A(I)-A(I)+C(I) C(I)-0. ENDIF IF (I.EQ.KI2) THEN A(I)=A(I)+B(I) B(I)-0. ENDIF 60 CONTINUE C KM-KI2-1 CALL CUGA(X,A.B,C,D,KM) DO 70 I12,KI2 VX(I, J)-X(I) 70 CONTINUE 50 CONTINUE C DO 80 J-2,KJ2 VX(1,J)-VX(2,J) 80 CONTINUE C C--Y-ADVECTION C COMPUTATION OF UY(100,100) C DO 110 I-2,KI2 DO 120 J=2,KJ2 A(J)-1. B(J)-VX(I,J)*DELT/(DELY(J)+DELY(J-1)) C(J)--VX(I,J)*DELT/(DELY(J)+DELY(J-1)) D(J)-UX(I, J) PAGE 110 102 IF(J.EQ.2) THEN A(J)-A(J)+C(J) C(J)-O. ENDIF IF(J.EQ.KJ2) THEN A(J)-A(J)+B(J) B(J)O0. ENDIF 120 CONTINUE C KM-KJ2-1 CALL CUGA(X,A,BC,D,KM) DO 130 J-2,KJ2 UY(I, J) X(J) 130 CONTINUE 110 CONTINUE C C COMPUTAITON OF VY(100,100) C DO 150 I-2,KI2 DO 160 J-2,KJ2 A(J)-1. B(J)-VX(IJ)*DELT/(DELY(J)+DELY(J-1)) C(J)--VX(I, J)*DELT/ (DELY(J) +DELY(J-1)) D(J)-VX(I.J) IF(J.EQ.2) THEN C(J)-0. ENDIF IF(J.EQ.KJ2) THEN B(J)-0. ENDIF 160 CONTINUE C KM-KJ2-1 CALL CUGA(XA,B,C,D,KM) DO 170 J-2,KJ2 VY(I, J)-X(J) 170 CONTINUE 150 CONTINUE C C--IMPOSE BOUNDARY CONDITIONS C DO 200 J=2,KJ2 UY(1, J)-UY(2,J) VY(1,J)-VY(2.J) UY(KI, J)0. VY(KIJ)-VY(KI2,J) 200 CONTINUE C DO 220 I-1,KIO UY(I, 1)-UY(I, 2) UY (I, KJO) UY(I, KJ2) VYC(I.1)-0. VY(I.KJO)-0. PAGE 111 103 220 CONTINUE C C--TRANSFER UY,VY INTO UAD,VAD C DO 240 I=1,KIO DO 240 J-1,KJO UAD(I,J)-UY(I,J)*HUV(I,J) VAD(I,J)-VY(I,J)*HUV(I,J) 240 CONTINUE C RETURN END C************* SUBROUTINE ADVT2 ************* C**THIS SUBROUTINE IS USED FOR THE CASE OF NO ADVECTION. C SUBROUTINE ADVT2 INCLUDE 'DIM.FOR' C DO 20 I-1,KIO DO 20 J-KY1(I),KY2(I) UAD(I,J)-U(I,J) VAD(I,J)-V(I,J) 20 CONTINUE C RETURN END C******.******** SUBROUTINE DIFN1 **************** C**PERFORM THE CALCULATION OF THE DIFFUSION C THE DIFFUSION IS SOLVED USING THE ADI METHOD. C SUBROUTINE DIFN1 C DIMENSION A(150),B(150),C(1O0), D(150),X(150) DIMENSION AA(150),BB(150),CC(150),DD(150) INCLUDE 'DIM.FOR' C C--COMPUTE UD(I,J) C KJ2=KJO-1 DO 50 J-2.KJ2 N1-KXI(J)+1 N2-KX2(J)-1 DO 60 I-N1,N2 IF (J.EQ.2) THEN TERMI-2*DC*(UAD(I,J+) -UAD(I,J)) k /DELY(J)/(DELY(J)+DELY(J-1)) ELSE IF (J.EQ.KJ2) THEN TERM1-2*DC*(UAD(IJ)-UAD(I,J-1)) 6 /DELY(J-1)/(DELY(J-1)+DELY(J)) PAGE 112 104 ELSE TM-(UAD(I.J+1)-UAD(I,J))/DELY(J) & -(UAD(I,J)-UAD(IJ-1))/DELY(J-1) TERM1=2*DC*TM/(DELY(J-1)+DELY(J)) ENDIF ENDIF C IF(I.EQ.N1) THEN A(I)-I.+DC*DELT/DELX(I)/DELX(I) B(I)--DC*DELT/DELX(I)/DELX(I) C(I)-0. D(I)-UAD(I, J)+DELT* (F*VAD(I,J)+TERM1) ELSE IF(I.EQ.N2) THEN A(I)-I.+2*DC*DELT/DELX(I)/DELX(I-1) B(I)=0. C(I)--2*DC*DELT/DELX(I)/DELX(I-1) D(I)-UAD(I,J)+DELT*(F*VAD(I,J)+TERM1) ELSE A(I)-1.+2*DC*DELT/DELX(I)/DELX(I-1) B(I)--2*DC*DELT/DELX(I)/(DELX(I)+DELX(I-1)) C(I)--2*DC*DELT/DELX(I-1)/(DELX(I)+DELX(I-1)) D(I)-UAD(I,J)+DELT*(F*VAD(I,J)+TERM1) ENDIF ENDIF 60 CONTINUE C N-N2-N1+1 DO 65 I-2,N+1 AA(I)-A(N1-2+I) BB(I)-B(N1-2+I) CC(I)-C(N1-2+I) DD(I)-D(N1-2+I) 65 CONTINUE C CALL CUGA(X.AA,BB.CC,DD,N) C DO 70 I-N1,N2 UD(I,J)-X(I-N1+2) 70 CONTINUE C 50 CONTINUE C C--COMPUTE VD(I,J) C KI2-KIO-1 DO 110 I-2,KI2 N1-KYI(I)+1 N2-KY2(I)-1 DO 120 J-N1,N2 IF (I.EQ.2) THEN TM-(VAD(I+1,J)-VAD(I,J))/DELX(I) TERMI-DC*TM/DELX(I) ELSE PAGE 113 105 IF(I.EQ.KI2) THEN TM=-(VAD(IJ)-VAD(I-1,J))/DELX(I-1) TERM1=2*DC*TM/(DELX(I)+DELX(I-1)) ELSE TM-(VAD(I+1,J)-VAD(I,J))/DELX(I) & -(VAD(I,J)-VAD(I-1,J))/DELX(I-1) TERMI=2SDC*TM/(DELX(I)+DELX(I-1)) ENDIF ENDIF C IF (J.EQ.N1) THEN A(J)-1.+2*DC*DELT/DELY(J)/DELY(J-1) B (J)--2*DC*DELT/DELY(J) /(DELY(J) +DELY(J-1)) C(J)-0. D(J)-VAD(I,J)+DELT*(-F*UAD(IJ)+TERMI) ELSE IF (J.EQ. N2) THEN A(J)1.+2*DC*DELT/DELY(J-1)/DELY(J-1) B(J)-0. C(J)--2*DC*DELT/DELY(J-1)/(DELY(J)+DELY(J-1)) D(J)-VAD(I.J)+DELT*(-F*UAD(I,J)+TERMI) ELSE A(J)-I.+2*DC*DELT/DELY(J)/DELY(J-1) B(J)--2.*DC*DELT/DELY(J)/(DELY(J)+DELY(J-1)) C(J)--2.*DC*DELT/DELY(J-1)/(DELY(J)+DELY(J-1)) D(J)-VAD(I,J)+DELT*(-F*UAD(I.J)+TERM1) ENDIF ENDIF 120 CONTINUE C N-N2-NI+1 DO 130 J-2,N+1 AA(J)-A(Ni-2+J) BB(J)-B(N1-2+J) CC(J)-C(N1-2+J) DD(J)-D(NI-2+J) 130 CONTINUE C CALL CUGA(XAABB,CC,DD,N) C DO 140 J-NI,N2 VD(I,J)-X(J-N1+2) 140 CONTINUE 110 CONTINUE C C IMPOSE THE CLOSE BOUNDARY CONDITIONS C DO 200 J=1,KJO UD(1,J)-0. VD(1,J)-0. UD(KIO,J)-0. VD(KIO,J)-0. 200 CONTINUE DO 210 I-2,KIO-1 PAGE 114 106 IF(I.LT.25) THEN Ni-KY1(I-1) N2-KY2(I-1) ELSE N1-KY1(I+1) N2-KY2(I+l) ENDIF DO 220 J=1,N1 UD(I,J)-0. VD(I,J)-0. 220 CONTINUE DO 230 J-N2,KJO UD(I,J)-0. VD(IJ)-0. 230 CONTINUE 210 CONTINUE C RETURN END C************** SUBROUTINE DIFN2 ********** C**THIS SUBROUTINE IS USED FOR THE CASE OF NO DIFFUSION. C SUBROUTINE DIFN2 INCLUDE 'DIM.FOR' C DO 10 J-2,KJO-1 DO 20 I-1,KIO-1 UD(I, J)-UAD(I, J)+DELT*F*VAD(I, J) VD(I, J)-VAD(I, J)-DELT*F*UD(I, J) 20 CONTINUE 10 CONTINUE C C--IMPOSE THE BOUNDARY CONDITIONS C DO 30 J-1,KJO-1 UD(KIO.J)-0. VD(KIO,J)-VD(KIO-1, J) 30 CONTINUE DO 40 I-1,KIO-1 UD(I, 1)-UD(I,2) UD(I,KJO)-UD(I,KJO-1) VD(I,1)-0. VD(I,KJO)-0. 40 CONTINUE C RETURN END PAGE 115 107 C************** SUBROUTINE PROP ************* C** PERFORM THE COMPUTATION OF THE PROPAGATION C THE PROPAGATION IS SOLVED USING THE CONJUGATE GRADIENT C ITERATION METHOD. C SUBROUTINE PROP(NCG,EPSLON) C INCLUDE 'DIM.FOR' DIMENSION Q(150,150),P(150,150),R(150,150) C C C--PERFORM THE CONJUGATE GRADIENT ITERATION C DO 50 N=-lNCG C IF (NT.EQ.1) THEN IF(N.EQ.1) THEN DO 60 I-2,KIO IF(I.LT.25) THEN NI-KYI(I-I)+I N2-KY2(I-1) ELSE N1-KYI(I)+1 N2-KY2(I) ENDIF DO 60 J=N1,N2 Q(IJ)-0. 60 CONTINUE ENDIF ENDIF C C COMPUTE DZI(IJ) AND DZ2(I,J) C IF(N.EQ.1) THEN CALL CIVSE(Q) ELSE CALL CIVSEl(Q) ENDIF C C CALCULATE THE SUM OF ERRORS: TE C TE-0. DO 80 I=2,KIO IF(I.LT.25) THEN Nl-KY1(I-1)+1 N2-KY2(I-1) ELSE NI-KY1(I)+l N2-KY2(I) ENDIF DO 80 J-N1,N2 TE-TE+(DZ(I, J)-DZ2 (I,J))**2 80 CONTINUE C PAGE 116 108 IF (TE.LT.EPSLON) GOTO 200 C C EVALUATE THE COEFFICIENT A2 C IF(N.NE.1) THEN A2-TE/TE1 ENDIF TE1-TE C C COMPUTE RESIDUAL VECTOR R(I,J) C DO 90 I-2,KIO IF(I.LT.25) THEN NI=KY1(I-1)+1 N2-KY2(I-1) ELSE NI-KY1(I)+1 N2-KY2 (I) ENDIF DO 90 J-N1,N2 R(I,J)-DZ1(I, J)-DZ2(IJ) 90 CONTINUE C C COMPUTE THE DIRECTION VECTOR P(I,J) C DO 100 I-2,KIO IF(I.LT.25) THEN Nl-KY1(I-1)+1 N2-KY2(I-1) ELSE N1-KY1(I)+1 N2-KY2(I) ENDIF DO 100 J-N1,N2 IF (N.EQ.1) THEN P(I.J)-R(I. J) ELSE P(I,J)R(I, J)+A2*P(I,J) ENDIF 100 CONTINUE C C COMPUTE YM1(I,J) AND YM2(I,J) CALL COYM1(P) C C EVALUATE THE COEFFICIENT Al C PDZO0. PYM=0. DO 120 I-2,KIO IF(I.LT.25) THEN NI-KYI(I-1)+1 N2-KY2(I-1) ELSE PAGE 117 109 Nl-KYI(I)+1 N2-KY2(I) ENDIF DO 120 J-N1,N2 PDZ-PDZ+(DZ1l(I,J)-DZ2(I,J))*P(,J) PYM-PYM+(YM1(I,J)+YM2(I,J))*P(I,J) 120 CONTINUE Al-PDZ/PYM C C UPDATE THE LAGRANGE COORDINATE Q(I,J) C DO 140 I-2,KIO IF(I.LT.25) THEN NI-KY1(I-1)+1 N2-KY2(I-1) ELSE N1-KY1(I)+1 N2-KY2(I) ENDIF DO 140 J-N1,N2 Q(IJ)-q(IJ)+Ai*P(I,J) 140 CONTINUE C 50 CONTINUE C C--EVALUATE DZ(I.J) C 200 DO 160 I-2,KIO IF(I.LT.25) THEN NI-KYI(I-1)+1 N2-KY2(I-1) ELSE NI-KY1(I)+1 N2-KY2(I) ENDIF DO 160 J-N1,N2 DZ(I,J)-(DZi(I,J)+DZ2(IJ))/2. 160 CONTINUE C C--IMPOSE THE BOUNDARY CONDITIONS C DO 170 I-2,KIO IF(I.LT.25) THEN DZ(I,KY1(I-1))=DZ(I,KY1(I-1)+1) DZ(IKY2(I-1)+1)=DZ(I.KY2(I-1)) ELSE DZ(I,KY1(I))-DZ(I.KYI(I)+1) DZ(I,KY2(I)+1)-DZ(IKY2(I)) ENDIF 170 CONTINUE DO 180 J-KYI(1),KY2(1) DZ(1,J)-DZ(2,J) !CLOSE BOUNDARY AT I1180 CONTINUE DO 190 J=KYI(KIp),KY2(KIO) PAGE 118 110 DZ(KIO+1,J)-DZ(KIO,J) !CLOSE BOUNDARY AT I-KIO 190 CONTINUE DO 195 J-2.KJO IF(J.LT.25) THEN DZ(I,KX1(J-1))-DZ(I.KX1(J-1)+1) DZ(I,KX2(J-1)+1)-DZ(IKX2(J-1)) ELSE DZ(I,KX1(J))-DZ(I.KX1(J)+1) DZ(I,KX2(J)+1)-DZ(IKX2(J)) ENDIF 195 CONTINUE C RETURN END C**a****"********* SUBROUTINE CIWSE **************** C** PERFORM THE COMPUTATION OF DZI(IJ) AND DZ2(I,J),GIVEN Q(I,J) C SUBROUTINE CIWSE(Q) C DIMENSION Q(150,150) DIMENSION A(150).B(150),C(150),D(150),X(150) INCLUDE 'DIM.FOR' C C--COMPUTATION OF DZI(IJ) C DO 20 J-2,KJO IF (J.LT.25) THEN N1-KXI(J-1)+1 N2-KX2(J-1) ELSE Ni-KX1(J)+1 N2-KX2(J) ENDIF DO 30 I-N1,N2 DX-(ZDELX(I)+ZDELX(I-1))/2. IF(I.EQ.N2) THEN TERM1-1./2./G/DELT/DELT TM-HX(I-1.J)+(Z(I.J)-Z(I-1,J))*GAMA(I-1,J) TERM2-ALFA**2*TM/DX/ZDELX(I-1) TM-UD(I,J)/HX(I,J)-GAMA(I-1,J)*UD(I-1.J)/HX(I-1,J) TERM3-ALFA*TM/G/DELT/DX AX(I.J)-TERM1+TERM2+TERM3 C BX(I.J)-0. C TM--HX(I-1,J)+(1-GAMA(I-1,J))*(Z(I.J)-Z(I-1,J)) TERM1-ALFA**2*TM/DX/ZDELX(I-1) TM-(1-GAMA(I-1,J))*UD(I-1,J) TERM2--TM*ALFA/G/DELT/HX(I-1,J)/DX CX(I.J)-TERM1+TERM2 C TM"(1-ALFA)*(U(I.J)-U(I-1,J))+ALFA*(UD(I.J)-UD(I-1.J)) PAGE 119 111 TERM1--TM/G/DELX (I-1)/DELT TM--HX(I-1,J)*(Z(I,J)-Z(I-1,J)) TERM2-ALFA*TM/DX/ZDELX (I-) FC=116.*CN**2/HR(I,J)**(1./3.) TMI-SQRT(U(I,J)**2+V(IJ)**2)*U(I.J)*FC/8./HX(I,J)**2 FC-116.*CN**2/HR(I-1,J)**(1./3.) TM2=SQRT(U(I-1,J)**2+V(I-1.J)**2)*U(I-1.J) & *FC/8./HX(I-1,J)**2 TM3-3.*0.000001*ABS(W(I,J))*W(I,J) TM4-3.*0.000001*ABS(W(I-1,J))*W(I-1,J) TERM3=ALFA*(TM1-TM2-TM3+TM4)/G/DX DDX(I,J) TERMi+TERM2+TERM3-Q(I,J) C ELSE C TERM1-1./2./G/DELT/DELT TM1-HX(I,J)-(Z(I+1,J)-Z(I,J))*(1-GAMA(I,J)) TM2-HX(I-1,J)+(Z(I,J)-Z(I-1,J))*GAMA(I-1,J) TM-TM1/ZDELX(I)+TM2/ZDELX(I-1) TERM2-ALFA**2*TM/DX TM1-UD(I.J)*(1-GAMA(I,J))/HX(I,J) TM2-UD(I-1, J) *GAMA(I-1,J)/HX(I-1, J) TERM3-ALFA*(TM1-TM2)/G/DELT/DX AX(I,J)-TERMI+TERM2+TERM3 C TERM1i-ALFA**2*HX(I.J)/DX/ZDELX(I) TERM2--ALFA**2*(Z(I+1,J)-Z(I,J))*GAMA(I,J)/DX/ZDELX(I) TERM3-ALFA*UD(I.J)*GAMA(I,J)/G/DELT/HX(I J)/DX BX(I,J)-TERM1+TERM2+TERM3 C TERMI--ALFA**2*HX(I-1,J)/DX/ZDELX(I-1) TERM2-ALFA**2*(Z(I,J)-Z (I-1,J)) & *(1-GAMA(I-1,J))/DX/ZDELX(I-1) TERM3=-ALFA*UD(I-1,J) *(1-GAMA(I-I J)) & /G/DELT/HX(I-1,J)/DX CX(I,J)-TERMi+TERM2+TERM3 C TERM1(1-ALFA)* (U(I, J)-U(I-1, J))/G/DELT/DELX(I-1) TERM2--ALFA* (UD(I, J)-UD(I-1, J))/G/DELT/DELX(I-1) TMi-HX(I.J)*(Z(I+1,J)-Z(I.J))/ZDELX(I) TM2-HX(I-1,J)*(Z(I,J)-Z(I-1,J))/ZDELX(I-1) TERM3=ALFA*(TM1-TM2)/DX FC-116.*CN**2/HR(I,J)**(1./3.) TM1-SQRT(U(I,J)**2+V(I,J)**2)*U(I,J)*FC/8./HX(I,J)**2 FC=116.*CN**2/HR(I-, J)**(1./3.) TM2-SQRT(U(I-1,J)**2+V(I-1, J)**2)*U(I-1,J) & *FC/8./HX(I-1,J)**2 TM3=3.*0.000001*ABS(W(I.J))*W(I,J) TM43. *0.000001*ABS(WVI-1,J))*W(I-1,J) TERM4-ALFA* (TM1-TM2-TM3+TM4)/G/ZDELX(I-1) DDX(I,J) TERMI+TERM2+TERM3+TERM4-Q(I,J) C IF(I.EQ.Ni) THEN AX(I,J)=AX(I.J)+CX(I,J) ICLOSE BOUNDARY AT 1-2 PAGE 120 112 CX(I,J)-0. ENDIF C ENDIF 30 CONTINUE C N=N2-NI+1 DO 35 I-2,N+1 A(I)-AX(N1-2+IJ) B(I)-BX(NI-2+I,J) C(I)-CX(N1-2+I,J) D(I)-DDX(N1-2+I,J) 35 CONTINUE C CALL CUGA(X,A,B,C,D,N) C DO 40 I-Ni,N2 DZ1(I.J)-X(I-NI+2) 40 CONTINUE C 20 CONTINUE C C--COMPUTATION OF DZ2(I,J) C DO 50 I-2,KIO IF(I.LT.25) THEN N1-KY1(I-1)+ N2-KY2(I-1) ELSE N1-KYI(I)+1 N2-KY2(I) ENDIF DO 60 J-N1,N2 DY-(ZDELY(J)+ZDELY(J-1))/2. IF (J.EQ.N1) THEN TERM-=1./2./G/DELT/DELT TM--HY(I,J)+(1-GAMA(I,J))*(Z(I,J+1)-Z(I,J)) TERM2--ALFA**2*TM/DY/ZDELY(J) TM-(1-GAMA(I,J))*VD(I,J) TERM3-ALFA*TM/G/DELT/DY/HY(IJ) AY(I,J)-TERM1+TERM2+TERM3 C TM-HY(I,J)+GAMA(I.J)*(Z(I,J+1)-Z(I,J)) TERM1--ALFA**2*TM/DY/ZDELY (J) TM-GAMA(I,J)*VD(I,J) TERM2-ALFA*TM/G/DELT/DY/HY(IJ) BY(I,J) TERM1+TERM2 C CY(I,J)-0. C TERM1-( (1-ALFA)*V(I,J)+ALFA*VD(IJ))/G/DELT/DELY(J-1) TERM2-ALFA*HY(I,J)*(Z(I,J+1)-Z(I,J))/ZDELY(J)/DY FC=116.*CN**2/HR(I,J)**(1./3.) TM1-FC*SQRT(U(I, J)**2+V(I J)**2)*V(I.J)/8./HY(I,J)**2 PAGE 121 113 TERM3-ALFA*TM1/G/DY DDY(I,J)-TERM1+TERM2+TERM3+Q(IJ) C ELSE IF(J.EQ.N2) THEN C TERMI11./2./G/DELT/DELT TM-HY(I.J-1)+(Z(I,J)-Z(I,J-1))*GAMA(I,J-1) TERM2-ALFA**2*TM/ZDELY(J-1)/DY TM-GAMA(I,J-1)*VD(I,J-1) TERM3--ALFA*TM/G/DELT/DY/HY(I,J-1) AY(IJ)-TERM1+TERM2+TERM3 C BY(I,J)-0. C TM--HY(I,J-1)+(Z(I,J)-Z(I,J-1)') *(1-GAMA(I,J-i)) TERMI-ALFA**2*TM/ZDELY(J-1)/DY TM-(1-GAMA(I.J-1))*VD(IJ-1) TERM2--ALFA*TM/G/DELT/DY/HY(I,J-1) CY(I,J)-TERMI+TERM2 C TMi1(1-ALFA)*V(I,J-l)+ALFA*VD(I,J-l) TERM1-TM1/G/DELT/DELY(J-1) TERM2--ALFA*HY(I,J-i)*(Z(I,J)-Z(I,J-1))/DY/ZDELY(J-1) FC-116.*CN**2/HR(I,J-1)**(1./3:) TM1IFC*SQRT(U(IJ-l)**2+V(I.J-l)**2)*V(I,J-l) /8./HY(IJ-i)**2 TERM3--ALFA*TM1/G/DY DDY(I,J)-TERMI+TERM2+TERM3+Q(I,J) C ELSE C TERM1-1./2./G/DELT/DELT TM1--HY(I.J)+(Z(I.J+1)-Z(I,J))*(1-GAMA(I.J)) TM2--HY(I,J-1)-(Z(I,J)-Z(I.J-1))*GAMA(I,J-1) TM-TM1/ZDELY(J)+TM2/ZDELY(J-1) TERM2--ALFA**2*TM/DY TM1-VD(I,J)*(1-GAMA(I.J))/HY(I.J) TM2-VD(I.J-1)*GAMA(I,J-1)/HY(IJ-1) TERM3-ALFA*(TM1-TM2)/G/DELT/DY AY(I.J)-TERM1+TERM2+TERM3 C TERMI--ALFA**2*(HY(I,J)+GAMA(I,J) & *(Z(I,J+1)-Z(I,J)))/DY/ZDELY(J) TERM2-ALFA*GAMA(I,J)*VD(I,J)/G/DELT/DY/HY(I,J) BY(I J)-TERMI+TERM2 C TERM1--ALFA**2*HY(I.J-1)/DY/ZDELY(J-1) TERM2-ALFA**2*(Z(I,J)-Z(I.J-1)) & *(1-GAMA(I,J-1))/DY/ZDELY(J-1) TERM3--ALFA*VD(I,J-1)*(1-GAMA(I,J-I)) & /G/DELT/DY/HY(I,J-1) CY(I,J) TER+TERM+E2+TERM3 C PAGE 122 114 TERM1=-(1-ALFA)*(V(I,J)-V(IJ-1))/G/DELT/DELY(J-1) TERM2=-ALFA*(VD(I.J)-VD(I J-1))/G/DELT/DELY(J-1) TERM3-ALFA*HY(I,J)*(Z(I,J+1)-Z(I,J))/DY/ZDELY(J) TERM4--ALFA*HY(I,J-1)*(Z(I,J)-Z(I,J-1))/DY/ZDELY(J-1) C -------COMPUTATION OF FRICTION----C FC=116.*CN**2/HR(I.J)**(1./3.) TM1-FC*SQRT(U(I,J)**2+V(I,J)**2)*V(I,J)/8./HY(IJ) *2 FC'116.*CN**2/HR(I,J-1)**(1./3.) TM2-FC*SQRT(U(I,J-1)**2+V(I.J-1)**2) A *V(I,J-l)/8./HY(I,J-1)**2 TERM5-ALFA/G/DY*(TMi-TM2) DDY(IJ)-TERMi+TERM2+TERM3+TERM4+TERM5+Q(I,J) C ENDIF ENDIF 60 CONTINUE C N-N2-Ni+2 DO 65 J-2,N+I A(J)=AY(I,NI-2+J) B(J)-BY(I,Ni-2+J) C(J)-CY(I.Ni-2+J) D(J)-DDY(I,N1-2+J) 65 CONTINUE C CALL CUGA(X,A.B,C.D.N) C DO 70 J-NI,N2 DZ2(I,J)-X(J-N1+2) 70 CONTINUE C 50 CONTINUE C RETURN END C************ **** SUBROUTINE CIWSE1 *********** C** PERFORM THE COMPUTATION OF DZ1(I,J) AND DZ2(I,J) C SUBROUTINE CIWSE1(Q) C DIMENSION Q(150,150) DIMENSION A(150),B(150),C(150),D(150),X(150) INCLUDE 'DIM.FOR' C C--COMPUTATION OF DZI(I,J) C DO 20 J-2,KJO IF(J.LT.25) THEN NI-KX1(J-1)+1 N2-KX2(J-l) PAGE 123 115 ELSE NI1KX1(J)+1 N2-KX2 (J) ENDIF N-N2-NI+1 DO 30 I-2,N+l A(I)-AX(Ni-2+IJ) B(I)-BX(N1-2+I.J) C(I)-CX(Ni-2+I, J) D(I)-DDX(N1-2+I,J) -Q (N1-2+I J) 30 CONTINUE C CALL CUGA(XA,B,C,D,N) C DO 40 I-N1,N2 DZ1(I,J)-X(I-Ni+2) 40 CONTINUE C 20 CONTINUE C C--COMPUTATION OF DZ2(I,J) C DO 50 I-2,KIO IF(I.LT.25) THEN NIlKYl (I-1)+1 N2-KY2 (I-1) ELSE Ni-KY1(I)+1 N2-KY2(I) ENDIF N-N2-N1+1 DO 60 J-2.N+1 A(J)-AY(I.N1-2+J) B(J)-BY(I,Ni-2+J) C(J)-CY(I,NI-2+J) D(J)-DDY(I,Ni-2+J)+Q(I.N1-2+J) 60 CONTINUE C CALL CUGA(X,A,B,C.D,N) C DO 70 J-N1.N2 DZ2(I,J)-X(J-N1+2) 70 CONTINUE C 50 CONTINUE C RETURN END C**************** SUBROUTINE COYM1 ************* C** PERFORM THE COMPUTATION OF YMl(I,J) AND YM2(I,J), C GIVEN P(I,J) C PAGE 124 116 SUBROUTINE COYMi(P) C DIMENSION P(150,150) DIMENSION A(150),B(150),C(150).D(150),X(150) INCLUDE 'DIM.FOR' C C--COMPUTATION OF YMI(I,J) C DO 20 J=2,KJO IF(J.LT.25) THEN NI-KX1(J-1)+1 N2-KX2(J-1) ELSE N1-KX1(J)+1 N2-KX2(J) ENDIF N-N2-NI+1 DO 30 I12,N+l A(I)=AX(N1-2+I,J) B(I)-BX(NI-2+I,J) C(I)-CX(N1-2+I,J) D(I)-P(Ni-2+I,J) 30 CONTINUE C CALL CUGA(X,A.B,C,D,N) C DO 40 I=N1.N2 YM1(I, J)-X(I-N1+2) 40 CONTINUE C 20 CONTINUE C C--COPUTATION OF YM2(I.J) C DO 50 I-2,KIO IF(I.LT.25) THEN NI-KY1(I-1)+l N2-KY2(I-1) ELSE N1-KY1(I)+1 N2-KY2(I) ENDIF N-N2-N1+1 DO 60 J-2,N+1 A(J)-AY(I.N1-2+J) B(J)-BY(I, N-2+J) C(J)-CY(IN1-2+J) D(J)-P(I, N-2+J) 60 CONTINUE C CALL CUGA(X,AB,C.D,N) C DO 70 J-N1,N2 YM2(I. J)-X(J-N1+2) PAGE 125 117 70 CONTINUE C 50 CONTINUE C RETURN END C****************** SUBROUTINE CUV ************** C** ADVANCE THE VELOCITIES C SUBROUTINE CUV DIMENSION ZP(150,150) ,HP(150,150),HUVP(150.150) INCLUDE 'DIM.FOR' C C--RESTORE THE OLD WATER SURFACE ELEVATION AND WATER DEPTH C TO ZP(I,J), HP(I,J) AND HUVP(I,J) C DO 20 I-1,KIO DO 20 J-1,KJO ZP(I.J)-Z(I.J) HP(I,J)-H(I,J) HUVP(I,J)-HUV(I,J) 20 CONTINUE C C--ADVANCE WATER SURFACE ELEVATION AND WATER DEPTH C CALL DEPT C C--ADVANCE UNIT-WIDTH DISCHARGES U(I,J) AND V(I,J) C DO 60 I-2,KIO-1 IF(I.LT.25) THEN Nl-KY1(I-1)+1 N2-KY2(I-1)-1 ELSE NI-KYI(I+I)+l N2-KY2(I+1)-1 ENDIF DO 60 J-N1,N2 UV-(U(I,J)**2+V(IJ)**2)**0.5 C C COMPUTATION OF U(I,J) C TMI-(Z(I+1,J)+Z(I+1.J+1))/2.-(Z(I.J)+Z(I,J+1))/2. TM2-(ZP(I+1,J)+ZP(I+IJ+1))/2.-(ZP(I.J)+ZP(I,J+l))/2. TERM1--DELT*G*(ALFA*HUV(I,J)*TM1 k +(I-ALFA)*HUVP(I,J)*TM2)/ZDELX(I) TM1-(Z(I+1,J)+Z(I+I,J+1)+Z(I,J)++Z(IJ+1))/4. TM2-(ZP(I+1,J)+ZP(I+1,J+1)+ZP(I,J)+ZP(I,J+1))/4. TERM2=UD(I,J)*(TMI-TM2)/HUVP(I,J) TM=1UV*U(I,J) FC-116.*CN**2/HR(I,J)**(1/3.) TRMI-FC*TM1/8./HUVP(I.J)**2 PAGE 126 118 TRM2=3.*0.000001*ABS(W(I,J))*W(I.J) TERM3=-DELT*(TRM1-TRM2) U(I,J)-UD(I,J) +TERMI+TERM2+TERM3 C C COMPUTAION OF V(I.J) C TM1-(Z(I,J+l)+Z(I+1,J+1))/2.-(Z(I,J)+Z(I+1,J))/2. TM2=(ZP(I,J+1)+ZP(I+1.J+1))/2.-(ZP(I,J)+ZP(I+1,J))/2. TERM1--DELT*G*(ALFA*HUV(I,J)*TM1 & +(1.-ALFA)*HUVP(I,J)*TM2)/ZDELY(J) TMI-(Z(I+1,J)+Z(I+1,J+1)+Z(I,J)+Z(I,J+1))/4. TM2-(ZP(I+1,J)+ZP(I+1,J+1))+ZP(I,J)+ZP(I,J+))/4. TERM2-VD(I,J)*(TM1-TM2)/HUVP(I,J) TM1-UV*V(I,J) TERM3--DELT*FC*TM1/8./HUVP(IJ)**2 V(I,J)-VD(I,J)+TERM1+TERM2+TERM3 60 CONTINUE C C---IMPOSE CLOSE BOUNDARY CONDITIONS C DO 70 J-1,KJO U(1.J)-0. V(1,J)=0. U(KIO,J)O0. V(KIO,J)-0. 70 CONTINUE C DO 80 I-2,KIO-1 IF(I.LT.25) THEN Ni-KYi(I-1) N2-KY2(I-1) ELSE N1-KY1(I+1) N2-KY2(I+1) ENDIF DO 85 J-1,Ni U(IJ)-0. V(I,J)-0. 85 CONTINUE DO 90 J-N2,KJO U(IJ)-0. V(I,J)=0. 90 CONTINUE 80 CONTINUE C C--ADVANCE VELOCITIES: PU(I,J) AND PV(I,J) C DO 100 I-1,KIO DO 100 J-KY1(I),KY2(I) PU(I,J)-U(I,J)/HUV(I,J) PV(I,J)-V(I,J)/HUV(I,J) 100 CONTINUE C RETURN PAGE 127 119 END C C C*.*********** SUBROUTINE DEPT ************** C** ADVANCE WATER DEPTH AND WATER SURFACE ELEVATION C SUBROUTINE DEPT INCLUDE 'DIM.FOR' C DO 30 I-2,KIO DO 30 J-2,KJO Z(I,J)-Z(IJ)+DZ(I,J) H(I,J)-Z(IJ)-ZB(I,J) 30 CONTINUE DO 29 I=,.KIO Z(I.1)-Z(I,2) H(I.1)-H(I.2) 29 CONTINUE DO 31 J-1KJO Z(1.J)-Z(2,J) H(1,J)-H(2,J) 31 CONTINUE C DO 32 I-1,KIO DO 32 J-1.KJO IF(H(I,J).LT.HMIN) THEN H(I,J)-HMIN Z(I,J)-ZB(I,J)+HMIN ENDIF 32 CONTINUE C DO 35 J-1,KJO DO 35 I-1,KIO-1 HX(I. J)-GAMA(I,J)*H(I+, J)+(-GAMA(I, J))*H(I, J) 35 CONTINUE DO 36 J-l.KJO HX(KIO, J)-H(KIO.J) 36 CONTINUE C DO 40 J=IKJO-1 DO 40 I-=.KIO HUV(IJ)-(HX(I,J)+HX(I, J+1))/2. 40 CONTINUE DO 45 I-1,KIO HUV(I,KJO)-H(I.KJO) 45 CONTINUE C RETURN END C*************** SUBROUTINE FMOVB ************* C** PERFORM THE COMPUTAION OF THE TIDAL FLODDING AND DRYING C PAGE 128 120 SUBROUTINE FMOVB C INCLUDE 'DIM.FOR' COMMON/W/BAT(150,150) C C--CHECK TO SEE IF THE FLOW COMPUTED FROM THE NAVIER-STOKES C EQUATIONS IS CONSISTENT WITH THE FLOW CONTROLLED BY THE C BOTTOM FRICTION C AFF-2. DO 5 J-2.KJO-1 NI-KX1(J) N2-KX2(J) DO 5. IN1,N2 C C EXAMINE U(I,J) C SGR-AFF*ABS(G*(Z(I+1,J)-Z(I,J))/DELX(I)) FC-116.*CN**2/HUV(IJ) *(1/3.) BFR-FC/8.*U(I,J)**2/HUV(I,J)**3. C IF(H(I,J).LT.50.) THEN IF(BFR.GT.SGR) THEN TM-G*ABS(Z(I+1,J)-Z(I.J))/DELX(I)*HUV(I,J)**3.*8/FC U(I,J)-AFF*SQRT(TM)*U(I,J)/ABS(U(I,J)) ELSE U(I,J)-U(IJ) ENDIF ENDIF C C EXAMINE V(I,J) C SGR1-AFF*ABS(G*(Z(I,J+1)-Z(I,J))/DELY(J)) FC-116.*CN**2/HUV(I.J)**(1/3.) BFR1-FC/8.*V(IJ)**2/HUV(I,J)**3. C IF(H(I,J).LT.50.) THEN IF(BFRI.GT.SGR1) THEN TM-G*ABS(Z(I, -Z(IJ I,J))/DELY(J)*HUV(IJ)**3.*8/FC V(I,J)-AFF*SQRT(TM)*V(I,J)/ABS(V(I.J)) ELSE V(IJ)-v(I,J) ENDIF ENDIF C 5 CONTINUE C C--UPDATE THE WEIGHTING COEFFICIENTS DUE TO THE MOTION C OF THE BOUNDARY C DO 10 J-2,KJO-1 IF(J.LT.25) THEN NI-KX1(J-1)+1 N2-KX2(J-1) PAGE 129 121 ELSE N1-KX1(J)+1 N2-KX2(J) ENDIF DO 10 I-N1.N2 IF(U(I,J).LT.O) THEN C C FLOODING C TM-3.*(Z(IJ)-Z(I+1,J)) TMM-(TM-H(I+1,J)+H(I,J)) IF(TMM.EQ.0) THEN BAT(I,J)-0.5 ELSE BAT(I, J)-(TM-H(I+1, J))/TMM ENDIF IF (BAT(IJ).LT.O.) THEN BAT(I,J) -.5 ELSE IF(BAT(I,J).GT.1.) THEN BAT(I,J)-1. ELSE BAT(I',J)-BAT(I,J) ENDIF ENDIF GAMA(I, J)-BAT(I, J) C ELSE C C DRYING C IF(Z(I,J).LT.ZB(I+1,J)) THEN BAT(I.J)-1. ELSE TM-3.*(Z(I+1,J)-Z(I,J)) TMM-(TM-H(I,J)+H(I+1,J)) IF(TMM.EQ.O) THEN BAT(I,J)-0.5 ELSE BAT(I, J)-(TM-H(I.J))/TMM ENDIF IF (BAT(I.J).LT.O.) THEN BAT(I,J)O0.2 ELSE IF(BAT(I,J).GT.1.) THEN BAT(I..J)-1. ELSE BAT(IJ)-BAT(I,J) ENDIF ENDIF ENDIF GAMA(I, J)-1-BAT(I, J) ENDIF ENDIF PAGE 130 122 10 CONTINUE C RETURN END C********t******s SUBROUTINE UVDEP ******t*** C** PERFORM CALCULATION OF HX(I,J), HY(I.J). HYDRAULIC RADIUS, AND C WATER DEPTH AT (U,V) GRID C SUBROUTINE UVDEP INCLUDE 'DIM.FOR' C C--EVALUATE HX(I,J) AND HY(I.J) C DO 10 I-1,KIO-1 DO 10 J-I.KJO HX(I,J)-GAMA(I,J)*H(I+1,J)+(1-GAMA(I,J))*H(I,J) 10 CONTINUE DO 20 J-l.KJO HX(KIO.J)-H(KIO,J) 20 CONTINUE C DO 30 I-1,KIO DO 30 J-1.KJO-1 HY(I,J)-BATAY(I.J)*H(I.J+1)+(1-BATAY(I.J))*H(I,J) 30 CONTINUE DO 40 I-1,KIO HY(I,KJO)-H(I,KJO) 40 CONTINUE C C--EVALUATE WATER DEPTH AT (UV) GRID C DO 50 I-l.KIO DO 50 J-1,KJO-1 HUV(I,J)-(HX(I.J)+HX(I,J+1))/2. 50 CONTINUE DO 60 I-1.KIO HUV(I.KJO)-H(I.KJO) 60 CONTINUE C C--EVALUATE THE HYDRAULIC RADIUS C DO 70 I-1,KIO DO 70 J-1,KJO HR(I.J)-HUV(I.J) 70 CONTINUE C RETURN END C**~************* SUBROUTINE T2DOT ******** ** C*** THIS SUBROUTINE IS USED TO OUTPUT THE RESULTS PAGE 131 123 C SUBROUTINE T2DOT INCLUDE 'DIM.FOR' C C--OUTPUT THE VELOCITIES C IF(NT.EQ.200) THEN OPEN(UNIT-1,NAME-'VELOCITY.DAT' ,STATUSNEW') DO 10 Il.,KIO DO 10 J-1,KJO TM1-U(I,J)/HUV(I,J) TM2-V(I,J)/HUV(I,J) WRITE(1,*) NTTMi,TM2 10 CONTINUE CLOSE(1) ENDIF C C--OUTPUT THE LOCATION OF THE BOUNDARY C IF(NT.EQ.100) THEN OPEN (UNIT-1, NAMEBOUNDARY.DAT', STATUS-'NEW') 40 DO 20 J=1,KJO DO 30 I-1.KIO IF(H(I+1, J) .GT.HMIN) THEN VRITE(1,*)J.I GOTO 40 ENDIF 30 CONTINUE 20 CONTINUE CLOSE(1) ENDIF C RETURN END PAGE 132 124 C CUGA----------MATRIX INVERSE C SUBROUTINE CUGA(X,A,B.CD,M) C C PERFORM TRIDIAGONAL MATRIX INVERSE C DIMENSION X(150),A(150),B(150),C(150),D(150) DIMENSION P(200). R(200),T(200),W(200) C MN-M-1 P(1)=A(2) DO 10 K-1.MN 10 R(K)-B(K+1) DO 20 K-2,M T(K)-C(K+1)/P(K-1) P(K)-A(K+1)-T(K)*R(K-1) 20 CONTINUE V(1)-D(2) DO 30 K-2,M V(K)=D(K+1)-T(K)*W(K-1) 30 CONTINUE X(M+1)W(M)/P(M) DO 40 KK-1,MN K-MN-KK+1 X(K+1)-(W(K)-R(K)*X(K+2))/P(K) 40 CONTINUE RETURN END C**************DIM.FOR***************** C** DIMENSION FILE C COMMON/A/Z(70.70),H(70,70) ,R(70,70),HUV(70,70),HX(70,70) S, HY(70,70),ZB(70,70) COMMON/B/U(70,70).V(70,70) ,UD(70,70),VD(70,70).UAD(70,70) A ,VAD(70,70) COMMON/C/DZ(70,70).DZ1(70,70),DZ2(70,70) COMMON/D/ZDELX(70) ,ZDELY(70),DELX(70),DELY(70) COMMON/E/GAMA(70,70),ALFA COMMON/F/CCNDELT,DC.F,WOMAGA,ATAHMIN COMMON/G/KIO.KJO,NT,KX1(70),KX2(70),KY1(70),KY2(70) COMMON/I/YM1(70.70).YM2(70,70) COMMON/J/PU(70,70),PV(70,70) COMMON/K/AX(70,70),BX(70,70),CX(70.70),DDX(70,70) COMMON/L/AY(70.70).BY(70,70),CY(70,70),DDY(70.70) COMMON/M/V(70.70) PAGE 133 BIBLIOGRAPHY Beckman, F. S., "The Solution of Linear Equations by Conjugate Gradient Method," Mathematical Methods for Digital Computers, Vol. 1, 1959, pp 62-72. Benque, J. P., J. A. Cunge, J. F. Hauguel, and F. M. A. Holly, "New Method for Tidal Current Circulation," Journal of the Waterway Port Coastal and Ocean Division, Vol. 108, No. WW3, ASCE, 1982a, pp 396-417. Benque, J.-P., A. Hauguel, and P.-L. Viollet, Engineering Applications of Computational Hydraulics, Vol. 1, Pitman Advanced Publishing Program, Boston, 1982b. Benque, J.-P., A. Hauguel, and P.-L. Viollet, Engineering Applications of Computational Hydraulics, Vol. 2, Pitman Advanced Publishing Program, Boston, 1982c. Garrett, J. R., "Review of Drag Coefficients Over Ocean and Continents," Monthly Weather Review, Vol. 7, 1977, pp 915-929. Carrier, G. F. and H. P. Greenspan, "Water Waves of Finite Amplitude on a Sloping Beach," Journal of Fluid Mechanics, Vol. 4, 1958, pp 97-109. Gopalakrishhan, T. C. and C. C. Tung, "Numerical Analysis of a Moving Boundary Problem in Coastal Hydrodynamics," International Journal for Numerical Methods in Fluids, Vol. 3, 1983, pp 179-200. Gray, W. G. (Editor), Physics-Based Modeling of Lakes, Reservoirs, and Impoundments, ASCE, New York, 1986. Hirt, C. W. and B. D. Nichols, "Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries," Journal of Computational Physics, Vol. 39, 1981, pp 201-225. Isaacson, E. and H. B. Keller, Analysis of Numerical Methods, John Wiley & Sons Inc., New York, 1966. Khosla, P. K. and S. G. Rubin, "A Conjugate Gradient Iterative Method," Computers and Fluids, Vol. 9, 1981, pp 109-121. Leendertse, J. J., "A Water-Quality Simulation Model for WellMixed Estuaries and Coastal Seas: Volume I," Principales of Computation, RM6230-RC, The Rand Corporation, Santa Monica, Calif., Feb., 1970. 125 PAGE 134 126 Lynch, D. R. and W. G. Gray, "Finite Element Simulation of Flow in Deforming Regions," Journal of Computational Physics, Vol. 36, No. 2, 1980, pp 135153. Lynch, D. R., "Moving Boundary Numerical Surge Model," Journal of the Waterway Port Coastal and Ocean Division, Vol. 106, No. WW3, ASCE, 1980, pp 425-428. Marchuk, G. I., Methods of Numerical Mathematics, Springer-Verlag, New York, 1975. Melkonian, S., Nonlinear Waves on Thin Films and Related Phenomena, Ph. D. Thesis, McGill University, Canada, 1987. Murty, T. S., Storm Surges, Canadian Government Publishing Centre, Canada, 1984. Peyret, R. and T. D. Taylor, Computational Methods for Fluid Flow, SpringerVerlag, New York, 1986. Rahman, M., "Analytical Solutions for Tidal Propagation in a Rectangular Basin," Advances in Water Resources, Vol. 6, 1983, pp 44-53. Reid, R. O. and B. R. Bodine, "Numerical Model for Storm Surges in Galveston Bay," Journal of the Waterways and Harbors Division, Vol. 94, No. WW1, ASCE, 1968, pp 33-57. Schmalz, R. A., "A Numerical Investigation of Hurricane Induced Water Level Fluctuations in Lake Okeechobee: Report 1, Forcecasting and Design," US Army Engineer District, Jacksonville, 1986. Sheng, Y. P., "Mathematical Modeling of Three-Dimensional Coastal Currents and Sediment Dispersion: Model Development and Application," A. R. A. P. Report CERC-83-2, Princeton, 1983. Smith, G. D., Numerical Solution of Partial Differential Equations, Oxford University Press, London, 1969. Weare, T. J., Errors Arising from Irregular boundaries in ADI Solutions of the Shallow-water Equations," International Journal for Numerical Methods in Engineering, Vol. 14, 1979, pp. 921-931. Yanenko, N. N., The Method of Fractional Steps, Springer-Verlag, New York, 1971. Yeh, G.-T. and F.-K. Chou, "Moving Boundary Numerical Surge Model," Journal of the Waterway Port Coastal and Ocean Division, Vol. 105, No. WW3, ASCE, 1979, pp 247-263. Yeh, G.-T. and F.-K. Chou, "Moving Boundary Numerical Surge Model," Journal of the Waterway Port Coastal and Ocean Division, Vol. 107, No. WW1, ASCE, 1981, pp 34-36. PAGE 135 127 Yeh, G.-T. and F.-F. Yeh, "Generalized Model for Storm Surges," Coastal Engineering, Proceedings of the Fifteenth Coastal Engineering Conference, Honolulu, Hawaii, 1976, Vol. 1, pp 921-940. Yih, C. S., "Stability of Liquid Flow Down an Inclined Plane," Phys. Fluids 6, 1963, pp 321-331. Zelt, J. A., "Tsunamis: the Response of Harbors with Sloping Boundary to Long Wave Excitation," Report No. KH-R-47, W. M. Keck Laboratory of Hydraulics and Water Resources, Division of Engineering and Applied Science, California Institute of Technology, 1986. xml version 1.0 encoding UTF-8 REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EHU5VZTJ0_KDLT1T INGEST_TIME 2017-12-01T19:43:46Z PACKAGE UF00076135_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES |