SELECTIVE WITHDRAWAL OF A LINEARLY STATIFIED FLUID IN A TRIANGULAR RESERVOIR By Stephen D. Hnidei B.Sc. University of Windsor 1988 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S T H E I N S T I T U T E O F A P P L I E D M A T H E M A T I C S We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A 1990 © Stephen D. Hnidei, 1990 In presenting this thesis in part ial fulfilment of the requirements for an advanced degree at the Universi ty of Br i t i sh Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of The Institute of Appl ied Mathematics The Universi ty of Br i t i sh Columbia 1956 M a i n M a l l Vancouver, Canada Date: In loving memory of my mother and grandfather. n List of Symbols C concentration of the impuri ty Cn wave speed of the nth mode cp specific heat D diffusion coefficient F Froude number ( Q / N # 2 ) g acceleration due to gravity G r Grashof number (N2L4v~2) H half height of the reservoir Ti Heaviside step function k Fourier heat conductivity coefficient L length of the reservoir m slope of the bot tom N Brunt-Vaisal la frequency (y/eg) p pressure P r P rand t l number (V/K) Q discharge per unit width R dimensionless parameter (F(Gr)s) t t ime T temperature i n —* u velocity vector u horizontal velocity component V vertical velocity component X horizontal coordinate y vertical cordinate 8 withdrawal layer thickness 6 measure of density stratification vort ici ty K molecular diffusivity A aspect ratio (L/2H = m _ 1 ) V- coefficient of viscosity V kinematic viscosity P perturbation in density Po mean density Pe ambient variation in density due to stratification PT total density streamfunction internal wave frequency IV Abstract The water in most reservoirs is density stratified wi th depth. This stratification leads to the inhibi t ion of vertical movement, consequently, when water is withdrawn from the reservoir it tends to move i n a jet-like layer called a withdrawal layer, towards the sink. B y placing the sink at a certain depth, one is able to selectively withdrawal water from a l imited range of depths and thus obtain water of a desired quality. M u c h work has been done in this field by considering a simplified boundary geometry, usually rectangular. However l i t t le attention has been given to the effects of accurate boundary geometry. For this thesis, five numerical experiments were conducted for the problem of a two-dimensional, viscous, incompressible, slightly-stratified flow towards a sink in a triangular reservoir. v Table of Contents List of Symbols iii Abstract v Acknowledgement xiv 1 Introduction 1 1.1 Density Stratified Flows 1 1.2 Previous Research 2 1.3 A Br ie f Summary of the Thesis 6 2 Equations of Mot ion 7 2.1 State Equations 8 2.2 Momen tum Equations 10 3 Numerical Calculations 17 3.1 Computat ional Details 17 3.2 Efficiency 20 3.3 Marker Particles and Interpolation 22 3.4 Updat ing Part icle Posit ion 23 3.5 Free Surface Considerations 24 3.6 Computat ional Stabil i ty and Accuracy 25 3.7 T i m e Step Dependance 26 v i 4 Computational Results 27 4.1 R u n l 31 4.2 Run2 35 4.3 Run3 36 4.4 Run4 37 4.5 Run5 38 5 Conclusions and Recomendation for Further Research 40 5.1 Conclusions 40 5.2 Further Research 41 Bibliography 43 A Stability Analysis 45 B The Program 48 C Figures 6 6 v i i List of Tables 1.1 F low regimes classified by Imberger, Thompson and Fandry 5 3.2 Over-relaxation coefficients 21 4.1 T ime taken for the first shear wave to pass selected points 28 4.2 Summary of varied quantities in the five experiments 29 v i n List of Figures 1.1 Part icle displacement and density profile 66 1.2 Withdrawal layer 66 1.3 Y i h ' s solution for F = 0.35 67 1.4 Kao's method for F < TT" 1 67 1.5 Potential flow 68 1.6 Withdrawal layer formation 68 2.1 Definition sketch of reservoir 69 2.2 Purely buoyant solution 69 3.1 Computat ional Domain 70 3.2 Staggered G r i d 71 3.3 Horizontal velocity interpolation 72 3.4 A Point Near The Sloping Boundary 72 3.5 Mirror-Image Technique 73 3.6 Updat ing Vor t ic i ty and Density 73 3.7 The location of the two points used for the grid testing 74 3.8 G r i d Dependance - coarse mesh at the point x — 0.2L and y = 0.2H . . . 74 3.9 G r i d Dependance - intermediate mesh at the point x = 0.2L and y = 0.2H 75 3.10 G r i d Dependance - fine mesh at the point x = 0.2L and y = 0.2H . . . . 75 3.11 G r i d Dependance - fine mesh at the point x = 0.2L and y = 0.6H . . . . 76 3.12 T i m e Step Dependance at x = Q.2L and y = 0.2H for a 60x120 grid . . . 76 4.1 The location of the nine points at y = 0.5H 77 4.2 The location of the seven points at y = 0.75H 77 i x 4.3 ip versus t at y = 0.5H for the points listed in Table 4.1 78 4.4 ip versus t at y = 0.75H for the points listed in Table 4.1 78 4.5 Propagation of the first shear wave 79 4.6 The location of the eleven points at x — 0.25L 79 4.7 tj) versus t at x = 0.25L, for equally spaced points starting at y = 0.283/7 (16,35) and ending at y = 0.7H (16,85) 80 4.8 A displaced spherical particle in a stratified fluid 80 4.9 Lines of constant density near the sloping bottom 81 4.10 ip versus t for points near the sloping bottom 81 4.11 Streamlines at t = 0.0 seconds, R u n l ; . . . . 82 4.12 Streamlines at t = 5.0 seconds, R u n l 82 4.13 Streamlines at t = 10.0 seconds, R u n l 83 4.14 Streamlines at t = 20.0 seconds, R u n l 83 4.15 Streamlines at t = 30.0 seconds, R u n l 84 4.16 Streamlines at t — 40.0 seconds, R u n l . . 84 4.17 Streamlines at t — 60.0 seconds, R u n l 85 4.18 Streamlines at t = 120.0 seconds, R u n l 85 4.19 Definition sketch for 6+ and S~ 86 4.20 Horizontal velocity profile at t = 15.0 seconds at x = 0.25L 86 4.21 Horizontal velocity (t = 20.0 s, x = 0.251), R u n l 87 4.22 Horizontal velocity (t = 35.0 s, x = 0.25L), R u n l 87 4.23 Horizontal velocity (t = 45.0 s, x = 0.25X), R u n l 88 4.24 Horizontal velocity (t = 50.0 s, x = 0.25L), R u n l 88 4.25 Horizontal velocity (t = 55.0 s, x = 0.25L), R u n l 89 4.26 Horizontal velocity (t = 85.0 s, x = 0.25L), R u n l 89 4.27 Horizontal velocity (t = 120.0 s, x = 0.251), R u n l 90 x 4.28 Horizontal velocity (t = 25.0 s, x = 0.5L), R u n l 90 4.29 Horizontal velocity ( i = 35.0 s, x = 0.5L), R u n l 91 4.30 Horizontal velocity (t = 40.0 s, x = 0.5L), R u n l 91 4.31 Horizontal velocity (t = 50.0 s, x = 0.51), R u n l 92 4.32 Horizontal velocity (t = 60.0 s, x = 0 .5L), R u n l 92 4.33 Horizontal velocity (t = 75.0 s, x = 0.51), R u n l 93 4.34 Horizontal velocity (t = 120.0 s, x = 0.57L), R u n l 93 4.35 Peak horizontal velocity at 3 stations, R u n l 94 4.36 Peak horizontal velocity at 3 stations, R u n l 94 4.37 Peak horizontal velocity at 3 stations, R u n l 95 4.38 S+ at 3 stations, R u n l 95 4.39 6+ at 3 stations, R u n l 96 4.40 6+ at 3 stations, R u n l 96 4.41 Average peak velocity, R u n l 97 4.42 Wi thdrawn particles in i t ia l position after 120 s, R u n l 97 4.43 Streamlines at t = 40.0 seconds, Run2 98 4.44 Streamlines at t = 80.0 seconds, Run2 98 4.45 Streamlines at t = 120.0 seconds, Run2 \ . 99 4.46 ip versus t at y = O.hH for the points listed in Table 4.1, Run2 99 4.47 Averaged peak horizontal velocity comparison between flows wi th TV = 0, 0.0764 and 0.764, R u n l and Run2 100 4.48 Streamlines at t = 5.0 seconds, Run3 101 4.49 Streamlines at t = 20.0 seconds, Run3 101 4.50 Streamlines at t = 120.0 seconds, Run3 102 4.51 Horizontal velocity (t = 15.0 s, x = 0.25L), Run3 102 4.52 Horizontal velocity (t = 35.0 s, x = 0.25L), Run3 103 xi 4.53 Horizontal velocity (t = 120.0 s, x = 0.251), Run3 103 4.54 Horizontal velocity (t = 35.0 s, x = 0.51), Run3 104 4.55 Horizontal velocity (t = 40.0 s, x = 0.51), Run3 104 4.56 Horizontal velocity (t = 120.0 s, x = 0.5X), Run3 105 4.57 Averaged peak horizontal velocity comparison between R u n l and Run3 . 105 4.58 Averaged 6+ comparison between R u n l and Run3 106 4.59 Streamlines at t = 0.0 seconds, Run4 106 4.60 Streamlines at t — 5.0 seconds, Run4 107 4.61 Horizontal velocity (t = 5.0 s, x = 0.25Z), Run4 107 4.62 Horizontal velocity (t = 5.0 s, x = 0.5L), Run4 108 4.63 Streamlines at t = 10.0 seconds, Run4 108 4.64 Horizontal velocity (t = 10.0 s, x = 0.25L), Run4 109 4.65 Horizontal velocity (t = 10.0 s, x = 0.5X), Run4 109 4.66 Streamlines at t = 15.0 seconds, Run4 110 4.67 Horizontal velocity (t = 15.0 s, x = 0.25Z), Run4 110 4.68 Horizontal velocity (t = 15.0 s, x = 0.5Z), Run4 I l l 4.69 Streamlines at t = 20.0 seconds, Run4 I l l 4.70 Horizontal velocity (t = 20.0 s, x = 0.25Z,), Run4 112 4.71 Horizontal velocity (t = 20.0 s, x = 0.5L), Run4 112 4.72 Horizontal velocity (t = 30.0 s, x = 0.25X), Run4 113 4.73 Horizontal velocity ( i = 120.0 s, x = 0.251), Run4 113 4.74 Horizontal velocity {t = 25.0 s, x = 0.5X), Run4 114 4.75 Horizontal velocity (t = 30.0 s, x = 0.5X), Run4 114 4.76 Horizontal velocity (t = 35.0 s, x = 0 .5L), Run4 115 4.77 Horizontal velocity (t = 45.0 s, x = 0.5L), Run4 . . . 115 4.78 Horizontal velocity (t = 55.0 s, x = 0.5L), Run4 116 x i i 4.79 Horizontal velocity (t = 70.0 s, x = Q.5L), Run4 116 4.80 Horizontal velocity (t = 75.0 s, x = 0.57:), Run4 117 4.81 Horizontal velocity (t = 95.0 s, x = 0.51), Run4 117 4.82 Horizontal velocity (t = 120.0 s, x = 0 .5L), Run4 118 4.83 Horizontal velocity (t = 40.0 s, x = 0.25L), Run5 118 4.84 Horizontal velocity (t = 50.0 s, x = 0.25L), Run5 119 4.85 Horizontal velocity (t = 65.0 s, x = 0.25L), Run5 119 4.86 Horizontal velocity (t = 120.0 s, x = 0.25X), Run5 120 4.87 Horizontal velocity (t = 25.0 s, x = 0.5X), Run5 120 4.88 Horizontal velocity (t = 50.0 s, x = 0 .5£) , Run5 121 4.89 Horizontal velocity (t = 65.0 s, x = 0.5Z), Run5 121 4.90 Horizontal velocity (t = 85.0 s, x = 0.5X), Run5 122 4.91 Horizontal velocity (t = 100.0 s, x = 0.5X), Run5 122 4.92 Horizontal velocity (t = 115.0 s, x = 0.5L), Run5 123 4.93 Horizontal velocity (t = 120.0 s, x = 0.5L), Run5 123 4.94 8+ at 3 stations, Run5 124 4.95 8+ at 3 stations, Run5 124 4.96 8+ at 3 stations, Run5 125 4.97 Average 8+ comparisons between Runs 1,4 and 5 125 4.98 8+ at t = 120.0 seconds for A = 1, 4 and 10 126 x i i i Acknowledgement I would like to express my appreciation to my supervisor Dr . Gregory Lawrence of the C i v i l Engineering department for introducing me to this topic and for his many helpful suggestions during my research on this thesis. I would also like to acknowledge the following: Dr . W i l l i a m Hsieh of the Oceanography department for his careful reading of my early manuscript, the secretaries in the Mathematics department for doing what they do, Anthony Weaver and Owen Walsh for their different insights into fluid mechanics and Loren Daltoe for her motherly support throughout. x i v Chapter 1 Introduction 1.1 Density Stratified Flows M a n y bodies of water exhibit a vertical variation of density due to changes in temperature or concentration of dissolved solids, usually salt. This density variation can often be approximated as linear over a significant portion of the depth. There are many naturally occuring examples of fluid flow in which density stratification plays a key role. Examples include the motion of cold and warm fronts in the atmosphere, which are density stratified due to temperature, and ocean currents that may have been density stratified due to both temperature and salinity. Another example of stratified flow is called Selective Withdrawal , which is the topic of this thesis. In a stratified flow which is statically stable, work must be done to displace a fluid particle vertically since there is a force due to gravity always tending to oppose this movement. The density gradient tends to inhibit vertical motion in favour of horizontal motion. We can see this i f we assume density decreases wi th height (dp/dy < 0) and therefore the stratification is stable. However if the stratification were to be unstable, in a l l practicali ty it would eventually turn over and become stable. Second, consider a fluid particle at a height y2, wi th density p(y2), volume V and experiencing gravitational acceleration g. Now if we move the fluid particle to a new height y\ the work done is which is positive regardless whether y 2 > V\ o r Vi < V\- See Figure 1.1. Now if the 1 Chapter 1. Introduction 2 fluid where not stratified, then p(y) would be a constant and thus the work done would be equal to zero. This is the main idea of Selective Withdrawal . Consider a reservoir wi th a sink (water intake) located at a certain depth below the surface. If the water were not stratified, water from different depths would move radially towards the sink since there is no work involved. In the case of a stratified fluid water w i l l tend to come mostly from one depth or layer, called a withdrawal layer because vertical motion tends to be inhibited as mentioned earlier. This is depicted in Figure 1.2. In this way water quality can be controlled. B y placing the sink at a particular depth, water of a certain temperature, dissolved oxygen content or foreign mineral content can be withdrawn and thus the name Selective Withdrawal . B y using Selective Withdrawal , one can withdraw water of a certain temperature or oxygen content for drinking purposes or for irrigation. One can also use Selective Withdrawal for the removal of salt water that has encroached upon a supply of fresh water by using a scour outlet located near the base of the dam. 1.2 Previous Research The first investigation of two-dimensional linearly density stratified flow towards a sink was done by Y i h (1958). He assumed the sink was located in the bottom corner of a semi-infinite channel and that the velocity and density distributions at infinity were undisturbed. He obtained the steady-state solution in terms of a Fourier series for values of the densimetric Froude number F > 7 r _ 1 , where F = Q/NH2. Q is the discharge rate, N is the Brunt-Vaisal la frequency and H is the height of the reservoir. The Brunt-Vaisal la frequency is defined as N — y/eg, where e = —p^^dpe/dy), where po is the mean density of the fluid and pe is the ambient variation in density due to the stratification of the fluid. The solution, shown in Figure 1.3, is almost uniform flow towards the sink with a corner eddy above the sink. For values of F decreasing towards 7 r - 1 , the corner eddy Chapter 1. Introduction 3 extends horizontally upstream unt i l F = 7 r _ 1 where we have a returning flow at infinity, violating the upstream condition. Debbler (1959), using an 18 ft. long flume tried to experimentally verify Y i h ' s results. He filled the flume wi th 15 different layers of saltwater of different densities. Each layer placed in the flume was progressively less dense. Al ternat ing layers where then dyed so as to observe the flow pattern during the emptying of the flume. Before conducting his experiment, he waited eighteen hours, so as to obtain a linear density distribution by diffusion. He qualitatively showed that the flow pattern near the cr i t ical Froude number was a region of stagnant water near the top wi th nearly uniform flow towards the sink i n the bot tom part of the flow field. He found the cri t ical Froude number to be 0.28. K a o (1965) obtained a solution for F < T T - 1 by introducing a uniform sink distribution as shown in Figure 1.4. The flow field exhibited a dividing streamline in which fluid below leaves through the bottom corner sink and fluid above leaves through the distribution of sinks. Then by making the pressure on the dividing streamline equal to the hydrostatic pressure of the fluid above i t , the dividing streamline separates the flow leaving through the corner sink from the now stagnant upper region. He obtained a cri t ical value of the Froude number to be 0.345 whereas 7 r _ 1 « 0.318. K o h (1966) obtained a steady-state solution which incorporated diffusive and viscous effects, but which neglected the nonlinear terms in the governing equations by making boundary-layer like assumptions. His solution predicted the withdrawal layer grows in thickness wi th distance x from the sink proportional to x^ and that the horizontal velocity profile was self-similar wi th the maximum horizontal velocity at any station proportional to x~z. Pao and K a o (1974) using a horizontal-duct model were able to identify the successive propagation of a wave-like disturbance. They called these wave- like disturbance modes "columnar disturbances" or "columnar waves of zero frequency" because they travelled Chapter 1. Introduction 4 upstream in horizontal jet-like columns wi th various modes travelling at a velocity Cn = Nd/rnr, where d is the channel depth. They identified these waves as the mechanism responsible for the flow concentrating into a withdrawal layer. K a o , Pao and Wei (1974) also conducted experiments using a 30.5 ft. long flume that was linearly stratified wi th a salt-water solution. They observed the successive arrival of "columnar disturbance modes". Theoretical and experimental work was conducted by M c E w a n and Baines (1974) in a rectangular reservoir that was filled wi th a continuously stratified fluid. The shear waves in their study were created by impulsively setting the end walls i n a uniform shear motion. They found that the shear waves travelled at a speed of Nd/mr wi th a front that widens wi th time. The front width was found to be (Nt)^dn~l, where t is time. Imberger, Thompson and Fandry (1976) obtained a solution for the flow in a finite two-dimensional tank, that was density stratified and had a free surface. They showed that such flows depend on the value of R = F(Gr)^, where Gr = N2L4i/~2. The flow dependence on R is summarized in Table 1.1. The shear waves they discovered are identical to those described by Pao and K a o (1974) and M c E w a n and Baines (1974). U p o n ini t ia t ion of the flow, vorticity must be zero everywhere in the flow field and we have potential flow as shown in Figure 1.5. Flow is radial towards the sink wi th large vertical velocity components near the sink. This causes much fluid to be moved from its place of neutral buoyancy. The restoring effect is the creation of these shear waves. We note that at mid-depth each shear wave that passes adds to the horizontal velocity towards the sink. Above and below the sink, different modes wi l l tend to negate each others horizontal velocity component and this w i l l leave the highest velocities at mid-depth and thus create a withdrawal layer. See Figure 1.6, from Silvester (1979). If the rear wall of the reservoir is vertical, the shear waves are reflected but continue to propogate horizontally. However, we seldom expect a real reservoir to be of a rectangular Chapter 1. Introduction 5 Time scale Comments Flow establishment N'1 Time scale for internal waves. # _ 1A _ 1 Time taken for the 1st shsar wave to reach end. TV - U - * Time taken for 1st mode to expand to a standing wave. Approach to steady state (a) R < Pr-i A , - 1 Gri Initial withdrawal-layer formation. Layer collapses owing to convection of in situ gradient. S = 0{LOr~i). N^GriPri Pseudo-steady-state withdrawal as described by §3. S = O(LGr-lPr-l). N-iQrlPr-ip-1 Further collapse due to self-induced convection. (b) Pr-i < P. < P r - i A7"1Gri Initial withdrawal formation. S — O(LQr-l). N~lGrlR-l Self-induced convection accelerates collapse. S = O(LGr-lRl). N^GriPr^ Steady state, as described by Imberger (1972), is reached; viscosity dominates. Ljxe > 1. (c) Pr-i < R < 1 N-iQri Initial withdrawal formation. 6 = 0(LGh—l). N-iQrlR-l Self-induced convection accelerates collapse. S = O(LGr-lRl). N-iF-iR-i Steady state, as described by Imberger (1972), is reached; inertia dominates. Ljxc < 1. (d) R > 1 N-ip-i Inertia dominates and solution behaves as described in §3. S = O(LFl). Large times N'WGrl Time taken for all wave motion to decay. LHQ-1 Time taken for tank to empty. The hierachy of time scales. 6 is the characteristic scale of the withdrawal-layer thickness. N~xOri is the time taken for a viscous-buoyancy layer to form, or the time taken for waves of wavelength smaller than LOr'i to decay, or the time taken for mode n = XOri to reach the end of the tank. N^GrlR is the time taken for the inertial zone to be set up. N-ijp-t i s the time needed for steady state to be achieved, or the time for wave n = XF~i to reach the end of the tank, or the time for fluid to fall through a distance LFi. Table 1.1: F low regimes classified by Imberger, Thompson and Fandry Chapter 1. Introduction 6 shape. We would expect the bottom of the reservoir to be gently sloping as shown in Figure 1.2. Cacchione and Wunsch (1974) have investigated the effects of a slope on internal wave motion i n a continuously stratified fluid. They showed the wave motion falls into one of three categories depending on the magnitude the of slope m compared to the parameter c, where c 2 = u2/ (TV2 — u2) and u is the internal wave frequency. When m < c, absorption of the internal waves into the acute corner was noted. When m > c, the internal waves tended to be reflected back to the interior of the fluid domain and when m = c, the internal waves are heavily damped. In light of this result, we would expect the sloping bot tom to somehow alter the selective withdrawal process. To examine this more closely, this thesis w i l l investigate the effects of the selective withdrawal process in a triangular reservoir. 1.3 A Brief Summary of the Thesis We shall now briefly summarize the scope of this thesis. In Chapter 2 the governing equations describing a density stratified flow are derived. The stratifying agent is assumed to be either temperature or salinity, but not both. Furthermore we w i l l assume the flow satisfies the Boussinesq approximation. For insight into the flow, the governing equations are solved exactly for a simplified l imi t ing steady-state case when both diffusive and convective terms are neglected. The computational procedure is described in Chapter 3. A n explanation of how the Marker- in-Cel l method and Seidel's iteration method are employed in the computer program is given. To increase the efficiency of the program, over-relaxation is used to help obtain quick convergence. Also to increase efficiency an Alternat ing Direct ion method is used to sweep the field. Chapter 4 deals wi th the results obtained by implementing the computer program. We first test the validity of the model. The rest of the chapter analyses the five numerical experiments. Chapter 2 Equations of Motion In dealing wi th a fluid, we must consider that the Principle of Mass Conservation not be violated and thus the velocity distribution of the fluid must satisfy what is known as the continuity equation. For a homogeneous fluid the Equat ion of Cont inui ty is dt + V • (pTu) = 0. (2.1) where px is the density of the fluid and u is the velocity vector. If we consider a density stratified fluid and the stratification is due to temperature then (2.1) remains val id. Now if the fluid is non-homogeneous ( there exist impurities in the fluid, such as salt in our case ) then we must consider molecular diffusion of this foreign substance because there wi l l be an additional transfer of mass. Fick ' s Law of Diffusion is: net outflux of mass = - j (DVC) • n dS (2.2) where C is the concentration of the impurity, D is the diffusion coefficient, S is the closed surface and n is the unit outward normal. 7 Chapter 2. Equations of Motion 8 B y the Divergence thereom we may write (2.2) as -Js(DVC)-n dS = - JvV • (DVC) dV. (2.3) Wri t ing (2.1) as a volume integral and using (2.2) and (2.3) we have that 'v Since V is arbitrary, we can write dpr dt This can be rewritten as DpT **L + V • (pTu) - V • (DVC) dV = 0. + V • (PTu) - V • (DVC) = 0. Dt + pTV • u = V • (DVC), (2.4) D 0 _ „ , where rTt = m + {u-v) and -j^ is called the material derivative. 2.1 State Equations For a homogeneous flow we know that p? — constant. However, for a stratified fluid we are not so lucky and this simple relation does not hold. We must consider the following relations. If the stratification is due to temperature, then from Fourier's Law of Heat Conduction we have that DT Dt = V - ( / c V T ) , (2.5) where T is temperature, k prep k is the Fourier heat conductivity coefficient, and c p is specific heat. Chapter 2. Equations of Motion 9 If the stratification is due to some foreign substance then from Conservation of Mass, we have that DC — = V - {DVC) • (2.6) Now if we assume that density variations are small (Apr/pr *C 1), so that K ~ constant and D ~ constant,then we can write (2.5) and (2.6) as DT Dt DC = «V 2r, (2.7) DVZC. (2.8) Dt To relate these equations to the continuity equation we need a suitable relation be-tween p and the two new variables C and T . If we assume a linear relation between them exists such that pT - po = a(T - T0) + b(C - C0) where p0, To, Co, a and b are constants. To simplify things even more we wi l l assume that the stratification wi l l be caused by only one of these two possible mechanisms and not both. In the case where stratification is caused by temperature, we have that T = ^ I _ ^ + T o . (2.9) a a Since (2.1) holds for this case, we can rewrite it as D p T + pT^ • u = 0. (2.10) Dt Also, if we substitute (2.9) into (2.7) we get DpT Dt r Chapter 2. Equations of Motion 10 Similar i ly, we see that when the stratification is caused by a foreign substance we have that C = P - T - T + C°- (2 - 1 1 ) o o W i t h a l l the new assumptions and upon substituting (2.11), (2.4) can be writ ten as ^ + PTV-U = JV2PT. (2.12) Substi tuting (2.11) into (2.8) we get D P T = DV'PT. (2.13) Dt Rewri t ing (2.12), we have that b r ± Dt U p o n substituting (2.13) into (2.14) we then have > r V - i l = ^ - % - (2-14) ™ = ( i-0^£. (2-15) .b J pT Dt From (2.10) and (2.15), we see that both equations are of the form V • u = (constant) — ( 2 . 1 6 ) pT Dt Since we have assumed that the density variations are small , we have for our Cont i -nuity Equat ion a good approximation to be V - u = 0. (2.17) 2.2 Momentum Equations In general we have that Du PTJ^ = -VPT + V-T + PTG, (2.18) Chapter 2. Equations of Motion 11 where px is the total pressure at a point, r is a stress tensor —• and G is an external body force. For a Newtonian fluid, Stoke's relation between stress and strain is valid. It is T = n Vt? + ( V u ) T l + A 7 V • u , (2.19) where p is the coefficient of viscosity, 2 A + —u is called the bulk or volume viscosity o and I is the identity tensor. Substituting (2.17) into (2.19), we have that T = / / [ V U + ( V U ) t ] , (2.20) Substi tuting (2.20) into (2.18) and assuming p to be constant gives Du In our case the external body force is that due to gravity, so G = (0, — g). So the governing equations of a density stratified non-homogeneous fluid are (2.13), (2.17) and (2.21). -Dt = D V PT' V • u = 0, Chapter 2. Equations of Motion 12 Du 2 -> p r — = -VpT + pW u + pTG. The constant D in the first equation above is either the diffusion coefficient or heat diffusivity depending on what mechanism caused the density stratification. For convenience we will decompose the total density PT, of the stratified fluid as PT (X, y, t) = po + pe (y) + p (x, y, t) where po is the mean density, pe is the ambient variation in density due to stratification and p is the perturbation in density due to the motion of the fluid. Similarly the total pressure at a point PT, can be decomposed as PT ( Z , y, t) = pe (y) + p (ar, y, t) where pe is the ambient pressure and p is the perturbation in pressure due to the motion of the fluid. We now can replace px and p? by their decompositions and substitute them into our governing equations after realizing ry Pe (y) = Pe (y) 9 dy Jo and we get ^ - ep0v = DV2P, Chapter 2. Equations of Motion 13 where u = (u, v) and V • u = 0, Po^ = -Vp + pV2u + pG, 1 dpe e = Po dy ' where e is a measure of the stratification of the fluid. Also PT~fr~ has been replaced by p0—-. This is called the Boussinesq approximation and is valid since the variation in the density throughout the flow is assumed small. The above three equations writ ten out explici t ly i n two-dimensional cartesian coor-dinates are: Pt + upx + vpy - ep0v = D [pxx + pyy), (2.22) ux + vy = 0, (2.23) p0 [ut + uux + vuy] = —px + u [uxx + uyy]. (2.24) Po [vt + uvx + vvy) — -py - pg + p [vxx + vyy] (2.25) The presence of (2.23) suggests that we introduce a stream function ip such that u = ipy and v = — ipx. Chapter 2. Equations of Motion 14 In this way (2.23) w i l l automatically be satisfied. Now we wi l l define vorticity as £ = V V = uy ~ vx. (2.26) Then differentiating (2.24) wi th respect to y and subtracting (2.25) after it has been differentiated wi th respect to x and then using (2.26), we get 6 + + 1>£v = — + " [U + tyy] (2.27) Po where v = ppo1 . Alternat ively (2.22) and (2.27) can be writ ten as ^ + £ft,V>. = £ V 2 p , (2.28) § = J L ^ + „VV. (2.29) Dt podx If we omit the diffusive and convective terms from (2.28) and (2.29) we obtain Pt + tpo^x = 0, (2.30) 6 = —px. (2.31) Po Differentiating (2.30) wi th respect to x and (2.31) wi th respect to t and then eliminating the pxt term yields the following equation We w i l l define the Brunt-Vaisal la frequency as N = (eg)^. Using (2.26) the equation becomes 1>xxtt + i>yytt + N2ipxx = 0. (2.32) If we now consider the steady-state case of (2.32), we have i>xx = 0. (2.33) Chapter 2. Equations of Motion 15 The solution to (2.33) subject to the appropriate boundary conditions w i l l be called the Pure ly Buoyant Solution. This solution is important because we expect after sufficient t ime that the flow w i l l adjust to this state. The following wi l l be the boundary conditions used for the triangular reservoir with mid-depth withdrawal as shown in Figure 2.1: ip = Q[l-H(y)] on x = 0, (2.34) xp = ^ on y = H, (2.35) xp = Q on y = mx — H, (2.36) where H.(y) is the Heaviside function and m = ^ . In view of (2.33), we can readily see the general solution is 1>(x,y) = f(v)x + 9(y), (2-37) where f(y) and g(y) are arbitrary functions of their arguments. Substituting (2.34) into (2.37) yields g(y) = Q[l - 7i(y)]. Therefore, ^(x, y) = f(y)x + Q[l - H(y)}. (2.38) Substi tuting (2.35) into (2.38) yields f = f(H). Substituting (2.36) into (2.38) yields QH(mx -H) = f(mx - H)x. Chapter 2. Equations of Motion 16 Therefore, f(y) must be inversely proportional to x when y > 0 and f(y) = 0 when y < 0. The previous observations yield /(„) = g L n y ) and therefore ^ y) = + - « ( y ) ] . (2-39) The purely buoyant solution is shown in Figure 2.2. Note that the potential flow solution, Figure 1.5, and the buoyant solution are solutions to the governing equations under two entirely different sets of conditions. We ini t ia l ly start wi th potential flow. This is successively modified by each shear wave that passes and we would expect the flow to adjust towards the flow pattern exhibited by the buoyant solution after some time. The t ime taken for the shear waves to decay, due to viscous effects maybe much longer than the time to empty the reservoir, so the flow may not obtain the flow pattern described by the purely buoyant solution. One must also recall that we have neglected diffusive and convective effects in obtaining the buoyant solution. Keeping this i n mind, one sees from the purely buoyant solution, that there exists everywhere above mid-depth a forward flow, i.e., a positive horizontal component of flow velocity. Whereas the purely buoyant solution for the rectangular case has a zero horizontal velocity component above mid-depth. Chapter 3 Numerical Calculations 3.1 Computational Details This chapter describes how the equations that govern the motion of a density stratified fluid in a triangular reservoir can be solved using a computer. The first step is to divide the computational domain into a series of small rectangular regions called cells. For our particular problem the computational domain was divided so that there were N X cells in the horizontal direction and N Y cells in the vertical direction. Each cell was D X long and D Y high, so that we had uniform grid spacing in both the x direction and the y direction. However, there were twice as many cells in the y direction as in the x direction, so as to facilitate better resolution of the withdrawal layer structure. See Figure 3.1. In each of these computational cells, physical quantities of the fluid such as velocity, vorticity, density, etc. are calculated at different locations wi th in the cell. This technique is commonly refered to as the staggered grid or staggered node method. To see the implimentat ion of the staggered grid see Figure 3.2. The relative positions of PSI(I,J) and U(I,J) were chosen so as to stagger U(I,J) between PSI(I ,J) and PSI(I ,J+1) since u — ipy. The positioning of V(I , J ) and Z(I,J) is analogous. This implementation is also convenient for imposing the boundary conditions. Now that the computational domain has been divided and the staggered grid has been formed, we must take our governing equations and discretize them using what is known as finite-difference approximation. To illustrate how we obtain our descretized 17 Chapter 3. Numerical Calculations 18 version of our governing equations, consider (2.26) at the cell numbered (I,J). We have that d2PSI(I,J) d2PSI(I,J) + = Z(I,J), dx2 ' dy2 where PSI and Z are our discrete representations of ip and £ for our computer program. Now let us consider the first term in this equation. The treatment of the second term wi l l be analogous. Consider the Taylor Series expansion of P S I ( X 0 + AX, Y0) about the point (X0, Y0). We have PSI(X0 + AX, Y0) = PSI(X0, Y0) + dPSI(X, Y) + d2PSI(X, Y) dX2 (X0,Y0) (Axy 2! + ...+ dX dnPSI{X, Y) dX* AX (Xo,Y0) (AX)a (r,Y0) n! where XQ < r < (Xo + AX). Upon rearrangement we have 8PSI(X, Y) dX PSI(X0 + AX, Y0) - PSI(X0, Y0) d2PSI(X, Y) AX dX2 AX (X0,Yo) 2! + 0(Ax), (Xo,Y0) The Forward Difference Approximat ion of ipx becomes dPSI(I,J) _ PSI(I+l,J)- PSI(I,J) dx Ax where O(Ax) = K \ Ax\ for some K G IR. The Backwards Difference Approximat ion of ^becomes dPSI(I, J) _ PSI(I,J)- PSI(I - 1, J) dx Ax If we add the Taylor Series expansions for PSI(X0 + AX, Y0) to PSI(X0 - AX, Y0), we get the following equation for the second derivative a2psi(i, j) PSI{I +1, j) - 2PSI{I, j) + PSI(I -1, J) + 0(Ax). dx2 (AXy + 0(Ax2). (3.1) Chapter 3. Numerical Calculations 19 Similarly, we can obtain an approximation for xpyy , d2PS1(1, J) _ PS 1(1, J + 1) - 2PSI(I, J) + PSI(I, J-l) dy2 ~ (Ay)2 Using (2.26), (3.1) and (3.2) we obtain (3.2) PSI(I,J)=± PSI(I+1,J) + PSI(I-1,J) DX2 PSI(I,J+1) + PSI(I,J -1) (3.3) DY2 where - Z(I, J) B2 = 2(-±-+ 1 .DX2 DY2 Now if we write (3.3) in an iterative form we have P 5 / ( / , j ) ( n ) = PSI(I,J)^ -HLH D2 1 PSI(I + 1, J ) ( n _ 1 ) - 2PSI(I, J)( n- a) + PSI(I - 1, J ) ^ - 1 ) DX2 PSI(I, J + l ) ^ " 1 ) - 2PSI(I, J)*"" 1) + PSI(I, J - l ) ^ " 1 ) , (3.4) DY2 where PSI(I, J)^ is an initial guess field and then iteratively we calculate PSI(I, J)(n\ After many iterations we expect Jim [PSI(I, J)W - PSI(I, J)^" 1)] 0. Therefore the term in brackets is the difference between iterations or can be thought of as the error term and can be denoted by 75i?(n_1\ So we have then that PSI(I, J ) ( n ) = PSI(I, J ) ( n _ 1 ) + ^ER(I, J){n-X\ D2 We can calculate PSI(I,J) to any desired accuracy by making sure that the error term ER(I,J) is less than some fixed tolerance before we finish the iteration process. Chapter 3. Numerical Calculations 20 We should now consider the case where PSI(I+1,J) is below the sloping boundary as shown in Figure 3.4. The value of the stream function at the point (I + | , J ) is known to be exactly Q. It therefore would be better to use that point in our differencing formula. Repeating the same argument as was done earlier we obtain the following d2PSI(I, J) 2 [lPSI(I - 1, J ) - 2PSI(I, J) + IQ] dx2 DX2 d2PSI(I,J) _ Q - 2PSI(I, J) + PSI(I, J + 1) dy2 ~ DY2 and an iterative formula (3.5) (3.6) PSI(I, J)<n> = PSI(I, J ) ( n _ 1 ) + -jr-ER(I, J)*"" 1 ) , (3.7) where D3 = 2 + ^-). 3.2 Efficiency Note we must calculate PS 1(1, J ) ( k + 1 ) for a l l (I,J) i n our computational domain and also note that PSI(I, J )< k + 1 ) is related to PSI(I + 1, J)^,PSI(I - 1, J )< k \ PSI(I, J + 1 ) « and PSI(I,J — l )( k ) through (3.4) and (3.7). For any of these quantities, i f we hap-pen to know its value at the k+1 iteration, we should use this new value in calculating PSI(I, J ) ( k + 1 ) . This method is called Seidel's Method. This w i l l increase the efficiency of the program. Since we must solve for PSI(I ,J) for a l l (I,J) in the computational domain, the order in which we solve (or sweep the field) is important. If we solve for PSI(I,J) in an order that uses as many of these new values as possible, Forsythe and Wasow (1960) have shown this can double efficiency. The method so employed is called the Alternat-ing Direct ion Method. The order i n which the PSI(I,J) 's are calculated is given by the following: a) Sweep the computational domain from the bottom to the surface and from left to Chapter 3. Numerical Calculations 21 m 1.0 0.5 0.2 0.1 o^pt 1.642 1.698 1.747 1.769 Table 3.2: Over-relaxation coefficients right. b) Sweep the computational domain from the bottom to the surface and from right to left. c) Sweep the computational domain from the surface to the bottom and from left to right. d) Sweep the computational domain from the surface to the bottom and from right to left. Another way to increase the efficiency of calculating PSI(I ,J) iteratively is to reduce the difference between iterates as quickly as possible. If we denote this residual by R(I, J ) ( n ) = ^r-ER(I, J ) ( n ) or ^-ER(I, J ) W depending i f we are near the the sloping boundary or not and consider PSI(I, J)( n + 1) = PSI(I, J ) ( n ) + uR(I, J) ( n ). Forsythe and Wasow (1960) have shown that for 0 < to < 2 the above scheme wi l l converge. This is called the Method of Relaxation. If cu = 1 relaxation is absent. If 0 < LO < 1 then we have under-relaxation. If 1 < co < 2 then we have over-relaxation. Under-relaxation is used when the iterates PSI(I, J)^ oscillate towards PSI(I, J)(n\ Over-relaxation is used when the iterates are monotone increasing or decreasing. Lawerence (1976) has shown that for different height-to-depth ratio's (m) the optimal choice of a; is as shown in Table 3.1. Chapter 3. Numerical Calculations 22 3.3 Marker Particles and Interpolation U p unt i l now, we have only discussed the Euler ian fields ( ie. velocity ( U , V ), density ( R ) and vort ici ty ( Z ) ). Also used is a set of Lagrangian marker particles. The marker particles move according to the velocity components in their vicinity. These particles are introduced to help visualize the flow of fluid parcels. They show which cells have fluid and which cells lie on the surface. The position of the k t h particle is denoted by it's center point X ( k ) , Y ( k ) . Each particle is considered to occupy a space D X by D Y . Ini t ial ly the particles are positioned uniformly through the reservoir at a spacing of D X / 2 and D Y / 2 . A s the flow is initiated, the particles w i l l move, but for s implici ty the motion wi l l be purely translational. Since there is no rotation, each particle is uniquely defined by it's center point X ( k ) , Y ( k ) . W i t h each particle there is a corresponding vort ici ty Z E T A ( k ) and density R H O ( k ) . The particles and the fields are associated by the use of bilinear interpolation. For this interpolation we need the velocities at four nearest grid points near the particle. To help understand this interpolation, consider the horizontal velocity component, U I N T , of the k t h particle as shown in Figure 3.3. The horizontal velocity components interpolated at the points A and B as shown are: rr DX XI T T / r -r\ XI RT, _ T. U A = DX ( ' } + ~DX + ' J ) ' DX — XT XT U B = DX U { H J P ) + D J C U ( I + 1 ; J P ) ' where X I and E T A are as shown in Figure 3.3. Now using the velocity components at A and B we can define U I N T as: If we define Cl = C2 = ^ r , C 3 = 1 - Cl and CA = 1 - C 2 then UINT = C72*C3*/7(7, JP)+C2*Cl*U(I+l, JP)+C4*C3*U(I, J)+C4*C1*U(I+1, J). Chapter 3. Numerical Calculations 23 Near the sloping boundary this presents problems. Consider a particle at the point A near the sloping boundary as shown in Figure 3.4. We can see from Figure 3.4 that V(I+1,J) and V ( I + l , J - f 1) are necessary for bilinear interpolation. To calculate the imaginary values for U and V below the boundary we wi l l assign imaginary values to P S I below the boundary. To ensure that no particles move through the boundary we wi l l use a Mirror-Image Technique. We w i l l define PSI for points below the boundary to have the value PSI(I ,J) = 2Q - PSI(J , I ) . See Figure 3.5. W i t h this accomplished we can define U and V . F rom the definition of the stream function, we have that and V = - £ . oy ox Using Finite-Differences, as discussed earlier, we obtain: • U(I ,J) = [ PSI(I ,J+1) - PSI(I ,J) ] / D Y , • V( I , J ) = [ PSI(I ,J) - PSI(I+1,J) ] / D X . The governing equation (2.29) is computer coded by the following: • ZT( I , J ) = - G V R O * R X + V I S C O S * ( Z X X + Z Z Z ). The terms R X , R X X , R Z Z , Z X X and Z Z Z represent the derivatives of R(I ,J) and Z(I,J) and are calculated similarly to (3.1). Equation (2.28) is computer coded likewise. 3.4 Updat ing Particle Position The kth marker particle and its properties are updated as follows: a) Calculate the velocity at X ( k ) , Y ( k ) . b) Calculate where the particle would move in half a t ime step. C a l l this point X P , Y P . Chapter 3. Numerical Calculations 24 c) Calculate velocities at X P , Y P . d) Using the velocities obtained in c), calculate where the particle would travel in a full t ime step starting at X ( k ) , Y ( k ) . e) If the particles new position is below the boundary, place it on the nearest point on the boundary. f) Interpolate R T and Z T at the new position. g) R H O and Z E T A are updated by • i2#0(Jfc)(new) = RHO{ky°lV + RT* DT, • ZETA(k)(new) = ZETA^ + ZT * DT. To understand how the Euler ian quantities R(I ,J) and Z(I,J) are updated, consider the k t h marker particle as shown in Figure 3.6. We assign the area ( D X - X I ) * ( D Y - E T A ) to A ( I , J ) . S imi lar ly the area ( D X - X I ) * ( E T A ) is assigned to A(I,J+1) and so on. The contribution of a l l the particles that affect Z(I,J) can be calculated by ZA(I,J)*ZETA(k) Z ( 7 ' J ) - z W J ) • 3.5 Free Surface Considerations In the previous chapter, we had said at y = H. Since our reservoir is to be considered finite in size, the free surface level should drop as l iquid is withdrawn through the sink. In the program, as the surface drops, there w i l l be cells that have no marker particles. A cell that has no marker particles is to be considered to have no fluid. The program keeps the reservoir full by adding marker Chapter 3. Numerical Calculations 25 particles to the surface and therefore being able to apply the said boundary condition at y = H. This is a fair assumption to make if the product of Q and the max imum time the program simulates is small compared to the cross-sectional area of the reservoir. 3.6 Computational Stability and Accuracy W h e n descretizing our governing equations, we have assumed that terms of O(Ax) and 0 ( A y ) i n the spacial derivatives were small and thus could be neglected. To check on this grid spacing dependance, eleven computational runs were carried out for different grid spacing. These eleven runs were carried out wi th a t ime step At = 0.025 seconds, a Froude Number of 0.143E-3 and a bottom slope wi th A = 4. For each run, we have plotted tb versus t at the point x = 0.2L and y = 0.2H because we expect there to be a large gradient in rb there ( thus the largest truncation error ). Note that in Figure 3.8, each curve deviates from each other rather quickly. As we expect wi th a finer grid, as shown i n Figure 3.9, the four curves exhibit the same behaviour for t < 15 seconds. However for larger t ime there is a rather large discrepancy between the four curves. Refining the mesh even further, as shown in Figure 3.10, reduces the discrepancy between the curves even further, to a tolerable level. For locations further away from the sink, these curves tend to agree a l i t t le more, see Figure 3.11. Further reduction of the mesh is s t i l l possible, however the trade off between C . P . U . time and improved performance doesn't warrant further reduction. In Figures 3.10 and 3.11 we have indicated at t = 120 seconds the value of ib calculated using the purely buoyant solution ( P B S ) .In Figure 3.10 the agreement between the buoyant solution and the computational solution is poor. This is expected at this location, near the sink, because the bouyant solution predicts an infinitesimally th in withdrawal layer. However for locations away from the sink, we would expect much better agreement between the buoyant solution and the computational solution. This is Chapter 3. Numerical Calculations 26 the case as shown in Figure 3.11 3.7 Time Step Dependance A similar problem to what happened in our spatial derivatives if Ax and Ay were not small enough will happen in our temporal derivates if At is not taken small enough. The two main restrictions on the time step At are the Courant condition, A Umax and the diffusion restriction 2v ( see Appendix 1 ). Where the constants A and B are defined as follows 2 , „ 1 A - TTZ—r~; TZ—r—rr and B = ((Ax)" 1 + (Ay ) - i ) ((Ax)" 2 + (Ay)" 2 ) ' For NX = 60, NY = 120, A = 4, the height of the tank = 0.6096 m a n d i / = 0.1E -3 m 2 s _ 1 we get Ati = 0.211 seconds and At2 = 0.127 seconds. Both conditions may be considered together by requiring that ,A* i At2, Therefore At « 0.08 seconds. As shown in figure 3.11, at a point near the sink, the time step of At — 0.05 seconds is suitable. Chapter 4 Computational Results In the previous chapter, we have shown that for NX = 60, NY = 120, At = 0.025 seconds and a slope of 0.25, our computer model is stable. To check on the programs accuracy, we check to see if the computed solutions exhibit proper physical properties. A s mentioned in Chapter 1, we know the withdrawal of fluid should produce shear waves. Therefore we should check to see i f the computational solution yields these shear waves. The propagation of shear waves along the length of the tank is shown clearly in Figures 4.3 and 4.4. Plots of xp versus t were made at various points in a horizontal line at two different heights, from the run wi th NX = 60, NY = 120, N - 0.764 and a slope of 0.25. The location of the points chosen is shown in Figures 4.1 and 4.2. The passage of the first shear front past each of these points is indicated by a point of inflection on the plots of xp versus t given in Figures 4.3 and 4.4. The passage of the second mode is not so easily identified. Two possible reasons for this are that the second mode interacts wi th the the widening front of the first mode, and wi th reflections from the sloping boundary. These effects w i l l be discussed below in more detail. Subsequent modes are even harder to see, but arguably are present. Before the first wave has encountered the sloping boundary, it 's properties should be the same as a similar wave travelling i n a rectangular reservoir. M c E w a n and Baines (1974), have shown that the n t h mode should propogate wi th a speed c = * £ . n7r 27 Chapter 4. Computational Results 28 slope = 0.25 Y =0.5H slope = 0.25 Y= =0.75H Pts . X t ( 8 ) 7T / N (s) Pts. X t ( s ) TT / N (s) 9 1.07H 6.1 1.472 12 1.46H 6.8 1.663 13 1.60H 7.6 1.856 17 2.13H 10.3 2.495 17 2.13H 10.8 2.62 22 2.79H 12.6 3.071 21 2.67H 12.6 3.072 27 3.46H 16.3 3.967 25 3.20H 15.3 3.711 32 4.13H 17.9 4.351 29 3.73H 16.6 4.03 37 4.79H 21.1 5.119 33 4.27H 17.6 4.287 42 5.46H 25.3 6.1 37 4.80H 21.1 5.11 41 5.3H 22.1 5.37 Table 4.1: T i m e taken for the first shear wave to pass selected points F rom Figure 4.3 and Figure 4.4 we can obtain how much time it took the first wave to pass the nine points. Since we know the distance to each point we can calculate the speed of the wave. This is presented in the Table 4.1. Since there are 60 possible grid points in the horizontal direction, we have labelled these points sequentially from 1 to 60, wi th 1 being at the vertical wall . We have also done this for the 120 points in the vertical direction. These values were obtained, as mentioned earlier, by t rying to interpret where the graph exhibited a point of inflection. This procedure is a l i t t le bit subjective, but nevertheless gave results comparable wi th the theoretical values. This is best illustrated by graphing the above data wi th the theoretical values as shown i n Figure 4.5. In Figures 4.3 and 4.4, we have indicated where the first and second waves should be using t = xirn/NH, where x is the position of each test point. Note the first wave appears to coincide wi th the points of inflection and the points lie in a straight line, however the second wave doesn't. This is not the fault of the program. Recal l the wave speed C = NH/nir was derived for a rectangular reservoir. We should expect reflection of each wave to effect each following wave. We also know, by M c E w a n and Baines (1974), that the ini t ia l ly Chapter 4. Computational Results 29 R u n l Run2 Run3 Run4 Run5 Height (m) 0.61 0.61 0.61 1.2 0.61 Length (m) 2.4 2.4 2.4 1.2 6.1 A 4 4 4 1 10 Q * 1 0 - 3 m 2 / s 0.65 0.65 0.35 0.65 0.65 N(s~l) 0.76 0.076 0.76 0.76 0.76 F * 10~ 4 1.4 14.0 0.77 5.7 0.234 USCL * l O " 1 0.22 0.071 0.16 0.22 0.22 R 0.18 0.39 0.098 0.29 0.010 P r 1000 1000 1000 1000 1000 Regime C C B / C C B / C Table 4.2: Summary of varied quantities in the five experiments sharp wave front w i l l expand due to frequency dispersion at a rate proportional to tz. Therefore we expect the points of inflection and frontal passage to agree better wi th the first wave then the second. We had mentioned earlier that the run was calculated on a 60x120 grid wi th N = 0.764 and a slope of 0.25. However not al l runs were calculated for the same mentioned parameters. A brief summary of a l l the pertinent parameters is shown i n Table 4.2, accurate to two significant figures. U S C L is used in non-dimensionalizing the velocity components and is defined by Q / T H K where THK = max ( L * y/FtL*(GrPr)-k). A l l the numerical experiments where conducted on a 60x120 grid wi th a fixed kinematic viscosity and molecular diffusivity of v = 0.1E — 3 m 2 s _ 1 and K = 0.1 .£7 — 6 m 2 s _ 1 respectively. Runs 1, 2 and 4 would be classified by Imberger, Thompson and Fandry as being i n regime C by Table 1.1 ( page 5 ). Runs 3 and 5 fall on the border between regime B and regime C flows. In R u n l we chose the height and length of the computational reservoir to match that of a reservoir constructed in the laboratory. Now that we have seen that our program accurately models shear waves, we should check to see i f these shear waves reflect after interacting wi th the sloping boundary. To Chapter 4. Computational Results 30 examine this we w i l l consider run4 which had a bottom slope of 1. A n argument on stabli l i ty like that shown in the last chapter wi l l show that this run was stable. The slope of 1 was chosen to be investigated first because we knew that a wave travelling in the horizontal direction leaving the sink would reflect off the sloping bottom and then travel in the vertical direction towards the free surface. Thus points ly ing in a vertical line should experience the incident first wave roughly simultaneously, however they should experience the reflected wave at different times. Points closer to the sloping bottom wi l l experience the reflected wave first, while points closer to the surface wi l l experience the reflected wave later. The location of the points chosen is in a horizontal line at x = 0.25L, starting at y = 0.283if ( the point (16,35) ) and ending at y = 0.7H ( the point (16,85) ). The points were evenly spaced five grid points apart in the vertical direction as shown in Figure 4.6. F rom Figure 4.7, we can see that these points a l l experience the ini t ia l encounter wi th the wave at roughly the same time. Notice that the top five streamlines in Figure 4.7, corresponding to the points nearest the sloping boundary, have crests at times that are progressively later. These points are labelled 1 through 5 in Figure 4.7. This is due to the fact that points closer to the bottom experience the reflected wave earlier than points directly above. The other points by this t ime are being affected by the second mode and the reflection is not as clear. Besides checking the behaviour of the shear waves, we should check to see if surging along the sloping bottom is present. To understand this, consider a small spherical fluid particle of volume 6V displaced a distance y from its e q u i l i b r i u m position. The buoyant restoring force induced by the displacement is — epoygSV while the resulting inertial force is p6V(d2y/dt2). See Figure 4.8. Disregarding other forces and balancing these two forces results i n the following equation of motion for the spherical fluid particle Chapter 4. Computational Results 31 This results in oscillations of frequency N = y/eg. The situation for the fluid particle will be different in the presence of the sloping bottom. Since this wall is impermeable, we require that the normal density gradient be zero there. This forces the lines of constant density, which for the most part are horizontal, to bend to meet the sloping bottom perpendicularly as shown in Figure 4.9. This in turn forces the fluid particle near the sloping bottom to oscillate parallel to the slope. Following a balance of force argument as before and recalling that only the vertical component of the displacement causes a buoyant restoring force yields the following equation of motion for the spherical fluid particle near the sloping bottom + egsin26r) = 0, where r] is taken to be along the sloping bottom increasing upwards and 6 is the angle between the horizontal and the sloping bottom. This results in oscillations of frequency Ns'mO. To check our program, consider run4, with N = 0.764 and 6 — ir/4 (A = 1). The calculated frequency of oscillation as shown in Figure 4.10 is in excellent agreement with iV sin6 = 0.540. Now that we have shown that the computer program is reliable, we will turn our attention to analyzing the output for the five numerical experiments. 4.1 Runl In trying to analyse the numerical output from the computer program it is necessary to determine in what graphical form to display the output so as to be able to interpret it. The first form choosen was to plot lines of constant ib. By looking at these streamlines we can easily determine the flow direction. The streamline patterns are shown in Figures 4.11-4.18. Note that the vertical scale has been stretched by a factor of 2, so as to be able to see more clearly. At time t = 0.0 seconds, we have potential flow, which is more or less Chapter 4. Computational Results 32 just radial flow towards the sink. This wi l l be altered by the shear waves as mentioned previously. A t t = 5.0 seconds, the first shear wave has moved the 0.1 streamline. A t t = 10.0 seconds, the shear wave has moved the 0.3 streamline. Note that eddying has occured just above and below the sink. A t t = 20.0 seconds the 0.5 streamline has been displaced, the 1.0 streamline has become almost horizontal ( eddying w i l l occur below this from now on ) and the second shear wave has displaced the 0.2 streamline. A t t = 30.0 seconds the first wave has altered the 0.7 streamline, it has also reflected off the sloping wall and influenced the 0.9 streamline. Comparing the streamlines at 40.0 and 60.0 seconds, we can see that the second wave has interacted wi th the reflected first wave and the streamlines are starting to follow the buoyant l imi t . B y the time t = 120.0 seconds, we can see the streamlines just above mid-height have assumed the buoyant l imi t as shown in Figure 4.18 wi th the buoyant l imit superimposed on the numerical result. This suggests we have reached a steady state much earlier than that predicted by Imberger, Thompson and Fandry for the rectangular reservoir wi th the same fluid. F r o m Table 1.1, steady state should be achieved in approximately 1403 seconds. One possible cause for the short period of time taken to achieve this would be that of shear wave absorption in the top corner. From Figures 4.3 and 4.4, we can determine LO to be approximately 0.49 and the parameter c, from Cacchione and Wunsch to be 0.83, where c 2 = LO2/(N2 — LO2) and LO is the internal wave frequency. Since c > to we would expect absorption ( see page 6 ). To investigate the withdrawal layer, it is much easier to analyse velocity profiles. Fi rs t we must precisely define what is meant by a withdrawal layer. In the case of the rectangular reservoir, we can use as an appropriate definition, the region between the two lines of zero horizontal velocity as shown in Figure 1.6b. However in the case of a triangular reservoir, we expect ( possibly everywhere above mid-height ) a residual forward flow above mid-height which is conceptually different from a withdrawal layer. Chapter 4. Computational Results 33 They are conceptually different because a withdrawal layer is a th in jet-like region of fluid located near the depth of of the sink, moving towards the sink. A residual forward flow is a region of fluid near the surface, moving in the same horizontal direction as the withdrawal layer, however the fluid in this region wi l l not be withdrawn. For our definition we wi l l decompose our withdrawal layer into two components 6+ and 6~. The first component 6+ is measured by taking the difference in the height at which the maximum horizontal velocity occurs and the height ( above mid-depth ) at which half the max imum horizontal velocity occurs as shown in Figure 4.19. We wi l l focus our attention on 6+. Horizontal velocity profiles are shown in Figures 4.20-4.27, at the station x = 0.25Z. In the first three profiles, above mid-height, we see the contribution of the first three shear waves as expected by Figure 1.6b. Note the small region of forward flow in the profile at t — 35.0 seconds. The profiles from then on seem fairly similar in the withdrawal region and the forward flow region oscillates as seen in the profiles at t = 45.0, 50.0 and 55.0 seconds but rapidly decays to reach a steady-state as seen in Figure 4.27. Turning our attention to a station further from the sink, at x = 0.5L, we can see the effect of the first and second shear waves as shown in Figures 4.28 and 4.29. Note the peak velocity is much less at this station compared wi th the previous station as expected. However the contribution of the second shear wave is accentuated due to the reflection of the first shear wave as seen at t = 40.0 seconds. A profile of interest is that at t = 50.0 seconds, which makes it difficult to differentiate the withdrawal region and the forward flow region. This profile is caused by the reflection of the second shear wave. The third shear wave's influence is shown in the profile at t = 60.0 seconds. The profile shown at t = 75.0 seconds is slowly modified unt i l it reaches a profile as shown at t = 120.0 seconds. Chapter 4. Computational Results 34 From looking at velocity profiles in Figures 4.20-4.34 we can see that the peak hori-zontal velocity and the withdrawal layer thickness 8+ seem to approach a steady state. For the peak velocity this is clearly shown in Figures 4.35-4.37. This is also clearly shown for 8+ in Figures 4.38 and 4.39. However, this is not so clear in Figure 4.40 because a small change i n the peak velocity can cause a rather significant change in 8+, when ve-locity profiles such as Figure 4.31 make it very hard to differentiate between withdrawal region and forward flow region. The withdrawal layer thickness 8 from Table 1.1 is la-belled I T F ( for Imberger, Thompson and Fandry ) in Figure 4.39. F r o m Figure 4.39 it appears as if the withdrawal layer calculated in Table 1.1 ( page 6 ) is larger than the withdrawal layer we have numerically calculated. This should not be the case when one considers what the differences in the geometry of the triangular reservoir as opposed to the rectangular reservoir means in terms of the withdrawal layer. One would expect the withdrawal layer to be larger in the case of the triangular reservoir. However when one considers the fact that we have only plotted one component ( 8+ ) of the withdrawal layer and only considered it based on half the peak velocity, it is readily apparent that the withdrawal layer calculated in Table 1.1 is much smaller than that of our calculations. Next we have averaged the peak horizontal velocity over the 120.0 seconds of the run. This is plotted i n Figure 4.41. B y simple curve fitting, we can see that Um&x oc x~n, where n « 7/8, as opposed to n = 1/3 for the channel as presented by K o h (1966). In Figure 4.42, we have plotted the original position of the fluid that has left the reservoir by use of the marker particles mentioned in Chapter 3. We see the fluid is not removed symmetrically about mid-height, but roughly 2/3 of the fluid comes from above mid-height. Note that the withdrawal layer does not extend al l the way to the rear wall . Chapter 4. Computational Results 35 4.2 Run2 In Run2 it is useful to estimate the number of shear waves that w i l l be created by opening the sink. F rom Figure 1.6 we expect 6 = H/2n. We also expect the maximum horizontal velocities in the flow to be Umax — Q/S. F rom the work of M c E w a n and Baines (1974) or Pao and K a o (1974), we know the propogation speed of the n t h wave to be Cn = NH/mt. Equat ing these two quantities and solving for n yields For R u n l we expect the number of waves to be approximately 8. For Run2 we expect 2 waves. Keeping this in mind we wi l l now analyse Run2. In Run2 we have kept a l l the parameters the same as in R u n l , except that we have decreased the Brunt-Vaisal la frequency by a factor of 10 to observe the effect of stratifi-cation. In this case the stratification wi l l be less than in R u n l and we therefore expect the flow to resemble that of potential flow, slightly altered by the passage of a couple of shear waves. Presented i n Figures 4.43-4.45, to see the streamlines at times 40, 80 and 120 seconds. Even after 40 seconds, the potential flow pattern has not changed much due to the fact that this shear wave wi l l travel 10 times slower than in R u n l . After 80 seconds the flow pattern almost resembles that of R u n l after 5 seconds except with a larger withdrawal region. Even after 120 seconds only one shear wave can be seen, even though we expect the presence of a second shear wave by now. In Figure 4.46 we have plotted ib versus t for the same nine points as shown in Figure 4.1. F rom this it is clear that the second shear wave has not propagated. To try and draw some comparisons between R u n l and Run2, we have plotted the average peak horizontal velocity component over the 120 seconds in both runs. The case when no stratification is present, (N — 0), is also shown in Figure 4.47. For N = 0, (4.1) Chapter 4. Computational Results 36 there are no shear waves present to modify the original flow. W i t h a slight stratification we observe at least one shear wave modifying this flow, so it no longer has a uniform horizontal velocity. W i t h increasing stratification, we wi l l have more shear waves present as seen by equation (4.1). This w i l l cause the horizontal velocity profile to become further peaked for stations close to the sink. 4.3 Run3 In Run3 we have kept a l l the parameters the same as in R u n l , except that we have decreased the outflow rate Q by a factor of approximately 0.54. W h e n examining the streamlines it was noted that there is no apparent difference between R u n l and Run3 as can be seen in Figures 4.48-4.50. When analysing the horizontal velocity profiles it is important to remember that Um&x oc Q and that the speed of each shear wave is independent of Q. Examin ing the stations x = 0.25L and x — 0.5L, we notice that they have similar horizontal velocity profiles as i n R u n l . For comparison we present profiles at t — 15.0, 35.0 and 120.0 seconds at x = 0.25L and t — 35.0, 40.0 and 120.0 seconds at x = 0.5L as seen in Figures 4.51-4.56. A more appropriate way to compare R u n l and Run3 is by examining the average peak velocity and average S+. As seen by Figure 4.57, C / m a x oc Q as expected. Using (4.1) and 6 = H/2n, we expect 6 oc y/Q or that S+ in R u n l should be ~ 1/3 larger than 6+ in Run3. This , however is not the case as shown in Figure 4.58, where 6+ in R u n l is about 1/9 larger than 6+ in Run3. However for Run3 we expect about 11 shear waves to occur as compared to 8 in R u n l . Therefore we would expect 6+ for Run3 to be reduced sufficiently to agree wi th the 1/3 factor i f steady-state were achieved. Chapter 4. Computational Results 37 4.4 Run4 In Run4 we have kept a l l the parameters the same as R u n l , except that we have changed the bot tom slope so that A = 1. To do this we have halved the length of the reservoir and doubled the height so as to have the same volume of fluid as in R u n l . U p o n analysing the the streamlines, we see that we start wi th potential flow, as shown in Figure 4.59. After 5 seconds, the first shear wave has travelled the length of the tank, as seen in Figure 4.60. This can also be seen by looking at the horizontal velocity profiles at stations x = 0.25L and x = 0.5L, as shown in Figures 4.61 and 4.62. From the streamlines at 10 seconds in Figure 4.63 we can see that the second wave has displaced the 0.4 streamline. This is apparent from the velocity profile at x = 0.25L. From the velocity profile at x = 0.5Z-we can see that the second wave has not influenced this station. F rom streamlines at 15 seconds, we can see that the second wave has altered the 0.9 streamline. In Figures 4.67-4.68 we can see the effect of the second and third wave at stations x = 0.5L and x — 0.26L respectively. Note the forward flow region in Figure 4.68. After this it becomes difficult to differentiate between the effect of a new wave and that of a previously reflected wave as shown in Figures 4.69-4.71 at t = 20.0 seconds. A t x — 0.25X at most one more wave can be seen, however this maybe a reflected wave. If we continue to examine further profiles at this station, we w i l l see a rather constant withdraw region and a forward flow region that oscillates ( as in R u n l but in a more complex manner ) and decays to reach a steady-state as shown in Figures 4.72 and 4.73. However at the station x = 0.5L things are not as tame. A t 20 seconds we can see the thi rd wave has influenced the profile and the forward flow region that we have seen at 15 seconds has begun to resemble a withdrawal layer. After 25 seconds the forward flow region is faster than that of the withdrawal region and this must be attributed to wave reflection off the sloping bottom. A t 30 seconds, another forward flow region has formed. Chapter 4. Computational Results 38 After this, the surface region wi l l oscillate and decay whereas the other two regions wi l l repeatedly jo in and separate unt i l this procedure decays to yield a profile as seen at 120 seconds. See Figures 4.74-4.82. 4.5 Run5 In Run5 we have kept a l l the parameters the same as in R u n l , except that we have changed the bottom slope to A = 10 ( keeping the height of the reservoir the same as i n R u n l ). Streamlines were not plotted because it was felt that a sufficient analysis could be made on the basis of the velocity profiles. Since the only parameter that has been changed is the bot tom slope ( v ia increasing L by a factor of 2.5 ) when compared to R u n l , we should expect the same number of shear waves as in R u n l and for them to travel at the same speed. We should expect to see similar profiles at the station x = 0.25X, at a later time by a factor of 2-3 times that of R u n l , any other differences being attributed to reflection or absorption of shear waves by the sloping bottom. After 40.0 seconds we can see the effect of the first shear wave. We can see the second wave is starting to affect this station after 50 seconds and more clearly after 65 seconds as shown in Figures 4.83-4.85. After this the withdrawal region remains unchanged. The forward flow region oscillates and is shown again after 120.0 seconds in Figure 4.86. A t the station x = 0.5L the first shear wave is just starting to be felt after 25.0 seconds and after 65 seconds the expected profile is obtained, as shown in Figures 4.87 and 4.88. After 85.0 seconds, the second shear wave has influenced this station as depicted in Figure 4.89. After that, we have oscillations in the forward flow region as has been seen previously. See Figures 4.90-4.93. Looking at any station over the whole 120.0 seconds of the run, that is between the sink and x = 0AL, it seems that the withdrawal layer thickness 6+ at these stations is Chapter 4. Computational Results 39 approaching a steady-state as seen in Figure 4.94 and 4.95. For the three stations further from the sink, shown in Figure 4.96, this doesn't seem to be the case. Note the increase in 8+ ( the three succesive peaks ) at the three stations between 50.0 and 75.0 seconds. This is not an effect of the second shear wave because as shown in the previous velocity profiles, the second wave is not felt unt i l after this time. In Figures 4.88 and 4.89 we clearly see that only one wave has passed. Notice that the height at which Umax occurs in the second figure is much closer to the free surface. This is caused by the first wave reflecting off the sloping bottom. To make some comparison between R u n l , Run4 and Run5, we have averaged S+ at each station and plotted the results in Figure 4.97. As mentioned earlier, it seems we approach a steady-state 6+ at the stations between the sink and x = OAL, where the forward flow is less prominant. As seen in Figure 4.97, by increasing the bot tom slope A, we are increasing S+ in the region between the sink and x — OAL. After this point the effect of forward flow, it seems, is to increase 8+ wi th decreasing A. Therefore if one is interested in measuring the flow of fluid actually in the withdrawal layer and not just the fluid moving forward, one should l imi t ones measurements between x = 0.0L and x = OAL. This could explain why when comparisons have been made between theoretical and experimental models ( a l l in rectangular domains ) the withdrawal layer thickness has been smaller than that measured i n actual reservoirs ( see K a o , Pao and Wei (1974) and Imberger (1980) ). In Figure 4.98, we have plotted 8+ at 120 seconds for six stations between x — 0.0L and x = OAL. Chapter 5 Conclusions and Recomendation for Further Research 5.1 Conclusions From analysing the five numerical experiments, we draw the following conclusions: 1. Shear waves were observed and at least the first mode was seen to travel at the speed N H / n x as suggested by Pao and K a o (1974) and also by M c E w a n and Baines (1974). 2. F rom analysis of the velocity plots, a withdrawal layer ( a th in jet-like region of fluid located near the depth of the sink that is moving towards the sink ) and a forward flow region ( a region of fluid near the surface that is moving wi th positive horizontal velocity ) were observed. In some instances the forward flow region travelled wi th a greater positive velocity than that of the withdrawal layer ( i.e. Figure 4.74 ). 3. There is evidence that shear waves are both reflected off the sloping bottom and absorped into the acute corner. The oscillating forward flow region indicates reflec-tion, while a steady-state reached long before that of a rectangular reservoir wi th the same properties indicates absorption. 4. A stable withdrawal layer, based on our definition of 8+ ( Figure 4.18 ) was observed to form between the sink and x = OAL. 40 Chapter 5. Conclusions and Recomendation for Further Research 41 5. The sloping bottom establishes a forward flow region almost everywhere above mid-height. In the case of a rectangular reservoir, the vertical endwall reflects waves horizontally and therefore l i t t le to no forward flow region exists because the main contribution of each wave occurs at the level of the sink tending to reduce the withdrawal layer thickness. 6. Because the fluid in the withdrawal layer w i l l be removed through the sink and very l i t t le of the fluid in the forward flow region w i l l , great care must be taken when making field measurements to ascertain what region one is actually measuring. 7. A quasi steady state streamline pattern, slightly above the withdrawal layer, closely resembled that of the purely buoyant solution by the completion of a l l runs (120 seconds ). 5.2 Further Research As noted i n i tem 7 in Chapter 5, Section 1, the computed streamline pattern eventually resembled the purely buoyant solution. In view of this, a new definition of 6+ should be used i n future studies. Define 8+ by taking the difference i n the height at which the numerically calculated maximum horizontal velocity occurs and the height ( above mid-depth ) at which the horizontal velocity calculated from the purely buoyant solution and the numerically calculated solution intersect. M a n y changes can be made to the computer model to make it more efficient and model the physical problem more effectively. The use of a F in i te Element method could be used so that a no-slip condition could be imposed at a l l computational points along the sloping boundary and the vertical wall . This would also permit the use of a non-uniform grid near the sink for better detail in investigating the withdrawal structure. A M u l t i g r i d technique should be employed so as to reduce the number of iterations Chapter 5. Conclusions and Recomendation for Further Research 42 necessary for convergence. As the model stands at present it can only be run for a small t ime period so that the free-surface has not fallen much. To eliminate this problem a moving coordinate frame which contracts wi th the falling surface should be used. See M c G u i r k and Islam (1986). Other physical problems should be looked into. One of them would be the effect of variable stratification, so as to model the epil imnion ( the region of water near the surface that is warm due to circulation by the wind ) and the hypolimnion ( the region of cold water which is unaffected by the sun ). This could be done by changing the in i t ia l density stratification. The effect of a sink located at a position other than mid-depth should be investigated, as well as the effect of multiple sink interactions. As far as analytic work goes, there is not much hope of solving the governing equations. However it would be interesting to see the results if someone could solve (2.32) without assuming steady-state, even if it is done asymptotically for small t ime as was done in Imberger, Thompson and Fandry. Bibliography [1] Cacchione, D . and Wunsch, C . (1974) Experimental study of internal waves over a slope. J . F l u i d Mech . , 66, 223 [2] Debler, W . R . (1959) Stratified flow into a line sink. Journal of the Engineering Mechanics Divis ion , A . S . C . E . , 85, 51 [3] Imberger, J . (1980) Selective Withdrawal: A Review. In "Proceedings of the Sec-ond International Symposium on Stratified Flows", International Association for Hydraul ic Research ( I A H R ) , Trodheim, Norway, 381 [4] Imberger, J . and Fandry, C . (1975) Withdrawal of a Stratified Fluid from a Vertical Two Dimensional Duct. J . F l u i d Mech. , 70, 321 [5] Imberger, J . , Thompson, R . and Fandry, C . (1976) Selective withdrawal from a finite rectangular tank. J . F l u i d Mech . , 78, 489 [6] K a o , T . W . (1965) A free-streamline solution for stratified flow into a line sink. J . F l u i d Mech . , 21, 535 [7] Kao , T . W . , Pao, H . P . and Wei , S .N . Dynamics of establishment of selective with-drawal of a stratified fluid from a line sink. Part2. Experiment. J . F l u i d Mech. 65, 689 [8] K o h , R . C . Y . (1966) Viscous stratified flow towards a sink. J . F l u i d Mech . , 24, 555 [9] Lawrence, G . L . (1976) Selective Withdrawal of a density stratified fluid from a trian-gular tank. Honours Thesis, Department of C i v i l Engineering, Universi ty of Western 43 Bibliography 44 Aust ra l ia [10] M c E w a n , A . D . and Baines, P . G . (1974) Shear fronts and an experimental stratified shear flow. J . F l u i d Mech. , 63, 257 [11] M c G u i r k , and Islam (1986) Selective withdrawal from a stratified flow into a line sink. In " Proceedings of the International Symposium on Buoyant Flows", Athens, Greece, 498. [12] Pao, H . P . and Kao , T . W . (1974) Dynamics of establishment of selective withdrawal of a stratified fluid from a line sink. Parti. Theory. J . F l u i d Mech . , 65, 657 [13] Silvester, R . S . (1978) An experimental study of end-wall effects on selective with-drawal of linearly stratified liquid from a reservoir. P h . D . Thesis, Department of C i v i l Engineering, Universi ty of western Aust ra l ia [14] Y i h , C .S . (1964) Dynamics of Non-Homogeneous Fluids . M c M i l l a n Company Appendix A Stability Analysis Consider a very simple analysis of a one-dimensional model equation to study numerical instability. The equation is as follows dt dx dx2' ( A . l ) For simplicity, we w i l l now only consider the diffusion term and we wi l l assume that a steady-state solution £™ = 0 has been reached for a l l points i. Introducing a small perturbation e in which may occur due to numerical roundoff etc., we obtain the following discrete equation ff+i _ (ff + e) A t { £ 5 7 [ & 1 + ff-i - m + e)] B y previous assumptions this yields •n+1 2ve At (AXy or C + 1 = c(l - 2d), where the diffusion number is defined as d = vAt ( A x ) 2 " For stability, this disturbance must die out. This requires q i < 1 (A.2) (A.3) (A.4) (A.5) 45 Appendix A. Stability Analysis 46 or —1 < 1 — 2J < 1. (A.6) The right-hand inequality is automatically satisfied by the fact that v and d are both positive. The left-hand inequality is satisfied when d < 1. Requiring that the computer program models the physics in that "no overshoot" be allowed, we have ff+1 and by A.4 we have that > 0 (A.7) d<\. (A.8) This places a restriction on our time step that In two dimensional problems, this restriction becomes At<^~ where B = ( ( A x ) ~ 2 + (Ay)"2)"1 . (A.10) We now consider both advection and diffusion and without loss of generality assume u > 0. Descretizing A. l yields ff+1 ~ (ff + 0 «ffk ~ "ff.! . + ff.! - 2(g + 6 ) A l - 2A~x + V (Axy • ( A - U j Note that the advection term does not include any perturbation term at i. However, if we apply the same discretization at the point i + 1, we get = + e ) ~ + TX^ fe1 + & + e ) " + ^ - J ( A - 1 2 ) 2Ax vAt (Ax) Appendix A. Stability Analysis 47 or after assuming steady-state as before & = ~je + de, (A.13) where C = uAt/Ax is called the "Courant Number". Requir ing zero overshoot at i — 1, we have that tn+l > 0. (A.14) Using A.13 yields -j + d>0 (A.15) or d v Using A . 8 yields the "Courant condition" C uAx n < 2. (A.16) Ax or alternatively C = ^ < 1 (A.17) A x At<—. (A.18) u In two dimensions this restriction becomes A t < A/u where A = 2 ( ( A x ) - 1 + ( A y ) - 1 ) " 1 . (A.19) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 C 18 c 19 c 20 c 21 c 22 c 23 c 24 c 25 26 27 c 28 c 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 c 50 c 51 c 52 53 54 55 56 57 58 Appendix B The Program PROGRAM SRF IMPLICIT REAL*8 (A-H,0-Z) PARAMETER (INX=61 , INZ=121, INP«20000) DIMENSION PSI(INX,INZ),U(INX,INZ),V(INX,IN2) DIMENSION Z(INX,INZ),ZT(INX,INZ),R(INX,INZ),RT(INX,INZ) DIMENSION ROUT(INX) DIMENSION ZETA(INP),RHO(INP),XORG(INP),YORG(INP) DIMENSION X(INP),Y(INP),INK(INP),21(INP) DIMENSION IOUT(50) DIMENSION PSIDEL(INX,INZ) COMMON NX,NZ,DX,DZ COMMON DXSQ,DZSQ,ORF,Q,SML,D2 COMMON TIME,WO,STRATF,HTO COMMON U,V,ZT SRF- STRATIFIED RESERVOIR FLOWS: This program has been written in order to determine the effects of geometry on the dynamics of a s t r a t i f i e d reservoir. I/O UNITS: 1 -- INPUT UNIT "*«*.DAT" 2 -- OUPPUT DATA (i n c l contor data) in "***.OUT" 3 -- TIME HISTORY OUTPUT DATA UNIT "***.SEC-OPEN (2,FILE*'SRF.OUT*,STATUS='UNKNOWN') OPEN(8,FILE='PPLOT',STATUS-'UNKNOWN') CALL INPUT(NX,NZ,HTO,BVN,DT,TIMELM,VISCOS,DIFFUS,Q,IP, $ IOUT,ORF,SLOPE,NCASE) WRITE12,*) NX,NZ,DT QA •= DABS (Q) H ' HTO/2. XWIDTH •= HTO/SLOPE DX « HTO/(SLOPE*(FLOAT(NX)-1.)) DZ = HTO/(FLOAT(NZ)-1.) WRITEI2,*) DX,DZ WRITE(2,*) (IOUT(I),1-1,IP) DXSQ » DX**2 DZSQ - DZ**2 D2 « 2.*(1./DXSQ • 1.0/DZSQ) 0O=NZ/2 JOP1 • JO + 1 JOM1-JO-1 GVRO * 9.81/1000. STRATF • BVN**2/GVRO WO « Q/XWIDTH CALCULATE FROUDE AND RAYLEIGH NUMBERS RA ' BVN**2*XWIDTH**4/(DIFFUS*VISCOS) F « QA/(BVN*XWIDTH**2) RATIO - F*(RA**(1./3.)) RAMI6 - RA**(-l./6.) WRITE(2,*) ' F, RA, GR', F, RA, RA*DIFFUS/VISCOS 48 Appendix B. The Program 49 175 ZETA(NP) » 0. 176 RHO(NP) « STRATF*(H-Y(NP) ) 177 INK(NP) « 0 178 I F ( I .EQ. 16) INK(NP) ' 4 179 I F ( J .EQ. K) INK(NP) « 5 180 1160 CONTINUE 181 182 C 183 c SET THE RESIDUAL WHICH WILL BE TOLERATED AFTER RELAXATION 184 c OF THE POISSON EQUATION 185 c 186 187 SML = OA*0.00002*D2 188 DTMAX •= 3.1415926/(SLOPE*BVN*(NX-1)) 189 WRITE(2,280) DTMAX 190 '' 280 FORMAT('DTMAX ',F9.3) 191 192 C 193 C BEGIN MAIN LOOP OF PROGRAM FROM HERE, TIME STEPS START HERE 194 . C 195 196 C 197 c I N T I A L I Z E PSIDEL TO 0.0 198 c 199 200 DO 1238 J=NZ,1,-1 201 DO 1238 1=1,NX 202 1238 P S I D E L ( I , J ) = 0.0 203 204 TIME = 0.0 205 JTSTEP = 0 206 IPP=1 207 WRITE(3,«) I P 208 1180 JTSTEP = JTSTEP+1 209 TIME •= (JTSTEP-1) *DT 210 C WRITE(*,290) NP 211 C290 FORMAT (' NUMBER OF PARTICLES = ',14) 212 213 C 214 C RELAX FOR STREAM FUNCTION PSI 215 c 216 217 CALL RELAX(PSI,Z,SLOPE,PSIDEL) 218 219 c 220 c CALCULATE VELOCITY USING THE STREAM FUNCTION PSI 221 c 222 223 DO 1200 1=1,NX 224 K = 2*1-4 225 I F (K .LT. 1) K - l 226 DO 1200 J=K,NZ-1 227 1200 U ( I , J ) • ( P S I ( I , J + 1 ) - P S I ( I , J ) ) / D Z 228 DO 1220 I«1,NX-1 229 K » 2*1-3 230 I F (K .LT. 1) K • 1 231 DO 1220 J«K,NZ-1 232 1220 V ( I , J ) « ( P S I ( I , J ) - P S I ( 1 + 1 , J ) ) / D X Appendix B. The Program 50 117 P S l ( N X , K 2 - 3 ) " ( 2 . 0 - ( F L O A T ( N X ) - 2 . 5 ) / ( F L O A T ( N X ) - 1 118 PSI(NX,NZ-4) • (2.0-(FLOAT(NX)-3.0) / (FLOAT (NX)-1 119 ELSE 120 PSI(NX,NZ-1) • Q 121 PSI(NX,NZ-2) • Q 122 PSI(NX,NZ-3) « Q 123 PSI(NX,NZ-4) ' Q 124 PSI(NX-1,NZ-5) • (FLOAT(NX)/FLOAT(NX-1))*Q 125 PSI(NX-1,NZ-6) * (FLOAT(NX)/FLOAT(NX-1))*Q 126 END I F 127 128 C 129 C ASSIGN I N I T I A L VALUES FOR STREAMFUNCTION 130 C 131 132 DO 1080,I«2,NX-1 133 K « 2*1 134 XX = ( I - 1 ) * D X 135 DO 1080,J«K,NZ 136 y y = ( j - i ) * D Z 137 1080 PSI(I,J)«Q*XX*SLOPE/YY 138 139 C 140 c WALL BOUNDARY CONDITIONS 141 c 142 143 DO 1100,J<=2,NZ-1 144 U ( 1 , J ) = 0.0 145 R T ( 1 , J ) = 0.0 146 1100 Z T ( 1 , J ) = 0.0 147 148 I F (NCASE .NE. 1) THEN 149 DO 1110,J=2,JO 150 1110 P S I ( 1 , J ) « 0 151 DO 1120,J=JOP1 ,NZ 152 1120 P S I ( 1 , J ) - 0.0 153 ELSE 154 DO 1140,J=2,NZ 155 1140 P S I ( 1 , J ) • 0.0 156 END I F 157 C 158 C ASSIGN I N I T I A L PARTICLE PARAMETERS 159 C 160 161 NP = 0 162 MEBL ' 2 163 M « MEBL*(NX-1)+1 164 MJ « MEBL*(NZ-1)*1 165 DO 1160, I-1,M 166 K « 2*1-1 167 DO 1160,J«K,MJ 168 NP • NP+1 169 I F (NP .GT.25000) STOP 'STOP FOR NP > 25000' 170 X(NP) - (1-1.)*DX/FLOAT(MEBL) 171 Y(NP) » (J - l . ) * D Z / F L O A T ( M E B L ) 172 XORG(NP)«X(NP) 173 YORG(NP)«Y(NP) 174 21(NP) • FLOAT(I) Appendix B. The Program 51 59 THKDF • XWIDTH*RAM16 60 THKIN « XWIDTH*DSQRT(F) 61 THK * THKDF 62 I F (THKIN .GT. THK) THK • THKIN 63 TINO • THKIN*XWIDTH/QA 64 USCL « QA/THK 65 VSCL * QA/XWIDTH 66 VRTSCL «= USCL/THK 67 68 C 69 C ESTABLISH BOUNDARY CONDITIONS 70 C 71 72 c 73 C I N I T I A L I Z E ALL ARRAYS TO ZERO 74 C 75 DO 900, I « 1,NX 76 DO 900 J=1,NZ 77 P S I ( I , J ) •= 0.0 78 Z T ( I , J ) - 0.0 79 U ( I , J ) = 0.0 B0 V ( I , J ) « 0.0 81 R ( I , J ) « 0.0 62 900 R T ( I , J ) = 0.0 83 DO 920 1=1,NX 84 K=2*I-1 85 DO 920 0=K,NZ 86 920 Z(I,J>=0.0 67 88 C 89 C ASSIGN I N I T I A L DENSITY STRATIFICATION 90 C 91 92 DO 1020 I«1,NX-1 93 DO 1020 J=1,NZ-1 94 1020 R ( I , J ) * STRATF*(-(FLOAT (J)-0.5)*DZ+H) 95 96 C 97 C SURFACE BOUNDARY CONDITIONS 98 C 99 100 DO 1040 I=1,NX-1 101 1040 V(I , N Z ) = -WO 102 103 C 104 C SLOPE BOUNDARY CONDITIONS 105 C 106 107 DO 1050 I * 1,NX 108 ZT(I,NZ) « 0. 109 K-=2*I-1 110 Z T ( I , K ) » 0. 111 P S I ( I , K ) « 0 112 1050 CONTINUE 113 114 I F (SLOPE .EQ. 1.) THEN 115 PSI(NX,NZ-1) • ( 2.0-(FLOAT(NX)-1.5)/(FLOAT(NX) -1.) ) *Q 116 PSI(NX,NZ - 2 ) - ( 2 . 0 -(FLOAT(NX) - 2 . 0 )/(FLOAT(NX)-1.))*Q Appendix B. The Program 52 233 234 C 235 C CALL OUTPUT SUBROUTINE 236 C 237 C WRITE(4,1) JTSTEP-1 238 C 1 FORMAT('JTSTEP « ' , 13) 239 IF(IPP.LE.IP.AND.JTSTEP-1.EQ.IOUT(IPP) ) THEN 240 WRITE(2,«) TIME 241 CALL OUTPUr<ROUT,RVN,RAM16,TINO,X,Y,21,Z,VSCL, 242 s INK,ZETA,RHO,NP,PSI,R,USCL,VRTSCL,PSIDEL) 243 WRITE(3,' (14F8.4) ') T I M E , P S I ( 2 , 3 5 ) / Q , P S I ( 3 , 3 5 ) /Q, 244 s P S I ( 4 , 3 5 ) / Q , P S I ( 5,35)/Q,PSI( 6,35)/Q,PSI< 7,35)/Q, 245 s PSI ( 8, 35) /Q, PSI (9, 35) /Q, PSI (10, 35) /Q, PSI (11, 35) /Q, 246 s P S I ( 1 2 , 3 5 ) / Q , P S I ( 1 3 , 3 5 ) / Q , P S I ( 1 4 , 35) /Q 247 C DO 1385 J=1,NZ 248 C1385 WRITE ( 2 , ' (41F9.4) ') ( P S I ( I , J ) / Q , I = 1,NX) 249- IPP-IPP+1 250 END I F 251 I F (TIME .GE. TIMELM) STOP 'TIME OUT...' 252 DO 1380,1*2,NX-2 253 K ' 2*1*1 254 DO 1380 J«K,N2-1 255 RP = ( R ( I , J - l ) + R ( I , J ) ) 1 2 . 256 RM ' ( R ( I - 1 , J - 1 ) • R ( I - l , J ) ) / 2 . 257 RX «= (RP-RM) /DX 258 ZXX = ( Z { 1 + 1 , J ) - 2 . * Z ( I , J ) + Z ( 1 - 1 , J ) ) / D X S Q 259 ZZZ = ( 2 ( 1 , J + l ) - 2 . * Z ( I , J ) + Z ( I , J - l ) ) / D Z S Q 260 1380 Z T ( I , J ) = GVRO*RX+VISCOS*(ZXX+222) 261 262 DO 1400,1=2,NX-2 263 K=2«I+1 264 DO 1400,J=K,NZ-2 265 RXX = ( R ( I + 1 , J ) - 2 . * R ( I , J ) + R ( I - 1 , J ) ) / D X S Q 266 RZZ = ( R ( I , J + 1 ) - 2 . * R ( I , J ) + R ( I , J - 1 ) ) / D Z S Q 267 1400 R T ( I , J ) •= DIFFUS* (RXX + RZ2) 268 269 C 270 c UPDATE PARTICLE POSITIONS 271 c 272 273 K = 1 274 1410 XP = X(K) + 0.5 * DT * U I N T ( X ( K ) , Y ( K ) ) 275 YP « Y(K) + 0.5 * DT * V1NT(X(K),Y(K)> 276 I F (XP .LT. 0. .OR. YP .LT. 0.) THEN 277 X(K)=XP 278 Y(K)=YP 279 ELSE 280 XPP « X(K) • DT * UINT(XP,YP) 281 YPP « Y(K) • DT * VINT(XP,YP) 282 283 C 284 C ENSURE THAT PARTICLES ON SLOPING BOUNDARY TRAVEL 285 C ALONG SLOPING BOUNDARY ONLY 286 C 287 288 I F ( X ( K ) " S L O P E . L T . Y ( K ) ) GOTO 1420 289 XPP • (XPP+YPP/SLOPE)/2. 290 YPP " XPP'SLOPE Appendix B. The Program 291 INK(K) • 5 292 293 C 294 C ASSIGN POSITIONS ON THE BOUNDARY TO PARTICLES WHICH TRY 295 C CROSS THE BOUNDARY 296 C 297 298 1420 I F (YPP .GE. XPP*SLOPE) GOTO 1430 299 XPP = <XPP*YPP/SLOPE)12. 300 YPP « XPP*SLOPE 301 I N K ( K ) - 6 302 1430 XB = (X(K)*XPP)12. 303 X(K) ' XPP 304 YB ' (Y(K)+YPP)12. 305 Y(K) = YPP 306 C 307 c INTERPOPLATE TIME CHANGE 308 c 309 310 I F (YB .GT. HTO) GOTO 1460 311 RHO(K) • RHO(K) * RTINT(XB,YB,RT)*DT 312 ZETA(K) « ZETA(K)+ZTINT(XB,YB)*DT 313 IF (INK(K) .LT. 0) GOTO 1460 314 END I F 315 316 I F (X(K) .GE. 0. .AND. X(K)"SLOPE .LE. Y ( K ) ) GOTO 1460 317 INK(K) * -9999 318 1460 K = K+l 319 I F (K .LE. NP) GOTO 1410 320 321 C 322 C ERASE POINTS OUTSIDE, EXCEPT WHERE OK 323 C 324 325 K = 1 326 1470 I F (X(K) .LT. 0. .OR. X(K) .GT. XWIDTH 327 .OR. Y(K) .LT. 0.) THEN 328 C WRITE (12,*) TIME, XORG(K ) , YORG(K ) 329 NP • NP-1 330 I F (NP .LE. O) STOP 'BREAK AT NP < 0' 331 I F (K .GT. NP) GOTO 1600 332 DO 1490 KK«K,NP 333 KK1 » KK + 1 334 X(KK) • X(KK1) 335 Y(KK) * Y(KK1) 336 XORG(KK)«XORG(KKl) 337 YORG(KK)'YORG(KK1) 338 Z l (KK) " Z K K K l ) 339 ZETA(KK) « Z E T A ( K K l ) 340 RHO(KK) « RHO(KKl) 341 1490 INK(KK) - I N K ( K K l ) 342 GOTO 1470 343 END I F 344 K « K + l 345 I F (K .LE. NP) GOTO 1470 346 347 C 348 c COMPUTE NEW R AND Z Appendix B. The Program 54 349 C 350 351 1600 CALL INTPR(X,Y,Zl,NP,RHO,R,ZETA,INK,Z,XORG,YORG) 352 CALL INTPZ(X,Y,NP,ZETA,Z) 353 GOTO 1180 354 END 355 356 357 FUNCTION UINT(X,Y) 358 I M P L I C I T REAL*8 (A-H,0-Z) 359 PARAMETER (IKX*61 , INZ=121, INP'20000) 360 DIMENSION U ( I N X , I N Z ) , V ( I N X , I N Z ) , Z T ( I N X , I N Z ) 361 COMMON NX,NZ,DX,DZ 362 COMMON DXSQ,DZSQ /ORF,Q,SML,D2 363 COMMON TIME,WO,STRATF,HTO 364 COMMON U,V,ZT 365 366 C 367 C INTERPOLATE U WITH U ( I , J ) AT X = ( I - 1 ) * D X , J - ( J - O . 5 ) 368 c 369 370 I = INT( X/DX • 1.5 ) 371 J *INT( Y/DZ • 1. ) 372 c WRITE ( 2 , * ) X , Y , I , J 373 UINT ' 0. 374 I F (X .LT. 0. .OR. Y.LT.0.) RETURN 375 I F (U(1,J) .EQ. 0. .AND. X .EQ. 0.) RETURN 376 X I ' X - F L O A T ( I - l ) * D X 377 ETA = Y - ( F L O A T ( J ) - 0 . 5 ) * D Z 378 IS " 1 379 I F (XI .LT. 0.) I S • -1 380 I P - I • I S 381 J S « 1 382 I F (ETA .LT. 0.) J S = -1 383 J P » J+JS 384 I F ( I P .GT. NX) I P - NX 385 I F ( J .GE. NZ) J • NZ-1 386 I F ( J P .GE. NZ) J P « NZ-1 387 I F ( J P .LT. 1) J P • 1 388 C l " DABS(XI/DX) 389 C2 " DABS(ETA/DZ) 390 C3 * 1. - C l 391 C4 « 1. - C2 392 UINT « U ( I , J ) * C 3 * C 4 + U ( I P , J ) * C 1 * C 4 393 $ + U ( I , J P ) * C 3 * C 2 + U ( I P , J P ) * C 1 * C 2 394 RETURN 395 END 396 397 FUNCTION VINT(X,Y) 398 I M P L I C I T REAL*8 (A-H,0-Z) 399 PARAMETER (INX»61 , INZ=121, INP«20000) 400 DIMENSION U ( I N X , I N Z ) , V ( I N X , I N Z ) , Z T ( I N X , I N Z ) 401 COMMON NX,NZ,DX,DZ 402 COMMON DXSQ,DZSQ,ORF,Q,SML,D2 403 COMMON TIME,WO,STRATF,HTO 404 COMMON U,V,ZT 405 406 c Appendix B. The Program 407 C V ( I , J ) IS DEF AT X«(1-0.5)*DX, Y»(J-1)«DZ 406 C 409 410 VINT ' -WO 411 I F (X .LT. 0. .OR. Y .LT. 0.) RETURN 412 I «INT( X/DX+1.) 413 J »INT( Y/DZ+1.5) 414 I F (Y .EQ. HTO) RETURN 415 X I = X - ( F L O A T ( I ) - 0 . 5 ) * D X 416 ETA « Y - F L O A T ( J - l ) * D Z 417 I S ' 1 416 I F (XI .LT. 0.) I S * -1 419 I P » I • I S 420 J S * 1 421 I F (ETA .LT. 0.) J S • -1 422 J P = J + JS 423 424 C 425 C THE REFLECTIVE B.C.S. 426 C 427 428 I F ( I .GE. NX) I-NX-1 429 I F ( I P .GE. NX) IP=NX-1 430 I F ( I P .LT. 1) IP=1 431 I F (JP.GT.NZ) JP'NZ-1 432 CI * DABS(XI/DX) 433 C2 = DABS(ETA/DZ) 434 C3 « l . - C l 435 C4 «= 1.-C2 436 VINT ' V ( I , J ) * C 3 * C 4 * V ( I P , J ) * C 1 * C 4 437 S + V ( I , J P ) * C 3 * C 2 + V ( I P , J P ) * C 1 * C 2 438 RETURN 439 END 440 441 442 FUNCTION ZTINT(X,Y) 443 I M P L I C I T REAL*8 (A-H,0-Z) 444 PARAMETER (INX-61 , INZ=121, INP*20000) 445 DIMENSION U ( I N X , I N Z ) , V ( I N X , I N Z ) , ZT(INX,INZ) 446 COMMON NX,NZ,DX,DZ 447 COMMON DXSQ,DZSO,ORF,Q,SML,D2 448 COMMON TIME,WO,STRATF, HTO 449 COMMON U,V,ZT 450 451 C 452 C INTERP ZT. Z I ( I , J ) I S DEF AT X«(I-1)*DX, Y«(J 453 C 454 455 I «INT( X/DX+1.5 ) 456 J «INT( Y/DZ+1.5 ) 457 ZTINT • 0. 456 I F (X .LT. 0. .OR. Y .LT. 0.) RETURN 459 X I « X - F L O A T ( I - l ) * D X 460 ETA • Y - F L O A T ( J - l ) * D Z 461 I S - 1 462 I F ( X I .LT. 0.) I S - -1 463 I P • I • I S 464 I F ( I P .LT. 1) I P • 1 Appendix B. The Program 465 JS • l 466 IF (ETA .LT. 0.) JS • -1 467 JP «= J + JS 468 IF (IP .GT. NX) IP«NX 469 IF (JP .GT. NZ) JP*NZ 470 CI ' DABS(XI/DX) 471 C2 ' DABS(ETA/DZ) 472 C3 - l . - C l 473 C4 = 1.-C2 474 ZTINT = ZT(I,J)*C3*C4 + ZT(IP,J) *C1*C4 475 S +ZT(I,JP)*C3*C2+ZT(IP,JP)*C1*C2 4 76 RETURN 477 END 478 479 480 FUNCTION RTINT(X,Y,RT) 481 IMPLICIT REAL*8 (A-H,0-Z) 482 PARAMETER (INX=61 , INZ=121, INP=20000) 483 DIMENSION RT(INX,INZ) 484 COMMON NX,NZ,DX,DZ 485 486 C 487 C INTERP RT, DEF AT X=(1-0.5)*DX, Y*(J-0.5)*DZ 488 C 489 490 RTINT •= 0. 491 IF (X .LT. 0. .OR. Y .LT. 0.) RETURN 492 I = INT(X/DX+1.) 493 J =INT( Y/DZ+1.) 494 XI = X-(FLOAT(I)-0.5)*DX 495 ETA « Y-(FLOAT(J)-0.5)*DZ 496 IS = 1 497 IF (XI .LT. 0.) IS = -1 498 IP = I+IS 499 IF (IP .LT. 1) IP * 1 500 JS = 1 501 IF (ETA .LT. 0.) JS « -1 502 JP ' J+JS 503 IF (IP .LT. 1) IP=1 504 IF (JP .LT. 1) JP=1 505 IF (I .GE. NX) I«NX-1 506 IF (J .GE. NZ) J-NZ-1 507 IF (IP .GE. NX) IP«NX-1 508 IF (JP .GE. NZ) JP-NZ-1 509 CI * DABS(XI/DX) 510 C2 •= DABS (ETA/DZ) 511 C3 « l . - C l 512 C4 = 1.-C2 513 RTINT * RT(I,J)*C3*C4*RT(IP,J)*C1*C4 514 S •RT(I,JP)*C3*C2*RT(IP,JP)*C1*C2 515 RETURN 516 END 517 518 519 520 SUBROUTINE INTPZ(X,Y,NP,ZETA,Z) 521 522 IMPLICIT REAL'S <A-H,0-Z) Appendix B. The Program 523 PARAMETER (INX«61 , INZ»121, INP'20000) 524 DIMENSION SUM(INX,INZ), ZSUM(INX,INZ), Z(INX,INZ) 525 DIMENSION X(INP), Y(INP), ZETAUNP) 526 COMMON NX,NZ,DX,DZ 527 528 DO 2000 I ' 1,NX 529 K * 2*1 - 3 530 IF(K.LT.I) K • 1 531 DO 2000, J * K,NZ 532 SUM(I,J) « 0. 533 2000 ZSUM(I,J) * 0. 534 535 DO 2020, K * 1,NP 536 I «INT( X(K)/DX + 1.5 ) 537 J «INT( Y(K)/DZ +1.5 ) 538 XI « X(K) -FLOAT(I-l) *DX 539 ETA •= Y(K) - FLOAT(J-l)«DZ 540 IS •= 1 541 IF (XI .LT. 0.) IS •= -1 542 IP - I • IS 543 JS * 1 544 IF (ETA .LT. 0.) JS ' -1 545 JP «= J • JS 546 ABX ' DABS(XI/DX) 547 ABY * DABS(ETA/DZ) 54 8 A - (l.-ABX) * (l.-ABY) 54 9 SUM(I,J) ' SUM(I,J) + A 550 ZSUMU.J) * ZSUM(I,J) • A*ZETA(K) 551 IF (IP .GT. NX) GOTO 2030 552 A = ABX * (l.-ABY) 553 SUM(IP,J) « SUM(IP,J) + A 554 ZSUM(IP,J) « ZSUM(IP,J) + A*ZETA(K) 555 IF (JP .GT. NZ) GOTO 2020 556 A " ABX * ABY 557 SUM(IP,JP) - SUM(IP,JP) • A 558 ZSUM(IP,JP) • ZSUM(IP,JP) + A*ZETA(K) 559 2030 IF (JP.GT.NZ) GOTO 2020 560 A = (l.-ABX) * ABY 561 SUM(I,JP) * SUM(I,JP) * A 562 ZSUM(I,JP) « ZSUM(I,JP) * A*ZETA(K) 563 2020 CONTINUE 564 565 DO 2040, I « 1,NX 566 K = 2*1 - 1 567 DO 2040 J « K,NZ 568 IF (SUM(I,J) .EQ. 0.0) GOTO 2040 569 Z<I,J) « ZSUM(I,J)/SUM(I,J) 570 2040 CONTINUE 571 RETURN 572 END 573 574 SUBROUTINE INTPR(X,Y,Zl,NP,RHO,R,ZETA,INK,Z,XORG,YORG) 575 576 IMPLICIT REAL*8 (A-H,0-Z) 577 PARAMETER (INX-61 , INZ«121, INP-20000) 578 DIMENSION SUM(INX,INZ), R(INX,INZ) 579 DIMENSION X(INP),Y(INP),Zl(INP),RSUM(INX,INZ) ,XORG(INP) 580 DIMENSION ZETA(INP),INK(INP),Z(INX,INZ),RHO(INP),YORG(INP) Appendix B. The Program 58 5B1 582 COMMON NX,NZ,DX,DZ 583 COMMON DXSQ,DZSQ,ORF,Q,SML,D2 584 COMMON TIME,WO,STRATF,HTO 585 586 DO 2100, 1=1,NX-1 587 K • 2*1-3 588 IF (K .LT. 1) K = l 589 DO 2100,J = K,NZ-1 590 SUM(I,J) ' 0. 591 2100 RSUM(I,J) * 0. 592 593 DO 2120,K=1,NP 594 I = INT(X(K) /DX+1.) 595 J = INT (Y (K) /DZ+1.) 596 XI ' X(K)-(FLOAT(I)-0.5)*DX 597 ETA = Y (K) -(FLOAT(J)-0.5)*DZ 598 IS = 1 599 IF (XI .LT. 0.) IS * -1 600 . • IP = I + IS 601 JS = 1 602 IF (ETA .LT. 0.) JS = -1 603 JP = J+JS 604 ABX = DABS(XI/DX) 605 ABY = DABS (ETA/DZ) 606 IF (J .EQ. NZ) GOTO 2140 607 IF (I .EQ. NX) GOTO 2130 60B A = (1.-ABX)*(1.-ABY) 609 SUM(I,J) = SUM(I,J)+A 610 RSUM(I,J) = RSUM(I,J)+A*RHO(K) 611 2130 IF (IP .LT. 1 .OR. IP .GT. NX-1) GOTO 2150 612 A = ABX*(1.-ABY) 613 SUM(IP,J) = SUM(IP,J) + A 614 RSUM(IP,J) « RSUM(IP,J) • A*RHO(K) 615 2140 IF (JP .LT. 1 .OR. JP .GE. NZ) GOTO 2120 616 A « ABX*ABY 617 SUM(IP,JP) « SUM(IP,JP)+A 618 RSUM(IP,JP) = RSUM(IP,JP)+A*RHO(K) 619 2150 IF (JP .LT. 1 .OR. JP .GE. NZ .OR. I .EQ. NX) GOTO 2120 620 A = (l.-ABX)*ABY 621 SUM(I,JP) • SUM(I,JP)*A 622 RSUM(I,JP) » RSUM(I,JP)*A*RHO(K) 623 2120 CONTINUE 624 625 DO 2180,1=1,NX-2 626 K = 2*1+1 627 DO 2180,J • K,NZ-2 628 IF (SUM(I,J) .GE. 2.5) GOTO 2180 629 NP « NP*1 630 IF (NP .GT. 20000) STOP 'STOP IN INTPR' 631 632 C 633 C BREAK IN PROGRAM 634 C 635 636 X(NP) • (FLOAT(I)-0.5)*DX 637 Y(NP) • (FLOAT(J)-0.5)*DZ 638 XORG(NP)-X(NP) Appendix B. The Program 639 YORG(NP)«Y(NP) 640 21(NP) « FLOAT(I) 641 ZETA (NP) «= Z(I , J ) 642 RHO(NP) = R(I , J ) 643 SUM(I,J) « SUM(1,0) +1. 644 RSUM(I,J) • RSUM(I,J)+RHO(NP) 645 INK(NP) = -TIMEM00. 646 2180 R(I , J ) • RSUMd, J ) /SUMd, J ) 647 648 I » NX-2 649 J * NZ-2 650 DR = R(I-1 , J)-R(I-1 , J-1) 651 R d , J - l ) * Rd , J)-DR 652 DO 2190 1=1,NX-2 653 IF (SUMd,NZ-1) .GT. 3.5) GOTO 2200 654 NP « NP+1 655 X(NP) = (FLOAT(I)-0.5) *DX 656 Y(NP) -FLOAT(NZ-1)*DZ 657 XORG(NP)=X(NP) 658 YORG(NP/=Y(NP) 659 ZHNP) = FLOAT (I) 660 ZETA(NP) = 0. 661 RHO(NP) = -STRATF*(Y(NP)/2.+WO* TIME) 662 SUMd,NZ-1) = SUMd,NZ-1)+0.5 663 INK(NP) = -TIME*100. 664 2200 J = NZ-2 665 K = J - l 666 DR = Rd , J)-R(I,K) 667 2190 R(I,NZ-1) « R(I,NZ-2)*DR 668 DO 2210,1=1,NX-2 669 J = 2*1+2 670 K = 2*1+1 671 L = 2*1 672 M = 2*1-1 673 N = 2*1-2 674 DR = R (I, J ) -R (I, K) 675 Rd,L) = R(I,K) -DR 676 R(I,M) = R(I,L)-DR 677 IF (N -LT. 1) GOTO 2210 678 R(I,N) = R(I,M)-DR 679 2210 CONTINUE 680 I = NX-1 681 K « NZ-2 682 DRX • R(I-1,K)-R(I-2 ,K) 683 R(I,K+1) « R(1-1,K+l)+DRX 684 R(I,K( • R(I-1,K)+DRX 685 • R(I ,K-1) « R(I-1,K-1)+DRX 686 R(I ,K-2) « R(I-l,K-2)+DRX 687 RETURN 688 END 689 690 691 SUBROUTINE RELAX(PSI,Z,SLOPE,PSIDEL) 692 IMPLICIT REAL*8 (A-H, O-Z) 693 PARAMETER (INX-61 , INZ«121, INP« 20000) 694 DIMENSION PSI(INX,INZ),Z(INX,INZ) 695 DIMENSION PSIDEL(INX,INZ) 696 Appendix B. The Program 60 697 COMMON NX,NZ,DX,DZ 698 COMMON DXSQ,DZSQ,ORF,Q,SML,D2 699 700 C 701 C RELAX THE POISSON EQUATION BY SUCCESSIVELY PASSING THROUGH 702 C THE EULERIAN GRID IN FOUR WAYS: 703 C 1) FROM THE SLOPING BOTTOM TO THE SURF. LEFT TO RIGHT 704 C 21 FROM THE SLOPING BOTTOM TO THE SURF. RIGHT TO LEFT 705 c 3) FROM THE SURFACE TO THE SLOPING BOTTOM LEFT TO RIGHT 706 c 4) FROM THE SURFACE TO THE SLOPING BOTTOM RIGHT TO LEFT 707 c 708 709 c 710 c STORE PSI BEFORE ANY CHANGES ARE MADE 711 c 712 713 DO 1236 J«NZ,1,-1 714 DO 1236 I«1,NX 715 1236 PSIDEL(I,J) * PSI(I,J) 716 717 D3 •= 2.* (2. /DXSQ • 1 ./DZSQ) 718 DO 1715 MSWP • 1,1000 719 ERRMX "0.0 720 KDIR = MSWP-2*(MSWP/2) 721 IF (KDIR .NE. 1) THEN 722 DO 1660 KSWP=1,2 723 DO 1670 IP=2,NX-1 724 1 « IP 725 IF (KSWP .EQ. 2) I=NX+1-IP 726 K«2*l 727 DO 1670 JP*K,NZ-1 728 J • NZ-1+K-JP 729 PIJ * PSI(I,J) 730 P2 = 2.«PIJ 731 IF (J .NE. 2*1) THEN 732 ER ' (PSI(I,J+1) - P2 + PSI(I,J-1)> /DZSQ-Z(I,J) 733 ER • ER + <PSI(I+1,J) - P2 • PSI(I-1,J)) / DXSQ 734 PSI(I,J) « ORF*ER/D2 + PIJ 735 ELSE 736 ER « (Q-P2+PSI(I,J+l))/DZSQ-Z(I,J) 737 ER • ER+4./3.*(PSI(1-1,J)-3.*PIJ*2.*Q)/DXSQ 738 PSI(I,J) • ORF*ER/D3*PIJ 739 END IF 740 IF (ER .LT. 0.) ER « -ER 741 IF (ER .GT. ERRMX) ERRMX * ER 742 3670 CONTINUE 743 1660 CONTINUE 744 ELSE 745 DO 1695 KSWP=1,2 746 DO 1700 IP«2,NX-1 747 I * IP 748 IF (KSWP .EQ. 2) I«NX+1-IP 749 K « 2*1 750 DO 1700 J-K,NZ-1 751 PIJ - PSI(I,J> 752 P2 « 2.*PIJ 753 IF (J .NE. K> THEN 754 ER - (PSI(I,J+l)-P2+PSI(I,J-l))/DZSQ Appendix B. The Program 755 ER • ER-Z(I,J)*(PSI(I+1,J)-P2+PSI(I-1,J))/DXSQ 756 PSI(I,J) • ORF*ER/D2+PIJ 757 ELSE 75B ER ' (Q-P2+PSI(I,J+l))/DZSQ -Z(I,J) 759 ER*ER+(4.0/3.0)*(PSI(1-1,J) -3.0*PIJ+2.0*Q)/DXSQ 760 PSI(I,J) * ORF*ER/D3+PIJ 761 END IF 762 IF (ER .LT. 0.0) ER ' -ER 763 IF (ER .GT. ERRMX) ERRMX * ER 764 1700 CONTINUE 765 1695 CONTINUE 766 END IF 767 IF (ERRMX .LT. SML) GOTO 1730 76B C WRITE**,400) ERRMX 769 400 FORMAT(' CONVERGENCE ERROR= ',E15.6) 770 1715 CONTINUE 771 STOP 772 773 C 774 C BREAK IN PROGRAM 775 C 776 777 1730 CONTINUE 778 C WRITE(*,410) ERRMX 779 410 FORMAT(' ERRMX = ",E15.4) 780 C WRITE(*,420) MSWP 781 420 FORMAT(' NUMBER OF ITERATIONS = ',13) 782 MSWP = 1 783 784 C 785 C STREAM FUNCTION IS ASSIGNED AT POINTS EELOW AND PARALLEL TO 786 C THE BOUNDARY SO THAT FLOW THROUGH THE BOUNDARY DOES NOT OCCUR 787 C 788 789 IF (SLOPE .EC- 1.0) THEN 790 K « 2 791 DO 1740 I=K,NX-1 792 J « 2*1-2 793 PSI(I,J)«1.5«Q-(2.*PSI(1-1,J+l)+PSI(I- 1, J+2)+PSI(I,J+2))/8 794 J « J-1 795 PSI(I,J) « 2.0*Q-PSI(1-1,J+2) 796 J « J-1 797 IF <J .LT. 1) GOTO 1740 798 PS3(I,J) ' 2.«Q-(2.*PSI(1-1,J+3)+2. •PSI(1-2,J+3) 799 $ +PSI(1-2,J+2) + PSKI-2, J+4)+PSI(1-1,J+2) 800 $ +PSI(1-1,J+4))/8. 801 J • J-1 802 IF (J .LT. 1) GOTO 1740 803 PSI(I,J) - 2.*Q-PSI(1-2,J+4) 804 1740 CONTINUE 805 GOTO 1800 806 ELSE 807 DO 1780,1-2,NX-1 808 J « 2*1-2 809 PSI(I,J) « 2.*Q-PSI(I,J+2) 810 PSI(I,J-1) « 2.*Q-PSI(I,J+3) 811 1780 CONTINUE 812 DO 1790,I«3,NX-2 Appendix B. The Program 613 0 • 2*1-4 814 PSI(I,J) • 2.*Q-PSI U,0*6> 815 PSI(I,J-1) * 2.*Q-PSI(I,J + 7) 816 1790 CONTINUE 617 END IF 618 819 C 820 C CALCULATE CHANGE IN PSI 821 C ---822 823 1800 DO 1240 J*NZ,1,-1 824 DO 1240 1*1,NX 825 1240 PSIDEL(I,J> « <(PSKI,J) - PSIDEL (I, J) ) ) /Q 826 RETURN 827 END 828 829 830 831 SUBROUTINE INPUT(NX,N2,HTO,BVN,DT,TIMELM,VISCOS,DIFFUS,Q,IP, 832 $ IOUT,ORF,SLOPE,NCASE) 833 IMPLICIT REAL*8 (A-H,0-Z) 834 DIMENSION IOUT(100) 835 CHARACTER*1 ANS 836 CHARACTER*20 DUMM 837 C OPEN(1,FILE«'SRF.DAT',STATUS*'OLD') 838 WRITE(2,*) 639 READ(1,*) DUMM 84 0 READ(1,*) NX 841 WRITE(2,*) ' NX:',NX 842 READ(1,*) DUMM 843 READ(1,*) NZ 644 WRITE(2,*) ' NZ:',NZ 645 READ(1,*) DUMM 846 READ(1, *) HTO 847 WRITE(2,*) ' HEIGHT OF TANK:',HTO 848 READ(1,*) DUMM 849 READ(1,*) BVN 850 WRITE(2,*) ' BRUNT-VAISALA FREQUENCY:',BVN 851 READ(1,*) DUMM 852 READ(1,*) DT 853 WRITE (2,*) ' TIME INCREMENT:',DT 854 READ(1,*) DUMM 855 READ(1,*) TIMELM 656 WRITE(2,*) ' TIME LIMIT:',TIMELM 857 READU,*) DUMM 858 READU,*) VISCOS 859 WRITE(2,*) ' VISCOSITY:*,VISCOS 860 READ(2,*) DUMM 861 READ(1,*) DIFFUS 862 WRITE(2,*) ' DIFFUSIVITY:',DIFFUS 663 READ(1,*) DUMM 864 READ(1,«) Q B65 WRITE(2,*) ' DISCHARGE RATE:',Q 866 READ(1,*) DUMM 867 READ(1,*) IP 668 READU,*) (IOUT (I) »1*1»IP) 869 READU,*) DUMM 870 READU,*) ORF Appendix B. The Program 871 WRITE(2,*) * OVER RELAXATION FACTOR:',ORF 872 READd,*) DUMM 873 READd,*) SLOPE 874 WRITE(2, *) ' SLOPE:',SLOPE 875 WRITE(2,•) 876 WRITE(2,*) * CASE: (1) OUTLET AT BOTTOM OF DAM WALL' 877 WRITE(2,*) ' (2) OUTLET MIDWAY UP DAM WALL: ' 878 READ(1,*) DUMM 879 READd,*) NCASE 880 WRITE(2,*) ' NCASE: ',NCASE 881 882 C 883 C RESTRAINTS TO ENSURE THAT INSTABILITIES DO NOT OCCUR 884 C 885 886 DZSQ - (HTO/(NZ-1.0))**2 887 RESTRA - VISCOS*DT/DZSQ 688 WRITE(2,*) 889 WRITE(2,150) RESTRA 890 150 FORMAT(' A - (VISCOSITY)(dT)/(DZSQ) = ',F5.3> 891 IF (RESTRA .GT. 0.5) THEN 892 WRITE (2,*) 893 WRITE(2,*) ' UNSTABLE CONDITONS - CHANGE DATA FILE ' 894 WRITE(2,155) 0.5*DZSQ/VISCOS 895 155 FORMAT(' TIME INCREMENT SHOULD BE LESS THAN *,F9.4) 896 STOP 897 ENDIF 898 RETURN 899 END 900 901 902 903 SUBROUTINE OUTPUT(ROUT,RVN,RAMI6,TINO,X,Y,Zl,Z,VSCL, 904 $ INK,ZETA,RHO,NP,PSI,R,USCL,VRTSCL,PSIDEL) 905 IMPLICIT REAL*8 (A-H,0-Z) 906 PARAMETER (INX-61 , INZ-121, INP=20000) 907 DIMENSION PSI(INX,INZ),U(INX,INZ),V(INX,INZ) 908 DIMENSION Z(INX,INZ),ZT(INX,INZ),R(INX,INZ),RT(INX,INZ) 909 DIMENSION ROUT(INX) 910 DIMENSION ZETAdNP) ,RHO(INP) 911 DIMENSION X(INP),Y(INP),INK(INP),Zl(INP) 912 DIMENSION PSIDEL(INX,INZ) 913 914 COMMON NXrNZ,DX,DZ 915 COMMON DXSQ,DZSQ,ORF,Q,SML,D2 916 COMMON TIME,WO,STRATF,HTO 917 COMMON U,V,ZT 916 919 TVS « TIME*RVN*RAM16 920 TIN • TIME/TINO 921 922 C WRITE (8,*) NP 923 c WRITE (8,*) 924 c N10-INT(NP/10> 925 c DO 3333 I«1,N10 926 c ID-(I-1)*10 927 c WRITE (8,'(10F8.4)*) (X(ID+J),J-l,10) 928 c WRITE (8,'(10F8.4)') (Y(ID+J>,J-l,10) Appendix B. The Program 64 929 C WRITE (8,*(10F8 .4) ') (21(ID*J),J«l,10) 930 C3333 CONTINUE 931 C DO 3334 I«N10«10+1,NP 932 C3334 WRITE (8,'(3F8.4)') X(I), Y (I),21(I) 933 C GOTO 1251 934 C 935 C PRINT STREAMFUNCTION 936 C 937 WRITE(2,«) 93B WRITE(2,*) 'PSI AT:', TIME ,'USCL=',USCL 939 DO 1235,J-N2,1,-1 940 WRITE (2, ' (41F11.7) ') (PSI(I,J) /Q,I'1,NX) 941 1235 CONTINUE 942 WRITE (2,*) 943 C GOTO 19 944 C 945 C PRINT CHANGE IN STREAM FUNCTION 946 c 947 948 c WRITE(4,*) 'PSI INCREMENT' 949 950 c DO 1236,J=N2,1,-1 951 c WRITE(4,'(41F6.3)') (PSIDEL(I,J),1=1,NX) 952 C1236 CONTINUE 953 954 C 955 C PRINT DENSITY 956 C 957 958 19 WRITE(2,320) TIME,TVS,TIN 959 320 FORMAT C DENSITY FIELD AT TIMES: '.3F10.4) 960 DO 1250, JP=1,N2-1 961 J = N2-JP 962 WRITE (2,' (41F6.2) ') (R(I,J)/STRATF,1 = 1,NX-1) 963 1250 CONTINUE 964 WRITE(2,«) 965 C GOTO 1251 966 967 C 968 C PRINT DENSITY GRADIENT 969 c 970 c WRITE(2,*) 'DR/D2 AT:', TIME 971 c DO 1270,JP=2,N2-1 972 c J * N2+1-JP 973 c DO 1280, I«1,NX-1 974 C1280 ROUT(I) « (R(I,J-l)-R(I» J))/(DZ*STRATF) 975 C270 WRITE(2,'(41F6.2)') (ROUT(I),1=1,NX-1) 976 977 c 978 c PRINT VORTICITY 979 c 980 981 c WRITE<2,«) 'VORTICITY',VRTSCL 982 c DO 1290,JP«1,N2 983 c J • N2+1-JP 984 c WRITE(2, ' ( 4 1 F 1 1 . 7 ) ' ) (2(I,J)/VRTSCL,I«l,NX) 985 C1290 CONTINUE 986 C GOTO 29 Appendix B. The Program 65 987 988 C 989 C PRINT RATE OF CHANGE OF VORTICITY 990 C 991 C535 WRITE(2, *) ' RATE OF CHANGE OF VORTICITY' 992 C DO 1295,JP«1,NZ 993 C J • NZ+l-JP 994 C WRITE(2,* U1F6.2) ') (2T(I,J),1*2,NX) 995 C1295 CONTINUE 996 997 C 998 C PRINT INK 999 C 1000 C WRITE(2,*) 'INKED ' 1001 C DO 1310,K=1,NP 1002 C IF (INK(K) .EQ. 0) GOTO 1315 1003 C WRITE(2,») 'INK(K) X(K) Y(K) ZETA(K) RHO(K): 1004 C WRITE(2,360> INK(K) ,X(K) ,Y(K) ,ZETA (K) ,RHO(K) 1005 C360 FORMAT (I5,2F10.6,F13-3,F10.4) 1006 C1310 CONTINUE 1007 C1315 CONTINUE 1008 1009 C 1010 C PRINT HORIZANTAL VELOCITY U 1011 C 1012 1013 C29 WRITE(2,*> 'U' 1014 1015 C DO 1320,JP=1,NZ-1 1016 C J « NZ-JP 1017 C DO 1330,1=1,NX 1018 C ROUT(I) = U (I,J)* 100.0/USCL 1019 C1330 CONTINUE 1020 c WRITE(2,'(41F6.1)') (ROUT(I),1=1,NX) 1021 C1320 CONTINUE 1022 1023 C 1024 c PRINT VERTICAL VELOCITY V 1025 c 1026 1027 c WRITEC2,*) *V 1028 1029 c DO 1340,JP=1,NZ 1030 c J * NZ+l-JP 1031 c WRITE(2,*(41F6.2)') (V(I, J)/VSCL,1=1,NX-1) 1032 • C1340 CONTINUE 1033 1251 RETURN 1034 END Appendix C Figures Figure 1.2: Withdrawal layer Appendix C. Figures and Tables 6 7 Figure 1.3: Y i h ' s solution for F = 0.35 SINK Figure 1.4: Kao 's method for F < IT Appendix C. Figures and Tables 68 Figure 1.5 : Potential flow ^ M o d * N i*Tto * r > ; v (Q) Shear Fronts (b) Velocity Field Figure 1.6: Wi thdrawal layer formation Appendix C. Figures and Tables 69 X=0 X=L Q Y= Figure 2.1: Definition sketch of reservoir Figure 2.2: Purely buoyant solution Appendix C. Figures and Tables 70 (1.NY) (2, NY) (3.NY) (4,NY) (5.NY) (NX.Nyf (1.NY-1) (2.NY-1) (3.NY-1) (4.NY-1) (5.NY-1) i^XNY-1) (1.10) (2.10) (3,10) (4,10) ( 5 , 1 0 ) / / (NX.NY-2) (1.9) (2,9) (3,9) (4,9) / 1 5 . 9 ) (NX.NY-3) (1.8) (2,8) (3,8) ( 4 , 8 ) / (5,8) (1.7) (2,7) (3,7) (5,7) (1.6) (2,6) ( 3 , 6 y ^ (4,6) (1.5) (2,5) (4.5) (1,4) ( 2 , 4 ) / (3,4) • (1.3) / ( 2 . 3 ) (3,3) ( 1 , 2 ) , / (2,2) (2,1) Figure 3.1: Computat ional Domain Appendix C. Figures and Tables 71 PSI(I,J+1) 2(I,J+1) V(I,J+1) U(I.J) PSI(I.J) 2(I,J) PSI(I+1,J+1) Z(I+1,J+1) U(I+1,J) V(I.J) PSI(I+1,J) 2(1+1 ,J) Figure 3.2: Staggered G r i d Appendix C. Figures and Tables 72 U(I,JP) B U(I+1,JP) DY - ETA ETA X(k),Y(k) U(I,J) XI A U(I+1,J) D X - X I Figure 3.3: Horizontal velocity interpolation V(I,J+1) PSI(I+1,J+ / S L O P I N G BOTTOM \Y V(I + 1,J+1) Y • m / A / I « — 1 i • PSI(I,J) V(I,J) PSI(I + 1,J) V (M ,J ) Figure 3.4: A Point Near The Sloping Boundary Appendix C. Figures and Tables 73 PSI-0 PSI-Q PSl-O.20 PSU0.5O PSl-O.80 PSI-Q PSl-1.20 PSI.1.50 PSl.1.80 PSI-Q PSI-2Q PSI(B)-2Q-PSI(A) Figure 3.5: Mirror-Image Technique DX / ETA DY-ETA XI XfM.Yflri ETA 1 GRID CELL J+1 DY AREA A(I,J) DX-XI XI +1 MARKER PARTICLE Figure 3.6: Updat ing Vort ic i ty and Density Appendix C. Figures and Tables Figure 3.7: The location of the two points used for the grid testing Figure 3.8: G r i d Dependance - coarse mesh at the point x = 0.2L and y = 0. Appendix C. Figures and Tables Figure 3.10: G r i d Dependance - fine mesh at the point x = 0.2L and y = 0.2H Appendix C. Figures and Tables a GRID 60X120 o GRID 55X110 , GRID 50X100 „ GRID 45X90 Figure 3.11: G r i d Dependance - fine mesh at the point x = 0.2L and y = 0.6H z • o < o TIME STEP=0.025 s # TIME STEP=0.05 s 0 0 10 0 «0 0 30 0 40 0 M O 60 0 70 0 80.0 90.0 100 0 110 0 180 0 TIME sec Figure 3.12: T ime Step Dependance at x = 0.2L and y = 0.2H for a 60x120 grid Appendix C. Figures and Tables Figure 4.1: The location of the nine points at y = 0.5H Figure 4.2: The location of the seven points at y = 0.75/f Appendix C. Figures and Tables 0 0 5.0 10 0 ISO 20.0 25 0 30 0 35 0 40 0 45 0 50 0 55 0 60 0 TIME sec (TIME STEP = 0 05 sec) Figure 4.3: t/> versus t at y = 0.5H for the points listed in Table 4.1 100 15.0 20 0 25 0 30 0 35 0 40 0 TIME sec (TIME STEP = 0 05 sec) SO 0 60 0 Figure 4.4: xp versus t at y = 0.75H for the points listed in Table 4.1 Appendix C. Figures and Tables 79 0.0 10 8.0 3 0 4.0 6 0 6.0 7.0 TIME n/N sec. Figure 4.5: Propagation of the first shear wave Figure 4.6: The location of the eleven points at x = 0.25L Appendix C. Figures and Tables 80 Figure 4.8: A displaced spherical particle in a stratified fluid Appendix C. Figures and Tables Figure 4.10: rp versus t for points near the sloping bottom Appendix C. Figures and Tables 82 Figure 4.12: Streamlines at t = 5.0 seconds, R u n l Appendix C. Figures and Tables Figure 4.13: Streamlines at t = 10.0 seconds, R u n l Figure 4.14: Streamlines at t = 20.0 seconds, R u n l Appendix C. Figures and Tables 84 Figure 4.15: Streamlines at t = 30.0 seconds, R u n l Figure 4.16: Streamlines at t = 40.0 seconds, R u n l Appendix C. Figures and Tables 85 Figure 4.18: Streamlines at t = 120.0 seconds, R u n l Appendix C. Figures and Tables Figure 4.19: Definition sketch for 6+ and 8~ -2.5 0.0 2.5 5 0 VELOCITY (U)M00/USCL (m/s). 7.5 X = 0.25L. 10 0 175 Figure 4.20: Horizontal velocity profile at t = 15.0 seconds at x = 0.25L Appendix C. Figures and Tables Figure 4.22: Horizontal velocity (t = 35.0 s, x = 0.251), Runl Figure 4.23: Horizontal velocity (t = 45.0 s, x = 0.251), R u n l o -2.0 -1 0 00 1 0 20 30 40 50 60 7.0 80 90 100 11.0 12 0 130 14 0 150 160 VELOCITY (U)*100/USCL (m/s). X = 0 25L Figure 4.24: Horizontal velocity (t = 50.0 s, x = 0.25L), R u n l Appendix C. Figures and Tables > x o J5 •» -20 2 0 4 0 tO 8 0 10 0 12 0 VELOCITY (U)M00/USCL (m/ s). X = 0 25L. is o Figure 4.26: Horizontal velocity (t = 85.0 s, x = 0.25L), Runl Appendix C. Figures and Tables e> o o I i t • , , , I -20 00 20 40 60 80 100 12 0 140 160 VELOCITY (U)')00/USCL (m/s). X = 0 25L Figure 4.27: Horizontal velocity (t = 120.0 s, x = 0.25L), R u n l o l \ l 6-o *•• I l , l l -3 0 -2 0 -1.0 0 0 1 0 2 0 3.0 4 0 SO 6 0 7 0 8 0 9 0 10 0 110 VELOCITY (U)*100/USCL (m/s), X = O.fSOL Figure 4.28: Horizontal velocity (t = 25.0 s, x = 0 .5L) , R u n l Appendix C. Figures and Tables Figure 4.30: Horizontal velocity (t = 40.0 s, x = 0.5L), R u n l Appendix C. Figures and Tables 0 5 1.0 15 2 0 2 5 3 0 3 5 4.0 4.5 5 0 5 5 6 0 6 5 7 0 VELOCITY (U)*100/USCL (m/s). X = 0.50L. Figure 4.31: Horizontal velocity (t = 50.0 s, x = 0.5X), R u n l i.o 2 0 3 0 4 0 5 0 6 0 VELOCITY (U)M00/USCL (m/s). X = 0 50L, Figure 4.32: Horizontal velocity (t = 60.0 s, x = 0 .5L) , R u n l Figure 4.34: Horizontal velocity (t = 120.0 s, x = 0.5L), Runl Appendix C. Figures and Tables 94 i i i i i i 10 0 20.0 3 0 0 40 0 60 0 6 0 0 70 0 (0 0 TIME sec (TIME STEP = 0.05 sec) 1200 Figure 4.35: Peak horizontal velocity at 3 stations, Runl • X=0.30L o X=035L o X=0.40L e o — e -o.o i i i i i 400 M>.0 60.0 70.0 60.0 TIME sec (TIME STEP = 0.05 sec) •o.o 100.0 110 0 120.0 Figure 4.36: Peak horizontal velocity at 3 stations, Runl Appendix C. Figures and Tables • X = 0 4 5 L o X=0.50L o X=0 55L o o i 10.0 I I I I I 40 0 60 0 60 0 70 0 80 0 TIME sec (TIME STEP = 0.05 sec) too.o 110.0 120 0 Figure 4.37: Peak horizontal velocity at 3 stations, Runl I X=015L o X=0.20L n X=0.25L S! N o DC o o z ; + o •O S. 30.0 400 600 600 70.0 60.0 TIME (s). TIME STEP=0.05 s •oo Figure 4.38: c5+ at 3 stations, Runl 100.0 n o o 120.0 Appendix C. Figures and Tables Q o <- n OS ° 1° =•2 4 o •o # X=0.30L o X=0 35L „ X=040L eo.o 30.0 40 0 60.0 60 0 TOO 60 0 TIME (s). TIME STEP=0 05 s •0.0 100.0 110.0 1200 Figure 4.39: S+ at 3 stations, Runl CO d ° N d J a* a: ° g o + d X=045L X=0.50L X=0.55L oo i 10.0 40 0 60.0 60 0 70 0 80.0 TIME (s). TIME STEP=0.05 s 100.0 110 0 120.0 Figure 4.40: c5+ at 3 stations, Runl Appendix C. Figures and Tables Figure 4.42: Wi thdrawn particles in i t ia l position after 120 s, R u n l Appendix C. Figures and Tables 98 Figure 4.43: Streamlines at t = 40.0 seconds, Run2 Figure 4.44: Streamlines at t = 80.0 seconds, Run2 Appendix C. Figures and Tables 99 Figure 4.46: rp versus t at y = 0.5// for the points listed in Table 4.1, Run2 Appendix C. Figures and Tables Figure 4.47: Averaged peak horizontal velocity comparison between flows with N = 0.0764 and 0.764, Runl and Run2 Appendix C. Figures and Tables Figure 4.49: Streamlines at t = 20.0 seconds, Run3 Appendix C. Figures and Tables 102 Figure 4.50: Streamlines at t = 120.0 seconds, Run3 p Figure 4.51: Horizontal velocity (t = 15.0 s, x = 0 .25L), Run3 Appendix C. Figures and Tables o > X o - 3 0 - 2 0 - 1 0 0 0 1 ^ 0 l o * 0 4 0 5 0 8 0 T O 8 0 S O 10.0 1 1 0 12 0 13.0 14.0 VELOCITY (U)*100/USCL (m/s). X = 0 25L. ' Figure 4.52: Horizontal velocity (i = 35.0 s, x - 0.25L), Run3 J j O - i o 00 1.0 20 30 4 0 50 60 70 80 9.0 100 11.0 12 0 13 0 14 0 15 0 VELOCITY (U)*100/USCL (m/s). X = 0.25L -Figure 4.53: Horizontal velocity (t = 120.0 s, x = 0.25L), Run3 Figure 4.54: Horizontal velocity (i = 35.0 s, x = 0.5L), Run3 VELOCITY (UrMOO/USCL (m/s). X = 0.5L.' Figure 4.55: Horizontal velocity (t = 40.0 s, x = 0.5L), Run3 Appendix C. Figures and Tables 105 h H X o u ' X 1 ' 2 5 3 0 3 5 «° 4 5 5 0 5 5 6 0 6 5 °° 0 5 1 0 'VELOCITY (U)'IOOAJSCL (m/s), X = 0.5L Figure 4.56: Horizontal velocity (t = 120.0 s, i = 0.5L), Run3 8 ° O n O O UUAX.Q=65E-3 # UUAX,Q=35E-3 .00 0 05 0.10 015 DISTANCE (0 Figure 4.57: Averaged peak horizontal velocity comparison between R u n l and Run3 Appendix C. Figures and Tables 106 Figure 4.59: Streamlines at t = 0.0 seconds, Run4 Appendix C. Figures and Tables Figure 4.60: Streamlines at t = 5.0 seconds, Run4 e VELOCITY (U)MO0/USCL (m/s). X = 0.25L. Figure 4.61: Horizontal velocity (t = 5.0 s, x = 0.25L), Run4 Appendix C. Figures and Tables 108 Figure 4.63: Streamlines at t = 10.0 seconds, Run4 Appendix C. Figures and Tables o - IM 6 o O ' i , i , , , , , , || -6.0 -4.0 -2.0 0 0 2 0 4 0 6.0 BO 10 0 12.0 14 0 16.0 18 0 VELOCITY (U)MO0/USCL (m/s). X = 0.25U Figure 4.64: Horizontal velocity (t = 10.0 s, x = 0.25L), Run4 o VEUXTTV (U)*100/USCL (m/s), X = 05L. Figure 4.65: Horizontal velocity (t = 10.0 s, x = 0 .5L) , Run4 Appendix C. Figures and Tables Figure 4.66: Streamlines at i = 15.0 seconds, Run4 o zoo Figure 4.67: Horizontal velocity [t = 15.0 s, x = 0.25L), Run4 Fig' ;ure 4.69: Streamlines at t = 20.0 seconds, Run4 Appendix C. Figures and Tables o ° l I 1 1 I I 1 1 -7.5 -5 0 -2 5 0 0 2 5 5 0 7.5 10 0 12 5 ISO 17 5 20 0 VEU>OTY (U)'100/USCL (m/s). X = 0 25L. ' Figure 4.70: Horizontal velocity (t = 20.0 s, x = 0.251), Run4 o - J , . , , , , , , , , , 1 -1.0 -0.5 0 0 0 5 1 0 1 5 2 0 2 5 3 0 3 5 4 0 4.5 6.0 5.5 6 0 6 5 VELOCrrY (U)*100/USCL (m/s). X « 0.5L, Figure 4.71: Horizontal velocity (t = 20.0 s, x = 0.5L) , Run4 Appendix C. Figures and Tables M 6-o O T , i i • -* 0 -2.0 0 0 2 0 4 0 6.0 8 0 10 0 12 0 14 0 18.0 18.0 VELOCITY (U)*I00/USCL (m/s). X = 0 25L. " Figure 4.72: Horizontal velocity (t - 30.0 s, i = 0.25L), Run4 o -4.0 -2 0 0.0 2 0 4 0 6 0 8 0 10 0 12.0 14 0 16 0 18.0 20 0 22.0 VELOCITY (U)M00/USCL ( m /s) . X - 0.25L. Figure 4.73: Horizontal velocity (t = 120.0 s, x = 0.25L), Run4 Appendix C. Figures and Tables Figure 4.74: Horizontal velocity (t = 25.0 s, x = 0 .5L), Run4 VELOCITY (U)MOOAISCL (m/s). X = 0.5L, Figure 4.75: Horizontal velocity (t = 30.0 s, x = 0 .5L) , Run4 Appendix C. Figures and Tables o CO 6 -o 6 ^ ; — i i i , , 1 , 1 i i i I -0 5 0 0 0-5 1.0 1.5 2 0 2 5 3 0 3 5 4 0 4 5 6 0 5 5 VELOCITY (U)*I00/USCL (m/s). X = 0 5L. ' Figure 4.76: Horizontal velocity (t = 35.0 s, z = 0.5L), Run4 6-o ° J 1 1 1 1 r — 1 , 1 1 1 I 0.00 025 0.50 075 1.00 1.25 1.50 175 2 00 2 25 2 50 275 3 00 3 25 3.50 3 75 4 00 VELOCITY (U)MO0/USCL (m/s). X = 0.5L Figure 4.77: Horizontal velocity (t = 45.0 s, z = 0.5X), Run4 Appendix C. Figures and Tables o.o 1.0 15 2 0 2 5 3.0 VELOCITY (U)M00/USCL (m/s). X = 0 5L. 3 5 « o Figure 4.78: Horizontal velocity (t = 55.0 s, x = 0 .5L), Run4 1 0 1 5 2 0 2.5 30 3 5 VELOCTTY (U)M00/USCL (m/s). X - 0 5L " Figure 4.79: Horizontal velocity (t = 70.0 s, z = 0.51), Run4 Figure 4.80: Horizontal velocity (t = 75.0 s, x = 0 .5L), Run4 1 0 I S 8 0 2 5 3 0 VELOCrrY (U)M00/USCL (m/s). X •= 0.5L.' Figure 4.81: Horizontal velocity (t = 95.0 s, x = 0.51), Run4 Appendix C. Figures and Tables O U -X ° " tu o-o ° 1 I I I I I I 1 I—-— I I 0 0 0 5 1.0 1.5 2 0 2 5 3 0 3 5 4 0 4.5 SO VELOCITY (U)*I00/USCL (m/s). X = 0 5L Figure 4.82: Horizontal velocity (t = 120.0 s, x = 0.5X), Run4 PJ b • e ^ f , , -6 0 -6 0 -4 0 -2 0 0 0 2 0 4 0 0 0 B0 10 0 12 0 14.0 16 0 VELOCITY (U)*100/USCL (m/s). X = 0.25L Figure 4.83: Horizontal velocity (t = 40.0 s, x = 0.25L), Run5 Appendix C. Figures and Tables o O " o -6 0 -4.0 -2.0 00 20 40 60 60 10 0 12 0 14 0 VELOCITY (U)*100/USCL (m/s). X = 0.25L. Figure 4.84: Horizontal velocity (t = 50.0 s, x = 0.25L), Run5 o 6 -o °1 I i » 1 1 i 1 r — . 1 r — i i 1 • ' -3 0 -Z 0 -1.0 0 0 1 0 2 0 3 0 4.0 5 0 6 0 7 0 60 9.0 10.0 11.0 12.0 13 0 14 0 VELOCITY (U)M00/USCL (m/s). X = 0.25L Figure 4.85: Horizontal velocity (t = 65.0 s, x = 0.25L), Run5 Appendix C. Figures and Tables 2 oL . , , . , , , , I -2.0 - 1 0 00 1 0 20 30 4 0 50 8.0 7.0 BO 9 0 10 0 11 0 12 0 13 0 14 0 VELOCITY (U)'100/USCL (m/s). X = 0.25L. " Figure 4.86: Horizontal velocity (t = 120.0 s, x = 0.251), Run5 SJ , , , , , , 1 0.0 0 5 1 0 1.5 2 0 2 5 3 0 3.5 4 0 4.5 5.0 VELOCITY (U)*100/USCL (m/s). X - 0.5L ' Figure 4.87: Horizontal velocity (t = 25.0 s, x = 0 .5L) , Run5 Appendix C. Figures and Tables 00 10 20 30 4 0 50 6 0 70 00 90 100 11.0 120 VELOCITY (U)'100/USCL (m/s). X = 0 5L ' Figure 4.88: Horizontal velocity (t — 50.0 s, z = 0.5X), Run5 -30 -20 0 0 10 2 0 3 0 4 0 50 «0 VELOCITY (U)*100/USCL (m/s). X = 0.5L Figure 4.89: Horizontal velocity (t = 65.0 s, z = 0.5L), Run5 Appendix C. Figures and Tables > t-x u S3 • X < 0 0 0 5 |7 ^ To 7& 30 35 4.0 45 5.0 55 60 65 VELOCITY (U)MOOAISCL (m/s). X = 0.5L Figure 4.90: Horizontal velocity (t = 85.0 s, x = 0.5L), Run5 o Figure 4.91: Horizontal velocity (t = 100.0 s, x = 0.5L), Run5 Appendix o C. Figures and Tables x o w • x c 00 0 5 10 1 5 20 2 5 30 35 4.0 4 5 50 55 €0 «.5 '0 75 VELOCITY (U)'I00/USCL (m/s), X = 0.5L Figure 4.92: Horizontal velocity (t = 115.0 s, x = 0.5L) , Run5 VELOCITY (U)MO0/USCL (m/s). X = 0.5L. " Figure 4.93: Horizontal velocity (t = 120.0 s, x = 0 .5L) , Run5 Appendix C. Figures and Tables I ". o I >-CD O U N „ 1° 5 CC o z-. X=0.15L o X=0.20L o X=0 25L IT i i i i i i i i i 30.0 4 0 0 6 0 0 6 0 0 70 0 60 0 00 0 100 0 110.0 120 TIME (s). TIME STEP=0.05 s i i 100 £00 Figure 4.94: 6+ at 3 stations, Run5 o 6 I I I 1 1 I ' I 1 I . . m 10.0 20 0 30.0 40.0 60.0 60 0 70.0 60.0 S0 0 100.0 110.0 120 TIME (s), TIME STEP=0.05 s Figure 4.95: <5+ at 3 stations, Run5 Appendix C. Figures and Tables 125 Figure 4.97: Average «5+ comparisons between Runs 1,4 and 5 Appendix C. Figures and Tables Figure 4.98: 6+ at t = 120.0 seconds for A = 1, 4 and 10
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Selective withdrawal of a linearly stratified fluid...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Selective withdrawal of a linearly stratified fluid in a triangular reservoir Hnidei, Stephen D. 1990
pdf
Page Metadata
Item Metadata
Title | Selective withdrawal of a linearly stratified fluid in a triangular reservoir |
Creator |
Hnidei, Stephen D. |
Publisher | University of British Columbia |
Date Issued | 1990 |
Description | The water in most reservoirs is density stratified with depth. This stratification leads to the inhibition of vertical movement, consequently, when water is withdrawn from the reservoir it tends to move in a jet-like layer called a withdrawal layer, towards the sink. By placing the sink at a certain depth, one is able to selectively withdrawal water from a limited range of depths and thus obtain water of a desired quality. Much work has been done in this field by considering a simplified boundary geometry, usually rectangular. However little attention has been given to the effects of accurate boundary geometry. For this thesis, five numerical experiments were conducted for the problem of a two-dimensional, viscous, incompressible, slightly-stratified flow towards a sink in a triangular reservoir. |
Subject |
Water withdrawals -- Mathematical models Stratified flow -- Mathematical models |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-09-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080424 |
URI | http://hdl.handle.net/2429/28834 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1990_A6_7 H55.pdf [ 5.55MB ]
- Metadata
- JSON: 831-1.0080424.json
- JSON-LD: 831-1.0080424-ld.json
- RDF/XML (Pretty): 831-1.0080424-rdf.xml
- RDF/JSON: 831-1.0080424-rdf.json
- Turtle: 831-1.0080424-turtle.txt
- N-Triples: 831-1.0080424-rdf-ntriples.txt
- Original Record: 831-1.0080424-source.json
- Full Text
- 831-1.0080424-fulltext.txt
- Citation
- 831-1.0080424.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080424/manifest