NUMERICAL SIMULATION OF CURRENTS MONTEREY BAY by Roland Albert Garcia DUDLEY KNOX LIBRARY KAVAL POSTGRADUATE SCHOOL BONTEREY. CALIFORNIA 93943 United ■--.-.. ^ avai Postgraduate School THESIS NUMERICAL SIMULATION OF CURRENTS IN MONTEREY BAY by Roland Albert Garcia Th ssis Advisor E. B. Thornton March 1971 Approved {^on. puhtlc IctzaAe.; distribution witimltzd. T .138224 Numerical Simulation of Currents in Monterey Bay by Roland Albert^jGarcia Lieutenant, United S rates Navy B.A., Temple University, 1966 Submitted in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE IN OCEANOGRAPHY from the NAVAL POSTGRADUATE SCHOOL March 1971 i 'OSTGEfAIXJATF '^rmnm ABSTRACT Recent interest in pollution control and the proximity of Monterey Bay to the Naval Postgraduate School prompted an investigation of the circulation in the bay. The first phase of the study consists of solving the simple cavity flow problem. A vorticity-stream function relationship is solved using an explicit, time dependent, finite dif- ference scheme. Solutions for selected Reynolds' numbers and length to width ratios of the cavity are obtained. Values are chosen to give an indication of the flow patterns occurring over a wide range of these parameters. Equations for a refined model are derived to include the effects of the bottom topography, frictional forces and the Coriolis force. A numerical procedure similar to the one applied to the simple cavity flow problem is used on the refined equations. The topography of Monterey Bay is used in this study. Results indicate that closed circulation patterns are more probable in Monterey Bay than in other embayments due to the presence of the submarine canyon. TABLE OF CONTENTS I. INTRODUCTION ] 5 A. BACKGROUND . . . ] 5 B. NUMERICAL MODELING ]9 C. SCOPE OF PRESENT PAPER 20 II. NUMERICAL METHODS 22 A. INTRODUCTION 22 B. NUMERICAL RELATIONSHIPS :-2 1. Taylor's Theorem 22 2. Representation of Derivatives 23 3. Representation of Arakawa Jacobian 25 C. RELAXATION METHODS 27 1. Residuals 27 2. Relaxation 28 3. Overrelaxation 29 4. Optimum Overrelaxation Factor, Ropt 29 a. Introduction 2" b. Theoretically Determined Value 2° c. Empirically Determined Value 31 D. CONSERVATIVE AND TRANSPORTATIVE PROPERTIES 34 E. CONVERGENCE AND STABILITY 38 III. GOVERNING EQUATIONS 40 A. BASIC EQUATIONS OF MOTION 40 1. Horizontal Momentum Equations 40 2. Continuity Equation 41 3. Integrated Equations of Motion 42 B. DEVELOPMENT OF VORTICITY EQUATION A3 1. Derivation 43 2. Nondimensionalization 45 3. Nondimensional Parameters 50 C. DEVELOPMENT OF STREAM FUNCTION EQUATION 52 1. Derivation 52 2. Nondimensionalization 53 IV. SIMPLE CAVITY FLOW 54 A. GOVERNING EQUATIONS 54 B. NUMERICAL APPROACH '. . 56 1 . Upwind Differencing for Solving the Vorticity Equation 56 2. Method for Solving Vorticity Equation Employing the Arakawa Jacobian 57 3. Optimum Overrelaxation to Solve Poisson's Equation for the Stream Function 58 4. Boundary Conditions 59 a. Inflow 60 (1) Velocity Profile 60 (2) Stream Function Profile 63 (3) Vorticity 66 b. Walls 67 c. Downstream Continuation 68 d. Lid 68 e. Corners 69 5. Convergence and Stability Criteria 70 6. Computational Sequence 71 7. Programming Optimization 72 C. RESULTS 73 D. CONCLUSIONS 78 V. REFINED MODEL 79 A. GOVERNING EQUATIONS 79 1. Bottom Friction Equations 79 ?.. Final Form for the Vorticity and Stream Function ... 81 3. Convergence and Stability 85 B. INVESTIGATION OF THE RELATIVE IMPORTANCE OF THE CORIOLIS FORCE AND BOTTOM FRICTION 86 1. Coriolis Force 86 2. Bottom Friction 89 C. NUMERICAL APPROACH 90 1. Method for Solving the Vorticity Equation 90 2. Optimum Overrelaxation for Solving the Stream Function Equation 90 3. Boundary Conditions 91 a. Inflow 91 b. Lateral Boundaries 5 3 c. Downstream Continuation 95 d. Lid 95 e. Bottom Topography 96 4. Computational Sequence 98 D. RESULTS 98 1. Neglecting the Coriolis Force and Bottom Friction . .104 2. Neglecting Bottom Friction 10/ E. CONCLUSIONS H5 APPENDIX A - Development of the Bottom Friction 117 APPENDIX B - Solutions for the Simple Cavity Flow Problem 123 REFERENCES 1*8 INITIAL DISTRIBUTION LIST 149 FORM DD 1473 153 LIST OF FIGURES Figure Page 1.1 Location of Monterey Bay 17 1.2 Monterey Bay 18 2.1 Numerical grid 31 2.2 Plot of overrelaxation factor vs. number of passes required for convergence for Re = 1 , 10, 50, and 100 and aspect ratio of 1 35 2.3 Plot of overrelaxation factor vs. number of passes required for convergence for Re = 1 , 10, 50, and 100 and aspect ratio of 2 36 2.4 Plot of overrelaxation factor vs. number of passes required for convergence for Re = 1, 10, 50, and 100 and aspect ratio of 3 37 4.1 Specification of boundaries 60 4.2 Velocity boundary layer profiles 61 4.3 Unsuitable boundary layer prof Lie 64 4.4 Stream function boundary layer profiles 65 4.5 Portion of numerical grid plotted 74 4.6 Comparison of four solutions for the stream function 76 4.7 Sample vorticity plot, Re = 10 , B = 2 77 5.1 Frictional coefficient, Rl 83 5.2 Frictional coefficient, R2 84 5.3 Grid for refined model 94 5.4 Hypothetical bottom topography 9 7 5.5 Bottom topography simulating Monterey Bay 99 5.6 Bottom topography with no smoothing passes 100 5.7 Bottom topography with one smoothing pass 101 5.8 Bottom topography with two smoothing passes 102 5.9 Bottom topography with five smoothing passes . . . .10 3 5.10 Refined model, hypothetical bottom topography, volume transport stream function, Re = .01-100, eF = 0 105 5.11 Refined model, topography simulating Monterey Bay, volume transport stream function, Re = .01-100 eF = 0 106 5.12 Simple cavity flow, stream function, Re = .01, B = 2 . 103 5.13 Refined model, hypothetical bottom topography, volume transport stream function, Re = .01-1.0, eF = .23 110 5.14 Refined model, hypothetical bottom topography, volume transport stream function, Re = .01-1.0 > eF = 23.0 Ill 5.15 Refined model, bottom topography simulating Monterey Bay, volume transport stream functions Re = .01, eF = .23 x 10-2 112 5.16 Refined model, bottom topography simulating Monterey Bay, volume transport stream function, Re = .01, eF = .23 x 102 113 5.17 Refined model, bottom topography simulating Monterey Bay, volume transport stream function, Re = .01-1.0, eF = 1.0 x 10-3 114 B.l Simple cavity flow, stream function, Re = 1, B = 1 124 B.2 Simple cavity flow, vorticity, Re = 1, B = 1 .... 125 B.3 Simple cavity flow, stream function, Re=l,B = 2 . 126 B.4 Simple cavity flow, vorticity, Re = 1 , B = 2 .... 127 B.5 Simple cavity flow, stream function, Re = 1, B = 3 128 B.6 Simple cavity flow, vorticity, Re = 1, B = 3 .... 129 B.7 Simple cavity flow, stream function, Re = 10, B = 1 130 B.8 Simple cavity flow, vorticity, Re = 10 , B = 1 .... 131 B.9 Simple cavity flow, stream function, Re = 10, B = 2 132 B.10 Simple cavity flow, vorticity, Re = 10 , B = 2 .... 133 B.ll Simple cavity flow, stream function, Re = 10, B = 3 ". ... 134 B.12 Simple cavity flow, vorticity, Re = 10 , B = 3 .... 13!) B.13 Simple cavity flow, stream function, Re = 100, B = 1 136 B.14 Simple cavity flow, vorticity, Re = 100, B = 1 ... 137 B.15 Simple cavity flow, stream function, Re = 100, B = 2 133 B.16 Simple cavity flow, vorticity, Re = 100, B = 2 ... 139 B.17 Simple cavity flow, stream function, Re = 100, B = 3 1*0 B.18 Simple cavity flow, vorticity, Re = 100, B = 3 ... 1*1 B.19 Simple cavity flow, stream function, Re = 1000, B = 1 142 B.20 Simple cavity flow, vorticity, Re = 1000, B = 1 ... 1*3 B.21 Simple cavity flow, stream function, Re = 1000, . . , 144 B = 2 B.22 Simple cavity flow, vorticity, Re = 1000, B = 2 ... 145 B.23 Simple cavity flow, stream function, Re = 1000, B = 3 l^6 B.24 Simple cavity flow, vorticity, Re = 1000, B = 3 ... 147 LIST OF TABLES Table Page I Tabulation of the rate of convergence of Ropt .... 33 II Frictional coefficients 82 10 TABLE OF SYMBOLS AND ABBREVIATIONS a Argument in Ekman solution for slope currents A^ Coefficient of horizontal eddy diffusivity A^. Coefficient of vertical eddy diffusivity B Aspect ratio, AX/AY D Characteristic depth e 2.7182818 f Coriolis parameter, 2fl sin 9 F Nondimensional coriolis parameter g Acceleration due to gravity h Depth H Nondimensional depth i /=! J Jacob ian In Natural logarithm log Common logarithm L Characteristic length 0 Order of magnitude of P Pressure r Radius of earth Rl, R2 Frictional coefficients U L oo Re Reynolds number, - Ropt Optimum overrelaxation factor t Time T Nondimensional time 11 Tr Transport u,v x,y components of velocity U,V x,y components of nondimersional velocity U Characteristic velocity CO J w u + iv x,y,z Space coordinates X,Y Nondimensional space coordinates Z Nondimensional vorticity a /ill 3 Beta plane approximation to coriolis force t Time increment t Critical time increment c Ax, Ay Grid spacing V Del operator V Nondimensional del operator V Laplacian 2 V Nondimensional Laplacian 8 Partial differentiation BL2 e Nondimensional parameter, — — CO £ Vorticity K Von Karman's constant ri Height of sea surface 9 Latitute o tt 3.1415926 p Density bx by j. , x , T x,y components of bottom stress t , f x,y components of wind stress 12 t Turbulent shearing stress \\> Stream function Y Nondimensional stream function ft Angular velocity of the earth Subscripts i,j x,y grid positions x,y Partial differentiation with respect to x,y w Wall value Superscripts k Iteration level n Time step level 13 ACKNOWLEDGEMENTS Dr. Edward B. Thornton proposed the topic for this thesis and guided me during my investigation. I am especially indebted to him for his assistance, encouragement and professional advice. Several other members of the Naval Postgraduate School staff were helpful in the development of this study. Dr. Jerry A. Gait assisted me in the physical interpretation of the problem and also suggested the approach used to incorporate the bottom friction into the refined model. Dr. Henry Crew and Dr. Frank D. Faulkner aided me in the pro- gramming and mathematical portions of this thesis. This study would have been impossible without the facilities and cooperation of the staff at the William R. Church Computer Center, which is part of the Naval Postgraduate School. All programming was done for the IBM 360/67 Computer, which is located at this computer center. 14 I . INTRODUCTION A. BACKGROUND Monterey Bay, located on the West Coast of the United States approximately one hundred and twenty miles south of San Francisco, presents a good opportunity to study currents within an embayment (Figure 1.1). This bay is well suited for circulation studies because of the following geographic and oceanographic features. The bay is a large one and, to a close approximation, it has symmetrical boundaries. It is semi-elliptical in shape with a major axis of twenty miles and a width of about eight miles (Figure 1.2). Bisecting the bay is a large submarine canyon. This canyon originates in the bay offshore of Moss Landing. At the mouth of the bay, eight miles from Moss Landing, the canyon attains a depth of six-thousand feet. This topographic feature is as large as the Grand Canyon in vertical relief. The effect of the Monterey Submarine Canyon on the circulation within the bay is one of the problems which is investigated in this study. A desirable oceanographic feature is the presence of currents which flow past the bay. The Davidson Current flows north along the West Coast of North America from November to February. Flowing southward during the spring and summer is the California Current. It is assumed in this study that oceanic currents, such as these, are the driving force for the circulation within an embayment. The currents drive the circulation by advecting momentum into the area in question. The effects of tidal forces and local winds are neglected. It is not meant to be implied that these forces are unimportant, but that only a single component of the total system is being considered in this analysis. 15 Monterey Bay is conveniently located for conducting circulation studies for several reasons. First, the Naval Postgraduate School is located on the shores of the bay in the city of Monterey. Both the Postgraduate School and the United States Navy are interested in circu- lation studies of this kind. The facilities at this institution are well adapted to carry out such studies. The Department of Oceanography at the Naval Postgraduate School is interested in circulation studies and has access to oceanographic research vessels. Also, located at the School is a large computer center with a well-trained staff and an IBM 360/67 Computer along with all the necessary peripheral equipment. Secondly, there are many communities located on the shores of Monterey Bay which are becoming increasingly aware of pollution problems. An important question that needs to be answered before pollution can be controlled is: what is the effect of currents on the discharged effluent? To answer this question the circulation regime of the area involved must be determined. It is hoped that the results of this study will provide some answers to the above question. Prior to this study, little had been done to ascertain the circu- lation patterns within Monterey Bay. This thesis is part of a larger research effort being conducted by the Naval Postgraduate School aimed at determining the forces which dominate in controlling current patterns within an embayment and the types of patterns which are generated. Pre- sently, there are three phases to the project. Two phases concern numerical modeling. One numerical model, the subject of this thesis, assumes oceanic currents drive the circulation within a bay. A second model simulates the tide and wind induced circulation. The third phase involves field studies conducted in order to assess the models which have been developed. 16 ...„ __ Figure 1.1 Location of Monterey Bay 17 Figure 1.2 Monterey Bay 18 B. NUMERICAL MODELING Recently, with new developments in partial differential equation theory and finite difference methods, the field of fluid mechanics has progressed tremendously. The advent of large-scale computers has spurred these developments and made it possible to solve partial dif- ferential equations numerically which have proven, so far, to be impossible to solve analytically. Due to the great speed of these computers it has become possible to solve large systems of simultaneous equations rapidly enough to put numerical prediction models in fields such as Meteorology and Oceanography on a real time basis. Most circulation problems, basically, involve the solution of a large number of simultaneous equations with specified boundary con- ditions. Practically all studies of this type are concerned with either small-scale flows with simple geometries such as: channel flow, flow around a cylinder, and simple cavity flow or with large-scale oceanic circulation studies. Little has been done to try to apply these equations and numerical techniques to medium-scale flow such as the circulation within an embayment. There are several possible reasons why more medium- scale flow problems have not been investigated. One reason is the high cost of studies of this type. Computer time is very expensive and not always readily available. Medium-scale problems are in general highly variable and each case must be dealt with individually. Another reason is that, before the increased awareness of pollution problems, there was little interest in the regional circulation beyond that caused by waves and tides. Lastly, the application of the equations of motion to large or medium-scale flow problems requires assumptions about the forces and terms involved which have not yet been universally accepted as being valid, 19 There are many numerical procedures that can be used to try to solve a system of simultaneous partial differential equations. An interesting result of applying different numerical techniques to the same problem is the possibility of obtaining different solutions. One of the problems of numerical modeling is determining whether the correct solution has been reached, if only one exists. There has been a lot of research into the problem of the solution of partial differential equations by Ames [Ref. 2], Considering the Navier-Stokes Equations, he proposes that a unique solution exists only below a certain Reynolds' number and that above this Reynolds' number several solutions exist; whereas, above a larger' critical number no solutions exist. The critical Reynolds' number marks the transition to turbulent flow. To illustrate the existence of multiple solutions to a partial differential equation, he gives an example of a quasi-linear, elliptical partial differential equation for which there is no unique solution. An extensive study of various numerical procedures used in the solution of compressible and incompressible laminar separated flow for simple geometries was made by Roach [Ref. 7]. His dissertation is quite thorough in its consideration of the problems encountered and the techniques used to resolve them. Some of Roach's techniques are extended in this study to the problem of motion within an embayment in the simple cavity flow section of this thesis. C. SCOPE OF PRESENT PAPER This thesis attempts to take the equations of motion and apply them to the problem of the circulation in an embayment and solve them numerically. A progression is made from the simple cavity flow problem to the refined model. The simple cavity flow problem considers: local 20 rate of change of vorticity with time, advective terms, and lateral shear stresses. In addition, the refined model considers: planetary vorticity tendency, topographic vorticity tendency and bottom friction. Procedures for numerically solving these equations are examined while economy of computer time and the validity of the solution are kept in mind. An attempt is made to generalize the computer program so that it can be used on any embayment. The generality is obtained by allowing the input of an arbitrary bottom topography. The model is tested on Monterey Bay in order to discover, not only the circulation patterns resulting from the interaction of the forces involved, but also the relative importance of the forces. To determine the importance of the forces involved, various computer runs are made with different combinations of the forces acting, and the effect on the circulation patterns is examined. It is hoped that the results of these investigations will shed some light on the processes which control circu- lation patterns and help provide answers to some of the problems which now exist. 21 II. NUMERICAL METHODS A. INTRODUCTION This section consists of a brief discussion of some of the aspects of numerical modeling which are used in this thesis. There are many excellent texts which contain thorough discussions of this material; therefore, the material is not covered in depth. All of the repre- sentations derived in the following section, II-B, will be used in the development of the numerical equations solved in later sections. The following paragraphs introduce these concepts and develop the necessary relationships . ♦ B. NUMERICAL RELATIONSHIPS To represent a partial differential equation in numerical form the starting point is usually Taylor's Theorem. From this theorem, the necessary expressions which represent derivatives in their various forms can be derived. Numerical relationships for terms such as the Laplacian and the Jacobian can then be developed from these expressions. Basically, these expressions determine the value of a derivative at a given point in terms of the values at surrounding points and the grid spacing. 1. Taylor's Theorem The finite difference representation for the derivatives of a function at a given point in terms of the values of the function at surrounding points can be derived from Taylor's Theorem. Taylor's Theorem in Two Dimensional Space: If u and its derivatives of order _< p+1 are single valued, finite and continuous at every point on the line segment 22 from (x,y) to (x+Ax,y+Ay) then there is a point (a,b) on that line segment such that P 3N+M u(x+Ax,y+Ay) = u(x,y) + E ^77 '"f^^l (Ax)N (Ay)M + R (2.1) (N+M)=l ' * 3x3 y where the remainder, R, is given by 9N+M t. ! v u(a,b) - .N ,. .M R = rar N+M +,~J^ ( } ( y) N+M=p+1 3x dy It can be seen that if Ax and Ay < 1 as N and M ■*■ « then R -*■ 0 ; and in the limit, the expansion equals the function. If y is held constant (i.e., Ay = 0) there is a contribution to the summation terms only when M = 0, therefore P N u(x+Ax,y) = u(x,y) + £ |t ^~ (Ax)N + 0 (xP+1) (2.2) N=l * 3x If x is held constant, (Ax = 0), similarly P M u(x,y+Ay) = u(x,y) + I jfr ^ (Ay)M + 0 (yP+1) (2.3) M=l 3y Using these equations, expressions for derivatives at a point are developed in Section II-B-2. 2. Representation of Derivatives The following grid illustrates, diagrammatical ly , the notation that is used to fix the grid position of a given point. 23 x-Ax,y+Ay x,y+Ay x+Ax , y+Ay Ax Ui-l,j+l x-Ax,y Ay l.j+l x,y i+l,j+l x,y+Ay u. 1-l.J x-Ax,y~Ay u x,y-Ay um,j x+Ax,y-Ay i-l.j-l Ui,j"l 1+l.J-l Using Taylor's Theorem with y held constant Equation 2.2 in the new notation becomes 3u, _LLL a2 0 u. . u;x1 , ■ u. . + — ^- Ax + , l+l, j i,j 3x 2 3x2 1 1,1 - .2 . . ,. .3 (Ax) + 0 (Ax) If - Ax is substituted for Ax the equation becomes 9u. 32u u. - . = u. . - — ±-^ Ax + ± i^l (Ax)Z + 0 (Ax) J l-l, j i,j 3x 2 3x2 Linear combinations of 2.4 and 2.5 with the higher order terms neglected yields Representation Mag, of Error Diff. Scheme a u .- . - u. . 3x Ax u. . - u . ax Ax 0 (Ax) 0 (Ax) Forward Backward (2.4) (2.5) » u, n . - u . n 3u = i+l, ,1 i-l,,1 9x 2Ax .2 u... . + u. 1 . 3 u _ i+l,, 1 i-l,, 1 - 2u, 3x (Ax) 0 (Ax) i*i 0 (Ax): Central 24 Similarly, derivatives with respect to Y can be found. Representation Mag, of Error Diff. Scheme a u- -.1.1 " u- • 3y Ay r. U. . " U. Ay 0 (Ay) Forward 0 (Ay) Backward 3y 2Ay 0 (Ay) Central .2 u. , + u. . _ - 2u. . A-" = -Uii ^J"1 i-J. n M,^2 2 2 3y (Ay)Z 0 (Ay)' A representation for the Laplacian, V u, can be developed from the above expressions. 2 2 „2 3 u . 3 u V u = — 7? + 2 2 3x 3y o u.,n , + u. . . - 2u. . u. . ,. + u. . - - 2u. y2u = 1+1>J 1~1>.1 JiX + 1?Jtl 1>.l-1 LA (2,6) (Ax)2 (Ay)2 3. Representation of the Arakawa Jacobian This study utilizes Arakawa' s technique for representing the Jacobian which is an average of three different expressions for the Jacobian [Ref. 3]. The method Arakawa developed is used for two reasons. First, it has the desirable property of conserving certain quantities two of which are the vorticity and the vorticity squared. A further discussion of the conservative property is contained in Section II-D. Secondly, the representation is a very stable one which helps to minimize convergence problems in solving the numerical equations, The price that is paid for this stability is the high degree of smoothing 25 inherent in this representation which may eliminate some small-scale features. Although these features may be real, it is not felt that their elimination is very important because this study is only inter- ested in the general flow patterns. Examination of small-scale features in the flow patterns with the simplifying assumptions which are made in this study would not be realistic. In order to develop the desired relationship for the Arakawa Jacobian, the following three expressions for the Jacobian are considered. ... _. 9A 3B dk 3B TA . . J(A,B) = — — — = J (2.7) 3x 3y dy dx J ■ h - k (A £ } = jB (2-8) JB> - i - jC 3"l 1,3+1 1+1,3+1 i-l,3+l 1,3-1 1+1,3-1 1-1,3-1 JC ■ ; \ [B. ..-(A..- ,.,-A..- . -) - B. . - (A. n ,..,-A. , , J 4AxAy L i,3+l 1+1,3+1 1+1,3-1 i,j-l 1-1,3+1 i-l,j-l - B... .(A.,- ,.,-A,., . .) + B. , .(A. . ,.,-A. . . ,)] 1+1,3 1+1,3+1 1+1,3-1 i-l,3 i-l,3+l 1-1,3-1 26 (2.10) Collecting terms and factoring yields the desired form. J(A,B) = l [B... . (A. . ,+A.^ . -A. ..--A..- ..-) + 12AxAy i+l,3 i,j-l 1+l.j-l i,j+l 1+1, j+1' B4 «X1 (A--l-1 -+A--i.1 -U-l-A- 1 -_A- 1 iA.0 + i,j+l i+l, J i+l, 3+1 i-l,J 1-1,3+1 B. . . (A. .,,+A. ,,,-A. . -A. . ) + i-l, J 1,3+1 1-1,3+1 1,3-1 1-1,3-1 B. . . (A. ,+A. . --A..- ,-A... . .) + 1,3-1 i-l,3 1-1,3-1 i+l, J 1+1,3-1 1+1,3+1 i+l,3 1,3+1 B. . .,- (A. ,,,-A. . .) + 1-1,3+1 1,3+1 1-1,3 B. . (A. .-A. . J + 1-1,3-1 i-l,3 1,3-1 B... . - (A. . ,-A.., .)] i+l, 3-1 1,3-1 1+1,3 C. RELAXATION METHODS One general approach to solving a set of simultaneous equations is the relaxation method. This method attacks the problem by determining the difference between an estimated solution to the equations and the actual solution. Each successive solution to the equations differs from the "true" solution by a determined amount. The aim of this method is to manipulate the equations so that the differences from the "true" solution approach zero. The manner in which this is achieved is explained thoroughly in Ref. 1 and discussed briefly in the following three sections, (II-C-1, II-C-2, II-C-3) . This discussion begins with the definition of a residual. 1. Residual One can write a set of simultaneous equations such that all the terms are on one side of the equality sign, for example: 27 X + Y - 3 = 0 Residuals are defined as quantities which represent the values of the equations when values are assigned to the variables, X and Y. The solution to the equations is reached when the residuals equal zero. X + Y - 3 = R (2.11) x 2X + 3Y + 2 = R (2.12) y The residuals, R and R , are a measure of the closeness of X and Y to x y the exact solution. The problem of solving simultaneous equations then becomes one of finding values for X and Y which make the residuals approach zero. One method for doing this is called relaxation. 2. Relaxation The relaxation method involves taking the largest residual and computing a value for its associated variable which makes the residual zero. Using the above two equations as an example, the variable X is arbitrarily associated with the residual in Equation 2.11. Likewise, the variable Y is associated with the residual in Equation 2.12. Assuming that R is the larger residual, then to make R = 0 only the value of X would be adjusted. Now that R = 0, R is the j x y largest residual, and the value of Y in Equation 2.12 would be adjusted to make R =0. This process is continued until the values of R and R y x y are within acceptable limits. It is not economical in terms of computer time to determine which equation has the largest residual and relax only it. Making many sweeps through all of the equations and relaxing them simultaneously is a more efficient method. The reason for this is that logical statements 28 take much more computer time to execute than arithmetic statements. The process which relaxes all of the equations simultaneously until all of the residuals are smaller than a predetermined value is called simultaneous relaxation. Modifying this procedure, it is found that even more computer time is saved if simultaneous overrelaxation is used. 3. Overrelaxation Overrelaxation is a process whereby the residuals are not reduced to zero but are overcompensated for in order to speed the rate of convergence. This is the usual process used because it minimizes the computer time needed to solve the equations. This method is approached by writing the numerical equations in such a way that the residuals are multiplied by a convergence factor which is called the optimum overrelaxation factor, Ropt. A critical part of the problem is determining the factor which will yield the fastest convergence rate. A . Optimum Overrelaxation Factor, Ropt a. Introduction A factor, Ropt, can be determined which maximizes the convergence rate of a set of simultaneous equations. The convergence factor must be between 1.0 and 2.0 and its choice may be critical, because small differences in the factor may make large differences in the convergence rate. There are two ways of determining Ropt. One way is through theoretical considerations and the other is to determine it empirically. These methods are discussed in the following two sections. b. Theoretical Ropt Application of theoretical considerations must be restricted to simple partial differential equations and simple boundaries because 29 of the complexity of the problem [Ref. 7]. Frankel applied the concept of overrelaxation to Leibmann's method for solving the following equation Z = V2^ (2.13) Leibmann's method for solving the above equation involves deriving its numerical representation and solving for H*. This method produces yk+1 m k _R££t [VFk + M + B2 k + ^k+1 i»j i,J 2(1+B ) 1+1»1 1-1»J ' i»J+l i,J-l 9 k ? k - (Ax) Z , - 2 (1 + B ) V .] (2.14) The bracketed term represents the residual, and the solution is being reached as it approaches zero. Frankel determined that for a rectangular grid with Dirichlet boundary conditions the optimum overrelaxation factor is given by where and 1 _ A-a Ropt = 2 [- v-±-±] (2.15) cos (£) + B cos (J) a = [ 3 —] 1 + B Ax B - Ay N = Number of rows (Ax) M = Number of columns (Ay) It is found that for a 41 x 80 grid: 30 B = 1 Ropt - 1.8574 B = 2 Ropt = 1.8916 B = 3 Ropt = 1.9064 c. Empirically Determined Value It can be seen from the last section that Ropt, theoretically, depends upon the grid size ratio as represented by B and the dimensions of the array being relaxed. A study was conducted to experimentally determine Ropt for the following grid which is the one used in the simple cavity flow problem. 41 GRID POINTS 80 10 21 21 10 Figure 2.1 Numerical grid Leibmann's equation for solving for the stream function is the one relaxed. This is Equation 2.14 with Ropt = 1. Evaluating the equation involves initializing the grid and relaxing the equations until every value on the grid differs from the 31 previous value by less than 1.0 x 10 . This is approximately the limit of accuracy on the IBM 360/67 Computer without resorting to double-precision accuracy which would significantly increase compu- tational time. To further optimize use of computer time the convergence is checked only once every ten passes through the grid. Each pass con- sists of a sweep through the grid from left to right, top to bottom followed by a sweep from right to left, bottom to top. A limit of two-hundred passes through the grid is placed on each test. Results of this study are summarized in Table I and pre- sented in graphical form in Figures 2.2, 2.3, and 2.4. Table I gives the number of passes through the grid required for every value on the grid to converge. The criterion for convergence is that the difference between two successive passes is less than 1.0 x 10 . Results are given for four different Reynolds' numbers and for aspect ratios of 1, 2, and 3 represented by Ax = .05, .10, and .15 respectively. In all cases Ay is held constant at .05. Figures 2.2, 2.3, and 2.4 represent the same data in graphical form. This illustrates the effect of Ropt on the convergence of the stream function equation in a more explicit manner . One conclusion which can immediately be drawn from these empirical results is that Ropt is also a function of the Reynolds' number. The most significant result of this study is that overestimating the value to be used for Ropt may significantly increase the convergence rate; whereas, underestimating Ropt is not nearly as critical. A brief comparison between the theoretical and empirical results is as follows: 32 w ►J 9 p. o Pi MH o CJ 0) 60 0) > d o o CO Pi c o •H 3 CD H m m o o o O O o o o m O m n O o m m m n o O rH VO o o o o C o o o CM o in CM o o m rH en en o O rH CM CM CM OJ CM CM CM CM rH CM rH rH OJ OJ rH rH rH rH Ol CM iH o o o O o o O o m o m m O m o o o O O o o O II rH o o o o O O o en o CM o o On o o o o O o o O OJ Pi • CM CM CM CM CM CM CM rH CM rH rH CnI OJ OJ CM OJ CM CM OJ CM m o O O o O C o o O O o O O o o O o O o o O o o o o o O o o o o o o O o o o o o O o o O CM CM CM CM CM OJ CM CM CM CM OJ CM CM CM OJ C-J OJ CM CM OJ CM m m Ln m m m m m m m m in m m n m m m n in m O rH CM rH rH rH o o o o on On on ON On CO CO 00 03 o- r^ CO O rH rH rH rH rH rH rH rH CM o !-H o m m in m m m m m m m m m m m m m Ln in m m O II rH rH r-l rH o O o o a\ On CT: on On CO 00 CO t-^ r~- r-- r^- vD o • Pi • rH rH rH rH rH rH rH CM m m LO U") m m m m m m m m m m m Ln m in m m un in o o o on On OJ 00 00 i^» r^. r^. r-^ CO no CO o on o o rH rH en rH rH rH rH rH rH rH . m in in in in m m m m m m m m in m m m m m m in o rH CM rH rH rH iH o o o o CIA o-n o> o-> On 00 CO 00 r^ r^ h» o rH rH rH rH rH rH rH rH rH CM o m o in LO in m n m m in m m m m n m m m in m in m in II iH CM rH rH rH rH rH rH rH rH rH rH o rH o rH o rH On G> o> ON gn CO 03 oo oo OD r- rH on OA <^ on on ON ON on on ON o o rH CM C-) . r^ r^ O • rH rH rH rH rH rH rH rH rH rH CM o o rH o in m un in n m m m in m m m m m m Ln m in m m in II rH CM CM rH rH rH rH o o o o o o> ON o-n On CO CO 00 OO r-- CM Pi • rH rH rH H rH rH rH rH rH rH rH rH m m m m m m in m m m m in in in m in m m m m in o o rH rH rH rH o o o o o o O rH rH rH CM CM CO -3- MD r-» O rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH CM K O n O m o m O m O m o m O LO o m O m O n O O O oj m r-- o CM in r-» o CM m r^ o OJ m r- O CM m r-^ o / *j r^ r*» r~» r- CO CO co CO ON p- o rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH rH Pi 33 jJ Theoretical Ropt Empirical Ropt Re=100 Re=50 Re=10 Re=l 1 1.8574 1.650 1.725 1.625 2 1.8916 1.875 1.875 1.875 1.700 3 1.9064 1.850 1.850 1.825 1.800 From the above summary, it can be seen that the theoretical value for Ropt is usually greater than the empirical value. This situation might lead to a poor choice for Ropt; therefore, trying to apply the theoreti- cal results from simple problems to more complex one might lead to a bad choice for Ropt. D. CONSERVATIVE AND TRANSPORTATIVE PROPERTIES One property which is desirable in a finite difference method is the conservative property. A finite difference scheme possesses the conservative property if it preserves a certain integral relation of the momentum equations [Ref. 7]. This integral relation states that over some region the time rate of change of a given quantity must equal the net flux of the quantity across the boundary plus the production rate of the quantity within the region. Whether a finite difference scheme preserves the conservative property depends upon the form of the momentum equations used and the numerical scheme employed. The possession of this property does not imply that the scheme is more accurate than one that does not possesses it unless one of the criteria of the problem is the conservation of a given property. Upwind differencing, used in the simple cavity flow problem, conserves vorticity. The Arakawa Jacobian conserves vorticity, vorticity squared, linear momentum and kinetic energy; however, all procedures that employ the Arakawa Jacobian do not necessarily conserve these properties. 34 200 CJ W O c*s S3 o c_> w CO CO o u Pi o Q Pi M cy CO w CO CO < o Pi w M 160 140 P3 120 100 1.4 1.5 16 1.7 18 VALUE ASSUMED FOR Ropt 1.9 Figure 2.4 Plot of overrelaxation factor vs. number of passes required for convergence for Re = 1, 10, 50, and 100 and aspect ratio of 3. 37 Another desirable property is the transportative property. The finite difference form of a flow equation possesses the transportative property if the effect of a perturbation on a transport property is advected only in the direction of the flow [Ref. 7], An intuitive approach to the problem leads one to agree that a perturbation in the flow field should only move in the direction of the flow. The method of upwind differencing possesses the transportative property. E. CONVERGENCE AND STABILITY The mathematical basis for stability and convergence lies in .the theory of linear partial differential equations. Linear theory is used as a guideline for analyzing non-linear equations such as those used in this study. A convergent finite difference scheme is one in which all values of the finite difference solution approach the values of the exact solution as the finite difference size approaches zero. All the terms in a finite difference representation can approach the corresponding terms in the analytic equation, and yet, it is possible that the entire equation will not approach the exact solution; therefore, the convergence criterion will be unfulfilled. Convergence, therefore, is concerned with the limit of the entire equation and not the individual terms. Stability is achieved when the cumulative effect of all round-off errors is negligible. This implies that the errors do not increase exponentially. In the time dependent approach used in both the simple cavity flow problem and the refined model, the stability of the numerical method employed to solve the vorticity equation depends upon the size of the time increment and the degree of convergence of the stream function 38 equation. A critical time step can be determined for each problem. Exceeding this time step causes the equations to become unstable and diverge. It is found that the critical time increment depends upon the grid spacing, Ax and Ay, the Reynolds' number and the depth (in the refined model) . The most stable solutions for the simple cavity flow case occur at Reynolds' numbers of approximately one hundred and stability decreases as the grid size ratio of Ax to Ay increases. For the refined model, the greatest stability is found around Reynolds' numbers of one- tenth. Two is the only aspect ratio that is used in the refined model. 39 Ill . GOVERNING EQUATIONS A. BASIC EQUATIONS OF MOTION The governing equations in this study are derived from the hori- zontal momentum equations and the continuity equation. Developments of the momentum and continuity equations can be found in many books on Physical Oceanography such as: The Principles of Physical Oceanography [Ref. 5]. Applying Newton's second law to a continuous volume of fluid, the horizontal momentum equations can be derived. The equation of con- tinuity can be derived from the principle of the conservation of mass. Starting with these equations, the complete development of the governing equations used in this thesis is derived below. 1 . Horizontal Momentum Equations The Eulerian representation of the momentum equations desired in this study is derived from Newton's second law utilizing the following simplifying assumptions: a. The curvature of the earth is negligible for the distances considered in this study; therefore, cartesian coordinates can be used for the coordinate system. b. Fluid is homogeneous. c. Fluid is incompressible. d. Hydrostatic approximation. e. Vertical component of Coriolis force is negligible. f. Turbulent stresses are proportional to the gradient of the mean flow. These assumptions imply that the only external forces acting are: the Coriolis force, gravity, wind stress, and bottom stress. Assuming that 40 primes denote non-integrated variables, the Horizontal Momentum Equations are 9u' m t 3"' . . 3uf f , 1 3P' A „2 , , 92u' ,. 1N vr— + u' t — + v' fv' = — + A 7 u' + A — - (3.1) 3t 3x 3y p 3x H v,, 2 o Z 3v* , 3v' , 3v' . , 1 3P1 . _2 , , 3"v' T— + u1 - — + v' - — + fu' = | — + A V v' + A — — (3.2) 3t 3x 3y p 3y H v „^2 3z: where in Equation 3.1 3u' , 3u' , 3u' u' - — + v' - — 3y u - 3x fv' 1 3P' p 3x vv A 9 U v ^ 2 = Local rate of change of velocity with time = Non-linear field accelerations = Coriolis term = Pressure gradient term = Lateral shear stress = Vertical shear stress 3z 2 . Continuity Equation The Equation of Continuity is needed as one equation in a system of three equations in three unknowns. Physically, this equation imposes the condition that no two fluid particles can occupy the same space at the same time. Continuing the convention that primes denote non- integrated variables and neglecting any changes in the vertical, the Equation of Continuity can be expressed mathematically as |e. + leu: + saL . 0 0.3) 3t 3x 9y For incompressible, homogeneous flow, the density, p, is constant; therefore, Equation 3.3 becomes: 41 |^ + |^=0 (3.A) dx 9y 3. Integrated Equation of Motion Equations 3.1, 3.2, and 3.4 represent three equations in three unknowns: u, v, and p. The solution to these equations is fully deter- mined when proper boundary and initial conditions are specified. Two approaches may be used in order to solve these equations for the desired 1 xntegrated form. a. The equations can be integrated over the entire water column, but this cannot be done without the non-linear terms being known explicitly as a function of depth. b. The real problem can be replaced by one which assumes a homogeneous fluid and the last term in Equations 3.1 and 3.2 can be replaced by a body force in which u and v are independent of depth. The second method is employed; therefore, the last term in both equations is replaced by a body force which results from the wind stress and bottom stress, defined as wx bx .2 » m zr~ = A — o (3.5) ph ph v 8z2 wy by „2 ? — — - A (3.6) ph ph v 3z2 Substituting these terms into the horizontal momentum equations and dropping the primes to denote that the variables are now integrated over the entire water column, the integrated equations of motion are 9u , 9u 9u _ 1 9P , . _2 tWX t ,- 7v — - + u -r— + v fv= T~ + A_TV u + — r— (3.7) 3t 9x 9y p 9x H pn ph The following discussion is taken from Ref. 4, pp 12. 42 3v , 3v , 3v . 19P A 2 t"7 Tby -T-+U — + V — + fu = — + A 7 v + — — (3.8) 3t 3x 3y p 3y H ph ph Integrating the equation of continuity over depth, and applying Leibnitz Rule of integration yields the desired integrated continuity equation. If the sea surface is assumed to be of constant height, arbitrarily selected as 0, and if h equals the depth of the ocean, the equations take the following form o o J£4" + /#*»-° -h -h Using Leibnitz Rule where u and v are independent of depth o o |- f udz + |- / udZ + ly" / UC -h -h The desired equation is 3 (uh) + jLMO . 0 3x 3y B. DEVELOPMENT OF VORTICITY EQUATION 1. Derivation Starting with the equations of motion (3.7 and 3.8) and the integrated continuity equation (3.9), the desired vorticity equation with the wind stress assumed negligible is derived below. The local wind stress is assumed negligible as this force is not being considered because the oceanic current is the only forcing function being examined in this study. Assuming the wind stress is negligible implies wx wy t = t = 0 Cross differentiating the equations of motion, 3(3.8) _ 3(3.7) 3x 3y 43 yields : 3 ,3v 3u. , 3u ,3v 3u. 3 ,3v 3u, 3v ,3V 3u. 3t 3x 3y 3x 3x 3y 3x 3x 3y 3y 3x 3y , 3 ,3v 3us , r 3u 3f 3v , 3f /0 ,_. + V — (t — ) + f — + u ^— + f -r- + V T— - (3.10) dy dx 3y 3x 3x 3y 3y a v- t^y. ^u\ r 3 rJL—A _ + v I7 (? + f) ■ V ? + F = h u ai (3-18) 3t 3y 3x h 3x 3y h H 3y ph 3x ph If the Jacobian, J, is denoted as „,. „N 3A 3B 3A 3B ,„ ,„. J(A,B) = r- — - — — (3.19) 3x 3y 3y 3x Equation 3.18 can be written H + j <♦. <**» - V2^ + 17 (f > - k (f > <3-20) 2. Nondimensionalization The equations are nondimensionalized to generalize the results, freeing them from the constraints of dimensionality. Nondimensional equations can be written in many different forms depending on the way in which the nondimensional parameters are formed. For example, non- dimensional time can be defined as: 45 I- l L/U or T = l2/Ah where L = Characteristic length U = Characteristic velocity A^ = Horizontal eddy diffusivity The choice of methods depends upon the desired result which in turn depends upon the way in which the equation is going to be solved. Physical interpretations of the parameters developed may also dictate the method to be used. This is because certain parameters may be desired which can be used to help interpret the results. In nondimensionalizing, reference quantities must be defined. The following reference quantities are the ones used in this study: L = Characteristic length U = Characteristic velocity oo J D = Characteristic depth also used in the development are: 0 = Beta-plane approximation to the Coriolis force A^ = Horizontal eddy diffusivity Using these quantities the following relationships between nondimensional and dimensional quantities are derived Length scale x-r '-? Depth scale H - ± H " D 46 Velocity scale U -2- V =2- u u 00 o Time scale T- ' L/U Nondimensional Coriolis parameter The remaining quantities to be nondimensionalized can be done so in terms of the above nondimensional quantities. They are Nondimensional partials of the stream function Knowing that |t = - uh = - UU HD 3y °° then II - r X \ M 9Y = lU D; 3y Nondimensional vorticity From the definition of the vorticity _ _3v _ _3u _ . _a_V 3U, j» C " 8x 9y = l3X 3YJ L thus z - $-> c Nondimensional bottom stress From Appendix A the following relationship is obtained pL h 9y HD ^U°°U; 3Y 47 therefore b - b — = ( J— ) — pH V3LU ' ph Nondimensional Del operator The operator is defined as 2 2 2 2 3x 3y thus 2 2 2 VZ = (L ) VZ Notation is changed at this point for clarity of presentation. Sub- scripts are now used to denote differentiation. This notation transforms Equation 3.20 into r+f r+f ? Tbx xbx cT -* + *x <*& ■ v ? + <^r> - (^ x y y x All of the nondimensional terms are now substituted into the dimensional 1 equation. U zt (r)2 " U~D^ (U /L)Z+3LF DH - + U DH'x x L °° (U /L)Z+BLF DH 1 y l «> 2 *h 3 v z + T 6LU bx by factoring yields U 2 U 2 00 CO Z+%-F H x U 2 Z+^F H bx by AHU°° 2 -V v z + eu i (V-) - i T 3 °° pH y pH x Y = x dividing by: U 2 CO For notational purposes, small x's and y's will be used to denote non- dimensional space coordinates. 48 and defining the two nondiraensional parameters LU 00 Re = — — = Reynolds number e = ■ = Reciprocal of Rossby number CO the equation becomes 74-cF 7+pF 1 ? xbx rby zt - *y «*¥>, + \ 'f1, ■ h v z + E [$r>, - §r>x] (3-21) Employing the definition of the Jacobian, 7+rF 1 ? rbx rby ZT + J (I, (*£££» - ^ V2Z + e [(^r)y - (^-)x] (3.22) This form is compact and neat for notational purposes but it is not suitable for solving numerically because, it is hard to separately control the individual forces. Expanding the Jacobian transforms the equation into one which is more suitable. Expansion of the Jacobian yields j /¥ (Z+e¥)) = w fZ+eF) - w (Z+cY) J ^' lH jj x K H ;y y K H ;X H(Z+eF) -(Z+eF)H H(Z+eF) -(Z+eF)H = y r 1 Y_i _ y r * X-i x l h2 ] y [ h2 ] Y (Z+eF) Y (Z+eF) ,7J_ _. H H 7 L y x x yJ rl (f 2 -f Z ) = X y„ y X + £ (Y F -¥ F ) + (Z+^F) [H ! -H ? ] H Hxyyx u2 xyyx rl 4J C^Z) + f J (V,F) + C^^) J (H,f) H H R2 Substituting this form into Equation 3.22 49 ZT + I J (¥,Z) + f J (f ,F) + (^) J (H,V) = H. (3.23) , 9 bx by vzz + e [(ttt-)„ - Cr5-)J Re pH y pH x" where Z = Local rate of change of vorticity with time — J (Y,Z) = Advective term H — J (y,F) = Planetary vorticity tendency H (Z+eF) J (H,H0 = Topographic vorticity tendency H 1 2 — V Z = Frictional term Re bx by T T e [ (~rr~) - (~^~) I = Bottom friction pH y pH x This form is much more convenient for deleting, weighting or evaluating the acting forces individually. The method for incorporating the bottom friction will be explained in the section on the refined model and the complete development of the bottom friction is contained in Appendix A. 3. Nondimensional Parameters Two nondimensional parameters occur in the vorticity equation. The first one, the Reynolds' number, is defined as LU CO where L, U , and A^ have been defined previously. Physically the Reynolds' number is a ratio of the viscous to the inertial forces. The second parameter, e, is given by *" U 50 This can be interpreted as the reciprocal of the Rossby number which is a measure of the effect of the Coriolis force. The solution of the equations depends only on the values of the parameters, not on the values of the individual terms in the param- eters. Nevertheless, estimates must be made of the individual terms in order to determine the range of values for the nondimensional parameter which should be investigated. Certain terms, in particular 3 and to some degree A^, are fairly well known. The horizontal eddy diffusivity, iL , 7 9 2 has been found for turbulent flow to be between 10 and 10 cm /sec. 7 2 The value of 10 cm /sec will be used in this development. $ can be developed from its definition. 2fi cos 9 o where fi = Angular rotation of the earth = 7.29 x 10 rad/sec 9 = Latitude = 37 (for this study) o r = Radius of earth = 6 . 38 x 10 cm hence, a - 2(7.29 x 10"5) x .8 . p_ _ -13 -1 -1 P = ~ = 1.83 x 10 cm sec 6.38 x 10 Determining representative values for L and U which are the nondimensional oo length and velocity scales is not as easy. Physical measurements can be relied upon to give an estimate for U which can be considered the free oo stream velocity or in other words the average velocity of the current along the coast which is unaffected by the presence of the cavity. The best estimate that can be made indicates that a value of approximately one-tenth to one centimeter/second is a reasonable one. Theoretical 51 considerations are relied upon to obtain a representative value for L. The value chosen for L is the most questionable of all the character- istic values. A value of one mile (1.61 x 10 cm) is chosen since it represents the smallest possible size for any feature to be distinguish- able on the grid employed. This is because one mile is the size of one grid space on the arrays used to solve the equations numerically. Using these values, it is found that approximate values for the nondimensional parameters used in this study are Re = 1.61 x 10~2 e = 4.74 x 10"3 C. DEVELOPMENT OF STREAM FUNCTION EQUATION The vorticity equation is one equation in two unknowns; the vorticity, C, and the stream function, ty. In order to solve the vorticity equation, another equation relating to one or both of these variables must be found. A second equation is derived from the definitions of the vorticity and stream function and equating the two. 1. Derivation Respectively, the vorticity and the derivatives of the stream function were previously defined as C = v - u x y il> = vh x \b = -uh y therefore x y < - <-f > + H> 52 thus , h\b - ib h hil; - i> h _ xx x x 2 2 h h finally, separating terras C = £ V2^ - — (if, h + 1(1 h ) (3.2*) h , 2 x x y y Equation 3.24 can be solved for the stream function if it is converted to its numerical representation by substituting the proper relationships for the derivatives. This process is carried out in the numerical approach sections which occur in both the simple cavity flow problem and the refined model. 2. Nondimensionalization Following the same procedure as in Section III-B-3 and using the same nondimensional terms as used for the vorticity equation, the stream function equation can be nondimensionalized. Substituting into the stream function equation yields Z (J=) - V2 £-> , (U.LD) | (i) - Jj (ij) [Tx(UJ» Hx <2) + L n D \ ] collecting terms U „2UJ ULD . U D2 CO V Y oo oo Z (— ) = ^ (-=— ) " ^ [? H + ? H ] (— — ) L H L2D H2 x x y y L])2 therefore, all characteristic values cancel giving Z = ^ V2¥ - =rr (V H + Y H ) (3.25) H R2 x x y y This is the desired nondimensional stream function equation which is used in conjunction with the vorticity equation to solve for both the stream function and the vorticity. 53 IV. SIMPLE CAVITY FLOW A. GOVERNING EQUATIONS The governing equations for the simple cavity flow problem derive directly from the equations developed in Section III. Two equations are necessary: The vorticity transport equation, Equation 3.23, and the stream function equation, Equation 3.25. ZT + | J (¥,Z) + | J 0 R i R \ ' Zi+1 " UR < ° U. + u. . Z, - Z. 'if U > 0 L l-l L ZT = Z. if UT < 0 L l L The formula for the y-direction term is analogous to the x-direction formula. Interpreting the situation physically, a flow reversal appears as an artificial source or sink depending on the sign (- to + is a source of vorticity) . Averaging eliminates the source or sink. 2. Method for Solving Vorticity Equation Employing the Arakawa Jacob ian This approach to the problem is very simple. Forward time differencing is used for local rate of change of vorticity with time, 57 the representation for the Laplacian is used for the frictional term and the representation for the Jacobian is the one developed by Arakawa. Originally only upwind differencing was going to be used to solve the simple cavity flow problem but when the aspect ratio was increased to three the upwind differencing method would not converge. Utilizing the Arakava method, solutions were obtained but at the expense of decreased detail. The following is the numerical equation which was derived for this procedure. zn+1 - zn At ^ = ir [ k k+1 k Zr,, . + ZK , . - 2Z. 1+1 > J 1"1?J Lj k+1 k + Z. . , - 2Z. Re (Ax) 1 + i,j+l i,.j-l (Ay)2 i*l] (4.7) 12AxAy 1+1,3 i,J-l i+l,j~l x,j+l 1+1, j+1 V. .,- (Z,,. ,+Z.,, .,,-Z. . .-Z. _ ,,_) + i,j+l i+l,3 1+1,3+1 i-l,3 1-1,3+1 Y. . (Z. .,,+z. ..,-z, . ,-z, . , .) + 1-1,3 1,3+1 1-1,3+1 1,3-1 1-1,3-1 v. . , (z. , ,+Z. . --Z..- ,-z,.- . ) + i,j-l i-l,3 1-1,3-1 i+l,3 1+1,3-1 Vl,j+1 (Zi41,j-Zi,j+l) + v. , ,., (z. ..,-z, . .) + 1-1,3+1 i,j+l i-l,3 y4 , . , (z. 1 ,-z, . j + 1-1,3-1 1-1,3 1,3-1 Yi+i,j-i (zi,j-rzi+i,j} 3. Optimum Overrelaxation to Solve Poisson's Equation for the Stream Function The numerical equation used to solve for the stream function is derived in the following manner. Z = V2f Using the representation for Laplacian n = i+1>J j-kii i»i + i»J+1 i>J"1 i'i i,j (Ax)2 (Ay) Letting B = Ax/Ay (Ax)2 Zn , = ^kin , + «^+J . - 2Tk+l + B2 [Yk .,, +Vk+] , - 24>k+1] this yields vk+l ~r- rf,i .+^ .+B2(vFk .11+vfk+1 J-(Ax)2zn .] 1,3 2(1+B ) i+1»J 1_1'J 1,3+1 1,3-1 1,3 This method for deriving an equation for the stream function is known as Richardson's method. If new values are used as soon as they are available, it becomes Leibmann's method. A further refinement is achieved by using optimum overrelaxation. The formula for this technique is yk+1 = yk + Ropt [H,k +H,k+1 +B2(lJ/k +4,k+l ) i,3 i,3 2(1+B2) i+1'J i-1'3 i,3+l i,3"l (4.8) (Ax)2 Zn - 2(1+B2)yk .] 1,3 i,3 Ropt is the optimum overrelaxation factor which in this study has been found to be approximately 1.72 for the best overall results. The bracketed term is the residual and as the solution converges to the exact solution it approaches zero. 4. Boundary Conditions There are many solutions to any given partial differential equation. Boundary and initial conditions must be specified in order to determine a unique solution. Figure 4.1 identifies the boundaries which must be specified in the simple cavity flow problem. 59 1. Inflow 2. Wall 3. Wall 4. Wall 5. Lid 6. Outflow Figure 4.1 Specifications of Boundaries Each boundary and its associated equation are discussed separately in the following sections, a. Inflow (1) Velocity Profile To be able to specify the inflow, a velocity boundary layer profile must be defined. Three velocity profiles are considered: a polynomial, logarithmic, and hyperbolic tangent. It is found that there is little difference in the results obtained by using these different profiles. The profiles are plotted for comparison in Figure 4.2. The equations are normalized so that the values range between 0 and the maximum velocity, 1. Over this range the polynomial and the logarithmic profiles agree well although the logarithmic profile does not approach the maximum value, 1, asymptotically as desired. The first profile considered is a polynomial derived experimentally for small scale flow in Ref. 6. It is expressed mathematically by the following fifth degree polynomial: U = .000493 + 2.408031y - 2.841989y2 + 2.796410y3 - 1.904102y4 + .533821y5 (4.9) 60 H M CJ O ►J w > o H CO :z; w M Q O P5 O a >•. i- x; o ,0 >. ■u •H 0 o r-i > U n •H imvQKnog JO HXQIM AwaNnoa Roua aoKVisia 61 The second profile derives from Prandtl's equation [Ref. 8]. According to Prandtl's theory it can be shown that the turbulent shearing stress becomes where K = Von Karman's constant u = average velocity Directly from Equation 4.10 du _ K / o 1_ dy V p y If it is assumed that the shearing stress, t , is constant and the density, p, is constant then K 'p ./ — = constant = c. V p 1 therefore du _ fl dy " y and u = c. In y + c To normalize this profile it is required that u = 0 at y = 0 u = 1 at y = 1 thus c2 = 0 Cl = 1/ln (e + 1.0) y = ey' + 1.0 where e = 2.7182818 62 The velocity equation becomes U = In (e\ 1.0) ln (£y + L-0) (4al) The last profile considered is the hyperbolic tangent. This profile is considered because it is well behaved and possesses two desirsble characteristics. First, it approaches 1 asymptotically and secondly near the boundary it is, to a close approximation, linear. In order to make the value of the profile vary from 0 to 1 over the values of y = 0 to y = 1, the argument Try must be used. The resulting equation is U = tan h (Try) (4.12) (2) Stream Function Profile The stream function profile is obtained by integrating the velocity profile from the boundary to y. It is found that unless- a well-smoothed stream function profile is obtained from the velocity profile, instabilities in the vorticity equation develop. For example, the following technique does not work well in deriving the stream function from the velocity profile. Given central differencing yields i,.i+i Lozi_ u 2Ay where in the preceding equation the subscript i = 1 denotes the first grid column which is the inflow boundary (Figure 4.1). Thus \i+i - 'i.j-1 + 2AyU Using this procedure, a profile like the one in Figure 4.3 is developed. 63 Stream Function Figure 4.3 Unsuitable boundary layer profile This profile causes the vorticity field to alternate in sign. For this reason, all three velocity profiles examined are integrated analytically in order to obtain the stream function profiles. The three stream function profiles are POLYNOMIAL mn/cn 2.408031 2 2.841989 3 "P = . 00049 3y H ~ y - - y 2.796410 4 1.904102 5 .533821 6 + T y = y + 7 y (4.13) LOGARITHMIC (ey + 1.0) In (e + 1.0) [In (ey + 1.0) - 1.0] (4.14) HYPERBOLIC Y = — log (cos h iry) (4.15) The three profiles are plotted in Figure 4.4. It is apparent how similar the three profiles are. The choice of the hyperbolic profile for the one used in this study is based primarily on the fact that it 64 p o pq o w o a H § O O PC H o RBOLIC .4 .6 .8 1.0 NONDIMENSIONAL STREAM FUNCTION Figure A. A Stream function boundary layer profiles 65 possesses the two theoretically desirable properties previously discussed. (3) Vorticity The development of the vorticity input equation is much simpler than the stream function equation. Only one assumption is necessary in this development , that is f € > - ° 3x dX i , this can be approximated by _3V, = _3_V| 3x' . 3x' . 1,3 2, j Remembering the relationship between the velocity and stream function then 3U. 32f, 8y'l,3 3y2 i,J 3 V, 3x' . = 2 ■s 2', 3 4? „ 2 l,j 3x l,j 3x 2,j Defining the vorticity as 3Ui _ _3Vi substitution yields 2 2 31, 91 . 7 = + ^ 3yZ l,j 3xZ 2,j Centered differencing gives Z = ,J ,:] '-1 I 3,;1 1,:1 }1 (4.13) 1,j (Ay)2 (Ax)2 66 b. Walls All of the walls can be handled in a similar manner. Two assumptions are necessary in order to find the value of the stream function at the walls. First, it is assumed that the no slip condition applies at all the walls; therefore, at the walls U = V = 0 V = Y = 0 y x hence, V = constant Secondly, assuming that the arbitrary constant of integration is zero, then for all the walls ¥ - 0 (4.14) The vorticity is handled in a slightly different manner. Consider first the walls parallel to the X direction, at the wall it is obvious that V = 0, therefore z m™L\ = a2yi wall 3y' 1n .2' n1 wall 3y wall Expanding ¥ into a Taylor's series ° J»° ? Y n^i=X i i + "77 Ay+^r ;r (Ay) wall+1 / wall isy 2 _ 2 . - / / wall 9y wall this gives : Z „ -— SSim (4.15) Wal1 (Ay)2 U equals zero for the walls parallel to the y direction, thus = av, _ A wall 8X1 ... '21 .. wall 3x wall Expanding again yields: 67 o + |?| Ax+^-^4| (Ax)2 r- wall+1 /wall ^x1 n. 2 „ 2 wall 9x in wall therefore wall /, \2 (Ax) c. Downstream Continuation The stream function at the outflow is determined by a simple extrapolation procedure. Letting IL be the outflow boundary and IL-1 and IL-2 the two preceeding grid points, then the desired relation is given by (J) + vj; _ YIL-2 *IL IL-1 2 thus Y = 2^ - H* (A. ] 7) IL IL-l,j IL-2,j V^.J/; The vorticity at the outflow is assumed to be transferred out of the grid. This yields ZIL " ZIL-1 <*•«> This can be interpreted as meaning that there is no production of vorticity between the last two grid points. d. Lid If the Lid is considered to be a frictionless impermeable wall then it is implied that it is a streamline. A streamline is a line along which ¥ is constant; since, the stream function is assigned a value at the inflow boundary then it follows that all of the other values along the Lid must have the same value, or Y = T (4.19) i,LID 1,LID v ' 68 Determining the vorticity at the Lid is based upon two assumptions about the velocity at the Lid. They are f| =0 3X LID _aUi = 3U. 3y LID = 3y LID-l from these Z. TTT. - Z. TTn (4.20) x,LID i,LID-l This can be interpreted as a linear extrapolation of the U velocity component up to the Lid. e. Corners It is desirable to retain the effect of sharp corners in order to approximate reality to as closs a degree as possible. The stream function is no problem since the value at the corner is the same as the rest of the wall, or zero. The approach to the vorticity is a little more complicated and it is handled by using different values for the corner depending upon the point being evaluated. The scheme which is used is as follows A 2(Ti,i+i - Ti,i) 2 (Ay)Z •l.J+1 wall • i+l,j Ax 20' . - v. .) B = 1+1>J LiJ_ 2 (Ax)Z If the vorticity at i,j+l is being evaluated the value A is used for Z. . when Z.,, . is being evaluated the value B is used. This scheme i,J i+l, 3 retains the effect of sharp corners. A similar method is employed for the other corner. 69 5 . Convergence and Stability Criteria As is discussed in Section II-D, the convergence of the vorticity transport equation depends upon the degree of convergence of the stream function equation and the time step taken. The maximum allowable time step in turn depends upon the grid spacing, the maximum velocities and the Reynolds' number. This critical time step can be determined by employing a discrete perturbation, stability analysis [Ref. 7]. Using this method it has been determined that the critical time step for the upwind differencing method is given by Atc = 1 (4.21) max + max 2_ , 1 1 ■, Ax Ay Re ' .2 ,. .2J (Ax) (Ay) The time step taken must be smaller than this critical step if the equations are to converge. It is found, by making many test computer runs, that it takes approximately 100 relaxations of the stream function in order to allow the vorticity equation to converge to within acceptable limits. The number of relaxations is less for the method employing the Arakawa Jacobian than for the upwind differencing method since the former method averages three values for the advective term, minimizing the growth of instabilities. The degree of convergence of the equations is not very strict since this is a preliminary study, although it is felt that the conver- gence obtained is great enough so that any further changes that might occur is less than the accuracy of the plotting capabilities of the Calcomp plotters used. In general, the equations are relaxed until the residual is one-tenth of one percent of its original value. The original value 70 is the one determined by the first guess field which is obtained by making the entire field outside of the cavity equal to the input boundary layer profile and all the values in the cavity equal to zero. Using these criteria and a time step equal to eight-tenths of the critical value, it takes between twenty and forty minutes of computer time on the IBM 360/67 to obtain the desired results for each set of input parameters, which are the Reynolds' number and the aspect ratio. 6 . Computational Sequence Simply stated, the procedure used to solve the numerical equations involves calculating the vorticity transport equation, relaxing the stream function equation, taking a time step and calcu- lating the new vorticity field. This procedure is repeated until the equations are converged to within acceptable limits. A more compre- hensive list of the steps involved in this procedure is as follows: STEP1 - Define all constants which specify the problem and/or minimize computational time. STEP2 - Specify the initial stream function by either reading in an initial guess field or extrapolating values from the inflow boundary condition. STEP3 - (Needed only for upwind differencing method.) Specify the initial velocity fields utilizing the initial stream function field. STEP4 - Specify the initial vorticity field by either reading in an initial guess field or determining the field by using Poisson's equation. STEP5 - Calculate all the boundary conditions for the vorticity. STEP6 - (Needed only for upwind differencing method.) Calculate the velocity fields from the stream function field. 71 k+1 STEP7 - Calculate Z " using the vorticity transport equation. STEP8 - Check the vorticity field for convergence. Exit to STEP11 if convergence criteria is met, continue if not. STEP9 - Calculate the new *F distribution using optimum overrelaxation holding the vorticity field constant. STEP10 - Advance one time step and return to STEP5. STEP11 - Plot stream function and vorticity field computed on Calcomp plotters . 7. Programming Optimization One of the objects of solving the. simple cavity flow problem is to gain experience in programming problems of this nature. An important aspect of programming is being as efficient as possible in order to minimize the computer time necessary to solve the problem. Two basic principles to follow in order to do this are: minimize operations and structure the program efficiently. Minimizing operations can be divided into two categories. First, one can minimize computations by manipu- lating the numerical equations until they are in the most efficient form to be programmed. This basically breaks down to collecting like terms and defining new constants (to be put in the initialization part of the program) in order to eliminate unnecessary calculations. Secondly, all unnecessary computations should be eliminated. An example will demonstrate how this is different from the first category. The critical time step in the simple cavity flow problem involves the maximum components of the velocity (Equation 4.21). This requires many time consuming computations in order to determine the maximum values on the grid. It was found that after five time increments the maximum values of the velocity components varied very little and, therefore, the time increment stabilized. It was 72 therefore unnecessary to continue these time consuming calculations to determine the critical time step since it remained virtually constant. The second basic principle of structuring the program efficiently can make a tremendous difference in the amount of computer time used. A summary of some of the things this author has found helpful is contained in the following list. 1. Avoid subroutines, especially in DO loops 2. Avoid function statements 3. Avoid unnecessary indexing and if an index appears more than once in a calculation, equivalence it to a nonsubscripted variable A. Avoid mixed modes 5. Avoid computed and assigned GO TO statements 6. Convert divisions to multiplications where possible 7. Initialize values in DATA statements Most of the things in the list, which it is advisable to avoid, were developed for the convenience of the programmer, but the price that is paid for this convenience is increased computational time. C. RESULTS The complete results of the simple cavity flow problem, consisting of twenty-four Calcomp plots of the vorticity and stream function, is con- tained in Appendix B. These plots, twelve stream function and twelve vorticity, are for Aspect Ratios of 1, 2, 3, and Reynolds' numbers of 1, 10, 100, and 1000. The plots are of only a portion of the total grid as shown in the following diagram in which the shaded area is that area which is plotted. Contours are not plotted at equal intervals; therefore, the relative magnitude of the circulation within a plot cannot be determined from the 73 stream function gradients. Only the direction of the flow is indicated by the streamlines since the velocity vectors are tangent to the stream- lines. All plots are plotted using the same contour levels so that a realistic comparison between plots can be made. 80 40 GRID POINTS ' - ' ' " ' *'X-X-vXv.XvX\\\ J_ Lu ..-»., ^- *..- _.. y ■^^C'-ltAft^'A^UA^ 41 50 . . .. PORTION OF TOTAL GRID PLOTTED Figure 4.5 Portion of numerical grid plotted The main purpose of this part of the investigation is an examination of the different circulation patterns which might occur. With this pur- pose in mind, four of the more interesting flow patterns are presented 74 in Figure 4.6. An indication of the effect of increasing the magnitude of the flow by increasing the Reynolds' number is shown in the top two plots. Not only is the magnitude of the flow increased when the Reynolds' number is increased, but the pattern becomes more circular and the center of the gyre is positioned closer to the mouth of the cavity. The lower diagrams represent patterns which result if the Aspect Ratio is increased from 2 to 3. It is shown that the pattern starts to split in the case where Re = 10 and the Aspect Ratio = 3. If the Reynolds' number is lowered to 1, the split becomes complete. This latter case is intar- esting because there is some evidence that this type of pattern occurs in Monterey Bay. It had been speculated that the splitting of the cur- rent might be caused by the presence of the submarine canyon. But as seen later, this type of circulation only occurs for the simple cavity flow model. All of the vorticity plots are similar. A typical plot, for Reynolds' number of 10 and Aspect Ratio of 2, is illustrated in Figure 4.7. It can be seen that the walls and especially the corners have a significant effect on the vorticity. It is doubtful that the effect of the sharp corners can be easily extrapolated to what occurs in the real world. Defining closed circulation as that which occurs in any area where the streamlines form a closed curve, it is interesting to note that in all cases closed circulation does occur. It will be seen that this differs markedly from the results for the refined model. The closed circulation varies from occurring in only the corners of the cavity to filling the entire cavity. If it is assumed that this simplistic model approximates the real world, then, important conclusions about a bay as a place within which to deposit effluent can be deduced. These conclusions are discussed in the following section. 75 r— f— I \ 111 ! 1 i i i \\\wm I I IIV iMllW MV i I !!}i I..U i i f !! i! i * i W V» «■ — ' ■« 'i»«w o •H ■U O s QJ U •u CO CO 1^ a> u Zi t>0 •H f»4 77 D. CONCLUSIONS As was seen in the previous section, some form of closed circu- lation occurs in all of the solutions of the stream function obtained. If this is the case in reality, then, an embayment could be a poor place within which to deposit waste. Should the effluent be deposited within one of the gyres, diffusion and tidal currents would have to be relied upon to disperse the pollutant. Diffusion and tidal currents might not disperse the pollutant at a rate fast enough to keep the pollution below a safe level. 78 V. REFINED MODEL A. GOVERNING EQUATIONS The development of the final equations from Equations 3.23 and 3.24 lacks only the formulation for the bottom friction. The bottom friction terms in the vorticity transport equation (3.23) are not in a form which can be numerically solved; therefore, seme other relationship must be found from which they can be computed. Ihe following is a reasonable approach to the problem which yields a form for the bottom friction which can be solved numerically. 1. Bottom Friction Equations This development, due to its length, is derived in detail in Appendix A. The short explanation contained here will introduce the method used to obtain a relationship for the bottom friction. Consider- ing Ekman's solution for slope currents, it can be shown that 2 d w _ ifw g_ cU^ ,2 A A 9n dz v v where the complex velocity field is given by w = u + iv The solution to this equation is ig /£Osh_a£ - > d^ f cosh ah 8n where a- (r) V Now, the bottom stress equals 79 b A . T _ _V QW I ph " h 3z' z = -h. therefore £ = ^ (f«) tanh ah .£ (5.1) By integrating the velocity over depth it follows that the transport, Tr, equals Tr = (|£) (|^) - (tanh ah - ah) (5.2) i dr\ a then equating 5.1 and 5.2 :L_ 0 (11) Tr r tanh ah i ph h tanh ah - ah Reducing this equation to its x,y components and nondimensionalizing yields I — = Y [- — V + — V ] (5.3) pH l H x H yJ K J xby R? R1 I — zzFr-^li}/ .illy] (5#4) pH l H x H yJ v J where Rl and R2 are coefficients which are functions of ax which in turn is equal to ? 1/2 /2 ,fh \ ax = ^ (— ) v Assuming that the depth at a given point, once specified, remains con- stant and f and A are constant, then Rl/H and R2/H need only be computed once at each grid point in the beginning of the numerical program and stored for future use in each iteration of the vorticity equation. A closer look at the equations shows that the bottom friction is represented 80 by depth-dependent coefficients multiplied by the components of the velocity. This follows directly from the fact that F [- Rl r^ + R2 t^-] = - F [Rl-V + R2'U] n H The other bottom stress component takes a similar form. Table II contains values for Rl and R2 which are computed from the full equations. The frictional coefficients and approximations used to calculate them for large ax are plotted in Figure 5.1 and Figure 5.2. The approximations are given in Appendix A. It can be seen from the plots that the approximations are excellent as long as ax is greater than 2.0. Below this value, the full equations must be used. Relating this to the depth, it is found that the approximations are good until h < 3/A (meters) At very small values for the argument ax, which is dealing with depths on the order of centimeters, the solution oscillates about zero and the validity of the solution is questionable. This is not investigated, because at depths shallower than five meters the surf zone is encountered and this regime is controlled by a different set of processes. 2. Final Form for the Vorticity and Stream Function Equations for the vorticity and stream function of the refined model were derived in Section III. These equations with slight modifi- cations, depending upon which forces are being neglected, are used in this section. Substituting for the bottom friction, the equations arrived at in Section III-l, Equations 3.23 and 3.24 become V Re" V^ " H J(H/'Z) " H J(VF'F) " (Z+2F) J(HjVF) H (5.5) U H x H V H x H V 3 y x 81 TABLE II Frictional Coefficients ax Rl R2 .5 -0.19995 6.00305 1.0 -0.19899 1.51133 1.5 -0.19509 0.69132 2.0 -0.18587 0.41530 2.5 -0.17090 0.29412 3.0 -0.15284 0.22951 3.5 . -0.13526 0.18872 4.0 -0.12014 0.15995 4.5 -0.10773 0.13849 5.0 -0.09756 0.12197 6.0 -0.0819 7 0.09836 7.0 -0.07059 0.08235 8.0 -0.06195 0.07080 9.0 -0.05517 0.06207 10.0 -0.04972 0.05525 15.0 -0.03325 0.03563 20.0 -0.02497 0.02628 25.0 -0.01998 0.02082 30.0 -0.01666 0.01723 35.0 -0.01428 0.01470 40.0 -0.01250 0.01282 45.0 -0.01111 0.01136 50.0 -0.01000 0.01020 100.0 -0.00500 0.00505 200.0 -0.00250 0.00251 300.0 -0.00167 0.00167 400.0 -0.00125 0.00125 500.0 -0.00100 0.00100 1000.0 -0.00050 0.00050 2000.0 -0.00025 0.00025 82 1.0 .6 QC APPROXIMATE EQUATION FULL EQUATION Figure 5.1 Frictional coefficient, Rl 83 1.4 1.2 1.0 F) " ^f^ J(H^) H (5.7) + eF [|^ V2^ + J(§^,¥) + (I2-) » + (f2-) fx] y x Tha equations for the refined model contain the volume transport stream function rather than the velocity stream function which is con- sidered in the simple cavity flow problem. This study is interested only in the direction of the flow and net the volume transported. The velocity stream function has the property that the direction of the flow is at all points tangent to the streamlines. Here, it will be proven that the volume transport stream function has the same property. Considering v-w = u|^+ v I- dx oy therefore V-VY = U (+ UH) + V(-UH) thus v-vy E o This means that the velocity is perpendicular to the gradient of the volume transport stream function which in turn means the velocity is everywhere tangent to the streamlines. 3. Convergence and Stability The rate of convergence for the refined model, as expected, is much slower than for the simple cavity flow problem. To obtain 85 convergence to within one percent of the original residual, it takes between two and three hours of computer time. On the average, the equations are converged until the residual is one-half of one percent of the original value. The proper time step is determined experimentally and is found to be about one-fiftieth of the time step used in the simple cavity flow section. Convergence to within one percent of the original residual required approximately 600 time steps. Although the Arakawa Jacobian is a stable representation, it is found that the bottom topography introduced must be smoothed, eliminating some detail, in order to obtain the desired degree of convergence. B. INVESTIGATION OF THE RELATIVE IMPORTANCE OF THE CORIOLIS FORCE AND BOTTOM FRICTION One question which occurred during the development of this study was: how important are the Coriolis force and the bottom friction? To inves- tigate this question the nondimensional parameter, e, and the non- dimensional. Coriolis force, F, is examined since both forces contain these terms. 1. Coriolis Force The magnitude of the Coriolis force and, to a degree, the bottom friction depend upon the magnitude of the nondimensional parameter, e. Recalling that e=f^ (5.8) where and 2ft cos 0 3= °- F = F + Y o 86 where Y - Direction of true north F - t tan 9 o L o ft ~ (7.29x10 rad/sec) angular rotation of earth o r - (6.38x10 cm) radius of earth L - (1.609x10 cm) characteristic length = one mesh length in x direction 0 = (37°) latitude for this study o Carrying out the indicated arithmetic operations yields -13 -] -1 $ = 1.8 x 10 cm sec and _ 6.38 x 108 . ,_o _ F = c tan 37 + Y 1.61 x 10 which reduces to F = 3 x 103 + Y It can be seen from 5.8 that once the characteristic length is deter- mined, then E depends only upon the characteristic velocity, U . Fron the definition for the Reynolds' number, a relationship for U can be determined. ReA u =T-E Substituting into 5.8 gives SL3 e = R6AH Assuming a value of 10 for A^ and substituting the values obtained for 3 and L, then e becomes 7.64 x 10~5 e = - Re 87 This is a convenient relationship because solutions for various Reynolds' numbers are desired; therefore, once the Reynolds' number is specified, the value for e is determined. Assuming, at first, that the Coriolis force is constant, the term appealing as a coefficient in front of both the bottom friction and topographic tendency terms is eF = eF o Direct substitution yields .23 eF = - — Re The inclusion in the model of the constant term, eF, is necessary and presents no problems. Now, the effect of the inclusion of the variation of the Coriolis term with Latitude is examined. The term being examined is given by | J(Y,F) = | [V F -fF] H H x y y x It has been determined that eF = 3 x 103e + Ye Thus, if it is assumed that the grid is aligned such that Y is in the direction of True North, then (eF) = 7.64 xlO"5 but therefore y Re (eF) = eF y y § [T F _ ^y0] = 7-64 x 10~5 v* H L x y jf x Re H 88 This terra is at least several orders of magnitude less than the other terras in the equation and is, therefore, negligible. Thus, although the variation of the Coriolis force is negligible due to the scale used in this study, the presence of a constanc Coriolis force will be shown to have an effect upon the solutions which are obtained. 2. Bottom Friction Ex£mining the bottom friction terms in Equation 5.7 it can be seen that the frictional coefficients, Rl and R2 , appear in all of the expressions. Looking at Table II it is seen that the frictional coeffi- cients rapidly become small and the change of the coefficients in the horizontal direction, due to the change in the bottom topography, becomes negligible. Considering eF, with Reynolds' numbers greater than one, this factor which is multiplied into the frictional term is less than .2. Thus at Reynolds' numbers greater than one or at depths greater than thirty meters the bottom friction becomes negligible. Since this model uses a grid spacing of 0.5 miles in the Y direction and 1.0 miles in the X direction, the area covered and the topographic changes are large. The only place where the bottom friction would have any signifi- cance is at the few grid points immediately surrounding the boundaries where conditions are only approximated. At very shallow depths, less than five meters , the equations no longer represent the current regime because the surf zone is encountered. Also, considering that the velocity profile, irrespective of bottom friction, is generally a maximum at the surface and decreases with depth, the effect of bottom friction is minimized in this situation compared to the situation where the velocity is constant with depth. For these reasons, it is concluded that it is unrealistic to assume that the inclusion of the bottom friction would increase the accuracy of this model. 89 C. NUMERICAL APPROACH 1. Method for Solving the Vorticity Equation The full vorticity equation, Equation 5.7, neglecting the variation of the Coriolis force with Latitude and the bottom friction. reduces to ZT - |^ V2Z - | J(Y,Z) - (Z~^F) J(H,¥) (5.9) H in finite difference form this becomes _n+l _n . Zk .+Zk+| .-2Zk . Zk .^+Zk+1 -2Zk . Z =Z__ = 1_ r i+l,j i-l , j i^ 1,3+1 1,3-1 i,-| (Ax) (Ay) (5.10) 1 (Zi i + eF) i,j H2 . i,J where the Jacobians in the above equation are solved by Arak^wa's technique which is outlined in Section II-B-3. For programming purposes, the above equation is manipulated to yield a similar equation which mini- mizes number of computations needed to obtain a solution. 2. Optimum Overrelaxation for Solving the Stream Function Equation Substituting the appropriate finite difference forms into Equation 5.6 yields z = -1 t 1+1>J i-i».i LJ + 1»-i+1 i».i-i x»-i] Hi,j (Ax)2 (Ay)2 _ _i r 1+1,3 1-1 ».1 1+1, .1 1-1, .1 . 1,J+1 1,1-1 1,3+1 1,1-1 , 2 L 2 2 HT . 4 (Ax) 4 (Ay) i > J 2 An expression for ¥ . .is desired; therefore, first multiply by (Ax) i » J and then collect terms obtaining: 90 2 2 Av k+1 Ay 2 [1 + (t21) ] V. . = (T» , + ¥ , . .) + (7^) (V. ... + T, . .) Ay i,j l+l, J 1-1,3 Ay 1,3+1 i,j-l - (Ax) H Zn - (Ax) [ 1+1'J 1"1>-1 1+1»J 1"1»J i,i i,j H. . / /. \2 1,3 4 (Ax) + i»3+l 1,3-1 i>3+l 1,3-1 j 4 (Ay)2 AX Letting B = — and employing Optimum Overrelaxation k+1 , . k RoDt k k+1 ? k k+1 rv. - (1.0 - Ropt) v. . + opu0 [H* ..+¥. . . + BZC1\ ,._ + Y. .) i»3 i,3 2(1+B ) ^J ,J 1,J 3-»J_1 / t k+l X , - (Ax)2 H. .Zn . - i+l,3 i-l, J i+l, J 1-1,3 i,3 i,3 4H -1- , J E k k+1 - 75 (*. ..1 ~ *. , -,) (H. ,,, - H. . )] 4H . i,3+l 1,3-1 i,3+l i,3-l J- , J Defining the constants and reorganizing :'.n order to minimize computations, this equation is used to calculate the stream function. 3. Boundary Conditions a. Inflow In this model the stream function is not related to the velocity, U, as in the simple cavity flow but to the volume transport, UH; therefore, the stream function is a function of both the depth and the velocity. It is assumed that the inflow boundary layer is completely dominated by the bottom profile with frictional effects being negligible. This being the case, the input velocity is assumed constant and the stream function boundary layer is determined directly from the bottom topography. Knowing that Y = -UH y 91 then one procedure for determining the input boundary layer is 2Ay l,j and vi • ,1 = ^i • t - 2Ay UH1 . This procedure for determining the input boundary layer is found to be unsatisfactory as it is in the case of the simple cavity flow. To solve the problem of specifying the inflow stream function boundary layer, the grid is extended five grid spaces upstream and four spaces downstream as transitional regions. The bottom is specified as linear at the artificial inflow boundary gradually changing into the "true" bottom profile. When this is done, a well behaved input profile can be determined by thus «F = -UH y yi Y - - I UH dy y0 but U is constant and H = CY, therefore finally yl Y = -UC J Y dy yo 2 yi «P = -ucf- yo where C is the slope of the artificial inflow bottom H - R\ C = Y-Y1 92 It follows directly from Equation 3.11 that the vorticity is zero at the inflow since the V component of the velocity is zero and the U component is constant. b. Lateral Boundaries The rectangular shape of the cavity is retained in the refined model. It is assumed that the effect of the bottom topography will simulate the effect of realistic boundaries. This is achieved by using the actual coded bottom but specifying a minimum depth to be used whenever the depth becomes too shallow or the boundary goes overland. This minimum depth is used until the wall is reached. Diagramatically , this is illustrated in Figure 5.3. Using this method, the problems of programming variable boundaries and its associated large increase in computer time is bypassed. In doing this the boundary condition at the wall still must be specified. The stream function at the walls is still equal to zero since no slip conditions are applied which implies that the velocity at the walls is zero. The vorticity can be determined in the following manner. First, consider walls alligned in the X direction. At the wall, V = 0 and wal1 * wall ^*■ *i*-*r*r**ii* * *'»->--• «*«*■- iJLjfcl — — - — ,. ., .-,.— c^ 1 5 rr— - i ■■ ""JU 2 1 Minimum depth region Transition region Walls 14 Figure 5.3 Grid for refined model. 94 wall+1 yKall y^\ in 2 yy , / y; ^ -^ 'wall 'wall therefore f _ wall+1 yy (Ay)2 Substitution yields 2T _ 1_ wall+1 wall "" H. . " . N2 i,J (Ay) Vertical walls are handled in a similar manner yielding 2V 1 wall+1 wall H. . ,. N2 i,j (Ax) Sharp corners are handled in the same way they are handled in the simple cavity flow problem. It is expected that the effects of the boundaries on the solution will be diminished due to the effect of the bottom topography. c. Downstream Continuation The same method as the one employed for the simple cavity flow problem is used. This method gives w = y = 2^ - ¥ outflow IL,J IL-1,J IL-2,J and outflow IL,J IL-1,J d. "Lid" It is assumed that the "Lid" is far enough from the cavity so that the assumption that there is no volume transport across this imaginary boundary will not seriously affect the flow patterns created. It is assumed that the current at the "Lid" is constant and parallel to the "Lid". This makes the "Lid" a stream line and as in the simple cavity flow problem: 95 Zi,LID " Zi,LID-l e. Bottom Topography Two basic bottom topographies are used in this study. The first bottom topography (Figure 5.4) used is a hypothetical one which is very simple in order to minimize any instability it might generate. This bottom is obtained by making the depth at the "Lid" a constant 1600 meters and assuming a linear decrease of depth along the inflow boundary to zero at the first grid after the wall. This makes the value at the wall the minimum non-zero depth. All of the other columns exterior to the cavity are equivalenced to the inflow column. The entire grid is then relaxed, holding the boundaries constant until a steady state is reached. The presence of the cavity craated the bottom topography in Figure 5. 4. To obtain a bottom topography which simulates that of Monterey Bay, a bathymetric chart of the bay is overlain by a 50 x 80 grid and the depths are coded at each grid intersection and then placed on computer cards. The field is then read into the computer and stored for use in the program as needed. The maximum depth of 1600 meters is chosen because it is approximately the maximum level of no motion used in geostrophic current calculations. The field is smoothed by averaging the surrounding four grid points in order to obtain the central grid value , i.e.: H... , + H. - . + H. . + H. H = i+l >.1 r-1, .1 i,.i+l 1..1-1 It is necessary to smooth the field in order to eliminate inconsistencies in the bottom topography. In order to retain the effect of the canyon, the boundaries of the 1600 meter contour are preserved. This boundary 96 ■a U t>0 O o 4-1 B o o 03 u •H ■P 0) o ex m a) M M •H Pn 97 is marked by the dark line in Figure 5.5. Figures 5.6 — 5.9 show a progression from the original field to one smoothed five times. It can be seen that after only five smoothing passes only the general outline of the canyon is still preserved. 4 . Computational Sequence The computational sequence for this portion of the study is very similar to the one used in the simple cavity flow section with a few additions. The first difference is that the bottom topography must either be generated within the program or must be read into the program from some external source: tape, disc, cards, etc. The other main difference is that the vorticity transport equation is much more complex and some of the dimensioned variables occur more than once in the equation. In order to minimize computer time, at each time step, the.se multiple-subscripted variables are equivalenced to non-subscripted variables before the vorticity transport equation is solved. D. RESULTS The vorticity equation used for the refined model is given by ZT = ^ V2Z - I J(^Z) - (Z+e2Y) J(H,Y) (5.11) H In the above equation the stream function is the volume transport stream function. There are two cases which develop naturally from this equation, Case I occurs when eF equals zero. Physically, this is the situation that arises when the Coriolis force is neglected. When eF is assumed to be a constant other than zero, then only the change of the Coriolis force with latitude is being neglected. This is case II. Solutions for both cases are obtained and discussed in the next two sections. A measure of 98 >. t-i a) ■p a o S W) d •H cd rH e •H t m M O o< o •u 6 o ■u o w 99 co 0 •H 100 CO to a to c •H Xl •u o o B CO QJ C o o- cd Vj 60 O P. o e o 4-1 o to o u •H P4 101 co cu CO CO cd P. GO C •H .d ■u o o e CO o X! 4J .d rt M O o •u E o •u +J o PQ oo in CD u 60 •H 102 W QJ W CO ctf a c 4J o o e w > •H 4-1 •H $ O a. o •u E o o m cu bO •H 103 the importance of the various terms in Ecuation 5.11 can be obtained by varying the Reynolds' number and eF. 1. Neglecting the Coriolis Force For this case, results are obtained for both bottom topographies discussed in Section V-3—e. In the solutions obtained here the flow is from the left to the right representing flow from the South to the North. For the hypothetical bottom topography illustrated in Figure 5.3, due to its symmetry, the solutions are the same for flow in either direction. This is true only because the input velocity is assumed constant and therefore the input vorticity is zero, "his is tested in several computer runs by reversing the position of the inflow and outflow and differencing in the opposite direction. No differences in the solutions develop from reversing the direction of the flow, but it is expected that if gyres occur which are not in the center of the cavity then the two solutions would have the gyre positioned on opposite sides of the cavity. If there is an input velocity boundary layer, different solutions would be obtained for flows in opposite directions because of the change in the sign of the input vorticity. Using the hypothetical bottom topography illustrated in Figure 5.3, the volume transport stream function plot in Figure 5.10 is the result obtained. Computer runs for Reynolds' numbers of .01, .1, 1, 10, 100 are made and all of the results are qualitatively identical to this plot. The plot of the volume transport stream function obtained when the bottom topography simulating Monterey Bay is used is contained in Figure 5.11. The results for this bottom topography are also identical over the range of Reynolds' numbers from .01 to 100. It can be seen from the results for the hypothetical bottom topography that within the bay, the streamlines follow the bottom contours. 104 II I It! I I ! i ! I i i i : ' I i ! ! ! I nun i ! I I i ' III I if m-A i \ \ \ \ \ \ \ I I nil I ' hii ! I / / // i I n i < / '7/ ' / / 1 •' / J r/ ' I / , r^ 1 1 I ii I ! !;; i i : ! , G o •H u o M-l s CO CU U M O CU CO h cO 'a, M o o a O iJ o iH «0 a •H JJ a) o o Oh >-. II ^ r » w iH - ^_ ■^\m\\l;fev,fe^ x-- N \ ! \ \\\\ V.' — x\x ^ vo \ \ \ a o ■H ■IJ a % !-t 4J W JJ O 0 ■H 4J iH 3 [3 ■H U) >i ,d P. nj >-i bO • o o tx O II 4J * (J rH CO O -O O O rH 6 I iH T3 O CD • C •H II m CU CD tf Pi in CD bO •H 106 This means that the direction of the flow within the bay is generally along the bottom contours; but at the larger depths encountered outside of the bay, the direction of the flow does not follow the contours. The fact that trie flow does not follow the contours outside of the bay while it does inside the bay is undoubtedly, not only, the result of shallower depths within the bay, but also, due to the fact that the volume trans- port in the bay is much slower than outside of the bay. For the bottom topography simulating Monterey Bay the entire cavity is filled with a gyre. Again, no significant difference is discerned between the solutions obtained over the range of values investigated. It is also noticed that the cavity seems to cause a disturbance downstream from the cavity. Although the field downstream is not long enough to be absolutely certain, it is reasonable to assume that this is a wave induced into the flow by the bottom topography. For comparison purposes, the simple cavity flow problem is solved for the case of Reynolds' number equal to .01 and Aspect Ratio of 2. The result is contained in Figure 5.12. It is evident what a significant difference the inclusion of the bottom topography has upon the resultant solution. 2. Assuming a Constant Coriolis Force After examining the case where the Coriolis force is neglected, the effect of including a constant Coriolis force is investigated. A range of values of the Reynolds' number and Coriolis force must be investigated because these parameters depend upon the characteristic scales used. Since some of the scales, such as the length scale, cannot be determined exactly, solutions must be obtained for a range of values of the nondimensional parameters. It is then assumed that the solution, which models the true circulation the best , falls within this range of values. 107 ' CT-^»» — ■■» "■-- <«JMH— \ \ \ '• V ' 1 U ; <• : : / j l 1 1 I fj ,111 / / / r/A w // / '■■-- ^^^> — . eg I! o 4J O CO •H % O CJ I CO CM m M •H 108 Using the hypothetical bottom topography shown in Figure 5.3, solutions are obtained for eF equal to .23 and 23. These are contained in Figures 5.13 and 5.14, respectively. Each figure represents solutions for Reynolds' numbers of .01 and 1 since the solutions do not differ over this range of values. As can be seen a change of two orders of magnitude of eF has some effect on the flow pattern generated; whereas, a change of two orders of magnitude of the Reynolds' number has a negligible effect upon the solution. This means that over this range of values the effect of the Coriolis force interacting with the bottom topography, ^ j(h,> U) H O 0) o 13 • 9 ,H e i iH id o OJ • ■H II m OJ OJ Pi pci CO m OJ u 60 •H Pn 110 '- —" ' II i ! I ! I I 1 ! i l I i i ' ! I I I I I ! ! ' I I ! i - ii — ., .■■„„ .-_ \i {W \\\ -A 11 \ ! I If) "/y ■ o •r! 4-1 O 4J CO 4-1 V< O CO to r4 4-' cu I-" 5 •a CO 60 O P. o E o 4-1 4-1 o rO rH CCj O ■H 4-1 • . II P^ •> to rH 0) O T3 . O rH e i rH Tj O i\»H\\\\\\\ \\\\\\\ \ \ \ \ r-r-w^.-™ \\\\\\\\W \\\\\\\\\\ \\\\\\\\\\ \ 1 \ \ \ \ I \\\\ -L'^S: c/V'^XW \ Cvp ll/^N \ \ 111! I l|| /I / / / ' ) II rJJ/// i l. . : - c ,—. * s cO i) t-i jj w jj !-J O A IS) a Cd !-i JJ (1) § iH O :> *i >> rd P3 >. QJ J-l CJ 4J ti o fcH K-l <>0 l-l •H J.J « «-H :i U •H • COCN ^ o fi H tx Cd X »J to CO o CM Ch • o J.) II e Pm o (J 4-J 4-1 iH o O x> • n II iH a) 0 •H 112 ll,,m\\\\\\\ \ \\\\\ \ \\ \ \ \\\\v\\\\\\\\ \ \ > w r i\\\\\\\\\V ' } - — -~ !! ! i Iijiji f 1 1 i Mr \\\\\\f^ 1 \ \ \ \) > i o CO a 1 rH 2 cd M >^ cu n 0) +J p o p H 4-1 cd iH P CO • c-g >>o CO X u to m O CN P. • o ■p II O CO +J H o o ,o • • it iH CU CD 'O Pi O E P T3 O 0) «H P U •H O 4-1 P «o 3 Pi <4-l vD CU p 00 fe 113 1 1 U ,\\ \ \ , \ \ \ \ \ All \\\\\\\\\\ ■ . I ' • ■ II III 1 mi III < - : hi I Mil I . - i'H I \ \ \ \ \ \ \ \ \ \\\\\\\\ \ \ 1 \ \ \ \ \ \ , : I ' if i \ \ \ \ \ \ i / ///ii l\ ! , ;',•':' if L ... .,..,.,..-„ , \ \ \ I \ \ \ \ \ i ! I i \ \ , 1 \ 11 1 \ \ i » i ! i I I I > 1 ! ! \ \\\ • i l/i V- v ^ ' \v> \ A \ ■7 , , \ . , \ I \ i ; i i \ v ii// \ i M i > ! i | ! ! | i i i i MM ! ! I i : ! I ] * ! I I till i i I III i I 1 1 ! » , i1 1 r/ i 'Jy/ J/iiw '7// / «W / / .— - A&y /Mil \ i ) ' i i I i I MM ! / 1 i f i i !i l II Mill I I IP i 1 I \\V i vA\- 0 ctj 0) M 4J CO 4-1 J-l o (X to c CO V-i ■U cu e 3 rH O > #\ >% cU « >> ^ ,G O (X • cr) rH ^ txO II O Cb Pm o (0 4-) O E • o rH 4J 1 U rH O o wo • r, II rH Q) QJ *U Pi O B ^ C T3 o cu •H C 4J •H CJ iw B CD ^ P4 M-l M •H P4 114 and eF. It can be concluded, then, that the closed circulation is the result of the presence of the Monterey Submarine Canyon and that the bottom topography has the most significant effect upon the circulaticn pattern generated. E. CONCLUSIONS The foundation for the development of a good model for investigating the current-driven circulation in an embayment has been established. Included in this model are the effects of advection, planetary vorticity tendency, topographic vorticity tendency and lateral shear stress. The relative order of magnitude of the bottom friction and the change of Coriolis force with latitude is investigated and determined to be negligible: compared to the other acting forces and assumptions that are made in the development of this first stage of the model. There are two areas in particular where the model could be improved. First, effects of the inclusion of the bottom friction should be investi- gated further. Effects of the bottom friction diminishes rapidly with increasing depth, but in shallow areas around the boundaries and within the bay, the bottom friction might have a significant influence upon the circulation patterns. Second, some of the boundary conditions should be refined. The walls representing the coastline might be represented in a more realistic manner. This would go hand in hand with the inclusion of the bottom friction. Also, the inflow and outflow should be represented by less restrictive equations. The representation for the bottom friction is an approach suggested by Dr. Jerry Gait and developed in this thesis by the author. It is believed to be a new approach and one worth investigating further. This approach 115 is of special interest because the bottom friction is a linear function of the velocity making it amenable to analytic solutions. It is assumed in this model that the input vorticity is zero; because of this, the solution does not depend upon the sign of the vorticity gen- erated at the boundary which is a function of the direction of the flow. If a velocity boundary layer exists, as many studies propose, then the vorticity will not be zero at the inflow and the solution will, in addition to the other variables, become a function of the direction of the flow. This is a case which should be investigated not only to find out the effect of a boundary layer upon the solution but to determine the effect on the solution of changing the direction of the flow. Results of this study indicate that the occurrence of closed circu- lation in Monterey Bay is a distinct possibility. This investigation indicates that the presence of the submarine canyon seems to increase the possibility of the occurrence of closed circulation. In conclusion, it can be stated that this study indicates that the bottom topography is the controlling factor in determining what type of circulation pattern is to be expected in a given bay. Over the ranges of values examined for Re and cF, the change of the vorticity from one case to the next had little effect upon the solution of the stream function. This is because in deeper water the volume transport becomes much larger than the vorticity and, therefore, the stream function equation becomes I V2Y _ I_^ (H/ H + y H ) = o H R2 x x y y which is, for all practical purposes, independent of the vorticity., This puts additional constraints upon the types of flow patterns which may occur, 116 APPENDIX A DEVELOPMENT OF BOTTOM FRICTION Assuming that the coefficient of vertical eddy diffusivity, A , is constant, the bottom stress is proportional to the change of the mean velocity gradient with depth. Using Ekman's solution for slope currents an expression for the velocity as a function of depth can be obtained, but in this relationship the velocity is also a function of the slope of the sea surface. In order to eliminate the sea surface slope, which is not easily obtainable, the velocity can be integrated over depth to give an expression for the volume transport. Equating these two relationships an expression for the bottom stress which does not depend upon the slope of the sea surface can be obtained. The first step in this development is to obtain the analytic solution for the velocity from the following three equations 2 fv + A 1-a . - 1 1£ (A-!) V 3Z2 p dx 2 -fu + A ^--if (A-2) v 8z2 p 9y P = -pgh (hydrostatic pressure) (A-3) These equations can be derived from equations II-l and II-2 with the addition of the following assumptions 1. Steady state conditions prevail 2. Non-linear terms are negligible 3. Lateral turbulent stress is negligible Ekman solved these equations analytically in the following manner: 117 and 3P 3h "3x~ = "Pg al 3P ah a? = ~pg ^ Substitution yields 2 8 u _ Jfv £_ _3h - 2 A A 8x 3Z V V 2 9 v _ j[u j*_ ah thus defining it follows . .2 A A ay 3z v v J —2 (u+xv) = - (u+iv) + — (— + i — ) 3z v v y w = u + iv 3ri 3x 3y A_w = ± fir +i_l£ (A_4) ,2 X A A 3n kA q; dz v v Ekman found that the solution to this equation is w i£ cos£_az _ |c f cosh ah 3n where a= (r) V Substituting A-5 into equation A-A readily shows that it is a solution. Given an expression for the velocity, expressions for the bottom stress and the transport can now be found. The bottom stress is given by 118 since b A - T _V 3w I ph = h 3z ' z=-h 3w 3? |5 = a (-ifi) tanh ah , 3 z f 3n therefore b A a . — = -t— (rf) tanh ah — - ph h f 3n (A-6) Determining the transport is as easy as the determination of the bottom friction. Knowing that o then thus T = J wdz -h •T - f (i&) (f ) (£2^i - 1) dZ r J f 3n cosh ah T - (^) (|5") [ ^— r -1- sinh az r f 3n cosh ah a - Z -h -h finally T = (§*) (|£) (-) (tanh ah - ah) r r dri a Equating A-6 and A- 7 gives — a T = 7— a (rf2-) (— ) tanh ah [1 - - — - — - h r h f 3n tanh ah (A-7) from this the desired expression for the bottom friction is obtained b A n T v 2 _ , tanh ah ] (A-8) This form is not convenient for programming. Transforming equation A-8 into a more convenient form involves separating it into its x,y components 119 This is accomplished in the following manner: first, the transport is defined in terms of the stream function but T = T + i T r x y T = UH = -Y x y T = VH = ¥ y x thus T = -¥ + i ¥ r y x secondly, the hyperbolic terms must be separated into its components. Defining 9 1/2 , .1/2 ,fhA .1/2 ah = i (— — ) =i x v therefore 1/2 T .. if „ r TANK i ' x where ph ' h Tr *„ .1/2 .1/2 ] (A 9) M TANH l x - x x .1/2 ._o . . ,_o ft . . Jl l = cos 45+i sin 45 = — r- + i — r = a+ia Substituting for the arguments of the hyperbolic terms yields 1/2 TANH i ' x TANH (ax+iax) _AMU .1/2 .1/2 TANH (ax+iax) - (ax+iax) TANH l x - l x through fundamental hyperbolic trigonometric identities and some simple complex variable manipulations the following relationship can be arrived at TANH (ax+iax) = Rl + i R2 TANH (ax +iax) - (ax+iax) in which : 120 R1 G-C G + 2 [A-C] R2 = 3 Z G + 2 [A-C] and A = ax [COSH 2ax + COS 2ax] C = SINH 2ax + SIN 2ax D = SINH 2ax - SIN 2ax _ (COSH 2ax - COS 2ax) G — (A-10) ax substituting into equation A-9 gives b V = 7T- (-V +iY ) (Rl + i R2) ph h y x expanding terms yields bx by f ■V- + i "tt- = r [-R1Y + R2^ + i (-R2Y - R1Y ) ] ph ph h x y x y equating real and imaginary terms leads to the desired bottom friction components, i.e.: bx ^7— = - (-R1Y + R2^F ) ph h x y and T^ f -^r— = - (-R2T - RTF ) ph h x y Approximations for the frictional coefficients, Rl and R2, can be derived when 2ax is greater than 4 by eliminating negligible terms in the equations A-10. The approximations can be shown to be 1 - ax Rl = 2 2 (ax) - 2ax + 1 121 R2 = f 2 (ax) - 2ctx + 1 After these equations were developed it was determined from the relative order of magnitude of the terms involved that the bottom friction is negligible; therefore, this representation was not tested. Plots of the frictional coefficients and tables of it's values are contained in Section V. By way of comparison, though, it can be shown that the results from this method can be interpreted to agree with the values for the bottom friction used by Hansen who represented the bottom friction by b t I 2a. 2 /TT^ t = (ru vu +v , rv vu +v ) where r = 3 x 10~3 This representation is the same as b 0 ,„-3 .,2 t = 3 x 10 x V where V is the velocity vector. The method described in this thesis gives approximately the same value for "r" over the depth range of 2000 - 3000 meters if it is assumed that A = 100. A big difference is that the v 2 coefficient in this study is multiplied by V instead of V , but if it is assumed that the mean currents over depth investigated are of the order of magnitude of 1 cm/sec then the results from both methods are similar. It is impossible to compare the two exactly because the methods for computing the bottom friction are dissimilar. This method does allow the bottom friction to be represented in a rational manner and has the added attraction of being linear in velocity making it amenable to inclusion in analytical solutions. Walter Hansen, Hydrodynamical Methods Applied to Oceanographic Problems, Proceedings of the Symposium on Mathematical-Hydrodynamical Methods of Physical Oceanography, p. 25-34. Institut fur Meereskunde der Universitat Hamburg (1962) . 122 APPENDIX B SOLUTIONS FOR THE SIMPLE CAVITY FLOW PROBLEM This appendix contains Calcomp plots of solutions to the simple cavity flow problem. There are twenty-four plots, twelve of the stream function and twelve of the vorticity. Two different techniques are used to obtain these solutions. Upwind differencing is used to obtain as many solutions as possible but when this technique failed to converge the second method, employing the Arakawa Jacobian for solving the advective terms, is used. The following is a breakdown of the technique used to obtain each solution. Reynolds ' Number 1,10,100,1000 1,10,100 1000 1,10,100,1000 Aspect Ratio 1 2 2 3 Technique Upwind differencing Upwind differencing Arakawa Arakawa 123 y*r»*' 1 I I ■ • • ; -. ■ - \\v> . ' ifr^z II rH II a rH t PQ 0) |X4 124 M 5 u •H O •H •P O > !^ •H O a) tH CM pq •H 125 II Si o •H O 4-1 Q) o '■W >> •U •H > o t CO « PQ 60 •H Pn 126 CM II & 4J o > cS tH •H a f CO 60 •H 127 CO II pq £ o •H U O s H M-l > CO O (1) I to O U bO 128 CO II PQ Pi P^ +J u o > o t-i •H > CO O •H C/3 PQ Ci) •r-l 129 II S a o •H U O 4-1 I CD w % U-l 4-> •H > a a) ex. a W a) •H 130 M & CJ •H ■U O > I r-l <4-4 >. 4J •H > cd O Q) I CO 00 CU u 00 •H 131 I I 11 I Mi ■ 11 \ ' ; \ J CM II PQ O o •H •U CJ § +J en R > O Q) rH I CO cr> a) 60 132 CM II •H O •H 4-t J-l O > > o a) iH §• •H CO 0) M •H 133 en il $ c o •H •U o % a) CO i 4-1 >^ •H o cu I PQ CD J-i M •H 134 CO n a) ■u •H O •H •P l-i O > o iH •4-1 > o a) CO CvJ PQ cu 60 •H p4 135 II it Mn in o o o •H o c QJ H o iH w CO QJ •rl 136 II PQ O O $ P4 >^ 4J •H O •H ■U !-i O > > CO O CU I CO « U bO •H 137 I . . ■ - ' 1 1 111 ■■■HMBBBI . :nx\\\\\ mi ii, I///// ., K J»//j ■z/y/j - ! I I -> t^Jfrv CM « o o o •H o J-l CO IM >^ •H > o a) iH d. E •H CO m 0) J-i 60 •H 138 ' ■ t ■ ' •— II O o rH Pi 4-1 O > 4-1 i o •H CO VO M u 60 •H 139 CO II O O JSJ d o •H 4-t o a e o 4-> w 5: O iH M-) >^ 4J •H > «d o a. •H 60 •H 140 CO II o o •H O •H 4J U m >> •H E o rH I oo rH •H 141 pq o o o &>&; Pi a o •H •U O § M-l E m CD CO +J > ctf CJ ex CO CT> PQ U 142 II PQ O O o Pi •H O •H 4J o > o H-i S o I o CM PQ 0) bO •H 143 hi .,*»***■ y.tnmjuw»Mw>wi.« ! I ! i ■ ■ i ■ \ ■ - Mk 11 \ I w , '• I ' I 1 / i ; | i ' i - I i I i :>'s 2/ "?t^! I ;5 I ; . ,J CN II M O O O iH £ c o •H o I •U CO o rH IW >^ •H > CO CJ 4) I iH CM PQ CD 60 p4 144 cm II to o o o & u u o > 8 fH B O t CM CM a) 3 •H 145 CO II o o o V d 0 O •H O I w i+-t •H & O 0) I CO C\| PQ 0) •H 146 •— ««-<-o—57" - v-.ii 4J •H > o a) e CM 60 •H 147 REFERENCES 1. Allen, D. N. deG. Relaxation Methods in Engineering and Science. McGraw-Hill, Inc., New York, 1954 (pp 1-74). 2. Ames, W. F. Nonlinear Partial Differential Equations in Engineering. Academic Press, New York, 1965. 3. Haltiner, G. J. Numerical Weather Prediction. Naval Weather Research Facility, Building R-48, Naval Air Station, Norfolk, Virginia, 1968. 4. Holland, W. R. Wind-Driven Circulation in an Ocean with Bottom Topography. Ph.D. Dissertation, Scripps Institution of Oceanography , University of California, San Diego, California, 1966. 5. Neumann, G. , and W. J. Pierson Jr. Principles of Physical Oceano- graphy . Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1966 (pp 197-219). 6. O'Leary, R. A., and T. J. Mueller. Correlation of Physical and Numerical Experiments for Incompressible Laminar Separated Flows. University of Notre Dame, College oi: Engineering, Project Themis UND-69-4, 1969. 7 . Roach , P . J . Numerical Solutions of Compressible and Incompressible Laminar Flow. Ph.D. Dissertation, University of Notre Dame, Ann Arbor, Michigan, 1968. 8. Schlicting, H. Boundary Layer Theory. Sixth Edition, McGraw-Hill, Inc. , New York, 1968. 148 INITIAL DISTRIBUTION LIST No. Copies 1. Defense Documentation Center 2 Cameron Station Alexandria, Virginia 22314 2. Library, Code 0212 2 Naval Postgraduate School Monterey, California 93940 3. Assistant Professor E. B. Thornton, Code 58 1 Department of Oceanography Naval Postgraduate School Monterey, California 93940 4. Lieutenant Roland A. Garcia 1 Fleet Weather Facility Keflavik, Iceland FPO New York 09571 5. Commanding Officer 1 U. S. Army, Corps of Engineers North Pacific Division Office San Francisco, California 94102 Attn: Mr. Orville T. Magoon 6. State of California 1 Department of Water Resources P. 0. Box 9137 Sacramento, California 95816 Attn: Mr. Robert Ford 7. Dr. John H. Phillips, Director . 1 Hopkins Marine Station Cabrillo Point Pacific Grove, California 9 3950 8. Dr. John Harville, Director 1 Moss Landing Marine Laboratory Moss Landing, California 95039 9. Director, Naval Research Laboratory 1 Attn: Technical Information Officer Washington, D. C. 20390 10. Director 1 Office of Naval Research Branch Office 1030 East Green Street Pasadena, California 91101 149 No. Copies 11. Commanding Officer 1 Office of Naval Research Branch Office Box 39 Fleet Post Office New York, New York 09510 12. Chief of Naval Research 1 Ocean Science and Technology Group Code 430 Office of Naval Research Washington, D. C. 20360 13. Chief of Naval Research 1 Code 416 Geophysics Branch Office of Naval Research Washington, D. C. 20360 14. Oceanographer of the Navy 1 Office of the Oceanographer of the Navy 732 North Washington Street Alexandria, Virginia 22314 15. Director 1 Coastal Engineering Research Center Corps of Engineers U . S . Army 5201 Little Falls Road, N. W. Washington, D. C. 20315 16 . Commandant 1 U. S. Coast Guard Headquarters Washington, D. C. 20591 17. Director 1 Coast and Geodetic Survey, ESSA Washington, D. C. 20230 18. Waterways Experiment Station 1 U. S. Corps of Engineers Vicksburg, Mississippi 39180 Attn: Head, Research Center Library 19. Commanding Officer 1 Fleet Numerical Weather Central Naval Postgraduate School Monterey, California 93940 150 No. Copies 20. Mr. Willard Branson, President 1 Associated Monterey Bay Area Governments 1107 Escalona Drive Santa Cruz, California 95060 21. Director 1 National Oceanic and Atmospheric Administration U. S. Department of Commerce Washington, D. C. 20235 22. Division of Oceanography 1 Environmental Science Service Administration Silver Spring, Maryland 20910 23. Dr. Henry Crew 1 Oregon State University Department of Oceanography Corvallis, Oregon 9 7331 24. Mr. Neil J. Wilding 1 Ebasco Services Incorporated Two Rector Street New York, New York 10006 25. Dr. C. S. Fang 1 Virginia Institute of Marine Science Gloucester Point, Virginia 23062 26. Dr. Donald R. F. Harleman 1 Ralph M. Parsons Laboratory Department of Civil Engineering, Bldg. 48-335 Massachusetts Institute of Technology Cambridge, Massachusetts 02139 27. Dr. Donald J. O'Connor 1 Manhattan College Bronx, New York 10471 28. Mr. Sheldon M. Lazanoff 1 NAVOCEANO Rep. Fleet Numerical Weather Central Monterey, California 93940 29. Mr. Steven M. Gertz 1 Drexel Institute of Technology 32nd and Chestnut Streets Philadelphia, Pennsylvania 19104 151 30. Mr. M. J. Doyle, Jr. Department of Engineering Research Pacific Gas and Electric Company 4245 Hollis Emeryville, California 94608 31. Associate Professor Warren Denner, Code 58Dw Department of Oceanography Naval Postgraduate School Monterey, California 93940 32. Assistant Professor Jerry A. Gait, Code 58G1 Department of Oceanography Naval Postgraduate School Monterey, California 93940 33. Department of Oceanography, Code 58 Naval Postgraduate School Monterey, California 93940 152 Security Classification DOCUMENT CONTROL DATA -R&D {Security c]as si fie ation ot title, body of abstract and indexing annotation must be entered when the overall report is classified) 1 originating activity (Corporate author) Naval Postgraduate School Monterey, California 93940 Zb. report security classification Unclassified 2b. GROUP 3 REPOR T TITLE Numerical Simulation of Currents in Monterey Bay *- DESCRIPTIVE NOTES (Type o( report and, inclusive dates) Masfpr's Thpsis: March 19 71 S AUTHORiSI (Fifsf name, middle initial, last name) Roland Albert Garcia 6 REPOR T DATE March 19 71 7a. TOTAL NO. OF PAGES 154 7b. NO. OF REFS 6a. CONTRACT OR GRANT NO b. PROJEC T NO 9a. ORIGINATOR'S REPORT NUMBER(S) 9b. OTHKR REPORT NO(S) (Any other numbers that may be assigned this report) 10 DISTRIBUTION STATEMENT Approved for public release; distribution unlimited, 11. SUPPLEMENTARY NOTES 12. SPONSORING MILITARY ACTIVITY Naval Postgraduate School Monterey, California 93940 13. ABSTR AC T Recent interest in pollution control and the proximity of Monterey Bay to the Naval Postgraduate School prompted an investigation of the circu- lation in the bay. The first phase of the study consists of solving the simple cavity flow problem. A vorticity-stream function relationship is solved using an explicit, time dependent, finite difference shceme. Solutions for selected Reynolds' numbers and length to width ratios of the cavity are obtained. Values are chosen to give an indication of the flow patterns occurring over a wide range of these parameters. Equations for a refined model are derived to include the effects of the bottom topography, frictional forces and the Coriolis force. A numerical procedure similar to the one applied to the simple cavity flow problem is used on the refined equations. The topography of Monterey Bay is used in this study. Results indicate that closed circulation patterns are more in Monterey Bay than in other embayments due to the presence of the submarine canyon . DD.F. NOV «8 I *T / O S/N 010) -807-681 1 (PAGE 1) 153 Security Classification A-31408 Security Classification Model Numerical Simulation Bay Embay men t Circulation KEY WO R DS DD ,?„"vu..1473 (back) S/N 0101-807-6821 154 Security Classification A- 3 I 409 Thesis I263S2 G182 Garcia c.l Numerical simulation of currents in Monte- rey ty^ay. thesG182 Numerical simulaion of currents in Monte 3 2768 002 01041 5 i DUDLEY KNOX LIBRARY "