Citation |

- Permanent Link:
- http://ufdc.ufl.edu/UF00091092/00001
## Material Information- Title:
- A numerical hydrodynamic model for inlet-river systems
- Series Title:
- UFLCOEL-97005
- Creator:
- Seidle, Peter N., 1972-
University of Florida -- Coastal and Oceanographic Engineering Dept - Place of Publication:
- Gainesville Fla
- Publisher:
- Coastal & Oceanographic Engineering Dept., University of Florida
- Publication Date:
- 1997
- Language:
- English
- Physical Description:
- x, 83 leaves : ill. ; 28 cm.
## Subjects- Subjects / Keywords:
- Inlets -- Mathematical models -- Florida ( lcsh )
Shore protection -- Mathematical models -- Florida ( lcsh ) Ponce de Leon Inlet (Fla.) ( lcsh ) - Genre:
- government publication (state, provincial, terriorial, dependent) ( marcgt )
bibliography ( marcgt ) non-fiction ( marcgt )
## Notes- Thesis:
- Thesis (M.S.)--University of Florida, 1997.
- Bibliography:
- Includes bibliographical references (leaves 81-82).
- Statement of Responsibility:
- by Peter N. Seidle.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- All applicable rights reserved by the source institution and holding location.
- Resource Identifier:
- 37856747 ( OCLC )
## UFDC Membership |

Full Text |

UFL/COEL-97/005
A NUMERICAL HYDRODYNAMIC MODEL FOR INLET-RIVER SYSTEMS by Peter N. Seidle Thesis 1997 A NUMERICAL HYDRODYNAMIC MODEL FOR INLET-RIVER SYSTEMS By PETER N. SEIDLE A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 1997 ACKNOWLEDGMENTS I would like to thank my advisor, Dr. Wang, for his guidance and time I have used in my research. I would also like to thank Professors Dean and Thieke for serving on my committee and reviewing this thesis. I would like to thank the Coastal and oceanographic Engineering (COE) Department for providing me with a research assistant position and funding for my studies. With the opportunity and resources the COE Department has provided me, I am able to fulfill my personal goals and can better contribute to the field of coastal engineering. I would like to thank each and every professor in the COE Department for piqued my interest and curiosity in their areas of research. I would like to especially thank Dr. Lihwa Lin for his patience, help, and guidance in my studies and research. I truly appreciate the help from Professor Sheng, Dr. Eduardo Yasuda, and Justin Davis for helping me with numerical modeling. Gratitude is owed to Wendy Smith for her help as well. I also thank each and every fellow student who has helped me with class work. Many thanks to Becky, Helen, John, Lucy, Sandra, and Sidney for the administrative and computer support. And saving the best for last, I would like to thank my parents for their love, support, and encouragement, as well as Micheal O'Shea and Pete's Wicked. TABLE OF CONTENTS ACKNOWLEDGMENTS.................. LIST OF FIGURES.................. ... . .. . .. . 1 ... . . . . . .vi ABSTRACT..................................................... ix 1. INTRODUCTION.................. 1.1 Background.............. 1.2 Objective............... 2. GOVERNING EQUATIONS........... 3. NUMERICAL METHODS............. 3.1 Linearized System...... 3.2 Grid System............. 3.3 Non-Linear Terms....... 3.3.1 Friction Term . 3.3.2 Convective Term 3.3.3 Diffusion Term 3.4 Numerical Stability.. .1.. .. . .. . . .1.. .. . .. . . ... . . .. . . .3 ... . . .. . . .4 ... . . . . . .10 ... . . . . . .10 ... . . . . . .15 ... . . . . . .17 ... . . . . . .17 ... . . . . . .19 ... . . . . . .23 ... . . . . . .24 4. BOUNDARY CONDITIONS...................................... 27 4.1 Open Boundary Conditions.......................... 29 4.1.1 Tide........................................ 29 4.1.2 Neumann..................................... 30 4.1.3 River and Canal............................ 31 4.2 Solid Boudaries.................................... 35 5. MODEL 5.1 5.2 5.3 VERIFICATION....................................... 36 Wind................................................ 36 Seiching........................................... 38 Massachisetts Bay.................................. 41 6. MODEL VALIDATION......................................... 49 6.1 Field Data......................................... 49 6.1.1 Bathymetry.................................. 50 6.1.2 Tide........................................ 51 6.1.3 Current..................................... 54 6.1.4 Drogues..................................... 57 6.1.5 Wind........................................ 60 6.2 Spectral Analysis.................................. 60 6.3 Grid................................................ 62 6.4 Open Boundary Conditions.......................... 64 6.4.1 Tidal OBC................................... 64 6.4.2 Neumann OBC................................. 65 6.4.3 Canal OBC................................... 66 6.5 Results............................................ 68 7. MODEL 7.1 7.2 APPLICATION........................................ 71 Single Jetty....................................... 72 Pair of Jetties.................................... 73 8. CONCLUSIONS AND RECOMMENDATIONS......................... 76 8.1 Conclusions........................................ 76 8.2 Recomminendations................................... 78 LIST OF REFERENCES.......................................... 81 BIOGRAPHICAL SKETCH......................................... 83 LIST OF FIGURES 2 .1 Relationship between drag coefficient (shear stress coefficient) and wind speed......................... 9 3.1 Rectilinear cell ................................... 16 3.2 NumericalCells applied to friction term ........... 18 3.3 Numerical cells applied to convective term ........ 20 3.4 Numerical cells applied to convective term ........ 22 4.1 Example of grid showing solid and open boundaries. ................................................28 4.2 Example tidal and Neumann open boundaries ......... 33 4.3 Example of a grid for a river or canal open boundary condition .......................................... 35 5.1 Comparison of analytical and numerically modeled set-up ............................................. 38 5.2 Numerically modeled seich wiith 1 node .............40 5.3 Numerically modeled seich with 2 nodes .............41 5.4 Massachusetts Bay .................................. 42 5.5 Geometry of 2D analytical solution................. 44 5.6 Analytical surface elevations at high tide and maximim ebb velocity vectors....................... 45 5.7 Maximum velocity vectors from numerically modeling Massachusetts Bay .................................. 47 5.8 Numerically modeled surface elevations at high tide in Massachusetts Bay ............................... 48 6.1 Bathymetric survey of Ponce de Leon Inlet ......... 51 6.2 Location of stations where tide and current data were collected ..................................... 52 6.3 Continuously observed tidal record and discrete current measuremets ................................ 53 6.4 Flood velocity data taken inside the canal ........ 55 6.5 Ebb velocity data taken in the canal ...............56 6.6 Flood current patterns studied using drogues......58 6.7 Ebb flow patterns studied using drogues ............59 6.8 One tidal ccyle extracted from tide data .......... 61 6.9 Discrete spectrum of the tidaal cycle ..............62 6.10 Area of Ponce de Leon Inlet which the grid is applied to ......................................... 63 6.11 Simulated tidal OC signal compared to the observed tidal signal ....................................... 66 6.12 Maximum flood velocity vectors for numerical model validation ......................................... 68 6.13 Maximum ebb velocity vetcors showing numerical model validation ......................................... 69 7.1 Flood velocity vecors of a single jetty applied to Punta Gorda, FL ................................... 72 7.2 Ebb velocity vecors of a single jetty applied to Punta Gorda, FL ................................... 73 7.3 Flood velocity vecors of a pair of jetties applied to Punta Gorda, FL ................................ 74 7.4 Ebb velocity vecors of a pair of jetties applied to Punta Gorda, FL ................................... 75 8.1 Bathymetric smoothing ............................. 68 8.2 An idealized channel/bay system .................. 80 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 NUMERICAL HYDRODYNAMIC MODEL FOR INLET-RIVER SYSTEMS By Peter N. Seidle May, 1997 Chairperson: Dr. Hsiang Wang Major Department: Coastal and Oceanographic Engineering A numerical model is developed for predicting tidal induced near shore hydrodynamics of an inlet-river system or inlet canal network. The model couples a two-dimensional model for the inlet region and a one-dimensional model representing the river and/or canal network. The model is based on depth averaged equations of motion including the convective acceleration, friction, turbulent diffusion, and wind stress terms. The numerical method used to solve the equations is the factorization scheme. The model can be driven using tides and river discharges. A canal system is treated as a river with zero discharge. The model is verified by comparing known analytical solutions to the numerical results for selected cases. The numerical model is then applied to Ponce de Leon Inlet which serves a canal system for the City of Punta Gorda, FL. The model is used to study the effects on the hydrodynamic processes when jetties are added to Ponce de Leon Inlet. CHAPTER 1 INTRODUCTION Coastal regions are important to communities for their leisure purposes and resources. Access to these resources may be through coastal waterways, residential canal system, harbor, etc. it is important to keep these resources accessible not only for leisure purposes, but also for commercial, fish and wildlife, and ecological reasons. There have been numerous studies and projects to learn how keep these resources open, accessible, and stable in all aspects. One area of study is numerical modeling. This thesis will develop a numerical model such that hydrodynamic studies of inlet improvement measures such as jetties, rip-rap fill, breakwaters, groins, etc., can be evaluated. 1.1 Backaround It is important to understand the hydrodynamics of water bodies, large and small, so that the benefits gained from the coastal resources can grow, or at least remain stable. Technology has given us the tools so that these hydrodynamic 2 studies can be conducted in a faster, more efficient manner. Numerical models, confirmed by analytical models and augmented by physical models, can be developed to learn more about the micro- and macro-scale processes which govern the hydrodynamics of coastal waterways, harbors, canals, and inlets. As an aid to the design of inlet and waterway improvements, the effects of the tidal circulation, a hydrodynamic process, near the inlet improvement measures can be studied using numerical modeling. Assuming that the tidal flow is not highly stratified and that the inlet width is large compared to its depth, 2D numerical modeling should be adequate to predict the tidal flows and velocities. As previously mentioned, researchers have conducted innumerable studies and published many papers about hydrodynamics of estuaries and numerically modeling tidal areas. Arcilla (1989) developed a general algorithm for generating numerical models. Estuaries and bays are a focus of study in two-dimensional circulation: Cheng (1982) studied Tampa Bay, FL for circulation; Schomaker (1983) developed a circulation model of Monterey Bay, CA; Reidel and Wilkinson (1976) studied Cockburn Sound, Western Australia; Reeve and Hiley (1992) studied the southern North Sea; Bathen (1976) developed a model to study Keehi Lagoon, HI; etc. But these 3 models were developed for larger basins than a navigational inlet and are only applicable to the water body for which they were designed, so a more generalized two-dimensional numerical model is needed which can be applied to different types and sizes of inlets. 1.2 Obetv The objective of this project is to develop an implicit two-dimensional numerical model which can easily be applied to a multitude of inlets with different geometrical configurations. The model will also be designed to aid in the placement of jetties, groins, rip-rap fill, or other such inlet improvements by simulating the tidal flows. Analysis of this tidal process will allow for the optimum placement of inlet improvement measures to ensure maximum benefits from the changes. This report will (1) present the equations of continuity and motion used; (2) derive the numerical method used; (3) show the differencing techniques applied to the governing equations; (3) verify the numerical models efficacy; (4) validate the model using data obtained from Punta Gorda, FL; and (5) show how the model is applied as an aid to inlet improvement design. CHAPTER 2 GOVERNING EQUATIONS The governing equations of motion used for this model were derived from the Navier-Stokes equations of fluid motion for a non-compressible, Newtonian fluid. The simplifying approximations applied were as follows: l)By order-of-magnitude analysis, the vertical accelerations are negligible compared to gravity. Hence, a hydrostatic approximation is applied. 2)Since the Rossby number, the ratio of Coriolis force to inertial forces, is small for inlets this model is intended for, the Coriolis force is neglected. 3)Assume the turbulent Reynolds stresses are proportional to velocity gradients by an eddy viscosity coefficient and that the viscous stresses are negligible compared to the turbulent 5 stresses. Also, for simplicity, assume the eddy viscosity coefficient is constant for the inlet. 4)Assume atmospheric pressure is constant and uniform for the region being modeled. 5)Assuming the inlet is well mixed and that there are no density driven currents, the baroclinic pressure gradient can be neglected. 6)Boussinesq approximation is applied to the barotropic pressure gradient. This assumes that the density gradient is small enough such that the density is assumed to be uniform. The continuity and resulting equations of motion are as follows: 89 BuD 8vD ----+uD+-vD=o (2.1) Bt ax By Bu Bu Bu 89 f 82U B2U x-wnd a+u-+v=-g- fuu+A,(-+2)+ (2.2) 8t 8x ay ax 8D 8x2 ay 2 pD 8V 8V 8V 84a f 8-%>2 B2V \ -wn -a +u-a+v--g l f V+A-v a2v +-)+- (2.3) at ax ay ay 8D ax2 y 2 pD where t=time u=depth averaged velocity in the x-direction, v=depth averaged velocity in the y-direction, u=depth averaged velocity vector, D=total depth, p=water surface elevation, g=gravity, x,y=right handed Cartesian coordinates, j=Darcy-Weisbach friction coefficient, A,=horizontal diffusion coefficient, tx=x-direction wind shear, Ty=y-direction wind shear, and p=density of sea water. There are two empirical coefficients in the above equations, the Darcy-Weisbach friction coefficient and the 7 eddy viscosity coefficient. In numerical modeling they are often adjusted to tune the model based upon field data. There are numerous studies concerning the friction coefficient for open channel flows. In engineering applications, the most common methods of developing a friction factor use the Manning coefficient or Chezy coefficient. The friction coefficient is a function of many types of environmental conditions, such as type and amount of vegetation on both the bottom and banks, and some of these conditions have been quantified in terms of the Manning coefficient (Henderson, 1966). For a sandy bottom along an open coast, a Darcy-Weisbach friction coefficient of approximately 0.01 is commonly used. The value, as well as the range of variability, of the eddy viscosity coefficient is far less well established. This is because eddy viscosity is a flow property which is related to turbulent mixing. Since flow turbulence is a complicated phenomenon and the scale of mixing varies greatly from case to case, it is not a simple matter to model its mixing effects using only one representative coefficient. In an oceanic environment the eddy viscosity coefficient is known to vary enormously. The order of magnitude for the coefficient suggested by Pedolsky (1979) ranges from 10 to 10000. For each scenario numerically modeled, a stability condition (discussed in chapter 3.4) must be considered in selecting the 8 eddy viscosity coefficient. The stability condition relates the grid size of the numerical model and time step of the numerical method to the coefficient. The wind shear, Twind is somewhat empirical as well and may be necessary to validate the numerical model. If field studies were conducted on a windy day, then the wind must be accounted for in order to validate the model. Once the numerical model has been validated, the wind speed and direction may be changed or set to zero depending on the long term predominant conditions. The wind shear forcing term can be written as: Z .id P.J Cd - i WjosinO (2.4) pD p D -wind Pair Cd (2.5) ___ -WjcoO 25 pD p D where Cd=drag coefficient, W10=wind speed at 10m, and O=wind bearing. Although the drag coefficient is empirical, researchers have e V Windsped Z Figure 2.1 Relationship between drag coefficient (shear stress coefficient) and wind speed. Figure was taken from Roll (1965). found that it is on the order of i0-. The drag coefficient could be found graphically from figure 2.1. Also, Wilson (1960) concluded that for light winds (typically less that 10 m/sec), Cd=1.49 x I0' and for stronger winds, Cd=2.37 x i0-3. For validating this numerical model, Wilson's findings were used. The governing equations (2.1, 2.1, and 2.3) are solved numerically for the problems posed in this thesis. The numerical scheme and the applicable boundary conditions are discussed in chapters 3 and 4. J iiiiiiii ii i!i CHAPTER 3 NUMERICAL METHODS The numerical method adopted for the model is the socalled "factorization" scheme. This method of implicitly solving a set of linearized equations was developed to study tidal currents in the vicinity of the Mississippi Sound in the Gulf of Mexico (Sheng, 1983). 3.1 Linearized System The governing equations (2.1, 2.2, and 2.3) are linearized in order to develop the "factorization" numerical method. The governing equations now become: a_+Dau +DA=0 (3.1) at ax ay u (3.2) at .-a, (3.3) a3t ay The above equations 3.1, 3.2, and 3.3 can be combined as a system of equations by letting w=[ u, V]T The system of equations can be written as: Iwt+Awx+BWy=O (3.4) where I is an identity matrix and A and B are coefficient matrices as follows: A=g 00 000 B= 0 0 0 g0. The subscripts t, x, and y indicate the first partial derivative respective to the subscript. Applying a forward finite difference approximation to the time derivative yields: n+l nn-i Wnt -+Awn +Bw i=0 (3.5) At x11 y Multiplying by At and factoring out wn"1 results in: wn+1(1+AtA8x+AtB8 )=w n (3. 6) where 6 is the finite difference respective to the subscript and the superscripts indicate the time step. The left hand side (LHS) of the equation can be expanded to: wnIl(1 +AtA8x+AtB86 +(At)2AB8x8 y)-02=w (3.7) where 02 denotes a 2nd order error term. This term will be resolved and removed in later steps of this derivation. A term of AtB6 is added and subtracted to the right hand side y (RHS) of the equation with no net change. The LHS of Equation 3.7 is then factored into two terms to give: wn+I(I +AtA6x)(1+AtBy)-02=(I -AtB8)w n+AtB8 w (3.8) An intermediate time step may be introduced, which is also referred to as time splitting. Letting w* represent the 1-AtB8Y intermediate time step, n+2 w*-n and equation 3.8 1+AtA8x can be rewritten as: w n+(1 +AtAx)(1 +AtB6y)-02=(1 +AtA6x)w *+AtB8 yw n Moving the error term to the RHS of the equation and reexpanding it the above equation then becomes: w n+(1 +AtA8)(1 +AtB6 )=(1 +AtA8x)w +AtB8 w n+(t)2AB88,w " (3.10) Factoring out AtB8y from the two last terms on the RHS yields: wn+(1 +AtA6x)(1 +AtB8y)=(l +AtA8 )w +AtB6y(l +AtA8x)w n (3.11) Lastly, removing the common term, 1+AtAx from the equation results in: (1 +AtB8 )wn.+l=w*+AtB8yWn (3.12) In summary, the time splitting equations used to solved the linearized system of governing equations (3.1, 3.2, and 3.3) are written as follows: (3.13) (3.14) (1 +AtA8x)w *=(1 -At )WB8 (I +AtBy)w"+1=w* +AtB yw n Expanding equation 3.13 gives the following: 7 +AtD8xu ="-AtD6yv " (3.15) (3.9) u *+Atg6x8' =un (3.16) v* =v "-Atg6 q]" (3.17) y and expanding equation 3.14 gives: "n +AtD8v""l=p*+AtD v n (3.18) un+I=u (3.19) v n+Atg6 pnI =v +Agyqn (3.20) Note that equations 3.15 and 3.16 are a system of simultaneous equations where T' and u* are the two unknowns, as well as equations 3.18 and 3.20, where ,l'n and vn"i are the unknowns. Since equation 3.17 is independent of r* and u* it can be substituted directly into equation 3.20. Likewise, equation 3.19 states that unl =+=u so equations 3.15, 3.16, 3.18, and 3.20 can be adjusted accordingly. The four resulting equations are as follows: ' +AtD8xun+=I"-AtD8 yv n (3.21) Atg8x'l*+un I=u n (3.22) rg n+1+AtD6 Vn+l=r*+AtD8V n (3.23) A t n+I+vn I=V (3.24) The above "factorization" method is similar to the Alternate Direction Implicit (ADI) method in that equations 3.21 and 3.25 are solved in the x-direction (x-sweep) and that equations 3.23 and 3.24 are solved in the y-direction (ysweep). It is also important to note that equations 3.22 and 3.24 are essentially the momentum equations in the x- and ydirections, respectively. This will allow the non-linear terms to be added. Moreover, when equation 3.21 is substituted into equation 3.23 for r* ,the original continuity equation (Equation 3.1) is recovered. 3.2 Grid System A rectilinear staggered grid system is used for this numerical model, shown in Figure 3.1. In applying finite differences to equations 3.21 through 3.24 the staggered grid system is examined more closely. A forward difference is vij+I Ui+Ij Figure 3.1 Rectilinear cell applied to the velocities for the continuity equations (2.21 and 3.23) and may be written as: (3.25) (3.26) l n+1n 1, +AtD(unj-u )=T-AtD(v,-vin) A g (l -T 7 1 + U 1= u n tg~ij-) l +uj =ij A backward difference is applied to the surface elevation for the momentum equations (3.22 and 3.24) and are as follows: Atg( n+l_ n+l+n+l-v n grij l) = n+l n+ n+I n rij +AtD( vjl-i )=il* +AtD( vij, - (3.27) (3.28) 3.3 Non-Linear Terms The non-linear terms in the governing x- and y- momentum equations (2.2 and 2.3) which have been neglected thus far can now be added to the numerical x- and y- momentum equations (3.26 and 3.28), respectively. It is necessary to examine the adjacent elements of the grid shown in figure 3.1 in order to properly examine each non-linear term of the momentum equations individually. 3.3.1 Friction Term: The frictional forcing term in the x-direction momentum equation is -Lulu This term is non-linear in that it is 8D quadratic in terms of flow velocity vector. Here the velocity vector term, Jul is expressed numerically in terms of the known velocity components in the previous time step, such as u and vn' Furthermore, since the x-direction momentum equation, eq. 3.26, is solving for un'at the grid center depicted in Figure 3.2, lul must be calculated at the proper location. The numerical equation used to calculate the Figure 3.2 Numerical cells applied to friction term. magnitude of the velocity vector for the x-direction momentum equation is as follows: n n n u (u )2 VI- I +1 I I V-1 ij)2 (3.29) lI l = uid)+ 4 The numerical expression of the frictional forcing term which is added to the LHS of x-direction momentum equation, eq. 3.26, is given as follows: n n n n f ,n 2 V I-Ii+I +Vi+il +Vi-lj+vifi2 n+l +v++ 4 ij 8D1 4 Similarly, the frictional forcing term added to the LHS of the numerical y-direction momentum equation, eq. 3.28, is: n n n n ( Id+ +J d ij-1 i+lj-12+(f. 2. + + 1 4 ) 2 v 7 3.3.2 Convective Terms: The convective terms have directional characteristics which must be appropriately accounted for in the numerical scheme. That is, the down stream behavior is affected by the upstream convection phenomenon (Hirsch, 1994). A preferred differencing technique is the upwind method, which is similar to a backward difference written with respect to the direction of the velocity vector. Because of this directional phenomenon, the direction of the velocity must be predetermined in the numerical expression of these convective terms. For the x-direction momentum equation, the convective du du acceleration terms are u-+v- taken from equation 2.2. In ax ay order to be included with the numerical equation of motion (equation 3.26), finite differencing must be applied to these convective acceleration terms such that they are congruous Uij Figure 3.3 Numerical cells applied to convective term. with the upstream convection phenomenon. Both convective au a terms, u- and v- need to be examined separately. ax ay To analyze the term u- first the velocity direction ax needs to be determined. An average of the velocities between the two cells (figure 3.3), ui-112j and u+112, are used to identify the direction of the velocity so that the proper n upwind differencing method may be applied. If ui 1/2j is positive, then the upwind difference is u'n-u The ij i-1j numerical equation used to evaluate the first convective term is: U U i+/2j n n u=u. uij-ui-f (3.30) -ax -1 /2 Ax However, if ui+1/2j is negative, then the upwind deference n n is uij-utl and the equivalent numerical difference is as follows: n n aU -Un Uij-Ui+lj (3.31) ax i12j Ax au The second convective term, v-- can now be examined ay with the grid shown in figure 3.4 using the same steps in analyzing the first convective term. If vin/2j is in the positive direction, then the numerical convection term is as follows: n n y pin Uli/ i_ ii-1 (3.32) Vay i-1/2' Ay And, if Vi1/2,+1 is negative, the numerical differencing is: n n u nH Uig-Uij+l (3.33) iV_ /2;i+1 ay Ay At every time step, the direction of the velocity is determined for each element so that either equations 3.30 or uij+1 A Vi-1/2)+l v i-1/2j uij-1 Figure 3.4 Numerical cells applied to the convective term. 3.31 and either equations 3.32 and 3.33 are added to the numerical x-direction momentum equation, eq. 3.26. The same process used to derive the numerical equations au Bu for u-+v- can be used to find the numerical equations for ax ay the convection terms in the y-direction momentum equation, av 8v u-+v- For a positive velocity direction, the two ax ay convective terms are as follows: n n av u- i vi-Ij (3.34) ax -Ui-112 Ax and n n v n Vi Vi;V-ay=E Vi-1/2 y (3.35) Likewise, if the velocity is negative then the terms are as: n n aV n Vij -ViIlj -ua-- u 1il/2 A 8x 1/2; Ax (3.36) and n n (9v n Vii'-Pid+1 ay fi-1/2 Ay (3.37) 3.3.3 Diffusion Term: The diffusion term is linear provided that the diffusion coefficient is not a function of position, time, velocity, or water surface fluctuation. The numerical expression for the x-direction diffusive term: A 2u + 2U aX2 By2 is obtained by simply applying a center difference to yield: n n N N U n+n A(ui+lzj i+u-lj uid+-jl ij u i-1) Ax2 A2 which is included on the RHS of the numerical x-direction momentum equation, eq. 3.26. The terms added to the RHS of the y-direction momentum equation, eq. 3.28 are: A n( 'n -2vyn" + v nzl -2nin +V- l) Av.i+lj 4j Vi-ij ti- I S V IJ ) Ax2 'y2 3.4 Numerical Stability Numerical stability criteria for the non-linear system presented here are not established at present. However, the stability criteria of a linear system are used instead. There are two criteria which need to be met for numerical stability: one related to tide propagation and another to diffusion. For propagation stability, the well-known Courant criterion is invoked. The Courant Number used in the present model is defined as: Cr= At ( +lu l+lvl) (3.38) AX2 2Ay2 where At=time step, Ax=x-direction grid spacing, Ay=y-direction grid spacing, h=maximum depth, g=gravity, lul=maximum x-direction velocity superimposed on the tidal velocities, and IvI=maximum y-direction velocity superimposed on the tidal velocities. For a chosen grid size the time step, At, is then constrained by the above equation for a critical Courant number, Cr. Generally, for implicit numerical techniques the critical Courant number can be greater than unity, whereas for explicit methods the Courant number must be less than unity to maintain numerical stability. For the (ADI) method, the Courant number is suggested to be less than 5 to flows that are not consistent with observed nature (Weiyan, 1992). The advantage of using the factorization method is that stability is maintained for Courant numbers of about 10 (Sheng, 1983). For diffusive stability the eddy viscosity coefficient is examined. The following condition has been suggested by Hirsch (1994): AV AX, +AY 2 < 1 (3.39) At 2 This criterion is used in the present model to guide the selection of the diffusion coefficient depending on the spatial grid sizes and the time step size. In summary, the grid size, Ax and Ay, are first selected. Next, the time step, At, and the eddy viscosity coefficient are then selected using the Courant number relationship, eq. 3.38, and the aforementioned diffusion criterion by Hirsch (1994), eq. 3.39. CHAPTER 4 BOUNDARY CONDITIONS This chapter addresses "geometrical" boundary conditions (BCs) used to define, drive, and control the numerical model. Since the rectilinear grid is used, the numerical model must be able to differentiate land and water cells. The land and water cells are what define the model's geometry, so that an inlet, lagoon, or harbor can be numerically simulated using a rectilinear grid. The bottom and surface BCs are addressed by the governing equations as friction and wind stress forcing terms, respectively, which were defined in chapter 2. The following sections in this chapter will address the details of applying the land and water boundary conditions. The "geometric" boundary conditions (ECs) used in the numerical model can be broken into two types, open and solid An open boundary (OB) is a cell which is located on the edge of the grid that covers the modeling area. The input to the model, such as tide, river discharge and length, or canal length, are applied to the OB cells. Inside the modeling j=J M I;:.. .. .. ... .......... 64 . .....63 .... =! i "i2 ii i iIi Figure 4.1 Example of grid showing solid and open boundaries. area, some cells are designated either as land or water to define the model's geometry. The land/water interface creates a solid boundary. Figure 4.1 shows an arbitrary example. Cells (6,2), (6,3), and (6,4) are defined as OBs, so input must be applied to these cells to drive the model. Whereas 29 cells (1,4), (1,5), (2,2), and (1,5) are defined as land or solid boundary cells. 4.1 OpTen Boundary Conditions Four types of open boundary conditions (OBCs) are incorporated in the model; tidal, Neumann, canal, and river. Though a Neumann OBC does not drive the numerical model as a tidal OBC would, it is used to control the hydrodynamic stability of the numerical model in certain tidal OBCs. There are some considerations to be made when selecting the model's grid and designating the OB cells as either a tide, river, or canal BC. In practice, one edge of a cell may be designated as a tidal OB and a different edge may be defined as a canal or river OB. These methods of application will be discussed in the following sections. 4.1.1 Tide The tidal OB is defined by water level changes such as a tidal fluctuation of the form r=Acos(ut+( ) where A is the amplitude, a is the period, and 0 is the phase angle. Though it is not necessary for all scenarios, the numerical model can accept a real tide as input using up to four Fourier series 30 harmonics. The Fourier series representing a mixed tide would be as follows: 4 il=E (A~cos(at+4,)) (4.1) := 1 When applying the tidal OB, the velocities remain as unknowns, otherwise, the model would be over-constrained and could not be validated. The edge of the tidal OB cells should be of sufficient distance offshore from the region of interest. 4.1.2 Neumann If the area that is going to be modeled requires placing tide OB cells on the corners, then the a Neumann BC must be incorporated. Figure 4.2 shows the grid which would be applied to such a situation. Where the cells marked with a T indicate where a tidal forcing condition is imposed. The cells marked with an N are where a Neumann BC is imposed. The Neumann BC requires that the surface elevation gradient equal zero between the OB cells and the cells in the J-1 row or the 2nd row of cells. j .: i i i .:'~ i i:.. . .. .. .. .. .. .. ..... . . . .... . 11.: :::::::: ...... ..... X ... . . . ..:::.. ............... J -1 :* ": . . .... .. ................ T T .. !..'.. :!.... ...: ............... T T T T T 2 ..:......:. T Figure 4.2 Example of tidal and Neumann open boundary conditions. 4.1.3 River and Canal In order to apply the river or canal OB, the velocities must be perpendicular to the edge of the grid. If there is circulation or lateral currents, with respect to the river or canal, then the model area must be extended out to where the velocity vectors become perpendicular to the OB. Also, in 32 applying a river OB, the head of the river must not be influenced by the tide. The discharge at the head of the river is used to drive the model and should be constant, though the model can be easily modified to allow for a timedependent river discharge. The model becomes unstable or may produce flows that are not indicative of nature if the tide propagates to the head of the river. A canal is treated as a river, except that the discharge at the head is zero. Moreover, since the discharge at the head of a canal is zero, the model will be stable and valid if the tide propagates to the head, or end, of the canal. The length, depth, and discharge at the head, as well as grid spacing and friction coefficient, are needed as input at the river/canal OB (figure 4.3). In addition, the user needs to define a grid, or cell, length for the river or canal boundaries. The cell can be significantly larger compared to that of cells in the model area. It is recommended that the river or canal be divided into 15 to 30 cells. If any more cells are used, then the computing time required to run the model is compromised. The governing equations used by the numerical code in the river or canal cells are as follows: L + auD=o (4.2) at ax River or Canal OB Cell Size Discharge at head N' Cell Figure 4.3 Example of a grid for a river or canal open boundary condition. au + au a,, f JuJu at ax ax 8D (4.3) where u=depth averaged velocity in the direction of the river or canal, t-time, I I I I I I D=total depth, Tl=water surface elevation, g=gravity, x=direction of the river or canal, and f=Darcy-Weisbach friction coefficient. In the previous chapter (chapter 3.3), statements regarding the governing equations dealing with the friction coefficient and convective term also apply to the above equations. The model uses the surface elevation from the previous time step and the river/canal discharge as boundary conditions to calculate the new velocity at the model's OB cells. The computed velocities at the OB cells are then used as a boundary condition to drive the model. The new surface elevations calculated for the cells at the OB are then used as input into the river/canal OB subroutine for the next time step. Since the flow at the OB is assumed to be perpendicular to the boundary, the transverse velocity at the OB is set to zero. For example, since the velocity vector at the OB in figure 4.3 should be in the x-direction, the v-velocity is set to zero. At all OB cells, the model becomes unstable if diffusive or convective terms of the momentum equations are included. This prevents eddies from occurring at the boundaries which 35 would cause instability of subsequent computations. If field data shows that there is circulation at the edges of the numerical model's grid, the grid area should be enlarged to an area where there is no circulation found by the field studies. 4.2 Solid Boundaries Finally, at land boundaries, the flow condition is satisfied by setting the normal velocity component to zero. There is also a non-slip condition being imposed so that the velocity parallel to the boundary is set to zero as well. CHAPTER 5 MODEL VERIFICATION The numerical code must be compared to analytical cases to verify the efficacy. The model's efficacy was tested using wind set-up and seiching in a closed basin, and by replicating tidal currents in Massachusetts Bay. Once the model is verified, it can be confidently applied to various inlets. 5.1 Wind The analytical set-up due to wind stress is (Choi, 1992): ll=Tind~x L) 51 pgD 2 where Twind =Shear stress by wind, p=density of water, x=distance from middle of basin, g=gravity, D=depth, and L=length of basin. A closed, flat, rectangular basin was used to develop the analytical solution. The length of the test basin which the wind was blowing across, was 450 mn long, with a grid size of 50 m. The depth of the basin was constant at 2 M and the wind speed was constant at 5 rn/s. The wind shearing force, calculated using the respective drag coefficient discussed previously in the governing equations in chapter 2, equals 0.0448 Pa. Figure 5.1 shows the surface elevation of a cross section from the numerical model compared to the analytical wind set-up. Only one cross section is presented since, as expected, every cross section has the same wind set-up profile. The reason there is a difference between the numerical and analytical set-ups is that the model's output is rounded to 0.1 nun, where the analytical line does not include any round off. For example, at x=25 m, the numerical set-up is -0.0004 m and the analytical set-up is -0.0004386 mn. If the numerical model's output included more significant digits for this test, then a closer comparison could be made. Moreover, if the analytical solution was rounded off to the 1/10 mm, then the two lines would correspond exactly. Wind Set-up 200 250 Distance, (m) Figure 5.1 Comparison set-up. of analytical and numerically modeled 5.2 Seichina The model was also tested for its ability to replicate a standing wave in a closed basin. The basin length required for seiching is: Ln 2 (5.2) x 10, where /=basin length, n=number of nodes, and L=wavelength of long wave. For a channel which is 2 m deep and has a tide period, T, of 12 hrs, the wavelength is 191,352 m. The basin lengths which would form a seiche with 1 and 2 nodes are 95,676 m and 191,352 m, respectively. The rectangular test basin used for both cases is 3 m deep and 450 m wide. Also, for both cases (1 and 2 nodes), the number of cells used for the width and length were 9 and 50, respectively. Which means that Ax for the two different cases are 1913.52 m and 3827.04 m, and Ay was 50 m for both cases as well. The Courant number for both cases was 2. The model was driven by specifying te water surface elavation at the right hand side of the basin as a function of time, i=Acos-) where A, the amplitude, was set to 0.1 m and T, period, was set to 12 hrs. Figures 5.2 and 5.3 present the surface elevation for a cross section of the basin for 1 and 2 nodal standing wave at times 1/4T, 1/2T, 3/4T, and T. Note that the numerical results are reasonably good, however, they did exceed the 40 I (U 1 .... ....I 0 10 20 30 40 50 60 70 80 90 100 0.1 S -0 .1 I I I I I I 0 10 20 30 40 50 60 70 80 90 100 .1 I I I I I I I -0.11 ... 0 10 20 30 40 50 60 70 80 go 1iO0 Distance, (kmn) Figure 5.2 Numerically modeled seiche with 1 node at time, from top to bottom; (a) =1/4T, (b) = 1/2T, (c) =3/4T, and (d) =T. input amplitude because friction was set to 0.0001 in the numerical model. The small friction coefficient did not allow enough energy to be damped from the system. Though it is not evident for times 1/4T and 3/4T, there were surface fluctuations which were less than 0.05 mm. The initial surface conditions were evaluated using a cosine function to represent a seiche, so the model quickly reached equilibrium and became accurate. I-I1 I I I I I I I o 10 20 30 40 50 60 70 80 90 10 -0.2 L o 0 i1 20 30 40 50 60 70 80 90 C .0 U -0.2 L 0 10 20 30 40 50 60 70 80 90 Ic Distance, (kmn) Figure 5.3 Numerically modeled seiche with 2 nodes at time, from top to bottom; (a) = 1/4T, (b) = 1/2T, (c) 3/4T, and (d) = T. 5.3 Massachusetts Bay The final verification of the model's efficacy was to reproduce the results obtained from an analytical solution of the 2D equations of continuity and momentum with geometry representative of Massachusetts Bay (shown in figure 5.4). The analytical solution for tide, derived by Briggs and Madsen elgure -.4 Massachusetts Bay. Figure taken from Briggs and Madsen (1973) 43 (1973), to the linearized continuity and momentum equations (eqs. 3.1, 3.2, and 3.3) is: = c 2m~sinmy(sinknX2 sinkjx1) coshmy =Ae i [Cos mYnE mk (x2-xI) sinhmyocosknx] (5.3) where Mn n x and x., x, X21, and y0 are the geometric dimensions of the bay as defined in figure 5.5. For simplicity, the depth, h, was kept constant. The tide at a single location anywhere in the region must be known and is used to solve for A.. Briggs and Madsen used the tide at the Boston Light House as the driving force in their analysis. However, Boston Light House is inside the bay, not at an open boundary, so the tide measurement there could not be used to drive this numerical model. From examining the tide (figure 5.6) at the open boundary, found analytically by Briggs and Madsen, a value of Figure 5.5 Geometry of 2D analytical solution. Figure taken from Briggs and Madsen (1973). 4.30 ft can be used to drive the numerical model. The depth that Briggs and Madsen used was 120 ft. They also used a sinusoidal tide as input at Boston Light House, so it will be assumed that the tide at the entrance to the bay will be a sinusoidal tide as well. Figure 5.6 also shows the analytically derived maximum ebb currents for Massachusetts Bay. J (D P- P-. O Q PO Q rt (D i' 0 L< - H- X (1 (D 0) (-- I- s n 0 HSH- ~0 r -Q (Dt H-0 0.0 1.0 2.0 L I I Velocity Scale CFt./sec.) I I I I (Nautical Miles) For the numerical model, the bay was divided into 30 increments in the x-direction and 10 increments in the ydirection. The Courant number used for the model was 2.2, the Darcy-Weisbach friction coefficient was 0.05, and the diffusion coefficient was set to zero since the analytical solutions neglected the diffusive terms. Note that the analytical solutions do not include the friction or convective terms. Briggs and Madsen used order-of-magnitude comparison to neglect the friction term so that the governing equations could be solved easier. Friction had to be include in the numerical model so that energy could be dissipated and the model would remain stable. A high friction coefficient was used because of the deepness of the Massachusetts Bay. This caused a discrepancy comparing the numerical results to the analytical solutions. Figure 5.7 shows the maximum ebb velocity vectors developed by the numerical model. By visual comparison the numerical model is in good agreement with the analytical results. Figure 5.8 shows the numerically modeled surface elevations at high tide in the bay. Again, the discrepancies could be due to the friction and convective terms which the numerical model includes, and the analytical solution neglects. 47 Velocity Vectors of Massachusetts Bay 0 10 20 30 40 50 6( Nautical Miles Figure 5.7 Maximum velocity vectors from numerically modeling Massachusetts Bay. L II I t i t tit I I I I t t I I i I % ' K.............................. . . . .. . .. =1 f/sec 48 Surface Elavation (ft) of Massachusetts Bay I I I I 30 25 25- 474.7 ,a15-4. 10- 4.8 4.85 5 w 465 z 5 .5 4.55 0 -5 -10 I I I I I 0 10 20 30 40 50 Nautical Miles Figure 5.8 Nummerically modeled surface elevations at high tide in Massachusetts Bay. CHAPTER 6 MODEL VALIDATION Since the model's efficacy was verified, the numerical model can now be applied to various inlets, but must be validated for each inlet. Validation can be an iterative process; the model's coefficients (friction, diffusion, and wind drag) need to be tuned until the numerical results simulate the actual hydrodynamic process (i.e. circulation and eddies) and the currents of the inlet being modeled. Once the model is validated, any inlet improvement measure, such as jetties, groins, or rip-rap fill, can be confidently simulated. 6.1 Field Data The Coastal and Oceanographic Engineering Department of the University of Florida was commissioned by the City of Punta Gorda, FL to examine the potential for shoreline erosion due to the navigational channel at Ponce de Leon Inlet in the vicinity of Punta Gorda City Park (Wang, et al., 1996). Ponce de Leon Inlet is located on the east side of Charlotte Harbor 50 at the mouth of the Peace River. The channel serves as the public boat ramp at Punta Gorda City Park as well as a residential canal system. The data which was collected and presented in a report to the City of Punta Gorda will be used to show the validation and application (chapter 7) of this numerical model. Data were collected from Ponce de Leon Inlet in November and December 1995 and March 1996. The data necessary to validate the numerical model are: a bathymetric survey; a tide record; current measurements; tidal flow pattern studies; and wind speed and direction measurements. However, other field data may be necessary for coastal structure design purposes. 6.1.1 Bathvmetry The bathymetric data collected should be complete enough such that bottom contours can be drawn accurately. if available, navigational charts may be used in lieu of, or to supplement the field data. No navigational charts were available for this inlet so a bathymetric field survey was conducted. The bottom contours obtained from this data are shown in figure 6.1. The depth for each cell in the numerical grid can be interpolated from these contours. Figure 6. 1 Bathymetric survey of Ponce de Leon Inlet, figure taken from Wang, et al. (1996). 6.1.2 Tide Tide gages were placed at channel markers 8 and 18 (figure 6.2) to measure tide height for the period of data collection, these data are shown in figure 6.3. The time lag of the tide between channel markers 8 and 18 is on the order of 5 to 10 minutes. The maximum difference in amplitude between the two tidal records is approximately 1.5 cm, though most of the tidal peaks show no difference in amplitude. Figure 6.2 Location of stations where tide and current data were collected, figure taken from Wang, et al. (1996). 0.5 0.4 0.3 0.2 . o.1 I * I p -IoIi I z -0.2 -0.3 "" - Tide at Channel Marker 8 S-lide at Channel Marker 18 -0.4- ... Current at inlet entrance ) Oddmeter data 0 2 3 4 5 6 December 1995 Figure 6.3 Continuously observed tidal record and discrete current measurements, figure taken from Wang, et al. (1996). 6.1.3 Current Two instruments were used to collect current data: an electro-magnetic meter, and a vane-type meter. The electromagnetic meter manufactured by Marsh-McBurney, was installed on a docking pile just inside the inlet (station 1 in figure 6.2). The data collected from this meter are also given in figure 6.3. Comparing the current data to the tide data, it may be concluded that the maximum currents occur when the tide is near zero. It is also easily seen that there is slack water when the tide changes either from flood to ebb or ebb to flood. Other current data were collected from an Ottmeter, a hand-held device which uses vanes driven by the flow to measure the current. One minute averages were taken at about 0.5 m below the surface and 0.5 m above the bottom. Figure 6.4 shows the data collected during the flood tide and figure 6.5 shows the ebb velocity data measured at station 1, 2, and 3 (see figure 6.2) in the inlet. Figure 6.3 shows the Ottmeter data taken at the same location as the electro-magnetic meter. The flood velocities between the two different methods are in agreement, though there is a large discrepancy at the ebb tide. This discrepancy can be attributed to the improper installation of Station 1 S0.44 e00.41 Tlme=16:43 S0.21 " 0-- 39 90.28 S0.37 11130/95 Flood current(m/s) 16:55 0 0.42 0.51 0.54 S0.39 11 *. U -2 .0 15 Horizontal Distance (m) Figure 6.4 Flood velocity data taken inside the canal using the Ottmeter, figure taken from Wang, et al. (1996). Time=17:09 S0.25 S0.32 Time=1 4:07 0 0.31 CM022o.26 . S0.27 1113(0/95 Ebb current(m/s) Station 2 Time=1 4:34 00.30 0.30 80.23 Station 3 Time=1 5:02 Horizontal Distance (m) Figure 6.5 Ebb velocity data taken in the canal using the Ottmeter, figure taken from Wang, et al. (1996). Qtiflnn 1 -2 -2 - 0 0.27 w 57 the electro-magnetic current meter. The meter should have been installed in deeper water, instead the meter was installed such that it was less than 0.2 m of water during low tide therefore giving the low, inaccurate measurements. Also, the current meter was about 0.75 m from the pile, there could have been some adverse affects on the current caused by the pile. 6.1.4 Drogues The ebb and flood current patterns were studied using drogues. As the drogues were placed in the water, their motions were video taped. From a video tape the velocities were estimated from using estimated distances and time elapsed on the video tape. When extracting the distances the drogues traveled from relation to channel markers of a known distance a problem is introduced with the perspective of the video camera, because the video camera was at fixed, stationary point on land. This is a crude method to obtain the current velocities; therefore, the values may vary greatly from those obtained through more non-Lagrangian methods. As a result, the quantitative velocity values should only be used as a relative reference between individual drogue path lines. The main purpose of the drogues is to qualitatively study the flow Ponce De Leo n Inlet O.5rm/sec CHARLOTTE HARBOR Figure 6.6 Flood current patterns estimated using drogues, Wang, et al. (1996). IL Figure 6.7 Ebb flow patterns studied using drogues, Wang, et al. (1996). 60 patterns. Figures 6.6 and 6.7 show the flow patterns from the drogue study for the flood and ebb tides. This is one of the more important field studies conducted for validation of the numerical model. 6.1.5 Wind The days that the tide and current measurements were taken, the atmospheric conditions were ideal, there was a high pressure system, a slight breeze, and no synoptic events had occurred in the previous days or were predicted. Hence, wind effects were negligible and not included for validating the numerical model. 6.2 Spectral Analysis Spectral analysis is performed on the real tide, so that the predominate frequencies and amplitudes can be used to simulate the real tide using Fourier series (eq. 4.1). The tide record observed at Punta Gorda (figure 6.3) shows that the tide is mixed diurnal and semi-diurnal. Moreover, if the record was taken for 30 days, a spring and neap tide should be evident. For the model validation, spectral analysis will be Observed Tidal Cycle -0.05 -0.1 -0.15 -0.2 Hours Figure 6.8 One tidal cycle extracted from the tide data. performed on only one period of the recorded tide in figure 6.8. The rational for selecting only one tidal period is: (1) excluding the scattered current measurements at the beginning of he current record (figure 6.3), the current data reflects only one tidal cycle; and (2) since the model is being used as only an aid to inlet morphological improvement measures, 62 residual currents and tidal flushing are not of interest. The continuous tide data taken were digitized to 1 hour increments for a cycle during the time that current measurements were being recorded. Twenty four data points were taken from the observed data record. Figure 6.9 shows the resulting spectrum of the tidal cycle. 6.3 Grid The grid chosen for the Ponce de Leon Inlet as shown in figure 6.10. The grid size Ax and Ay, is 6 m x 6 m, which results in 59 rows and 29 columns. This resolution will be sufficient to allow jetties to be inserted and show the currents close to the shore. A higher resolution grid could be used, but the computation time required will increase. The information gained by using a higher resolution grid does not justify the excess time needed to run the numerical model. Both the northern and southern boundaries were extended, so that the eddies which were studied using the drogues (figure 6.7) do not cause a numerical instability at the open boundaries. The drogue studies also show that the tidal flow 0) E 150 100 50 01 0 50 100 150 meters Figure 6. 10 Area of Ponce de Leon Inlet which a grid will be applied to. becomes one-dimensional in the canal; henceforth, the grid chosen should allow the model to be validated against the data collected in the field studies. Using the stability criteria discussed in chapter 3.4, the time step and eddy viscosity coefficient can be determined. The Courant number used was 10, which resulted in 64 a time step of 18 sec and an eddy viscosity coefficient of 0.12. The Darcy-Weisbach friction coefficient was set to 0.01, and, again, wind was not used as a parameter for this model, since field studies were conducted on a calm day. 6.4 0pen Boundary Conditions The open boundary conditions (OBCs) needed to drive the model for the case of Punta Gorda are tidal, Neumann, and canal open boundaries (OBs). The tidal OBs will be placed on the western edge of the grid, the Neumann OBs will be placed on the northern and southern edges of the grid which extend into Charlotte Harbor, and canal OBs will be placed on the cells on the eastern edge of the grid. The following sections will discuss the parameters used for these OBCs. 6.4.1 Tidal OBC The amplitudes, circular frequencies, and phase shifts for the most dominant harmonics, determined by spectral analysis, are shown in table 6.1. These are the parameters which will be used to drive the numerical model at the open boundary using equation 4.1. The simulated tidal signal is n-th Amplitude, Ai Circular Phase Angle, i Harmonic (M) Frequency, a1 (rad) (rad/sec) 1 0.0847 7.2722 x 10-5 0.8178 2 0.1654 1.4544 x 10-4 2.7241 3 0.0119 2.1817 x 10-4 0.9671 6 0.0079 4.3633 x 10-4 2.5830 Table 6.1 Parameters used for the tidal OBC. compared against the measured signal from figure 6.8 in figure 6.11. 6.4.2 Neumann OBC Since the grid extends into Charlotte Harbor where the tidal OBC is, the appropriate measures discussed in chapter 4.1 must be applied here. In comparing the area the numerical model's grid will cover (figure 6.10) to the bathymetric data (figure 6.1), it is shown that the boundaries were extended roughly 60 m to the north and 60 m to the south. This adds a buffer zone so that any circulation in the area of interest does not cause an instability in the OBC and so that any effects caused by the Neumann OBC do not adversely affect the flow patterns in the area of interest. During the model validation or its application, if the eddies occur near the Simulated Tidal Cycle Hours Figure 6.11 Simulated tidal OBC signal compared to the observed tidal signal in fig. 6.8. Neumann OB cells, the model's grid must be extended farther to the north and south. 6.4.3 Canal OBC The parameters necessary for this OBC are depth at the end of the canal, the Darcy-Weisbach friction coefficient, the size of the increment (AXcana), and the number of increments (Icanal) These parameters are somewhat empirical. The 67 friction coefficient may be greater in the canal than in the inlet due to the different bottom sediments characteristics, excess submerged vegetation (e.g. sea grass) or overgrown foliage on the banks (e.g. shrubs and mangroves) In the interest of computational time, as previously discussed, the total number of increments should not exceed 15. The measured or estimated length of the canal can be divided evenly into a number of increments. The right side of figure 6.2 shows the beginnings of a complex finger canal system. This may add damping, or friction, to the long wave, so the friction coefficient should be increased. The finger canal system also makes it very difficult to obtain an "effective" length if the canal were straight. For the Punta Gorda Isles finger canal system, the surface area of the system was estimated from maps. Since there are five canal OBCs used on the numerical model and 10 increments were used for each OBC, the surface area was divided by 50 to obtain the size of increment, Ax,,,,,,,, which was 1500 km. The friction factor used in the numerical model was 0.05. The canal system was assumed have uniform depth, hence the depth at the end of the canal was set to the depth at the OBC, or beginning of the canal. E 8O , I VI I -~ ~ -0. -mlse----c ----20 0 20 40 60 80 100 120 140 160 180 meters Figure 6.12 Maximum flood velocity vectors for numerical model validation. 6.5 Results The numerical model was ran for 5 tidal cycles, a sufficient time for the model to reach steady state from the initial conditions. Figure 6.12 shows the maximum flood velocities and figure 6.13 shows that maximum ebb velocities. The maximum velocities, both flood and ebb, are about 0.3 m/s, this value is a little low compared to the electro-magnetic current meter readings of about 0.4 m/s and the estimated 1100 80 . . . . . . 60 40 20 0II I I -20 0 20 40 60 Figure 6.13 Maximum ebb velocity vec model validation. drogue velocities of 0.5 m/s. 80 100 120 140 160 180 meters tors showing numerical The difference in velocities may be as a result of the canal OBC. One could adjust the parameters used in the canal OBC to produce better results, though that could take many iterations. Comparing figures 6.13 and 6.7, the numerical model does not show the circulation patterns that the drogue studies show. The eddy viscosity coefficient was increased to the point that the stability criteria (eq 3.44) was slightly violated, though no instabilities became apparent. A possible 70 causes for this difference are that the bathymetric data may have errors. Another possible cause is that the bathymetric relief near the channel may be too great for the model's resolution. This will be addressed later in chapter 8. Though a more likely reason that the numerical model was not able to produce the circulation is because the circulation may be initiated by superimposing the wave conditions on the ebb tidal flow. The predominant waves were found to come from the north which could cause a circulation pattern on the northern side of the channel during the ebb cycle. Also, the Peace River empties into Charlotte Harbor just north of Punta Gorda. The current from the river, a north to south current, or an ebb current throughout Charlotte Harbor can initiate a circulation pattern which was observed. CHAPTER 7 MODEL APPLICATION Recall, the objective of this research is to develop a numerical model to aid in the placement of coastal structures as an inlet improvement. When applying the numerical model, any differences or inconsistencies between the field data and the model validation results must be kept in mind. The effects of these discrepancies on the numerical output must be realized when scrutinizing the inlet improvement measures, such as where a jetty or groin is placed and its length. Ponce de Leon inlet in Punta Gorda, FL will be used to show how the numerical model can be used. Two cases will be demonstrated where a single jetty is installed and a pair of jetties are installed. The results presented in this section are a "zoomed in" portion of the model's grid. This allows the currents in the area of interest to be shown in better detail. E 80 . ................ .. o 40 .1 I I *~ ~ ~ .. .. . . -20 0 20 40 60 80 100 120 140 meters Figure 7.1 Flood velocity vectors of a single jetty to Punta Gorda, FL. 160 180 applied 7.1 Single Jetty Figures 7.1 and 7.2 show the maximum flood and ebb velocities, respectively, for a single jetty case. From comparing figures 7.1 and 6.12, the ebb currents become greater along the southern bank, the bank with no jetty. The jetty also causes return flow during the ebb current on the southern shoreline, though the velocities may not be large enough to initiate shoreline erosion. The improvement measures recommended to the City of Punta Gorda were to instal E 8060 4020 0 -20 0 Figure 7.2 Ebb velocity Punta Gorda, 80 100 120 140 160 180 meters a single jetty applied to vectors of FL. the single jetty as shown and to add rip-rap fill to the southern shoreline, Wang et al. (1996). 7.2 Pair of Jetties Figures 7.3 and 7.4 show the maximum flood and ebb velocities, respectively, for the case with a pair of jetties. The pair of jetties virtually eliminates the erosional currents in the proximity of the shoreline. Again, the discrepancies between -20 0 20 40 60 80 100 120 140 meters Figure 7.3 Flood velocity vectors of a pair applied to Punta Gorda, FL. of jetties the field data and the validation results, namely the drogue study and circulations eddies, must be kept in mind. It is possible that the pair of jetties may cause significant eddies outside of the channel near the shoreline which may result in higher, erosional velocities. iA. 160 140 100- . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . - - - - - - ------------------------------------------------------------------------------- ----------------------.......... ................ ............ ............ . . . . . . . . -- - - ............ 0 ., -.0.3m/sec ................ ..... ......... '. 7q, 160 180 l U I I I I I I I 160 140 120 & 100 E 80 . 0 6 0 .. .. 40 20 - 0.3 m/sec -20 0 20 40 60 80 100 120 140 160 180 meters Figure 7.4 Ebb velocity vectors of a pair of jetties applied to Punta Gorda, FL. CHAPTER 8 CONCLUSIONS AND RECOMMENDATIONS A two-dimensional numerical model was developed using a factorization method for numerical solution. The model developed accounts for constant friction, convective acceleration, turbulent diffusion, and wind stress. The numerical model can be adopted to a variety of inlet-river and -canal systems to study the tidal currents. This model can be applied to coastal engineering design projects involving the installation of jetties, breakwaters, groins, or rip-rap fill where it is necessary to study the feasibility of such inlet improvement measures. 8.1 Conclusions The numerical model's efficacy was verified reasonably well. There is a discrepancy when comparing the numerical model to the analytical solutions of Massachusetts Bay. A possible cause of the difference could be the friction factor. The analytical solutions did not account for friction, where as if friction is not included in the numerical model, then 77 the model may become unstable or produce erroneous results. The analytical solutions also neglected the inertial terms. Perhaps, if the convective terms were set to zero in the numerical model, then the results might change. Overall, the numerical model that was developed can be applied to a multitude of tidal inlets and areas. The model was shown that it could produce reasonable results. The numerical model developed was not able to produce ebb tidal patterns as measured in the field when applied to Ponce de Leon inlet at Punta Gorda, FL. Possible reasons for the differences between the numerical model and the field data could be bathynunetry and the canal OBC. Most navigable inlets consist of a deep channel dredged across a relatively shallow area. When numerically modeling these inlets, the model may produce better results if the channel is "numerically smoothed". Figure 8.1a shows the channel at Punta Gorda and figure 8.1b how one might "smooth" the relief of the bathymmetry. The sudden relief in bathyinmetry may cause a large gradient in the friction and convective terms with profound affects on the numerical model. 78 a 0 -0.5 -1.5 -2 40 -30 -20 -10 0 10 20 30 40 b 0 -0.5 -1.5 -2i 40 -30 -20 -10 0 10 20 30 40 Distance (m) Figures 8.1 a & b (a) shows the original bathymetry of channel and (b) shows the "smoothed" bathymmetry of channel. 8.2 Recommendations The following studies and improvements to the numerical model are recommended: l)The affects of "numerically smoothing" a channel should be closely studied. From testing the numerical model as it was being developed, it is reasonable that a "smoothed" channel will produce 79 more realistic results. Information was not presented in this thesis on the grounds that no analytical solutions or data were available at the time to support this hypothesis. 2)Develop anther OBC which allows for a lagoon connected by a small canal as shown in figure 8.2. The governing equation presented by Bruun (1978), but developed by earlier researchers, for this case are as follows: d2 'B F AB dB d'IB gAc gA( + + p dt 2 2Lc Ac di dt LCAB LCAB where 'IB=surface elevation of bay, t=time, F=frictional losses, Lc=length of channel, AB=surface area of bay, Ac=cross sectional area of channel, g=gravity, and To=surface elevation at entrance. aAREA- A- .. BAY FRESH PMM. c'--TH.h WATER OCEAN OLUME A' ELEVATION Lc SURFACE AREA AO ELE TIN 70 ELEVATION Figure 8.2 An idealized channel/bay system. Figure was taken from Bruun (1978). This lagoon/canal system can be applied to the numerical model as an OBC the same as the canal and river are discussed in chapter 4.1, except that the ocean elevation shown in figure 8.2 can be the surface elevation can be used to drive the model as at an OB. 3)Morphological changes and sediment transport should be added to the numerical model. This will allow the long term sediment transport, shoaling, and bank erosion/accretion to be evaluated for a given inlet improvement measure. LIST OF REFERENCES Arcilla, A.S. "An Integrated Numerical Approach for Coastal Engineering Problems." Journal of Coastal Research. Vol. 5, 1989: 603-616. Bathan, Karl H. "A Tidal and Wind Driven Numerical Model of a Lagoon Before and After a Major Modification." J.K.K. Look Laboratory of Oceanographic Engineering, University of Hawaii, Vol. 6:2, 1976. Briggs, Douglas A. and Madsen, Ole S. "Analytical Models for One- and Two- Layer Systems in Rectangular Basins." Report No. 172, School of Engineering, Massachusetts Institue of Technology, Cambridge, MA, 1973. Bruun, Per. Stability of Tidal Inlets. Elsevier Scientific Publishing Complany, Amsterdam, 1978. Cheng, Ralph T. "Modeling of Tidal Residual Circulation in San Francisco Bay, California." Two-Dimensional Flow Modeling. U.S. Army Corps of Engineers, The Hydrographic Engineering Center, Davis, California, 1982. Choi, Jei-Kook. "Three-Dimensional Curvilinear-Grid Modeling of Baroclinic Circulation and Mixing in a Partially Mixed Estuary." Dissertation, University of Florida, Gainesville, FL, 1992. Henderson, F. M. Open Channel Flow. Macmillan Publishing Co., Inc., New York, 1966. Hirsch, Charles. Numerical Computation of Internal and External Flows. Volume 1, John Wiley & Sons Ltd., Chichester, 1994. Pedlosky, Joseph. Geophysical Fluid Dynamics. SpringerVerlag, New York, 1979. Reeve, D. E. and Hiley, R. A. "Numerical Prediction of Tidal Flow in Shallow Water." Computer Modelling of Seas and Coastal Regions. Elsevier Science Publishers Ltd., New York, 1992: 197-208. Reidel, H. P. and Wilkinson, F. L. "Numerical Modeling An Aid to Assessing Field Data." Coastal Engineering. 1976: 3243-3262. Roll, H. U. Physics of the Marine Atmosphere. Academic Press, New York, 1965. Schomaker, Christine W. "A Model for Tidal Circulation Adapted to Monterey Bay, California." Thesis, Naval Postgraduate School, Monterey, California, 1983. Sheng, Y. P. Mathematical Modeling of Three-Dimensional Coastal Currents and Sediment Dispersion: Model Development and Application. Technical Report CERC-832, Aeronautical Research Associates of Princeton, Princeton, New Jersey, 1983. Wang, Hsiang; Lin, Lihwa; and Wang Xu. "Ponce de Leon Channel Improvement Studies, Punta Gorda, Florida." UFL/COEL-96/004, Coastal and Oceanographic Engineering Department, University of Florida, Gainesville, Florida, 1996. Weiyan, Tan. Shallow Water Hydrodynamics. Elsevier Science Publishing Company, Inc., New York, 1992. Wilson, Basil W. "Note on Surface Wind Stress over Water at Low and High Wind Speeds." Journal of Geophysical Research. Vol. 65, No. 10, 1960: 3377-3382. BIOGRAPHICAL SKETCH Peter N. Seidle was born in Binghamton, NY on 28 January, 1972. His family moved to Maryland in suburban Washington DC where he started kindergarten and graduated from Thomas S. Wootton High School, Rockville, MD, in 1990. In 1994 he graduated from Widener University, Chester, PA, with a Bachelor of Science in Mechanical Engineering. Shortly after graduating from college, he became employed by the ARIST Corporation to work as a program analyst in the Ballistic Missile Defense Organization, which is part of the Department of Defense. In August 1995 he began his graduate studies at the University of Florida in the Coastal and Oceanographic Engineering Department. In May 1997 he graduated with a Master of Science in the field of study he has and will always enjoy and be fascinated with. |