NUMERICAL MODELING OF RESONANT INTERNAL SEICHETRIGGERED BY UNSTEADY WITHDRAWALbyLI GUB.Eng., Hohai University, China, 1984M.Eng., Hohai University, China, 1989A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE STUDIES(Department of Civil Engineering)We accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAOctober 1993© Li Gu, 1993In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature) Department of ^C 4; / Ey ineiwin jThe University of British ColumbiaVancouver, CanadaDate /2 Octo her, #993DE-6 (2/88)AbstractMost lakes and reservoirs are temperature stratified for some period during the year, atwhich time a typical fresh water body can be divided into a well-mixed, warm, surface layercalled the epilimnion, which is separated from the colder and relatively stagnant bottomwater, called the hypolimnion, by an intermediate region with strong temperature gradientknown as the metalimnion, which in conditions of temperature stratification only, issometimes called the thermocline. Internal standing waves or internal temperature seicheswithin a stratified lake or reservoir play an important role in the physical and biologicaldynamics.Most previous studies have considered wind-generated internal seiches. Little attentionhas been given to the problem of triggering internal seiche by unsteady withdrawal, whichwas first proposed by Imberger (1980) as a possible pump-back scheme to optimize thewater quality in a lake or reservoir. The resonant internal seiche triggered by unsteadywithdrawal is studied numerically in this thesis using the C/C (Cloud-In-Cell) method. Theinvestigation is concerned with internal seiches in an viscous, incompressible and weaklystratified fluid which is confined within a rectangular tank with a horizontal line sink in anend wall. Two-layered fluid with a sharp interface and two-layered fluid with a diffuseinterface are investigated. Two-layered fluid with a sharp interface approximates conditionsof lakes and reservoirs in late summer or early autumn when the thermocline is quite thin,while two-layered fluid with a diffuse interface approximates conditions of lakes andreservoirs in spring and late autumn, when the thermocline can be quite thick.The CIC method used in present study combines the best features of both Lagrangianand Eulerian approaches, and is most advantageously applied to fluid dynamics problemswith density stratification. The numerical model solves, in finite difference form, thegoverning equations in terms of stream function and vorticity, and therefore avoids thecomplexity of solving the pressure field. The free surface boundary conditions in this caseare somewhat difficult to apply. However, for the case of small density variation, simplifiedfree surface boundary conditions may be chosen without affecting the main features of theflow.A two-layered fluid with a sharp interface supports only one vertical internal wavemode. With the fluids in the upper and lower layers moving in the opposite directions, largevelocity shear develops at the density interface. The introduction of a continuously stratifiedinterfacial region between two homogeneous fluids enables more than one vertical internalwave mode to exist, although the energy is usually partitioned into the first few internalwave modes. It is found that the natural frequency of internal waves decreases withincreasing interface thickness and viscosity. The natural periods of a two-layered fluidsystem with an exponentially stratified interfacial layer are well predicted by the numericalmodel, in comparison with the theoretical solutions of Fieldstad (1933). By discharging fluidfrom the tank, two distinct classes of internal waves are created, namely, forced and freeinternal waves. The forced internal wave follows the forcing discharge and sustains as longas the discharge is retained, while the free internal wave oscillates according to its naturalfrequency with decaying amplitude due to viscous damping. Periodic beating is observed asa result of the interaction between the forced and free internal waves. The resonant internalseiche is created when the frequency of unsteady forcing discharge corresponds to thenatural frequency of internal waves.ifiTable of ContentsAbstract^ iiTable of Contents^ ivList of Symbols viiList of Tables^ xList of Figures xiAcknowledgments^ xviChapter 1. Introduction 1Chapter 2. Literature Review^ 52.1 Introduction^ 52.2 Free Internal Seiches^ 62.2.1 Linear Inviscid Internal Seiches^ 72.2.2 Finite Amplitude Internal Seiches 92.2.3 Viscous Effect^ 102.2.4 The Effect of Finite Interface Thickness^ 112.3 Forced Internal Seiches^ 162.3.1 Mass-Spring-Dashpot System Analogy^ 162.3.2 Forced Internal Seiches^ 19Chapter 3. Governing Equations 223.1 Introduction^ 22iv3.2 Governing Equations in terms of Primitive Variables^ 233.3 Governing Equations in terms of Stream Function and Voracity^263.4 Boundary and Initial Conditions^ 28Chapter 4. Numerical Method^ 314.1 Introduction^ 314.2 Computational Procedures^ 334.2.1 The Eulerian Calculation 344.2.2 The Lagrangian Calculation^ 384.2.3 Updating Eulerian Field Variables pAij) and g(ij)^404.3 Accuracy^ 414.3.1 Truncation Error^ 424.3.2 Round-off Error 434.3.3 Artificial Viscosity and Diffusivity Errors^ 434.3.4 Aliasing Error^ 444.3.5 Free Surface Boundary Treatment^ 454.4 Stability Consideration^ 45Chapter 5. Results and Discussion 485.1 Free Internal Seiche in a Two-layered Fluid with a Sharp Interface^495.2 Forced Internal Seiche Triggered by Constant Withdrawal^525.2.1 Experiment 1^ 535.2.2 Experiment 2 565.3 Forced Internal Seiche Triggered by Unsteady Withdrawal^575.3.1 Experiment 1; Comparison with Physical Experiments 585.3.2 Experiment 2^ 605.3.3 Experiment 3 615.3.4 Experiment 4^ 63Chapter 6. Conclusions and Recommendations^ 66Bibliography^ 69Tables 73Figures^ 76viList of Symbolsa^amplitude of the internal seicheC^Courant numberC/C^cloud in cell methodd diffusion numberF^densimetric Froude numberg^acceleration due to gravityg' reduced gravityhd^elevation of the discharge outlethi^depth of layer i, i=1=top layer, i=2=bottom layerh' distance between discharge and density interfaceH total depth of the fluidK molecular diffusivityk^wave numberL length of the tankMAC marker and cell methodN Brunt-Vdisallã frequency/VX^number of divisions in x directionNY^number of divisions in y directionn^wave mode numberP pressurePIC^particle in cell methodQ^discharge rate per unit widthviiQave^average discharge rate of unsteady withdrawal processSOR^successive over relaxationt^timeT^wave periodTI^forcing periodu horizontal velocity component^ vertical velocity componentx^horizontal coordinatey^vertical coordinate45^interfacial layer thicknesszip^density difference between top and bottom layersdx^mesh size in x directionAy^mesh size in y directionAt^time stepAQ^amplitude of unsteady withdrawal processe convergence criterion of SORlli^displacement of the internal interfaceris^displacement of the free surfaceA^wave lengthAt^dynamic amplification factorPi^density of layer i, i=1=top layer, i=2=bottom layerPo^mean densitya^wave frequencyye^eddy viscosity^ molecular viscosityto^forcing frequency of the mass-spring-dashpot systemviii0) R^over-relaxation factorV'^stream functiong^vorticityixList of TablesTable (5-1)Table (5-2)Table (5-3)Table (5-4)Table (5-5)Free internal seiche in a two-layered fluid with a sharpinterface; experimental parameters and results^73Forced internal seiche triggered by constant withdrawal,experiment 1; experimental parameters 73Forced internal seiche triggered by constant withdrawal,experiment 1, run 1; comparison of observed and predictednatural periods of fundamental internal wave modes 74Forced internal seiche triggered by constant withdrawal,experiment 2; experimental parameters and results^74Forced internal seiche triggered by unsteady withdrawal;experimental parameters^ 75xList of FiguresFigure (1-1)Figure (2-2)Figure (2-3)Figure (2-6)Figure (2-7)Figure (2-8)Temperature profile from Windermere northern basin showingthe typical temperature profile during summer (redrawn fromMortimer, 1952) 76Defining sketch of a two-layered reservoir model^77Location map for the Kootenay river drainage basin (fromCarmack & Gray, 1982)^ 78(a) Hourly positions of the thermocline bounded by the 9° and110c isotherms on a longitudinal section of Loch Earn,redrawn from Wedderburn (1912), (b) general flow pattern offree internal seiches in a two-layered lake, from Wedderburn& 'Williams (1911)Interface profiles of internal seiches in a two-layered fluidwith a sharp interface; (a) first or fundamental horizontal mode(n=1), (b) second horizontal mode (n=2), (c) third horizontalmode (n 3) 80The variation of tanhkhi and tanhkh2 for which the variation ofnatural frequency with wave amplitude changes. In region A,the frequency decreases with increasing wave amplitude; Inregion B, the frequency increases with increasing waveamplitude (from Thorpe, 1968a) 81Viscous effect on the natural frequency of internal seiche in atwo-layered deep fluid^ 81Schematic plot of a three-layered fluid system; (a) the densityprofile, (b) the associated two vertical internal wave modes,(c) the oscillation mechanism of the first vertical wave mode,(d) the oscillation mechanism of the second vertical wavemode (redrawn from Stevens & Imberger, 1993)^82Density profile of a two-layered fluid system with aexponentially stratified interfacial layer^ 83(a) Hyperbolic tangent density profile in deep fluid, (b)associated dispersion relationship 84Linear damped vibration of mass-spring-dashpot system: (a)dynamic amplification factor g and (b) phase angle a, asfunctions of frequency ratio ova° and damping factor 2/3 (fromWilson, 1972)Figure (1-2)Figure (1-3)Figure (2-1)Figure (2-4)Figure (2-5)7985xiFigure (5-3)Figure (5-4)Response curves for internal seiches in deep fluids with 4p/p2= 4.7*10-3, where 2all is the wave steepness. The symbolsrepresent different total plunger amplitudes (from Thorpe,1968a)Response curves for internal seiches in shallow lower fluid,kh2=0.48, and deep upper fluid, and with dp/p2 = 4.7*10-3(from Thorpe, 1968a)Boundary conditions for stream function yrStaggered gridBilinear interpolationA sketch showing locations of points A and BGrid size dependence tests on zir ; stream function versus timeat point AGrid size dependence tests on dy, , coarse mesh; steam functionversus time at point BGrid size dependence tests on dy, fine mesh; stream functionversus time at point BTime step et dependence tests; (a) stream function versus timeat point A, (b) stream function versus time at point BTime history of interface elevation at xr,) showing aliasingerrors increase with decreasing viscosityFree internal seiche in a two-layered fluid with a sharpinterface, run 1; normalized stream function versus time atthree different stationsFree internal seiche in a two-layered fluid with a sharpinterface; comparison of observed fundamental natural periodswith the theoretical values based on the linear inviscid theoryof equation (2-3) 94Free internal seiche in a two-layered fluid with a sharpinterface, run 3 and run 9; time series of interface elevation atx=L showing the damping process of the free internal seiches 95Free internal seiche in a two-layered fluid with a sharpinterface, run 3 and run 9; surface plot of interface movementshowing the damping process of the free internal seiches 96Figure (2-9)Figure (2-10)Figure (3-1)Figure (4-1)Figure (4-2)Figure (4-3a)Figure (4-3b)Figure (4-3c)Figure (4-3d)Figure (4-4)Figure (4-5)Figure (5-1)Figure (5-2)868687888889899090919293xiiFigure (5-5)Figure (5-6)Free internal seiche in a two-layered fluid with a sharpinterface, run 3; (a) and (b) Time series of density interfacemovements at x=0 and x=L , (c) and (d) time series ofhorizontal velocity at central point of upper and lower layers 97Free internal seiche in a two-layered fluid with a sharpinterface, run 3; horizontal velocity profiles at a time intervalof 0.5Ti showing the large velocity shear at the densityinterface 98Free internal seiche in a two-layered fluid with a sharpinterface, run 3; velocity vector plots over one complete cycle^99Forced internal seiche triggered by constant withdrawal,experiment 1, run 1; (a) piecewise linear density profile, (b)First four vertical internal wave modes associated with thedensity profile of (a) 100Forced internal seiche triggered by constant withdrawal,experiment 1, run 1; horizontal velocity profiles at x=0.5L^101Forced internal seiche triggered by constant withdrawal,experiment 1, run 1; (a) spectrum from the oscillation data ofthe center of interfacial layer in the range of 0.0-0.16 Hz, (b)spectra from the oscillation data of both upper interface andthe center of interfacial layer in the range of 0.0-0.015 Hz 102Forced internal seiche triggered by constant withdrawal,experiment 1, run 2; horizontal velocity profiles at x=0.5L^103Forced internal seiche triggered by constant withdrawal,experiment 2; comparison of observed fundamental firstvertical mode natural periods with the theoretical values, (a)g'bulk = 0.06 m/s2, h1/L=h2/L=0.04-0.25, and (b) g'bullc = 0.06rn/s2, h1/L=h2/L4.04-0. 5 Forced internal seiche triggered by unsteady withdrawal,experiment 1, comparison with physical experiments;schematic of experimental setup (redrawn from Allan, 1993) 105Forced internal seiche triggered by unsteady withdrawal,experiment 1, comparison with physical experiments;interfacial oscillation data at x=0 106Forced internal seiche triggered by unsteady withdrawal,experiment 2; time series of interfacial oscillations at x4^108Forced internal seiche triggered by unsteady withdrawal,experiment 2, run 3; horizontal velocity profiles at x4.5L^109Figure (5-7)Figure (5-8)Figure (5-9)Figure (5-10)Figure (5-11)Figure (5-12)Figure (5-13)Figure (5-14)Figure (5-15)Figure (5-16)104Forced internal seiche triggered by unsteady withdrawal,experiment 2, run 3; (a) time history of interface elevation atx=0, (b) time history of horizontal velocity at central points ofboth upper and lower layers, (c) forcing discharge processForced internal seiche triggered by unsteady withdrawal,experiment 2, run 3; surface plot of density interfacemovement of resonant internal seicheForced internal seiche triggered by unsteady withdrawal,experiment 3; interfacial oscillation data at x3Forced internal seiche triggered by unsteady withdrawal,experiment 3; spectra of interfacial oscillation data at x:1Forced internal seiche triggered by unsteady withdrawal,experiment 3, run 5; horizontal velocity profiles at x4:15LForced internal seiche triggered by unsteady withdrawal,experiment 4, run 1; interfacial oscillation data of (a) thecenter of interfacial layer and (b) upper and lower interfaces atx:), (c) time series of interfacial layer thickness at x=0Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 1; spectra of interfacial oscillation data atx=0Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 2; interfacial oscillation data of (a) thecenter of interfacial layer and (b) upper and lower interfaces atx=0, (c) time series of interfacial layer thickness at x;)Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 2; spectra of interfacial oscillation data atx::1Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 3; interfacial oscillation data of (a) thecenter of interfacial layer and (b) upper and lower interfaces atx=0, (c) time series of interfacial layer thickness at x::1Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 3; spectra of interfacial oscillation data atx=0Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 4; interfacial oscillation data of (a) thecenter of interfacial layer and (b) upper and lower interfaces atx:31, (c) time series of interfacial layer thickness at x:)Figure (5-17)Figure (5-18)Figure (5-19)Figure (5-20)Figure (5-21)Figure (5-22)Figure (5-22d)Figure (5-23)Figure (5-23d)Figure (5-24)Figure (5-24d)Figure (5-25)110111112114116117118119120121122123xivForced internal seiche triggered by unsteady withdrawal,experiment 4, run 4; spectra of interfacial oscillation data at=0 124Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 5; interfacial oscillation data of (a) thecenter of interfacial layer and (b) upper and lower interfaces atx=0, (c) time series of interfacial layer thickness at x:1 125Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 5; spectra of interfacial oscillation data at=0 126Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 5; horizontal velocity profiles at x4■ .5L^127Figure (5-25d)Figure (5-26)Figure (5-26d)Figure (5-27)XVAcknowledgmentsI would like to express my sincere appreciation and gratitude to my thesis supervisor, Dr.Gregory Lawrence, for his advice and guidance throughout the course of the research workand in the preparation of this thesis. Special thanks also goes to Dr. Craig Stevens for hisgenerous help and many valuable suggestions, and to Gary Allan for providing theexperimental data.I would also like to acknowledge Dr. Michael Isaacson for his careful reading of myearly manuscript of this thesis. Finally the financial support from BC Hydro and theMinistry of Environment in the form of a Research Assistantship is gratefullyacknowledged.xviChapter 1IntroductionThe dynamics of a lake or reservoir are largely determined by the combined action of thesun, the wind, inflow, outflow and the density stratification. During summer, a typical freshwater body can be divided into a well-mixed, warm, surface layer called the epilimnion,which is separated from the colder and relatively stagnant bottom water, called thehypolimnion, by an intermediate region with strong temperature gradient known as themetalimnion, which in conditions of temperature stratification only, is sometimes called thethermocline. Figure (1-1) is redrawn after Mortimer (1952) showing the typical temperatureprofile of a lake during summer. The depth of the epilimnion is a measure of the verticalextent to which wind energy and night-time cooling can maintain mixing against the densitygradient set up by solar surface heating. The thermocline between the epilimnion and thehypolimnion acts as a barrier to the mixing and transport of heat, dissolved oxygen, andplant nutrients which sustain life throughout the water body.The vertical transport of nutrients and dissolved oxygen can be explained satisfactorilyby turbulent diffusion resulting from the predominantly horizontal mean flows. Animportant mechanism for generating this turbulent diffusion is provided by the internalstanding waves or internal temperature seiches, occurring not only as isolated events, butmore or less continuously during summer stratification.It has been clearly recognized that the wind is the most important driving force for theinternal seiches. The steady state response of a two-layered lake to a surface wind stress can1Chapter 1 Introductionbe illustrated with the follow simple argument. When a wind starts up, the shear stress at theair-water surface causes a surface tilt by forcing the fluid against the downwind end of thelake. At the same time, to maintain equilibrium of the whole system, an interfacial tiltappears to apply a balance force at the base of epilimnion. The resulting interface tilt is suchthat the thermocline rises at the upwind end of the lake and is depressed downwind. Themagnitude of this interface tilt is much larger than that of the air-water surface due to thereduced density difference. Figure (1-2) is a schematic plot showing the interfacial tilts of atwo-layered lake. When the wind drops, the hydrostatic forces, which formerly helped tobalance the wind drag, are now out of balance and both surface and internal free seiches areset up with most of the energy concentrated in the latter. The slope on the water surfaceresolves into a barotropic surface seiche with a period, to a first approximation for the firstmode or fundamental mode seiche, of2LT =gHwhere g is the acceleration due to gravity, L is the length of the lake along the water surfacewhich is much greater than H, the mean depth of the lake. Similarly, the very much largerthermocline slope resolves into a series of interfacial oscillations with decreasing amplitude.These oscillations have periods that are very much longer than that of the surface seiche. Forthe simplified case of a rectangular lake of length L much greater than depth H, and with ahomogeneous epilimnion and hypolimnion of thickness h1 and h2 and density pi and P2respectively, the period of the fundamental mode internal seiche, the baroclinic componentof oscillations, is given byT =2Lil 131h2 M 802— Pdhlik(1-2)2Chapter 1 IntroductionThis is the formula used by Watson (1904) and Schmidt (1908), and is usually referred to asWatson's formula.Another possible mechanism of triggering the internal seiche, even a resonant internalseiche, is unsteady withdrawal. This mechanism was first proposed by Imberger (1980) as apossible pump-back scheme to optimize the water quality in a reservoir or lake. Bywithdrawing fluid from one side of a reservoir or lake, the thermocline can be drawn up inthe vicinity of the discharge outlet and depressed at the far end of the discharge side, asshown in figure (1-2). Thus a forced internal seiche can be generated by the unsteadywithdrawal process, and if the discharge or forcing period is close to the natural period ofthe system a resonant internal seiche will be excited. Unsteady withdrawal has beensuggested as a possible means of enhancing the export of Mysis from Kootenay Lake,British Columbia, Canada, by establishing a near-resonant transverse internal seiche(Lawrence & Allan, 1991). Figure (1-3) is from Carmack & Gray (1982) showing thelocation for the Kootenay river drainage basin. The proposed mechanism for inducing such aresonant internal seiche is to vary the discharge rate sinusoidally at, or near, the naturalfrequency from the Corra-Linn Dam at the downstream end of the West Arm of the lake. Toevaluate the feasibility of triggering resonant internal seiches by unsteady withdrawal, thedynamic response of the thermocline to a transient forcing through the outlet needs to beinvestigated. This can be done, in an abstract way, through the combination of numericaland physical modeling.The simplest approximation of the epilimnion/thermocline/hypolimnion structure ofdensity stratification is the two-layered fluid with a sharp interface which neglects thethickness of the thermocline. This situation best approximates conditions in the late summerand early autumn, when, as shown in the measurements of Mortimer (1952), the thermocline3Chapter 1 Introductioncan be quite thin. However, it is a poor approximation of spring conditions or conditionsoccurring in late fall. At these times, the thermocline region can be quite thick.The resonant internal seiche triggered by unsteady withdrawal from a rectangular tank isinvestigated numerically in this thesis. Parallel physical experiments are carried out anddescribed elsewhere (Allan, 1993). The scope of this thesis is now briefly summarized.Chapter 2 seeks to review some relevant field, laboratory and numerical experiments, andanalytic solutions of both free and forced internal seiches. In Chapter 3 the governingequations describing a density stratified flow in rectangular coordinates are presented.Following the statement of equations written in the primitive variables of velocitycomponents and pressure, the Boussinesq approximation is introduced. Furthermore thegoverning equations based on the stream function and vorticity will be derived to avoid thecomplexity of solving the pressure field. Boundary and initial conditions are also presentedto define the problem. The numerical method is described in Chapter 4. The computerprogram employed in this research is based on that used by Imberger, Thompson & Fandry(1976), originally written to simulate selective withdrawal from a rectangular tank filledwith a linearly stratified fluid. It has been modified so that the non-slip boundary conditionscan be imposed. The numerical scheme involved in the computer program is so-called C/C(Cloud-In-Cell) method used by Birdsall & Fuss (1969) which actually is the combination ofEulerian and Lagrangian approaches. To increase the efficiency of the program, Gauss-Seidel iteration, SOR (Successive Over-Relaxation) and the ADE (Alternating DirectionExplicit) method are employed to help obtain quick convergence. Some general discussionson numerical accuracy are also be presented. Chapter 5 discusses the results obtained byimplementing the computer program. Conclusions and recommendations as to furtherresearch are presented in chapter 6.4Chapter 2Literature Review2.1 IntroductionRhythmic oscillations in the surface levels of lakes have long been known as standing wavesor seiches. Amplitudes of more than a few centimeters are only found on long lakes, andthen infrequently (Mortimer, 1952). In themselves, surface seiches play a minor part in thephysical and biological dynamics of a lake, but interest in them led directly to the study of amore important class of oscillatory motion, the internal seiche. When investigations of thevertical distribution of temperature became common, oscillations of much larger amplitudein the level of thermocline were occasionally noticed.The first correct interpretation of these oscillations was given by Watson (1904). Heobserved oscillations in the level of the thermocline in Loch Ness, Scotland, with amplitudesof the order of 30 meters and a period of 3 days. He interpreted them as standing waves onthe interface of two fluids of differing densities, and found fair agreement between theobserved period and that calculated from equation (1-2). Schmidt (1908) carried outexperiments on standing waves in a two-layered fluid with a sharp interface. His apparatuswas a rectangular tank which, when filled, was gently rocked and then held stationary whilethe period of the resulting interfacial wave was measured. More detailed field observationswere made at Loch Earn, Scotland, by Wedderburn (1912) and Wedderburm & Young(1915). Figure (2-1a), redrawn from Wedderburn (1912), shows the general behaviour of aninternal seiche. The great achievement of Wedderbum and his collaborators was the5Chapter 2 Literature Reviewillustration, in considerable detail, of the course of the isotherms during these oscillations,and the way in which they were affected by wind. Less success attended their efforts tomeasure the currents. However, the main features of flow due to the seiche alone werepredicted for an idealized lake by Wedderburn & Williams (1911), and are indicated infigure (2-1b). Except at the ends of the basin, where vertical components of flow becomeimportant, the seiche is characterized by a horizontal drift, which changes sign at each half-cycle and is opposite in phase in the epilimnion and hypolimnion.Some relevant field, laboratory and numerical experiments, and analytic solutions willbe reviewed in the following two sections. In the section on free internal seiches, the reviewwill mainly concentrate on the influence of finite amplitude, viscous damping and finiteinterface thickness on the natural period. The section on forced internal seiches deals withresponse of a density interface to external forcing mechanisms.2.2 Free Internal SeichesAny system in nature that can invoke a built-in restoring force for reestablishing itsequilibrium position, once it has been displaced from it, will perform free oscillationsprovided the disturbing forces responsible for the initial displacement are not sustained. Theresulting oscillations are characteristic of the system only and are independent of theexciting force, except, of course, as to initial magnitude. They must eventually die out underthe effect of friction as the system returns to rest. A free internal seiche is such a freeinterface oscillation of a water body with different densities. The restoring force is providedby the action of gravity tending always to return the interface to its horizontal equilibriumposition.6Chapter 2 Literature Review2.2.1 Linear Inviscid Internal SeichesConsider the internal standing waves, or internal seiches, in an incompressible and inviscidfluid confined within a rectangular tank. The fluid consists of two homogeneous layers ofslightly different densities separated by a sharp interface. The equations of the motion arelinearized by assuming the amplitude of the motion is small. These linearized equations arethen applied to each layer separately. Within each layer the motion is irrotational and sopotential theory can be used to describe the flows above and below the density interface (seeLamb, 1932). The interface profile is found to be sinusoidal, i.e.(rh).= a sin(at)cos(kx)^ (2-1)where ni is the displacement of the internal interface,n is the horizontal mode number of the internal seiche,a is the amplitude of internal seiche,k = 2ir /A, is the wave number,A = 2L/n, is the wave length,L is the length of the tank,27r I 7', is the angular frequency, andT is the period of internal seiche.Figure (2-2) shows density interface profiles for the first three horizontal modes ofinternal seiches in a two-layered fluid with a sharp interface, where hi, p1 and h2, p2 are theupper and lower fluid depths and densities respectively. It is noted that at half wave lengthintervals (x = )/4, 3A/4, 5A/4, ...), density interface displacement is continuously zero withtime. These points are nodes and the intermediate points of rise and fall of density interfaceare antinodes of the internal seiche.7Chapter 2 Literature ReviewThe general dispersion relationship is found to be(72 . gk(p2— pi)tanhkhi tanh kh2 tanh kh2 + p2 tanh khiThe period of internal seiche T is1/2T=21r/a= 2J tanh kh2 + p2 tanh ( gk(p2— tanh khi tanh kh2Various special cases can be recovered using the limits tanhkhi —> 1 as khi^, andtanh khi^khi as khi—> 0.When both layers are deep, tanhkhi —> 1, this leads to(2-2)(2-3)T=24 Pl P2 gk(p2—)1/2(2-4)When both layers are shallow, tanhkhi —> khi, this leads toT 2L p1h2 + p2hin g(P2 VIAwhich is the result arrived at by Watson (1904). In the natural environment the densitydifference (p2-pi) is much smaller than either pi or p2 so that equation (2-3) yields a periodof oscillation much greater than that of the surface seiche. Along with the period, theamplitude of the internal seiche is also greatly increased over that of the free surface.The above estimates of the natural period of internal seiches are subject to the affects offinite amplitude, viscous damping, and finite interface thickness. Furthermore, a betterestimation of the natural periods to the conditions in natural lakes can be made by applyingmore complex theory which takes the shape of the basin and the rotation of the earth into)1/2(2-5)8Chapter 2 Literature Reviewaccount (Mortimer, 1979). Hamblin (1978) pointed out that inclusion of the Coriolis forceand real bathymetry can lengthen the above estimates by about 10%.2.2.2 Finite Amplitude Internal SeichesBy using an expansion technique first used by Stokes for surface waves, Hunt (1961)obtained the theoretical profile and dispersion relation for finite amplitude standinginterfacial waves in deep fluids. Thorpe (1968a) made an extensive study, both theoreticaland experimental, of the internal seiche. In his theoretical study, Hunt's analysis wasextended to include finite fluid depths. It was found that, to third order, the dispersionrelation iscf2 . gk(p2— pi)tanhkhitanhkh2 [1+ ^a2k2Sh p1tanhkh2 + p2 tanh khi^32 tanh2 khitanh2 kh21where Sh is defined asSh = (9 tanh2 kh, + 9 tanh2kh2 —18 tanh khi tanh kh2 —6 tanh3khi tanh kh2—6 tanh3kh2 tanh kh, + 8 tanh2 khltanh2kh2)Apparently if Sh> 0 then the frequency increases with increasing wave amplitude, whilst ifSh < 0 then the frequency decreases with increasing wave amplitude.If both khi and kh2 are large then equation (2-6) reduces toa2 . gk(p2— pOr 1 a2k2 IP1 +P2 L^8 (2-8)Thus when both upper and lower fluids are deep, the frequency decreases with increasingwave amplitude.(2-6)(2-7)9Chapter 2 Literature ReviewIf khi is small and khi is large (i j), then equation (2-6) reduces to0.2 gk(p2 — pi) tanh khi[i+ ^a2k2 + pi tanh khi^32 tanh2 khi (9 24 tanh khi + 17 tanh2 khi — 6 tanh3 khi )] (2-9)Thus the frequency increases with increasing wave amplitude if tanhkhi <0.54, i.e. if thedepth to wave length ratio is less than about 0.096, whilst the trend of frequency increase isreversed if tanhkhi > 0.54.The general behaviour of the frequency a for various values of tanhkhi, tanhkh2 isindicated in figure (2-3). In region A, the frequency of internal seiches decreases withincrease of wave amplitude; in region B, the frequency increases with increase of waveamplitude.2.2.3 Viscous EffectThe damping effect will start acting as soon as the seiche begins. If such damping issufficiently great it will not merely cause the oscillation to die out but will have anappreciable effect in increasing the period. Harrison (1908) considered the friction at thedensity interface of a two-layered fluid with a sharp interface and showed that the effect ofviscosity at the interface between two deep fluids, if Lip p2, is to decrease the frequencya, from the inviscid theoretical frequency a-0, by an amountAr _kcr. 2(2a0)1/2(2-10)In figure (2-4), the frequency variation Aa/cro is plotted against the eddy viscosity ve forvarious values of k/Vcro. It is obvious that higher mode of oscillations are more severely10Chapter 2 Literature Reviewdamped. The effect of viscosity is also to decay the amplitude rii(t), such thatni(t) = , where i(0) is the initial amplitude.Equation (2-10) is only valid for deep water and is based on the assumption that nomixing occurs and that the interface is sharp. In the case of finite water depth, internal wavesare damped more when the interface is near the surface or near the bottom (Thorpe, 1968a,Ting, 1992). Furthermore a transition layer between upper layer and lower layer may welldecrease this estimation substantially (Thorpe, 1968a, Hyden, 1974).Heaps & Ramsbottom (1966) investigated theoretically the dynamic response of a two-layered lake to wind stress. They assumed that the fluid was initially at rest and thatsubsequent motions were unaffected by the earth's rotation and interfacial mixing. They alsoassumed internal friction at the interface between the layers is zero. They considered bothfrictionless and frictional bottom boundary conditions. The solutions obtained by usingfrictional boundary conditions showed that the internal seiches would be damped by bottomboundary layer dissipation. If the natural period is longer than the decaying time, themotions will be overdamped. This overdamping is most likely to occur in the late fall andspring when small density differences would cause longer periods than in the summer andearly fall.2.2.4 The Effect of Finite Interface ThicknessA diffuse interface of finite thickness enables more than one vertical internal wave mode toexist. The vertical internal modes are represented as vertical profiles of horizontal velocity.An n -layered fluid system supports n-1 internal vertical wave modes. Therefore theintroduction of a continuously stratified interfacial layer between two homogenous layersallows an infinite number of vertical internal wave modes. Furthermore the dispersion11Chapter 2 Literature Reviewrelation for the internal wave changes as the stratification changes. Hence, Thorpe (1968a)and Hyden (1974) observed, in their physical experiments, fundamental natural periods upto 12% larger than those calculated from linear inviscid theory of a two-layered fluid with asharp interface (equation 2-3). This is because the restoring force represented by gndpIdz,where 17 is the vertical displacement of a fluid element from its equilibrium position,decreases as the interfacial layer thickness increases.The current profile of Loch Ness measured by Thorpe (1977) indicated that the strongestcurrents frequently appear in the thermocline itself, not in the epilimnion. Also, currents inthe epilimnion and hypolimnion are often in phase with each other and out of phase withcurrents in the thermocline. This velocity distribution can not be explained by the two-layered theory. Mortimer (1952) applied a three-layered model to field data. He presented ananalysis of internal wave motions in a system of three homogenous fluid layers and foundnot one, but two possible vertical internal wave modes. The shorter-period vertical wavemode is similar to the baroclinic wave in a two-layered fluid with a sharp interface. It has acomparable period. The top and bottom layer fluids flow in the opposite directions, creatingthe largest shear at the middle layer. With this wave mode, the middle layer oscillates as awhole with little deformation. The second vertical wave mode is generally 2 to 3 timesslower than the first vertical wave mode. The top and bottom layer fluids flow in the samedirection with an opposite flow in the middle layer. This leads to the horizontal variations inthe middle layer thickness. See figure (2-5) for a schematic of the vertical internal wavestructures as well as the first and the second vertical modes oscillating mechanism in a three-layered fluid system.The analysis of internal waves uses the linearized vertical momentum equation of anisiviscid and Boussinesq fluid (Roberts, 1975). In the one-dimensional coordinate systemwith y-axis pointing up vertically and the origin at the bottom, which may be written as12Chapter 2 Literature Reviewd2v„ ( N2 — G2 jk202^v. = 0+dy2 (2-11)which is subject to the boundary conditions imposed at the free surface and the bottomv„ = 0 at y = Hvn=0 at y=0(2-12)where vn is the vertical velocity component for nth vertical internal wave mode, and H is thetotal fluid depth. The no-flow boundary condition at the free surface is usually satisfied,since an internal wave produces very small vertical displacements of the free surface if thedensity variation is very small (Phillips, 1966). N is the Brunt-Viiisiihd, or buoyancy,frequencyN^gaP—MY)(2-13)where p = p(y) is the density distribution and Po is a reference density which may be takenas the mean density as the changes in density are very small compared to the absolute value.At a given depth y, the buoyancy frequency N(y) determines the maximum frequency forinternal waves occurring at that depth.For a prescribed distribution of N(y), equations (2-11) and (2-12) pose an eigenvalueproblem in which the possible vertical wave modes and the dispersion relation an(k)associated with each wave mode are to be found. Although general analytic solutions are notpossible with arbitrary N(y), with particular density profiles some analytic solutions areavailable. Phillips (1966) considered the case of two homogenous fluids separated by ainterface of thickness 5 which is much smaller than the wave length. It was assumed that forthe first vertical mode of internal wave motion, the vertical velocity v(y) changes little across13Chapter 2 Literature Reviewthe interfacial layer. equation (2-11) was then integrated over the range 6 and the dispersionrelation for the first vertical wave mode was found to be0_2 _ gk(p2 — pi)(ko+ tanh khi + tanh kh2 IP.^tanh khi tanh kh2(2-14)Apparently, the above dispersion relation is just a modified version of that for a two-layeredfluid with a sharp interface, with k3 being the correction term. With this particular solution,no analytic expression is specified within the interfacial layer and the only dispersionrelation found is for the first vertical internal wave mode. Furthermore the interfacethickness 8 should be much smaller than the wave length to keep equation (2-14) valid.Fjeldstad (1933) considered the case of two homogenous fluids separated by aexponentially stratified interfacial layer (figure 2-6). The density profile is given by,PiAexpr -1-sif(h2+8- y)]L gP2 = p1 exp(N25/g)h2+8<y5_Hh2.yh2+3Oy<h2(2-15) where h2, 8 and h1 are the individual layer thickness, from top to bottom, respectively.Within the interface the buoyancy frequency N(y) is constant, and outside N(y) = 0. Thedispersion relation for the first vertical internal wave mode isiN2_021, tanh kk il + tanh kh2)ltan[( N202— o-21/2k51 +^o-2 = 0( N2 — a21^) tanh khi tanh kh20-2(2-16)The above exponential interface profile can be approximated, to the first approximation, bya linearly stratified interface with p(y)= A(1+ N2 (h2 + 8 — y)/g), because the densityvariation Zip/po is small. For the case of small k8, equation (2-16) reduces to14Chapter 2 Literature Review0.2 _ gk(p2 — p1)(ks 4. tanh IA + tanh kh2 +1 TP.^tanh khi tanh kh2which is different from the dispersion relation of equation (2-14) (Phillips, 1966) for thecase of a small k3 . This is because of the Boussinesq approximation used in the derivationof equation (2-16) (Keller & Munk, 1970).Krauss (1966) considered the case of a hyperbolic tangent density profile in deep water(figure 2-7a). With the y-axis being vertically upward and the origin at the center of theinterface, the density profile is given byp = po exp[— 3tanh(Y-)]^--. y 5_ +00^(2-17)P.^3where 3 is the effective interface thickness, Po = p (0) is the mean density andAp = p(—.) — p(+.) is the total density difference. The dispersion relation is found to beAP gk28P. G2 = k282 + 2(2n —1)k8 + n(n —1)(2-18)where n is the mode number of vertical internal wave. The variations of a2 / (kgAp/pc, )against the relative interface thickness k3 for the first three vertical internal wave modes areshown in figure (2-7b). In the case of 8 = 0, then only one vertical wave mode is possibleand the dispersion relation of equation (2-18) reduces to that of a two-layered fluid with asharp interface (equation 2-4). For given density profile and wave number k the naturalfrequency of the first vertical internal wave mode decreases monotonically with increasinginterface thickness. Also noticed is that the natural frequency gets progressively smaller forhigher vertical internal wave modes.15d2X ^P--X + c2X = F(ty49t2^. at^. m (2-19)Chapter 2 Literature ReviewFor an arbitrary density distribution N(y), equation (2-11) may be solved numerically.This was first done by Fieldstad (1933). Krauss (1966) discussed the Runge-Kuua numericalintegration method applied to this eigenvalue problem. In his method, after assuming aninitial value vn = 0 at y = 0, the dispersion relation an(k) is obtained, by trial and error, sothat the boundary condition vn =0 at y = H holds.While a seiche is usually a free oscillation of a fluid, it is possible to envision a forcedseiche, activated by a source of excitation such as unsteady withdrawal which continues tofeed energy to the system.2.3 Forced Internal Seiches2.3.1 Mass-Spring-Dashpot System AnalogyIt is instructive to examine the purely mechanical analogy of a vibrating system beforeinvestigating the resonant internal seiche. The equation of motion for the displacement X ofa linear vibrating mass-spring-dashpot system subject a disturbing force F(t), which is afunction of the time t, may be written aswhere mo:2 is the spring constant for restoring force,/3 is the non-dimensional damping coefficient, andao is an angular frequency.16Chapter 2 Literature ReviewThe solution to this equation comprises two parts, the transient and steady state, or thefree and forced oscillations. The free oscillation X° is obtained by solving equation (2-19)for the case of F(t) = 0, from which= e-13"[asin(at)+ bcos(at)]^(2-20)where a and b are amplitudes of the motion determined by initial conditions, and the naturalfrequency a = a0(1—/32)1/2. The friction damping in the system, prescribed by the value ofJ3, tends to cause the free oscillation to decay exponentially.The forced oscillation Xf must satisfy equation (2-19) in full. In the case of a periodicexcitation of frequency co and the formF(tY = F cos(am + 7)With arbitrary phase angle y, the forced oscillation is found to be= (FP/ 2)cos(cot + y — a)in which= 1[1 - / )2]2 + [21340 / 0A2r1/2(2-21)(2-22)(2-23)tanh a = 213(col a.)[1 — (a)I o-.)2]^(2-24)The quantity u is the dynamic amplification factor of the oscillation and a is the phaseangle by which the forced oscillation lags the excitation.The complete oscillation then isX=X.+Xf,^ (2-25)17Chapter 2 Literature Reviewof which X0 decays and Xf persists as long as the excitation is applied.The dynamic amplification factor p , from equation (2-23), obviously has the meaningp — Xmax — X/X.(F/a!)(2-26)where Xi is the amplitude of the input displacement. The well-known forms of functionsp(wArro) and a(akao) in terms of variable frequency ratio 0o/a0 are shown in figure (2-8).Certain properties of figure (2-8) are worthy of noting: (a) The natural frequency of thesystem is a= ao if the damping 2,3 « 1. Therefore the ratio acro is effectively the ratio ofthe forcing frequency to the natural frequency. The dynamic amplitude p in figure (2-7) isseen to approach its peak when acro 1. (b) When resonance occurs, the magnitude of themotion may be several times that of the excitation amplitude. (c) When the frequency issmall, co/o-0 « 1, there is little or no magnification, p 1, and the phase relationship tendsto ensure that the resulting motion follows the excitation, a —> 0. On the other hand, if thefrequency ratio is large, oVao» 1, the induced motion becomes only a small fraction of theinput motion, p —> 0, and the motions tend to become out of phase, a —> 1800. (d) Thedegree of damping 2fl has an important effect on the severity of resonance. The smaller thedamping, the sharper the resonance. Therefore for small damping, a small difference in comay resulting in large difference in the amplitude of resonance.The purely mechanical oscillator considered above differs in detail from the nonlinearoscillations of internal seiches which may invoke many modes of oscillations, however ascan be seen in the following sections the analogy is helpful to illustrate some importantfeatures of a forced internal seiche.18Chapter 2 Literature Review2.3.2 Forced Internal SeichesObservations of temperature structure, such as those in Loch Ness reported by Thorpe(1971), frequently indicate the asymmetrical shape of internal seiches. This asymmetry ofthe internal seiche can be interpreted as being due to a non-linear internal surge or bore,which is the form taken by large amplitude internal seiches (Thorpe, 1968a). This internalsurge is sometimes observed to retain its magnitude over several periods, suggesting thepossibility of resonance with the wind. The internal surge in Loch Ness was later modeledby Thorpe (1974) as a near-resonant internal seiche in a shallow two-layered fluid with asharp interface. The response of a two-layered fluid with a sharp interface to various windpatterns, and the characteristics of wind-driven circulation in the upper layer were examinedtheoretically by Heaps & Ramsbottom (1966). Their analysis indicates a large scale wind-driven circulation in the upper layer with a downwind surface flow and a upwind return flowin the lower part of the upper layer.The generation of internal seiche in a submarine rectangular trench by normally incidentsurface waves was investigated by Ting (1992). The density profile used consisted of twohomogeneous fluids of slightly different densities separated by a continuously stratifiedtransition layer. It was observed that, when the frequency of the surface wave corresponds tothe natural frequency of internal waves, the amplitude of internal seiche becomes largecompared to the amplitude of surface waves. The natural frequency of the first verticalinternal wave mode decreases as the thickness of the density interface increases and depth ofthe bottom layer decreases. It was also found that the internal wave was strongly dampedwhen the depth of the bottom fluid was small compared to the wavelength of the internalwaves.19Chapter 2 Literature ReviewThorpe (1968a) carried out physical experiments on forced internal seiche in arectangular tank to test some of the theoretical conclusions made in the same paper (seeSection 2.2.2). The apparatus consisted of a rectangular tank and two plungers which werefitted at the side walls. The forced internal seiche was then generated by moving theplungers in and out of the tank in single harmonic motion. The density profiles used in hisexperiments initially consisted of two homogeneous fluids with a sharp interface, howeverdue to the mixing of wave motion and molecular diffusion this interface tended to thicken.Figures (2-9) and (2-10) are response curves obtained showing the amplitude of the forcedinternal seiche subject to different forcing frequencies. Figure (2-9) is the response curve indeep fluid, indicating that in deep water the frequency does decrease with increasing waveamplitude, as predicted by the theory (equation 2-9). However the theoretical increase infrequency with increasing amplitude, when the lower layer is sufficiently shallow, is notobserved at kh2=0.48, as indicated in figure (2-10). This is due to the extra damping effectfrom the bottom when the lower fluid is shallow. The arrows on the a axis of figures (2-9)and (2-10) indicate the measured fundamental natural frequencies of small amplitudeinternal seiches, which are about 2.6% and 5.8% smaller than those calculated from thelinear inviscid theory of a two-layered fluid with a sharp interface (equation 2-2)respectively.Hyden (1974) studied the forced internal seiche of a two-layered fluid with a sharp bothnumerically and experimentally. The experiments were carried out in a flume. At one endthe flume ended with a fixed wall, while at the opposite end flows were induced in differentdirections in two layers with a piston. When the forcing period matched that of thefundamental internal seiches, resonance ensued. It was observed that, when at resonance, themovement of the density interface follows the forcing process. This is different from thelinear mass-spring-dashpot system. Also observed was a difference as large as 12% between20Chapter 2 Literature Reviewthe measured fundamental natural frequencies and those based on the linear inviscid theoryof a two-layered fluid with a sharp interface. This large discrepancy was mainly due to thethickening of the interfacial layer. The numerical solution was based on the method ofcharacteristics and the boundary conditions were modeled by a sinusoidal discharge processdetermined by the flows induced by the piston. Hyden's resonant internal seiche experimentwas later compared with the numerical model by Hodgins (1979), who used an implicitfinite difference scheme, the Crank-Nicolson scheme, to solve the equations of motion of atwo-layered fluid with a sharp interface in one dimension. Both numerical models used byHyden (1974) and Hodgins (1979) were formulated in terms of layer transport and thicknessin order to optimize the efficiency of the computer program. However, these two numericalmodels have some drawbacks. Among them the most serious one is the breakdown when thedensity profile departs from the two-layered fluid with a sharp interface due to thethickening of the interface, and the continuously stratified density profile is impossible tohandle.21Chapter 3Governing Equations3.1 IntroductionThe set of equations which describe the physical phenomenon of fluid motions are non-linear partial differential equations. The essential nonlinearity of these partial differentialequations renders them insolvable analytically except in special cases. One approach to suchproblems is to model the phenomenon on a laboratory scale. This approach is, however,often limited by the impossibility of properly reproducing simultaneously all the parameterswhich are important. It has the advantage of sometimes revealing physical processes whichhave been overlooked in the mathematical formulation, or may, once the importantparameters have been discovered, give some useful extension of the solution beyond therange for which solutions are impossible by other means.An alternative approach to such problems is to find numerical solutions using variousnumerical methods such as the finite-difference method. This approach has the advantage ofhaving the complete control over fluid properties such as density, viscosity, etc. It hasenormous flexibility in the choice of flow parameters and information can be easilycollected without disturbing the flow. But in no sense can numerical experiments everreplace physical experiments and theoretical analysis. The ability of numerical models tosimulate and predict physical processes is largely limited by the ability to represent theseprocesses in terms of equations and coefficients. Finally, numerical models are as limited asphysical experiments in that they only gives discrete data for a particular parametric22Chapter 3 Governing Equationscombination. They do not provide any functional relationship, beyond those obtained bydimensional analysis. It is therefore no substitute for even the simplest theoretical analysis.The present investigation is concerned with resonant internal seiches in anincompressible fluid triggered by unsteady withdrawal from a rectangular tank. Thedischarge occurs through a line sink at one side wall of the tank. The density profiles usedinclude the two-layered fluid with a sharp interface and the two-layered fluid with a diffuseinterface. The effect of the rotation of the earth is ignored and there is no motion of water intransverse direction, i.e., motion is two-dimensional in the longitudinal section of the tank.The fluid is viscous and diffusive. The viscous effects, including those from boundaries andinternal interface, may be represented by an eddy viscosity ve which describes themomentum flux transfer of the fluid. This eddy viscosity is assumed constant and uniformthroughout the fluid. While a typical eddy viscosity in real lakes is ye =10-4 m2/s (Thorpe,1971, Heaps & Ramsbottom, 1966), whereas, in laboratory experiments a value close to themolecular viscosity, v = 10-6 m2/s is appropriate (Thorpe, 1968a, Allan, 1993).In the following sections, the governing equations for two-dimensional incompressibleflow in rectangular coordinates will be presented. Following the statement of governingequations written in the primitive variables of velocity components and pressure terms, andthe introduction of the Boussinesq approximation, the governing equations based on thestream function iy and voracity ç will be derived. Furthermore boundary and initialconditions are presented to define the problem.3.2 Governing Equations in terms of Primitive VariablesThe governing equations of motion relating to the flow of a viscous incompressible fluid arethe Navier-Stokes equations. In a two-dimensional rectangular coordinate system, with the x-23Chapter 3 Governing Equationsaxis horizontal, the y-axis vertical and the origin at the bottom left-hand corner of the region,these may be writtenDu = 1 aPT v v2uDt pT dxP--v = — dPT — g + veV2vDt pT dywhere —D = + u + v denotes differentiation followin^Dt^g the motion of a particle,dr dx^2^192^(92V2 ax2 19y2u,v are velocity components in the x and y direction respectively,g is the acceleration due to gravity, andve is the eddy viscosity.Pr is the total pressure which can be decomposed such thatPT(Xly,r) = P(y)+ P(x,Y,t)^ (3-3)where P(y) is the hydrostatic pressure and P'(x,y,t) is the perturbation in pressure due to thefluid motion.Similarly, pr is the total density which can be decomposed such thatPT(x,y,r)= P(Y)+ P'(x,y,t)^ (3-4)where p(y) is the density profile when the flow is in hydrostatic equilibrium, and p'(x,y,t) isthe density variations induced by the motion which is much smaller than p(y) due to the factthat the density difference is very small.Substituting equation (3-4) into equation (3-1), (3-2) and rewriting them in the form(3-1)(3-2)24(1 -aPChapter 3 Governing Equations(3-5)121D u _^1-B7---p--Lax+^4-(1+q),V2upP( i 4. elDv 1 ir )1' 7. _ ii 4. 12:..)g + (1 + Ljv y 2 v^(3-6)1 Dt^p dy p pWhen the density ratio pilp is small, it produces only a small correction to the inertia andviscous terms compared to a fluid of density p, but it is of primary importance in buoyancyterm. Thus the Boussinesq approximation is introduced which consists essentially ofneglecting variations of density so far as they affect the inertia and viscous terms, butretaining them in the buoyancy term. Hence equations (3-5) and (3-6) becomeDu1 dP____ = __ T + v V2uDt^p dx^•Dv 1 dP D____ = ___ T _,Tg+ 1) eV2VDt p dy pFor completeness, the additional equations are the continuity equation for incompressiblefluid expressing the conservation of massdu av—+—=vdx dy(3-9)and the advection-diffusion equation for total density prDPT = Kv2prDt(3-10)where K is the molecular diffusivity.Although it is possible to obtain numerical solutions for the above equations by using theso-called the method of primitive variables, some complexities are encountered in thecalculation of the pressure field which are not expressed explicitly in the continuity equation(3-7)(3-8)25Chapter 3 Governing Equationsfor incompressible fluids. The central idea of the method of primitive variables is to transferthe indirect pressure information in the continuity equation, either into the direct pressurecalculating scheme, or into the correction of the velocity field. Many numerical schemeswith the pressure and velocity components being the primary variables have been developed.Among them the most famous one is the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) method of Patankar & Spalding (1971). In this method a firstapproximation to the velocity field is obtained based on an initial guess for the pressurefield, and then a series of successive corrections to the pressure field are made so as to bringthe velocity field into conformity with the continuity equation. An alternative method toavoid the pressure-solving complexities is based on stream function and vorticity equationsand will be referred to as yi—g method hereafter (Dix, 1963).3.3 Governing Equations in terms of Stream Function and VorticityThe pressure terms are eliminated from equation (3-7) and equation (3-8) by differentiatingequation (3-7) with respect to y and subtracting equation (3-8) after it has been differentiatedwith respect to x. Defining the vorticity g asdu dvC.= -This leads to the parabolic vorticity transport equationB.; g aPT + %VCDt p dx(3-11)(3-12)The presence of equation (3-9) suggests that a stream function vi exists such thatu= ay,^dyrand v = (3-13)26Chapter 3 Governing EquationsIn this way equation (3-9) will automatically be satisfied. Substituting equation (3-13) intoequation (3-11), equation (3-11) is recast as an elliptic Poisson equation(3-14)Note that no time derivative is specified in equation (3-14), the implication being that allvariables are prescribed at the same time. This property of the Poisson equation is directlyrelated to incompressibility of the fluid. That is, if a disturbance exists in one part of thefluid it affects all other parts of the fluid.Now that the governing equations are the Poisson equation (equation, 3-14), thetransport equation for vorticity g (equation, 3-12) as well as the transport equation for totaldensity pi, (equation, 3-10). If the transient pressure solution is required, take the partialderivative of equation (3-7) with respect to x, the derivative of equation (3-8) with respect toy, and then sum the resulting equations. Together with equation (3-9), the continuityequation, this leads towhere Sp is defined(32(PT/P) + d2(PT/P) = sax2^ay2^Pd'u^d2(uv) a2v .^a(PT /P)Sp = ():2 + 2 axdy + ay2) g ay (3-15)(3-16)If the transient pressure solution is not required, the ty—c method is preferable to themethod of primitive variables because it contains one less variable and avoids theunnecessary complexity of solving the pressure field. The v-g method has some drawbacks,however, in that it is impossible to define the stream function iv in most three-dimensionalflows, and the free surface boundary conditions are more difficult to apply. The v—g method27Chapter 3 Governing Equationsis used here because simplified free surface boundary conditions may be chosen withouteffecting the main features of the flow.3.4 Boundary and Initial Conditions3.4.1 Solid Boundary ConditionsThe solid boundaries including the bottom and side walls are impermeable, rigid and non-slip. Since the solid boundaries are stream lines, the stream function y/ may be specified byconstant values. The solid boundary conditions for stream function yi are illustrated in figure(3-1), i.e.on 0<x<L,y=0andx=L,0<y<Hon x=0,0<y5hdon x=0,hd<y5..H(3-17)where Q is the discharge rate per unit width of the tank, and hd is the elevation of thedischarge outletThe solid boundary vorticity is an extremely important evaluation. The vorticitytransport equation determines how the vorticity is advected and diffused, but the totalvorticity is conserved at interior points. At a non-slip boundary, however, the vorticity isproduced. The vorticity is obtained from the non-slip conditions, using the bottom boundaryas an example, y/(i, 2) is expanded by a Taylor series out from the bottom values y/(i, 1)a2,il.yr(i,2)= tif(i,1)+^IiiAy +^0(Ay3)dy '^2 dy2 ' (3-18)28Chapter 3 Governing EquationsBut —av I. = u(i,l) =0 by non-slip conditionsayavNow g = —au — —av , and Along the bottom,^i= 0aY ax^ax 'Therefore, -a2vj dui,1 =^= 091).(..- v(i,1)= const = 0)Substituting this into equation (3-18) and solving for g(i, 1) givesg(i,l) = 2( v(i, 2) — v(i,1)) / Ay2 + 0(Ay)Regardless of the wall orientation, this can be written as= 2( vn, —^/ An2 + 0(An),where A n is the distance from (w+1) to (w) normal to the solid boundary.(3-19)(3-20)It may also be possible to model a slip boundary condition with the viscous governingequations. This implies that the boundary layer is less than An thick. The v values can bespecified, and the vorticity is evaluated from a Neumann condition Cw = gw-f-i. However thevalidity of this approach is unproven. It cannot provide a mathematically consistent systemof equations, since the boundary layer thickness cannot remain less than An as An—)0(Roache, 1972).3.4.2 Free Surface Boundary ConditionsFor numerical schemes formulated in terms of stream function v and vorticity g, the freesurface stress condition is very hard to apply (Harlow & Welch, 1965). Hence the freesurface itself is considered as a static and slippery open boundary, on which, however, a29Chapter 3 Governing Equationsuniformly distributed vertical inflow is applied. This small inflow v ly.H = QIL is introducedhere to compensate for the discharge from the side wall. The free surface boundaryconditions written in terms of stream function viand voracity g arevfly,__H = Qx I LA^du av ft=H=---=vY^ay axdu(... —, = 0 on slippery boundary and, v = QIL= constant)0Y(3-21)(3-22)3.4.3 Initial ConditionsInitially the stream function yr and voracity g at interior points are assumed to be equal tozero, i.e.vci,D1,.0 = 0^ (3-23)s(i,i)11=0 = 0 (3-24)Obviously, the initial conditions will affect the initial transient solutions, but this effectoften decays exponentially in time (Roache, 1972). Furthermore, the initial conditionsusually make little difference on the computer time required to obtain a completelyconverged solution. The initial conditions make little difference because the error of theinitial guess of 'if and g is usually bounded and many orders of magnitude larger than theconvergence criteria c. For example, if the withdrawal rate is normalized to vfmax=0(1), thena guess of 0.0= 0 at all the internal points will only give errors of 0(1). Any improvementsin the initial guess will be insignificant when convergence criteria like e = 10-4 are beingused.30Chapter 4Numerical Method4.1 IntroductionThere are two approaches to the numerical representation of fluid motion. In one, called theLagrangian approach, the fluid is considered to be subdivided into a large number of finitezones each of which characterizes an element of the fluid. This mesh of cells follows thefluid motion in a manner which is determined by an approximation to the governingequations. By following the motion of fluid elements, the artificial viscosity and diffusivityerrors are largely reduced and fine details can be resolved more easily and accurately. At thesame time, however, the approach suffers from certain drawbacks. In particular, the mostserious is the breakdown when the fluid becomes strongly distorted or large slippage occurs.The other is the Eulerian approach. In this, the fluid is again subdivided into numerousdifference cells, but in this case the mesh of cells is fixed in a frame of reference at rest withthe observer and the fluid flows through the cells. The tremendous advantage of the Eulerianapproach is that the calculation proceeds without difficulty when there are large distortionsor slippage in the fluid. The Eulerian method has several disadvantages, however, in that itis very difficult to resolve small features of the flow moving with the fluid through a coarseEulerian grid mesh system, and interfaces are not easily handled.On the other hand, the PIG (Particle-in-Cell) method originally devised by Evans &Harlow (1957) seems to combine the best features of both Eulerian and Lagrangianapproaches and is most advantageously applied to interface problems, because the discrete31Chapter 4 Numerical Methodparticles may be assigned different densities and vorticities to represent different fluids. Inthe PIC method the fluid is considered to be subdivided into a large number of Eulerian cellswhich remain at rest. There is, in addition, however, the Lagrangian mesh of particlesrepresenting elements of the fluid which move through the Eulerian cell system. TheEulerian mesh is used for the purpose of characterizing the field variables, while theLagrangian particles are used to characterize the fluid itself.The Eulerian grid consists of a set of rectangular cells which cover the computationaldomain. Each cell is Az long and Ay high. The resolution of the grid is determined byNX*NY, where NX is the number of divisions in the horizontal direction and NY is thenumber of divisions in the vertical direction. In each of these computational cells, physicalquantities of the fluid such as velocity, vorticity, density and stream function are calculatedat different locations. This treatment is commonly referred to as the staggered grid whichwas first introduced by Harlow & Welch (1965). Figure (4-1) shows how the Eulerian fieldvariables u(i, j), v(i, D, 'xi, D, go, D and p-r(i, j) are arrayed in the xy plane. The relativepositions of v/(i,D and u(i, j) were chosen so as to stagger u(i, j) between w(i, j) and v(iji-1)since u = vry. The positions of v(i, j) and g(i, j) are analogous. This implementation is alsoconvenient for imposing the boundary conditions.The set of Lagrangian particles are initially located uniformly throughout the numericaldomain. When the average number of particles per cell is small, say of the order of four orfive, then the fluctuations caused by the motion of particles across cell boundaries are verylarge, and it should be expected that such fluctuations would cause large effects on theinstantaneous values of local field quantities. In order to simulate correctly internal gravitywaves which have the property that small displacements cause small, continuous changes,particles will not be thought of as being point-like, but rather as small volumes and may passthrough one another. Therefore this method is similar to the CIC (Cloud-in-Cell) method32Chapter 4 Numerical Methodused by Birdsall & Fuss (1969). In this method, each particle is considered to occupy aspace dx by Ay and initially the particles are positioned uniformly through the numericaldomain at a spacing Axa and Ay/2. As the flow is initiated, the particles will move, forsimplicity the motion will be purely translational. Since there is no rotation, the position ofthe K th particle can therefore be defined by its center point X(k),Y(k). Unlike the MAC(Marker and Cell) method (Harlow & Welch, 1965) in which marker particles are only forthe purpose of indicating fluid configuration, in the CIC method, associated with eachparticle, there are corresponding vorticity g(k) and density p7(k) which are assumed to beevenly distributed throughout each particular particle. The Lagrangian particles are movedaccording to the velocity components obtained using bilinear interpolation from the valuesat nearby Eulerian grid points.4.2 Computational ProceduresOnce the computational domain has been defined and the staggered grid has been formed,the governing equations are discretized and solved using the CIC method. The detailedderivation of the finite difference equations is based upon the following sequence of eventsby which the calculation is marched from one time step to the next:(a) The stream function field ty(i, j) is calculated as an Eulerian field variable at thebeginning of each cycle. This requires the solution of the Poisson equation (equation, 3-14),which may be accomplished by a relaxation technique.(b) The Lagrangian particles are moved according to the velocity components in theirvicinities.33Chapter 4 Numerical Method(c) The Lagrangian particle properties (g(k) and pi(k)) are updated. The two materialderivatives are calculated; the products of these with the time increment per cycle then givesthe changes to be added to the old values.(d) The Eulerian field quantities g(i, j) and p7(i, j) are updated.This, then, completes the advancement of one entire calculation cycle. Actually, severalcrucial points have been omitted in this brief summary; these are more easily discussed inthe detailed description to follow. The computational details will be outlined according tothe order of the calculation.4.2.1 The Eulerian CalculationAt the beginning of each time step, the complete field of vorticity^j) is known either as aresult of the previous time step of calculation or from the prescribed initial conditions, thestream function ip(i, j) is then calculated by solving the Poisson equation.As the first step the Poisson equation will be expressed in basic fmite-difference formsusing the Taylor series expansion. Consider the Taylor series expansion of w(i+1, j) and'Ai –1, j) about the point (i, j):dv/(i +1, j) vfo^ig^1 .92 ivax I ij•Ax2 + 1 Cr^1 42 dx2 44^6 ax3 I idax3 + —19vi—I . .Ax4 + HOT24 ax4 4.'(4-1)^= 1/10j)– 11-1^+ 1 (92/P I AXi,j2 1 -(92111 AX + (341If I idAr4 + HOT6 thC3 .j^24 ax4(4-2)34Chapter 4 Numerical Methodwhere HOT is an abbreviation for high-order terms. A centered finite-differenceapproximation of tv,z, can be derived by adding the equation (4-1) and (4-2)1/1(i ± 1,j) + Ilf(i -1,j) = 21g(i,j)+—a21111Ax2 +-1---.411f I. i Ax4 + HOTax2 i.^12 ax4 'Solving for yi.„, givesaz vi _ vf(i +1,D - 2 Vf(i,j)+ Vf(i -1,j) + c(1\x2)dx2 ii.; - Ax2where 0(4x2) is the order of truncation error.The centered finite-difference approximation of vi,„ is thena2 tg 1 _ yr(i +1, j) - 2 ty(i, j) + yr(i -1, j) ax2 li,j -^Ax2which has a truncation error of order 4x2, i.e., a second order accuracy.similarly, the finite-difference approximation for yfyy can be derived:a2 Vf I= Vf(i, .i + 1) - 2 V/(i,j)+ 1/1(i, j - 1) dy2 ii,J^Ay2which has a truncation error of order 4y2.Therefore the discretized form of the Poisson equation, using second order difference, is the5- point formulaty(i +1, j)- 2 vi(i, j) + v(i -1, j) ± vi(i, j +1) - 2 vt(i, j) + vi(i, j -1) = CV, .i)^(4-7)Ax2^ 42where g(i, j) is known.(4-3)(4-4)(4-5)(4-6)35Chapter 4 Numerical MethodIn a rectangular domain where the maximum i = /VX and the maximum j = NY, equation(4-7) and boundary equations give a system of N=(NX-1)*(NY-1) simultaneous linearalgebraic equations. These equations are generally solved by various iterative methods,because: (a) they are very easy to understand and program while remaining flexible, (b) theyare also faster than most of the direct methods because they make use of the spareness of thematrix formed by equation (4-7), i.e., all terms more than two positions off the diagonal arezero.Now if define the mesh aspect ratio /3 = dr/Ay and denote index m the number ofiterations, equation (4-7) can be rewritten in an iterative form1 ^rIr+1(i,../)= Ilim(i,./)+ 2(1+ /32) [r(j+1,./)+ Vfm(i-1,./)+P2r(i,./ +1)+02 Vim (i, i — 1) — 2(1 + 02)1//m(i,i) — Ax2C(i,./)]^(4-8)which is referred to either as the Jacobi method or as Richardson iteration.If the iterative process is convergent, it is expected thatLiml tifm+1(i, j) — ye' (i, j)I --> 0^ (4-9)M -1**Therefore the bracketed term in equation (4-8) is the difference between successiveiterations and can be thought of as the error term. The tv(ij) can be calculated to any desiredaccuracy by stopping the iterative process once the error terms for all grid points are lessthan some tolerance e, i.e.Max vf—i(i,f)— rci,DI e.^(4-10)Two methods have been adopted in the numerical scheme to speed up convergence ofthe iterative process. Noting that equation (4-8) is a two-level equation, requiring storage of36Chapter 4 Numerical MethodIr+1 and vn. If, at the (m+l)th iteration, the most recently computed values of v are happento be known, these new values should be used in calculating Ir+1. This method is referredto either as the Gauss-Seidel iteration or as the Liebmann method. Asymptotically, m Gauss-Seidel iterations are worth 2m Jacobi iterations, and only require half the core storage. Sincethe stream function field tp(i, j) for all (i, j) in the computational domain is needed to besolved, the order in which the Eulerian field is solved or swept is important. If w(i, j) issolved in an order that makes most use of previously determined values, Forsythe & Wasow(1960) have shown that this can further double the efficiency of the iterative process. Themethod so employed is usually referred to as the alternating direction method. The Eulerianfield is swept in following four successive passes:(a) Sweep the computational domain from the bottom to the surface and from left to 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 method to speed up convergence is so-called SOR (Successive Over-Relaxation) method. As mentioned early, the bracketed term in equation (4-8) is thedifference between iterations and will be denoted as residual Rm(i, j). In an iterative processthe aim is to reduce the residual Rm(i, j) to an acceptable level as quickly as possible. In theSOR method this residual term is multiplied by a relaxation factor 0)R i.e., equation (4-8)becomesvfm+1(i,i ) = vim^+ wRiv"(i,i)^(4-11)For convergence, it is required that 0 < coR < 2. If 0.)= 1 relaxation is absent. If 0 < coR <1the under-relaxation is applied. If 1 < coR <2 the over-relaxation is applied. Under-relaxation37Chapter 4 Numerical Methodis used when the residuals oscillate towards zero, while over-relaxation is used when theresiduals are decreasing monotonically.After the stream function field tv(i, j) is solved by successive over-relaxation of thePoisson equation, the velocity field u(i, j) and v(i, j) can be calculated from the definition ofthe stream function, i.e.du(1., D = —111 l• • and v(i, j)= —13-11aY '''^ax iijUsing centered finite-difference approximations, as discussed earlier, this leads tou(i, j) = lgi,.i + 1) — Vf(i,./) v(i, j) = 1110, D - 10j)D AxAy(4-12)(4-13)4.2.2 The Lagrangian CalculationsUp until now, only the Eulerian field quantities were changed and the fluid has beenassumed to be momentarily at rest. The next step is to move the Lagrangian particles andupdate their properties.The Lagrangian particles move according to the velocity components in their vicinity. Abilinear interpolation is performed to calculate the velocity with which a particle shouldmove. The interpolation weighting depends upon the distance of the particle center from thefour nearest velocity points in the Eulerian mesh of cells. This is equivalent to an areavelocity weighting in proportion to the area of overlap of a particle onto the four adjacentEulerian cells. To illustrate this bilinear interpolation, consider the horizontal velocitycomponent u(k) of the Kth particle as shown in figure (4-2). The area A(i+1, j+1) shaded38Chapter 4 Numerical Method(EZ3) is assigned to grid point (1+1, j+1), the area A(1+1, j) shaded MEI is assigned to gridpoint (1+1, j), the area A(i, j+1) shaded (la) is assigned to grid point (i, j+1) and the areaA(i, j) shaded (E2) is assigned to grid point (i, j). The horizontal velocity u(k) may beobtained fromu(k)= A(i +1, j +1)*u(i +1, j +1) + A(i +1, j)* u(i +1, j)+ A(i, j +1)*u(i, j +1) + A(i, j)*u(i, j) Ax* Ay(4-14)which is of first order accuracy.The movement of the Kth Lagrangian particle is accomplished as follows:(a) Interpolate velocities at Xn(k), Yn (k). Here the index n=tIAt specifying the time stepnumber.(b) Using velocities obtained in (a), calculate where the particle would move in half a timestep 1/2At. Denote this point as Xn+112(k), yn+112M.(c) Interpolate velocities at Xn+112(k), yn+112(k).(d) Using the velocities obtained in (c), calculate where the particle would travel in a fulltime step At starting at Xn(k),Yn (k)X'1+1(k)= X' (k) + Un+112 (k)At^ (4-15)andr+1(k)=T1 (k)+ v".112(k)At.^ (4-16)The Lagrangian particle properties ( g(k) and pr(k)) are updated from equations (3-9) and(3-11) by calculating the total rate of change following the particle motion.D pT . Kv2prDt^ (3-9)39)x+1/2 KV 2pCPr 1 1:74 ( DA.^+1 /2At^Dt(4-17)Chapter 4 Numerical MethodDC = g °Pr + v.V2g.Dt pax(3-11)The finite-difference expressions for above equations, using centered finite-differenceapproximations for time derivatives, are:n+1/2gn+1 — gn .^—(pc)ri+1/2 ( g dprAt+ v.V2g"112LDt )^po dx(4-18)Hence the properties of the Kth particle can be modified by usingpr (k)= p;.(k)+ (DA. I Dt)"+112 At^(4-19)andç1(k)= c."(k)+ (DO Dt)n+112 At^(4-20)where (Dpr/Dt)n+1/2 and (DgOtr+1/2 are calculated as Eulerian field quantities by usingcentered finite-difference approximations and then are interpolated back to the Kth particleat X^yn+1/2(k)n+112(k),^by using the same bilinear interpolation technique as used in thevelocity interpolations.4.2.3 Updating Eulerian Field Variables pi(ij) and goAs the third essential part of the CIC method, the Eulerian field variables pl(i, j) and go, Dare updated from Lagrangian particle properties p7(k) and g(k). To illustrate how p7(ij) andg(i, j) are updated, consider again the km particle as shown in figure (4-2). The influence ofthe Kth particle on grid point (i, j) is assumed proportional to the overlapping area A(ij).40Chapter 4 Numerical MethodTherefore for a large number of particles contributing to point (i, j), the density fieldp7m-1(i,j), for example, may be obtained by summing over the particles as(i, j). AU, DP1(k)A(i, j)where EA(i, j) is the sum of area of particles which overlap the cell centered at (i, j). Alsogn+1(i, j) can be updated by following the same procedures.So far the entire cycle of calculation has been advanced from one time step to the next. Anew cycle of calculation is ready to begin starting from Section 4.2.1, if the time limit is notreached.4.3 AccuracySome common errors which may affect the accuracy are discussed here with somenumerical experiments. While more physically meaningful and quantitative estimations ofaccuracy will be carried out later when investigating the physical properties of internalseiches. All numerical runs carried out in this section have two-layered density profiles, withthe reduced gravity g'= 0.06 m2/s. The relative fluid depths are kh1=kh2= 0.39. The viscosityused is 10-4 m2/s and the time step At = 0.025 seconds. The convergence criterion e for theSOR process is equal to 10-4. The discharge process is sinusoidal with a forcing periodT = 55.0 seconds, which is close to the fundamental natural period of the internal seiche.The average rate of the unsteady discharge Qa„ = 1.3*10-3 m2/s and the amplitude of thedischarge AQ = 1.2*10-3 m2/s.(4-21)41Chapter 4 Numerical Method4.3.1 Truncation ErrorStrictly speaking, truncation error is the error of the difference equation itself rather than theerror from the solving process. Some additional truncation errors will occur if the boundaryconditions are Dirichlet conditions (e.g. the solid boundary condition for vorticity c),because some discretizations are needed. The interpolation process involved in the CICscheme may also introduce some truncation errors. The overall order of truncation error isdetermined by the lowest order of the above errors. After Taylor expansion, the totaltruncation error in the numerical model is evaluated to be of 0(4x) and 0(4y) in space andof 0(At) in time domain.To check on the grid spacing dependence, five computational runs were carried out fordifferent grid size. The time histories of normalized stream function v/Q at two differentstations, point A and B, are shown in figure (4-3b), (4-3c) and (4-3d). See figure (4-3a) forthe locations of point A and B. Since point A and B are located near the antinode and nodepoints of the fundamental internal seiche respectively, it is expected that the largest gradientof y/ in x direction, thus the largest truncation error, occurs near point A and the largestgradient of tir in y direction occurs near point B. The yi/Q variations at point A are shown infigure (4-3b) which serves to illustrate the grid size dependence on Ax. It can be seen that thecurve for 11*41 grid deviates from other curves significantly, while the curve for 21*41 gridand that for 41*41 grid agree with each other quite well. The ty/Q variations at point B areshown in figure (4-3c) and (4-3d) which serve to illustrate the grid size dependence on Ay.In figure (4-3c) there is a rather large discrepancy between the curves for the 21*21 and21*41 grids. Further refining the grid size in y direction, as shown in figure (4-3d), results ina much better agreement. It seems that a 41*41 grid is suitable as far as truncation errors areconcerned.42Chapter 4 Numerical MethodThe same kinds of experiments were carried out for time step dependence test on the41*41 grid (see figure 4-4). It seems that a time step At = 0.025 seconds is a acceptablechoice. Another limitation on dt is the stability criteria which will be discussed later. Furtherreduction on grid size and time step is still possible, however the trade off betweenimproved performance and increased C.P.U. time as well as round-off error doesn't warrantfurther reduction. Unless this reduction is essential for other reasons.4.3.2 Round-off ErrorRound-off error is caused by the limitation of word-length of digital electrical computers.This error may systematically accumulate during the solution of difference equations,especially during iteration. By reducing the grid size and/or time step, descretizingtruncation errors can be effectively reduced while on the other hand this reduction couldmean a substantial increasing of round-off errors. Generally it is preferred not to use a veryfine grid and let truncation error control the round-off error.4.3.3 Artificial Viscosity and Diffusivity ErrorsArtificial viscosity and diffusivity are special kinds of truncation errors exhibited by somefinite difference schemes of advection equations. The first use of the term was by VonNeumann & Richtmyer (1950), who explicitly added a viscosity-like term to the inviscid gasdynamic equations in order to allow the calculation of shock waves. Their explicit artificialviscosity term was deliberately made proportional to Ax2, so as to assume mathematicalconsistency; that is, their explicit artificial viscosity term was indeed a second-ordertruncation error.43Chapter 4 Numerical MethodInstead of adding explicit artificial viscosity terms to the equations, artificial dampingmay be added implicitly, just from the form of the difference equations in the sense of a non-physical coefficient of second space derivatives. Noh & Protter (1963) first presented ananalysis of the implicit artificial viscosity of the upwind differencing scheme applied to thelinear model advection equation gt = -ucz. However by following the movements ofLagrangian particles in the CIC method, the advective terms are removed from the transportequations which reduces artificial viscosity and diffusivity.4.3.4 Aliasing ErrorAliasing errors are those errors which occur because finite-difference method cannot resolvewaves with wavelength shorter than, or equal to, twice the grid length. It is helpful to have alook at the turbulent flow analogy before discussing the aliasing errors of numerical model.Physically the energy of turbulence cascades down from large eddies to small ones bynonlinear interactions. Beyond small eddies, it is dissipated, or degraded into internal energyvia friction. The scale at which the energy is dissipated is determined by the Reynoldsnumber Re. The larger Re, the smaller the scale will be.Numerically aliasing errors are associated with energy exchange between Fouriercomponents (here the energy is in the general sense of the first moment of the transportproperties such as g2). As already mentioned, the shortest wavelength component which canbe discriminated by the computational mesh is 24x. Generally accuracy of the long waves isof interest, so this lack of short-wavelength discrimination might seem to be of littleconsequence. Indeed, if the velocity field is everywhere uniform, no aliasing occurs. But, innonlinear problems, it is well known that the components interact in such a way that energymoves down from the long wavelengths to the short. If there is not enough dissipation to44Chapter 4 Numerical Methodremove energy from the short wavelengths, and if the grid size is not small enough todistinguish these short wavelengths, the energy may reappear in the long wavelengths,distorting those very components of the greatest interest, and even resulting in a kind ofcomputational instability. Three sets of numerical experiments were carried out on a 41*41grid with the same parameters except the viscosity which varies from 10-4 to 10-6 m2/s. Thetime histories of the interface elevation at the discharging wall are shown in figure (4-5).Note that aliasing errors increase with decreasing viscosity.The presence of any damping, physical or numerical, in the computation will alleviatealiasing by providing a dissipative mechanism. Another way of reducing aliasing is toreduce the grid size. Since viscosity cannot adjusted arbitrarily, the grid size is furtherreduced to diminish the aliasing error.4.3.5 Free Surface Boundary TreatmentIn order to implement the yt-g method, the free surface boundary conditions are simplified sothat no free surface wave is allowed. In reality, however, both internal and free surfaceseiches can be set up as the result of the forcing discharge. However, due to the fact that themagnitude of the free surface wave is smaller than that of the internal wave by a factor of theorder Ap/p2 and travels faster than the internal wave by a factor of the order 4i2TAT), forthe case of small density variation, this free surface wave is very small and travels muchfaster than the internal wave. Hence the free surface movements are negligible and a staticfree surface will not affect the main features of the flow as far as the internal waves areconcerned.45Chapter 4 Numerical Method4.4 Stability ConsiderationThe most commonly used method of stability analysis was originated about 1944 by J. VonNeumann. In this method, a finite Fourier series expansion of the solution to a modeladvective-diffusion equation is made, and the decay or amplification of each mode isconsidered separately to determine stability or instability. After the Von Neumann stabilityanalysis, two restrictions are recognized in the case of a linear, constant coefficientadvection-diffusion equation in an infmite domain (see Roache, 1972).The first restriction is Courant conditionC.+Cy 1^ (4-22)where Cx = ltht/AX, Cy = vAt/Ay are Courant numbers in x and y directions respectively.The second restriction is the diffusion restrictiond x + dy _1 1 2^ (4-23)Where dx = veAtiza2, dy = veAtlAy2 are diffusion numbers in x and y directions respectively.Considering the tank of height H = 0.6 m and L = 1.2 m used by Allan (1993) in the parallelphysical experiments, for a grid system of 41*81, velocity components u = v = 0.01 m/s anda viscosity v = 10-6 m2/s, equation (4-31) and (4-32) give a stability criterion of Ati = 0.06seconds and At2 = 0.26 seconds respectively. Both criteria may be considered together byrequiring that1^1At (— + --) 1At2this gives a stability criterion of At 0.18 seconds.46Chapter 4 Numerical MethodDividing equation (4-31) by equation (4-32) yieldsAthy(wly + vAr) < 2 (4-24)ve(dx2 + Ay2) —For the special case of Ax = Ay = 4, this leads toR = (u + v)A 5.. 4^ (4-25)8^V,Where Rg is the grid Reynolds number based on local velocities and a characteristic lengthscale equal to the grid size A. Equation (4-25) presents a stability criterion which isindependent of time step At.It should be mentioned that above restrictions are necessary, but not sufficient, and onlyprovide some guidelines to practical stability, because the non-linearity, or even non-constant coefficient may also cause some numerical instabilities. Therefore it is advisable touse trial and error procedure when actually testing the program.47Chapter 5Results and DiscussionThis chapter discusses the numerical results obtained by implementing the computerprogram. The density profiles used vary from two-layered fluid with a sharp interface totwo-layered fluid with a diffuse interface. The two-layered fluid with a sharp interfaceconsists of two homogeneous fluids of slightly different densities separated by a sharpinterface of density discontinuity. This density profile approximates conditions of real lakesor reservoirs in late summer or early autumn when the thermocline is quite thin, as shown inthe measurements of Mortimer (1952) (figure 1-1). The two-layered fluid with a diffuseinterface consists of two homogenous fluids of slightly different densities separated by acontinuously stratified interfacial region. This density profile approximates conditions ofreal lakes and reservoirs in spring and late autumn, at these times the thermocline can bequite thick.The scope of this chapter is briefly outlined here. Section 5.1 seeks to illustrate someflow features of free internal seiches in a two-layered fluid with a sharp interface. Theobserved fundamental natural periods are compared with the theoretical values as well asthose observed in some physical experiments. Section 5.2 discusses two sets of numericalexperiments of forced internal seiche triggered by constant withdrawal. The first set ofexperiments, experiment 1, has a piecewise linear density profile, which consists of twohomogeneous fluids separated by a linearly stratified interfacial layer. This experimentserves to illustrate the flow features of internal seiche in a continuously stratified fluidsystem. Furthermore the observed natural periods are compared with predictions based on48Chapter 5 Results and Discussionthe eigenvalue analysis. The density profile used in experiment 2 consists of twohomogeneous fluids separated by a exponentially stratified interfacial layer. The observednatural periods are compared with the theoretical solutions of Fjeldstad (1933). The forcedinternal seiche triggered by unsteady withdrawal is studied in section 5.3, which consists offour sets of numerical experiments. Six numerical runs are carried out in experiment 1 tocompare with the parallel physical experiments of Allan (1993). Experiment 2 discusses thedynamic response of a two-layered fluid with a sharp interface to the unsteady withdrawalprocess. The dynamic responses of a piecewise linear fluid system to the unsteadywithdrawal are studied in experiments 3 and 4.5.1 Free Internal Seiche in a Two-layered Fluid with a Sharp InterfaceThis section discusses nine numerical runs of free internal seiche in a two-layered fluid witha sharp interface. The free internal seiche in a two-layered fluid with a sharp interface isrelatively simple in the sense that the external exciting force is absent and only one verticalinternal mode of oscillation can exist. Two different values of viscosity are used: an eddyviscosity ve = 10-4 m2/s representing the typical values in real lakes (Thorpe, 1971), and aviscosity close to the molecular viscosity v = 10-6 m2/s representing the typical values inlaboratory scale (Allan, 1993, Thorpe, 1968a). Detailed experimental parameters and resultsfor each experimental run are listed in table (5-1). The free internal seiche is started bytriggering the interface movement gently using a sinusoidal discharge process with thefundamental natural period calculated from the linear inviscid theory. This discharge processis only retained for one cycle, after that time the effect of initial discharge is negligible,except as to the initial magnitude and the resulting oscillation should be characteristic of thesystem only.49Chapter 5 Results and DiscussionTo check the accuracy of the computer program, it is necessary to check if the computedsolutions exhibit proper physical properties of free internal seiches. It has been shown inChapter 2 that, based on the linear inviscid theory, the fundamental natural period T1 for atwo-layered fluid with a sharp interface can be estimated from equation (2-3). To checkwhether the numerical model can give accurate estimations of this natural period, thefundamental natural periods of the nine numerical runs are measured from time histories ofthe stream function. As an example figure (5-1) shows the time series of the stream functionat three different stations. The fundamental natural periods can be measured from this plotaccording to the points of inflection. Comparison of the fundamental natural periodsbetween the numerical experiments and linear inviscid theory shows a fairly good agreementwith the observed values being always larger than the theoretical ones. The discrepanciesrange from 1.1% to 6.8%, as shown in table (5-1). This increasing of natural period ismainly due to following two factors:(a) The presence of finite thickness density interface. Numerically this interfacethickening is caused by some interfacial mixing and numerical truncation errors. Thus theprogram does not model the stratification exactly as a two-layered with a sharp interface.The average thickness of the interface is estimated, from the time series of the densityprofile, to be in the order of 2Ay, where .Ay is the grid size in the y direction. Hence for a gridsystem with 80 divisions in y direction, this represents a interfacial thickness of about 1/40of the total fluid depth H. The increasing of natural periods due to finite interface thicknesswas also noticed in the physical experiments of Allan (1993), Thorpe (1968a) and Hyden(1974) who observed a discrepancy as large as 12%.(b) The viscous damping. This effect is manifested by the fact that internal waves withlarger viscosity always have longer natural periods than those with smaller viscosity.50Chapter 5 Results and DiscussionIn figure(5-2) the observed fundamental natural periods together with those from theexperimental work of Allan (1993) and Hyden (1974) are plotted against those calculatedfrom the linear inviscid theory of a two-layered fluid with a sharp interface.The effect of viscosity is also to dissipate the wave energy thus decrease the amplitudeof internal seiche. Ting (1992) observed in physical experiments that internal seiches aredamped more when lower layer depth is small compared to the wave length. The decay ofthe free internal seiches of run 3 and run 9 are compared to illustrate whether the numericalmodel exhibits this kind of behaviour. Figure (5-3) shows the time history of densityinterface elevations for run 3 and run 9 at x = L. It can be seen that the free internal seichewith relatively shallow lower layer do decay much more rapidly than that with relativelydeep lower layer. The damping process can even be better demonstrated in the surface plotsof figure (5-4a) and figure (5-4b). These two graphs show the time history of the densityinterface movement throughout the whole length of the tank rather than at a single point. Itshould also be mentioned that the enhancement of damping due to shallow lower layer willalso have an important effect on the severity of the forced internal seiche. On the other hand,the upper layer is generally thinner in lakes or reservoirs, this thin upper layer has the similareffect of enhancing damping to a shallow lower layer, as shown by Thorpe (1968b).Some flow patterns of the free internal seiche, from run 3, are discussed in detail here.Figure (5-5) shows the time history of the horizontal velocity components at the centralpoints of both upper layer and lower layer. Also shown in the figure are the time history ofthe density interface elevations at both x = 0 and x =L. Except during flow establishmentwhen there is a withdrawal existing to set up the initial interface displacement of thesubsequent free internal seiche, the two velocity components are out of phase with eachother and change their signs twice a cycle. A 90° phase angle difference can be foundbetween the interface movements and two velocity components, indicating that the velocity51Chapter 5 Results and Discussionreaches its maximum value when interface is in its equilibrium position, while zero velocityoccurs when interface is at its maximum slope. Figure (5-6) shows, at a time interval of0.5T1, four horizontal velocity profiles at x = 0.5L, where the dotted line indicates theequilibrium position of density interface. With the fluid in the upper and lower layersmoving in the opposite directions, large velocity shear develops at the density interface. Forinternal seiches, the maximum shear occurs at the nodes, i.e. the position of the maximuminterface slope instead of at the crests and troughs as for progressive internal waves. Hence,Thorpe (1968a) observed, for sufficiently large wave amplitude, instability occurring at thenodes of the internal seiche with the form of a vortex which changes direction of rotationevery half cycle. Moreover due to the fact that the non-slip boundary conditions are imposedin the numerical model, there is another evident velocity shear near the bottom indicating theexistence of the bottom boundary layer. The velocity distribution is best illustrated in thevelocity vector plots of figure (5-7) which closely matches the general flow patterns asdepicted in figure (2-1b).5.2 Forced Internal Seiche Triggered by Constant WithdrawalThis section discusses two sets of numerical experiments of forced internal seiche subject toa forcing withdrawal with a constant discharge rate. This experiment serves to illustratesome basic flow features of forced internal seiches in a two-layered fluid with a diffuseinterface. The forced internal seiche in a two-layered fluid with a diffuse interface differsfrom the free internal seiche in a two-layered fluid with a sharp interface in two respects: (a)The introduction of a continuously stratified interfacial layer between two homogenouslayers allows an infinity of vertical internal wave modes to exist, although the energy isusually partitioned into the first few modes (LaZerte, 1980). These vertical wave modes may52Chapter 5 Results and Discussionbe combined with any horizontal wave mode to form a wide spectrum of oscillations, (b)The presence of an external exciting force allows both forced and free oscillations to exist.The free oscillation is caused by the restoring force, gravity, which always tried toreestablish the equilibrium state. This free oscillation will eventually die out due to friction.The external exciting force is responsible for the forced oscillation which persists as long asthe exciting force is applied.5.2.1 Experiment 1This experiment consists of two numerical runs both having piecewise linear densityprofiles, as shown in figure (5-8a) where the two dotted lines indicate the interfaces betweenthe homogenous fluids and the interfacial region. Detailed experimental parameters for thesetwo runs are tabulated in Table (5-2). Where in the table h1, S and h2 are the individual layerthickness, from top to bottom, respectively. p1 and p2 are top and bottom layer densitiesrespectively. Tii are the natural periods of different internal wave modes calculated from theeigenvalue analysis, where i and j are horizontal and vertical internal wave mode numberrespectively. Therefore T11 represents the natural period of the fundamental first verticalwave mode which is of first mode in both horizontal and vertical directions, and T12represents the natural period of the fundamental second vertical wave mode which is of firstmode in horizontal direction and of second mode in vertical direction. The numericalprogram used in the eigenvalue analysis uses the shooting method to solve, in finitedifference form, the one-dimensional wave equation (equation 2-11), which is subject to theno-flow boundary conditions at both free surface and the bottom (equation 2-22).F =Q110-13 is densimetric Froude number, where h' is the vertical distance between theline sink and the density interface. Craya (1949) showed from dimensional analysis that thedensimetric Froude number is of primary importance in determining the nature of the flow53Chapter 5 Results and Discussioninduced by withdrawal through a line sink from a two-layered fluid with a sharp interface.The densimetric Froude Number F is redefined here as a bulk parameter such that theoverall density difference from the top to the bottom is used to calculate the reduced gravityg', and h' = kr(h2-1-0.5 8) is the vertical distance between the line sink and the center of theinterfacial layer.Figure (5-8) shows, from the eigenvalue analysis, the vertical profiles of horizontalvelocity components for the first four internal vertical wave modes associated with thedensity profile of run 1. The velocity profile for the first vertical wave mode appears to beproportional to the density profile. With top and bottom layers moving in the oppositedirections, the interfacial layer oscillates as a whole with little deformation. However, thesecond vertical wave mode is dominant in the interfacial layer causing the advection of fluidwithin this layer. With top and bottom layers now moving in the same direction as eachother and in the opposite direction to the interfacial layer, this leads to a horizontal variationof interfacial layer thickness. See figures (2-5c) and (2-5d) for a schematic of oscillatingpatters for the first two vertical wave modes. The third and the fourth vertical wave modeshave the similar oscillating mechanisms to the first and the second vertical internal wavemodes respectively, but with different magnitudes and periods. The natural periods of thesevertical wave modes get progressively larger. As a result it is expected that the set-up timefor these wave modes also get progressively larger. Heaps & Ramsbottom (1966) showedthat for a particular wave mode to set up, the forcing must persist at least one quarter that ofthe wave period.The development of the horizontal velocity at x = 0.5L is shown in figure (5-9), wherethe dotted lines indicate the upper and lower density interfaces. The first profile (t = 0.43T11)shows dominant first vertical internal wave mode oscillation with top layer velocity beinglarger than that of the bottom layer due to the constant forcing discharge. By t = 0.86T1154Chapter 5 Results and Discussion(t = 0.40T12) the jet-like flow is developed in the interfacial layer showing thecommencement of the second vertical wave mode. The second jet-like return flow shows upwithin the interfacial layer by the time of subsequent profile at t = 1.29T11 (t = 0.36T13)indicating the development of the third vertical wave mode. The fourth vertical wave modeis visible at t = 1.72T11 (t 0.33T14). Visual inspection of subsequent velocity profiles showsno higher vertical modes to appear, since the energy is usually partitioned into the first fewwave modes and higher wave modes are more severely damped, as shown by LaZerte(1980) and Crampin & Dore (1970). The reflections of the first and the second vertical wavemodes can be identified, while the movements of the third and the fourth vertical wavemodes are relative hard to follow.The spectral analysis of interfacial oscillations are shown in figure (5-10). Figure (5-10a)is the spectrum from the oscillation data of interfacial layer's center at x = 0. The dominanceof low frequency band is apparent with most of the energy concentrated in the range of0.0-0.1 Hz. Figure (5-10b) shows the spectra calculated from the oscillation data of bothupper interface and interfacial layer's center in the range of 0.0-0.015 Hz, where the upperinterface is the interface between the top and interfacial layers. The four peaks with periodsof 35.7, 76.0, 133.3 and 192.4 seconds, examined in comparison with the predicted periodsbased on the eigenvalue analysis, correspond to the fundamental first, second, third andfourth free vertical internal mode oscillations (see Table 5-3). The Observed resonantfrequencies are reasonably close to predictions with the maximum discrepancy of 8%. Thisdiscrepancy is mainly due to the viscous damping, as the predictions are based on thesolutions a linearized vertical momentum equation for an inviscid fluid. As mentioned thesecond vertical wave mode is responsible for the advection of fluid within the interfaciallayer thus causing the expansion and shrinking of the interfacial layer, therefore theoscillation data of upper interface has stronger second vertical mode signal than that from55Chapter 5 Results and Discussionthe center of the interfacial layer, as indicated in figure (5-10b). The significant peak withperiod of 20.4 seconds corresponds the free internal wave mode which is of second mode inhorizontal direction and of first mode in vertical direction. Noting that the natural period forthis mode T21 is more than half that of the fundamental first vertical mode T11, this is due tothe fact that higher modes are more severely damped. An additional feature of note in figure(5-10b) is the large peak at zero frequency for upper interface oscillation data. This peakcorresponds to the steady-state tilting of the upper interface due to the constant forcingdischarge. The steady-state tilting of the center of the interfacial layer is much smallercompared to that of the upper interface.Run 2 has a significantly smaller interfacial layer thickness than run 1. While bothexperimental runs have comparable fundamental first vertical mode periods T11, run 2 hasmuch larger higher vertical mode periods than run 1 due to the much thinner interfaciallayer, as shown in table (5-2). Figure (5-11) shows the development of horizontal velocitycomponent of run 2 at x = 0.5L. Again the vertical first wave mode is dominant at the firstprofile. By t = 0.84T11 (t = 0.25T12) the second wave mode has developed but with muchsmaller magnitude compared to that in run 1. Subsequent velocity profiles shows noevidence of higher vertical wave modes. With natural periods of higher vertical wave modesbeing much larger than that of the first vertical wave mode and little energy beingpartitioned into higher vertical wave modes, the first vertical wave mode is dominant and theinternal wave in this case effectively behaves as a two-layered fluid system, as shown byMortimer (1953) and Wiegand & Carmack (1986).5.2.2 Experiment 2This experiment consists of fourteen numerical runs. The density profiles used consist oftwo homogeneous fluids separated by a exponentially stratified interfacial layer, as shown in56Chapter 5 Results and Discussionequation (2-15). Two different values of reduced gravity are used with ,g; = 0.03 m/s2 and= 0.06 m/s2, here g' is defined as a bulk parameter. The relative top and bottom layerthicknesses, hi/L and h2/L, vary from 0.13 to 0.23, and the relative interfacial layer thickness&L varies from 0.04 to 0.25. The detailed experimental parameters are listed in table (5-4).Some comparisons as to the fundamental first vertical mode period T11 can be made herebetween numerical results and theoretical predictions, since for this particular case thetheoretical solution of Fjeldstad (1933) is available (equation 2-16). The comparison isshown in figure (5-12), where T11 is normalized with respect to T1, the fundamental periodin a two-layered fluid with a sharp interface when 6 = 0. The agreement is excellent withmaximum discrepancy being 1.2%. Also plotted in figure (5-12a) is the theoretical solutionof Phillips (1966) (equation, 2-14). It is not surprising to see the large discrepancy betweenthe numerical results and theoretical predictions of Phillips (1966) for large interfacialthickness 61L, Since equation (2-14) is based on the assumption that k6 << 1.5.3 Forced Internal Seiche Triggered by Unsteady WithdrawalThis section discusses four sets of numerical experiments, with different density profiles,varying from two-layered fluid with a sharp interface to piecewise linear density profile. Thevalue of viscosity used ranges from 10-6 m2/s to 10-5 m2/s. The discharge process issinusoidal and determined byQ = Q + Vsin(2,7rtIT f — 0.5x)^ (5 - 1)where Q aye = 4.610-4 m2/s is the average discharge rate through the line sink andAQ0.6Q„,,e is the amplitude of the discharge process. T1 is the discharge or forcing period,which is stepped through a series of values for each numerical experiment to test the57Chapter 5 Results and Discussionresponse properties of forced internal seiche. Detailed experimental parameters for all thenumerical experiments carried out in this section are listed in Table (5-5). The introductionof a sinusoidal discharge process will set up a forced internal seiche with the same period asthe forcing. At the same time a set of free internal wave modes are also set up as soon as thedensity interface is displaced from its equilibrium position. These free internal wave modesoscillate according to their own natural periods. The resulting oscillation is the combinationof both forced and free internal seiches.5.3.1 Experiment 1; Comparison with Physical ExperimentsThe physical experiments are carried out by Allan (1993) in the Hydraulic Laboratory of theUniversity of British Columbia. The experimental set up consists of a tank 2.40 m long, 0.60m high and 0.59 m wide, a programmable pump under computer control, a conductivitysensor. The tank is subdivided into a experimental section and a reservoir section by a centerwall 0.525 m high, 1.235 m from the discharge sidewall. Pumping occurs through a line sink0.30 m above the bottom of the tank. The water pumped out is returned to the reservoirsection of the tank and flows over the center wall to maintain the overall water depth in thetank. The conductivity sensor is mounted 10 cm from the discharge side on a traverse suchthat it may be raised and lowered to desired depths within the tank. For each experiment thevariation of the interface elevation is interpolated from the varying conductivity data of theprobe, which is located at the center of the density interfacial region. Figure (5-13) isredrawn after Allan (1993) showing the schematic of the experimental set up.Six numerical runs are carried out for comparison with the corresponding physicalexperiments. The density stratification is approximated by a piecewise linear density profile,which closely resembles that determined by the conductivity probe in the physical58Chapter 5 Results and Discussionexperiment. The viscosity used in the numerical experiments is 10-6 m2/s which isconsidered a typical value in the laboratory scale (Thorpe, 1968a, Allan 1993). During theexperiments the forcing period Tf is stepped through a series of values ranging above andbelow the natural period of the fundamental first vertical internal wave mode Tn. Thenormalized forcing period T1/T11 = 0.92, 0.94, 0.98, 1.00, 1.01 and 1.04 from run 1 to run 6respectively.The responses of density interface to the sinusoidal discharge process are plotted infigure (5-14) showing the time series of interfacial movement 10 cm from the discharge side(x = 0). Visual inspection of these plots shows a fairly good agreement, both in magnitudesand in phase angles, between numerical and physical experiments. Two major features worthof note are: (a) The large phase angle differences between numerical and physicalexperiments in the first few cycles, (b) The beating phenomena. The large phase angledifferences are caused by the free oscillations of previous physical experimental runs whichare repeated one after another as a sequence of tests. As it may take as long as 30 to 40cycles for these free oscillations to die out, therefore when the new experimental run isstarted these free oscillations may be dominant in the first few cycles. As time passes theforced oscillation gradually becomes dominate and those free oscillations turn into highhorizontal mode oscillations with decaying amplitude. The interaction between the forcedinternal seiche and the fundamental first vertical mode free internal seiche is responsible forthe beating phenomena which will be discussed later. The close matching of the beatingenvelopes between the numerical and physical experiments reflects the high accuracy ofnumerical model in modeling the forcing discharge process and in predicting the naturalperiod of internal seiche.59Chapter 5 Results and Discussion5.3.2 Experiment 2Experiment 2 consists of six numerical runs, which serves to investigate the dynamicresponse of a two-layered fluid with a sharp interface to the unsteady withdrawal. Thenormalized forcing period Tfai is stepped from 0.83 to 1.20, where T1 = 38.7 seconds is thefundamental natural period calculated from the linear inviscid theory of a two-layered fluidwith a sharp interface (equation 2-3).The time series of density interface movement at x = 0 for all six numerical runs arepresented in figure (5-15) which shows the predominance of the forced oscillation withforcing periods Tp A near-resonance occurs at run 3 with Tf = 40.0 seconds, which isdiscussed in detail here. The development of horizontal velocity component at x = 0.5L isshown in figure (5-16), where the dotted line indicates the equilibrium position of thedensity interface. The first three profiles of the flow establishment show that both upper andlower layer fluids move in the same direction towards the sink side wall. This flow patternof the initial set-up stage is different from that of a wind-driven internal seiche where thereis a wind-driven large scale circulation in the upper layer with a downwind surface flow anda upwind return flow in the lower part of the upper layer (Heaps and Ramsbottom, 1966).Subsequent velocity profiles show the typical flow pattern of the internal seiche in a two-layered fluid with a sharp interface, with progressively larger oscillating amplitudes as theforcing discharge continuously feed energy into this particular wave mode. Figure (5-17)shows the time series of the density interface movement at x = 0 as well as thecorresponding discharge process. Also shown in figure (5-17) are the time series of thehorizontal velocities at the central points of both upper and lower layers. The importantfeature of note is that the forced oscillation follows the sinusoidal discharge process, thusdoes not exhibits the common phase angle lag of a linear mass-spring-dashpot system asdepicted in figure (2-7b). This agree with the results of Hyden (1974) and Hodgins (1979).60Chapter 5 Results and DiscussionThe time history of the resonant internal seiche are better illustrated by a surface plot whichshows the interface movement throughout the whole length of the tank rather than at a singlepoint, this is shown in the colour photograph of figure (5-18).5.3.3 Experiment 3Experiment 3 differs from experiment 2 in several respects, most significantly, is theintroduction of a linearly stratified interfacial layer. This experiment consists of eightnumerical runs with normalized forcing period Tian ranging from 0.25 to 1.68, whereT11 = 32.2 seconds is the fundamental first vertical mode natural period calculated from theeigenvalue analysis. As mentioned in section 2.2.4, a linearly stratified interface can beapproximated, to the first order approximation, by an exponentially stratified interface if thedensity variation is small. A fundamental first vertical mode natural period T11 thuscalculated is 31.9 seconds (see equation 2-16). The responses of the interfacial layer to thesinusoidal discharge process are presented in figure (5-19), showing the variations of thedisplacement of interfacial layer's center from its equilibrium position at x= 0. The spectralanalysis of run 1, 2, 5 and 8 are presented in figure (5-20).Visual inspection of the interfacial oscillation data of run 1 with Tf = 8.0 seconds( figure 5 -19a ) shows two dominant oscillations, the forced oscillation with a forcing periodTf and the fundamental first vertical mode free oscillation with a period of about 30 seconds.It is seen that the forced oscillation reaches its steady state very quickly and persists as longas the discharge is sustained, while the free oscillation decays due to viscous damping. Theinternal wave modes are more clearly illustrated in the spectrum of figure (5-20a), whichshows six apparent energy peaks. The two peaks with periods of 33.3 and 8.0 secondscorrespond to the fundamental first vertical mode free oscillation and the forced oscillation61Chapter 5 Results and Discussionrespectively. While the remaining four apparent peaks with progressively less partitioningenergy correspond to the second, the third, the fourth and the fifth horizontal free internalwave modes respectively. Again noticed is the increasing of natural periods due to viscousdamping. An additional feature in figure (5-20a) is the small peak with the period of 122.0seconds, which corresponds to the fundamental second vertical mode free oscillation, withthe predicted value of T12 being equal to 118.2 seconds. The point raised here is that, withlittle energy being partitioned into the second vertical mode and T12 >> T11, the internalwave will effectively act as a two-layered fluid system with little deformation of theinterfacial layer. The time history of interfacial movement for run 8 with T1= 54.0 seconds ispresented in figure (5-19h), which shows the predominance of the forced oscillation. Thespectral analysis of this run is presented in figure (5-20d) showing a dominant energy peakof forced oscillation as well as three much smaller peaks corresponding to higher horizontalfree oscillation modes.A further important feature of figure (5-19) is the periodic beating phenomena. This iscaused by the interaction between the forced internal seiche and the fundamental firstvertical mode free internal seiche. As mentioned earlier, the forced oscillation has the sameperiod as the forcing Tf, while the fundamental first vertical mode free internal seicheoscillates according to its own natural period T11. If both oscillations have comparableperiods and magnitudes, the periodic beating occurs. Figure (5-20b) shows the spectralanalysis result of run 2 with Tf = 28.2 seconds. Noticing in this figure the two comparabledominant energy peaks with periods of 28.2 and 33.3 seconds, these are the exact internalwave modes responsible for the beating. However as the free oscillation is damped by theviscosity, the beating will get progressively smaller.A near-resonance occurs at run 5 with Tf=33.0 seconds, as indicated in figures (5-19e)and (5-20c). The development of the horizontal velocity components of this run at x = 0.5L62Chapter 5 Results and Discussionis shown in figure (5-21), where the dotted lines indicate the upper and lower interfaces. Thepredominance of a forced first vertical wave mode is apparent in the first three profiles. Byt = 1.39Tii ( t = 0.357'12) the jet-like flow pattern is developed in the interfacial layerindicating the commencement of the second vertical free internal wave mode. However dueto the fact that less energy is partitioned into the higher vertical wave modes for the case ofrelatively small interfacial layer thickness, and higher vertical wave modes are dampedmore, Subsequent profiles show no higher vertical wave modes to appear. As the forcingdischarge process continuously feeds energy to the force oscillation, the second vertical freewave mode gets progressively harder to identify.53.4 Experiment 4This experiment consists of five numerical runs. While both experiments 3 and 4 have alinearly stratified interfacial layer, experiment 4 differs from experiment 3 in that thethickness of the interfacial layer 8 is significantly increased, as shown in table (5-5). Thefundamental first vertical mode periods Tu of both experiments are comparable, however,due to the increased interfacial layer thickness, experiment 4 has a fundamental secondvertical mode period T12 of 83.4 second comparing with a T12 = 118.2 seconds in experiment3. Hence it is expected that in this experiment higher vertical internal wave modes areallowed to play a role in the interfacial oscillation.The interfacial movement of run 1 (Tf = 10.0 seconds) is presented in figure (5-22a),which shows the time series of the displacement of interfacial layer's center from itsequilibrium position at x = 0. Visual inspection of this figure shows a forced oscillation,which is superimposed on a fundamental first vertical mode free oscillation with decayingamplitude. The oscillation data of both upper and lower interfaces are shown in63Chapter 5 Results and Discussionfigure (5- 22b) , which shows the time series of the displacements of both upper and lowerinterfaces from their equilibrium positions at x = 0. Note the expansion and shrinking of theinterfacial layer, this is due to the second vertical wave mode. The second vertical modeoscillation is better illustrated in figure (5-22c), which shows the time series of interfaciallayer thickness at x = 0. Two apparent oscillations can be observed. The first one with aperiod of about 10 seconds is the forced second vertical wave mode caused by the forcingdischarge, while the second one with a period of about 80 seconds corresponds to the freesecond vertical wave mode. It can be seen, after comparing figure (5-22a) with (5-22c), thatthe second vertical free wave mode is more severely damped compared to the first verticalfree wave mode. Figure (5-22d) shows the spectra calculated from the oscillation data ofboth upper interface and the center of the interfacial layer. The three apparent energy peakswith periods of 34.6, 83.3 and 156.3 seconds correspond to the fundamental first, second andthird vertical free internal wave modes respectively. The predicted natural periods for abovewave modes, from the eigenvalue analysis, are 34.2, 83.4 and 143.4 seconds respectively.The energy peak with a period of 20.0 seconds corresponds to the internal wave mode whichis of second mode in horizontal and of first mode in vertical, and the energy peak with aperiod of 10.0 seconds corresponds to the forced oscillation. An important feature of figure(5-22d) is that, owing to the increasing of the interfacial layer thickness, the first threevertical free internal wave modes are now in the same order of magnitude.Run 2, run 3 and run 4 have similar oscillation patterns, as shown in figures (5-23),(5-24) and (5-25). Run 3 with Tf = 33.0 seconds is discussed in detail here. With forcingperiod Tf being close to T11, the forced first vertical wave mode oscillation is selectivelyamplified, as shown in figure (5-24a). The predominance of the first vertical wave modeoscillation can also been found in figure (5-24b), where the lower interface almost followsthe movement of the upper interface indicating little deformation of interfacial layer, thus64Chapter 5 Results and Discussionless second vertical wave mode oscillation. Figure (5-24c) shows a dominant forced secondvertical wave mode with much smaller magnitude compared with that of the forced firstvertical wave mode. The predominance of the forced first vertical wave mode is also foundin the spectra analysis of figure (5-24d), where the energy peak with a period of Ti- isdominant.Run 5 has a forcing period T1 of 80.0 seconds, which is close to the natural period offundamental second vertical wave mode T12. With Tf being close to T12, a near-resonantsecond vertical mode is created. With upper interface and lower interface being out of phase,as shown in figure (5-26b), the second vertical wave mode is dominant in the interfacialoscillations, as shown in figure (5-26c). The dominance of the second vertical mode is alsoshown in the spectrum of figure (5-26d). The development of horizontal velocitycomponents at x = 0.5L is shown in figure (5-27). The first profile (t = 0.24T12) shows thecommencement of the first vertical wave mode. By t = 0.48T12 the jet-like return flow isdeveloped in the interfacial layer showing the development of the forced second verticalwave mode. The third vertical wave mode shows up at the t = 0.96T12. Subsequent profilesshow the predominance of the forced second vertical mode oscillation, as the forcingdischarge continuously feeds energy into this particular mode.65Chapter 6Conclusions and RecommendationsThe resonant internal seiche in a viscous and incompressible fluid triggered by unsteadywithdrawal from a rectangular tank is investigated numerically in this thesis using the CICmethod. In the numerical model, two density profiles are used: a two-layered fluid systemwith a sharp interface and a two-layered fluid system with a diffuse interface. The viscouseffects are represented by an eddy viscosity which is assumed constant and uniformthroughout the fluid. The following conclusions can be drawn from this investigation:(a) The CIC method used here combines the best features of both Lagrangian andEulerian approaches, and is most advantageously applied to the fluid dynamics problemswith density stratification. The numerical model solves, in finite difference form, thegoverning equations in terms of stream function and vorticity, therefore avoids thecomplexity of solving the pressure field. The free surface boundary conditions in this caseare more difficult to apply, however, for the case of small density variation, simplified firesurface boundary conditions may be chosen without affecting the main features of the flow,as shown in the comparison with the parallel physical experiments of Allan (1993).(b) The two-layered fluid system with a sharp interface allows only one vertical internalwave mode to exist. With the fluids in the upper and lower layers moving in the oppositedirections, large velocity shear develops at the density interface. The observed fundamentalnatural periods T1 agree well with the predictions based on the linear inviscid theory of atwo-layered fluid with a sharp interface (equation, 2-3). The discrepancies range from 1.1%to 6.8% with observed values being always larger than the predictions. This increasing of66Chapter 6 Conclusions and Recommendationsnatural periods is mainly due to the effect viscous damping and finite interface thickness, asshown in the physical experiments of Thorpe (1968a), Hyden (1974) and Allan (1993). Thenumerical experiments also show that higher internal wave modes are more heavily damped,this is supported by the analytical solution of Harrison (1908) and the numerical solutions ofCrampin & Dore (1970).(c) The introduction of a continuously stratified interfacial layer between twohomogeneous fluids enables more than one vertical internal wave modes to exist,furthermore the dispersion relation changes as the density profile changes. The naturalperiods of a two-layered fluid system with an exponentially stratified interfacial layer arewell predicted by the numerical model, in comparison with the theoretical solutions ofFieldstad (1933) (equation, 2-16). However large discrepancies are found between thenumerical results and the theoretical solutions of Phillips (1966) (equation, 2-14) for largeka because equation (2-14) is based on the assumption that the interface thickness is smallcompared to the wave length, i.e. a << 1. The numerical experiments with relatively smallinterface thickness show only two vertical internal wave modes. With less energy beingpartitioned into the second vertical wave mode and T12 >> T11, the internal wave in this caseeffectively acts as in a two-layered fluid system with a sharp interface, with littledeformation of the interfacial layer. This is supported by the observations of Wiegand &Carmack (1986). The numerical experiments with large interfacial layer thickness shows upto five vertical internal wave modes.(d) Both free and forced internal wave modes are created by unsteady withdrawal. Thefree internal wave modes oscillate according to their own natural periods, while the forcedinternal wave mode follows the forcing of the discharge. Periodic beating phenomena isobserved as a result of the interaction between the forced internal wave mode and thefundamental first vertical free internal wave mode. This periodic beating gets progressivelysmaller, as the free internal wave mode decays due to viscosity.67Chapter 6 Conclusions and RecommendationsSome recommendations as to further research are also made here. More modificationscan be made to the numerical model to make it more efficient and accurate. The use of anon-uniform grid system allows finer mesh in the density interfacial region without largelyincreasing C.P.0 time. Furthermore the multigrid technique can be employed so as to speedup the iteration process. The development of a more versatile numerical model with velocitycomponents and pressure being the primitive variables is recommended. This model wouldallow more accurate free surface boundary conditions and both free surface and internalwaves to exist for the case of large density variation. Other physical problems can be lookedinto. One of them would be the effect of the Coriolis force and real bathymetry on thenatural periods of internal seiches, as most analytical solutions are based on the assumptionthat the lake or reservoir basin is box-shaped and the flow is two-dimensional. Mortimer(1979) presented a one-dimensional numerical model for a two-layered fluid system with asharp interface which accounts for the cross-section variation of a lake, however a moregeneral three-dimensional model which may account for arbitrary density distributions andreal basin topographies still needs to be developed. Finally the effects of viscous dampingand finite amplitude on internal waves need to be investigated in more depth.68BibliographyAllan, G. (1993). "Experiments with forced resonant internal waves." M.A.Sc. Thesis, Univ.of British Columbia, Canada.Birdsall, C.K. & Fuss, D. (1969). Clouds-in-Clouds, Clouds-in-Cells Physics for Many-Body Plasma Simulation. J. of Computational Physics. 3,494-511.Carmack, E.C. & Gray, C.B.J. (1982). Patterns of Circulation and Nutrient Supply in AMedium Residence-time Reservoir Kootenay Lake, British Columbia. Canadian WaterResources Journal. 7(1), 51-69.Crampin, D.J. & Dore, B.D. (1970). Numerical Comparison of the Damping of InternalGravity Waves in Stratified Fluids. Pure Appl. Geophys. 79(2), 53-65.Craya, A. (1949). Recherches Th6oriques Sur L'Ecoulement De Couches Superposdes DeFluides De Densites. La Houille Blanche. Jan-Feb., 56-64.Defant, A. (1932). Beitrage zur Theoretischen Limnologic. Bertr. Freien Atmos. 29, 143-150.Dix, D.M. (1963). The Magnetohydrodynamic Flow Past A Non- conducting Flat Plate inthe Presence of a transverse Magnetic Field. J. of Fluid Mech. 15, 449-476.Evans, M.E. & Harlow, F.H. (1957). The Particles-in-cell Method for HydrodynamicCalculations. Los Alamos Scientific Lab. Rept. La-2139, Los Alamos, New Mexico.Fjeldstad, J.E. (1933). Interne Wellen. Geofys. Publi. 10(6), 1-35.Forsythe, G.E. & Wasow, W.R. (1960). "Finite Difference Method for Partial DifferentialEquations." New York, Wiley.Hamblin, P.F. (1978). Internal Kelvin Waves in a Fjord Lake J. of Geophysical Research.83, 2409-2418.69Harlow, F.H. & Welch, J.E. (1965). Numerical Calculation of Time-Dependent ViscousIncompressible Flow of the Fluid with Free Surface. Physics of Fluids. 8(12), 2183-2189.Harrison, W.J. (1908). The Influence of Viscosity and Capillarity on Waves of FiniteAmplitude. Proc. London Math. Soc. Ser. 2. 7, 107-121.Heaps, N.S. & Ramsbottom, A.E. (1966). Wind Effects on the Water in a Narrow Two-layered Lake. Phil. Trans. Roy. Soc., London. A259, 391-430.Hunt, J.N. (1961). Interfacial Waves of Finite Amplitude. La Houille Blanche. 16, 515-531.Imberger, J. (1980). Selective Withdrawal: A Review. Proc. 2nd Intl. Symp. on StratifiedFlow. IAHR, Trondheim, Norway. 381-398.Imberger, J., Thompson, R. & Fandry, C. (1976). Selective Withdrawal from a FiniteRectangular Tank. J. of Fluid Mech. 78(3), 489-512.Keller, J.B. & Munk, W.H. (1970). Internal Wave Wakes of a Body Moving in a StratifiedFluid. Phys. of Fluids. 13, 1425-1431.Krauss, W. (1966). "Methoden und Ergebnisse der Theoretischen Ozeanographie, II-InterneWellen." Berlin, Gebriider Borntraeger.Lamb, N. (1932). "Hydrodynamics." Dover Publications, New York.Lawrence, G.A. & Allan, G. (1991). "Resonant Seiche Enchanced Mysis Export fromKootenay Lake." A Proposal to Investigate, Univ. of British Columbia, Canada.LaZerte, B.D. (1980). The Dominating Higher Order Vertical Modes of the Internal Seichein a Small Lake. Limnol. Oceanogr. 25, 846-854.Mortimer, C.H. (1952). Water Movements in Lakes During Summer Stratification; Evidencefrom the Distribution of Temperature in Windermere. Phil. Trans. Roy. Soc., London. B236,355-404.Mortimer, C.H. (1979). Strategies for Coupling Data Collection and Analysis with DynamicModeling of the Lake Motions. Hydrodynamics of Lakes. Elsevier Sci. Publishing Co., NewYork. 183-222.70Noh, W.F. & Rotter, M.H. (1963). Equations of Hydrodynamics. J. of Math & Mech. 12,149-191.Patankar, S.V. & Spalding, D.B. (1971). A Calculation Procedure for Heat, Mass andMomentum Transfer in Three-dimensional Parabolic Flows. Int. J. of Heat & MassTransfer. 15, 1787-11806.Phillips, N.A. (1959). An Example of Non-linear Computational Instability. Atmosphere andSea in Motion. Rossby Memorial Volume, Rockefeller Institute Press, New York.Phillips, O.M. (1966). "The Dynamics of the Upper Ocean." Cambridge Univ. Press,Cambridge.Roache, P.J. (1972). "Computational Fluid Dynamics." Hermosa Publishers, Albuquerque.Roberts, J. (1975). "Internal Gravity Waves in the Ocean." Marine Science, Marcel DekkerInc., New York.Schmidt, W. (1908). Stechende Schwingungen in der Grenzschicht Zweeier Fliissigkeiyen.S.B. Akad. Wiss. Wien, Ha, 117,91-102.Stevens, C. & Imberger, J. (1993). The Initial Response of a Stratified Lake to a SurfaceShear Stress. Submitted to J. of Fluid Mech.Thorpe, S.A. (1968a). On Standing Internal Gravity Waves of Finite Amplitude. J. of FluidMech. 32(3), 489-528.Thorpe, S.A. (1968b). On the Shape of Progress Internal Waves. Phil. Trans. Roy. Soc.,London. A263, 563-614.Thorpe, S.A. (1971). Near-resonant Forcing in A Shallow Two-layer Fluid: A Model for theInternal surge in Loch Ness? J. of Fluid Mech. 63(3), 509-527.Thorpe, S.A. (1977). Turbulence and Mixing in A Scotish Loch. Phil. Trans. Roy. Soc.,London. A286, 125-181.Ting, F.C.K. (1992). On Forced Internal Waves in a Rectangular Trench. J. of Fluid Mech.235, 255-288.71Von Neumann, J. & Richtmyer, RD. (1950). A Method for the Numerical Calculations ofHydrodynamic Shocks. J. of App!. Phys. 21, 232-237.Watson, E.R. (1904). Movement of the Water of Loch Ness As Indicated by TemperatureObservations. Geogr. J. 24, 430-437.Wedderburn, E.M. (1912). Temperature Observations in Loch Earn; with a FurtherContribution to the Hydrodynamical Theory of the Temperature Seiche. Trans. Roy. Soc.,Edinburgh. 48, 629-695.Wedderburn, E.M. & Williams, A.M. (1911). The Temperature Seiche, Part II,Hydrodynamical Theory of Temperature Oscillations. Trans. Roy. Soc., Edinburgh. 47, 628-634.Wedderburn, E.M. & Young, A.W. (1915). Temperature Observations in Loch Earn. Trans.Roy. Soc., Edinburgh. 50, 741-767.Wilson, B.W. (1972). Seiches. Advances in Hydroscience. 8, 1-94. Academic Press, NewYork.72Run 1 2 3 4 5 6 7 8 9L (m) 2.40 2.40 1.20 1.20 2.40 2.40 1.20 1.20 1.20h1 (n) 0.30 0.30 0.30 0.30 0.30 0.30 0.30 0.30 0.45h2 (m) 0.30 0.30 0.30 0.30 0.30 0.30 0.30 0.30 0.15g' (m/s2) 0.06 0.06 0.06 0.06 0.02 0.02 0.08 0.08 0.06Ve (m2/S) 10-4 10-6 10-4 10-6 10-4 10-6 10-4 10-6 10-4Ti (theoretical) (s) 52.3 52.3 26.2 26.2 90.7 90.7 24.2 24.2 31.6Ti (numerical) (s) 53.9 53.3 26.9 26.5 93.9 92.8 25.5 25.1 33.8Table (5-1). Free internal seiche in a two-layered fluid with a sharp interface;experimental parameters and results.Parameter Run 1 Run 2L (m) 1.20 1.20h1 (m) 0.15 0.258 (m) 0.30 0.10h2 (m) 0.15 0.25P1 (kg/m3) 1000.0 1000.0N (Hz) 0.44 0.63P2 (kg/m3) 1006.0 1004.0v (n2/s) 10-6 10-6F 0.036 0.043T11 (s) 34.9 38.8T12 (s) 75.5 130.1T13 (s) 125.8 244.2T14 (s) 178.2 357.2Table (5-2). Forced internal seiche triggered by constantwithdrawal, experiment 1; experimental parameters.73Mode Predicted Value Observed ValueTil ,^34.9 s 35.7 sT12 75.5 s 76.0 sTi3 125.8 s 133.3 sT14 178.2 s 192.4 sTable (5-3). Forced internal seiche triggered by constantwithdrawal, experiment 1, run 1; comparison of observed andpredicted natural periods of fundamental internal wave modes.Run 1 2 3 4 5 6 7 8 9 10 11 12 13 14L (m) 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2h1 (C111) 27.5 25.0 22.5 20.0 17.5 15.0 20.0 15.0 27.5 25.0 22.5 20.0 17.5 15.08 (Cm) 5.0 10.0 15.0 20.0 25.0 30.0 20.0 30.0 5.0 10.0 15.0 20.0 25.0 30.0(CIII),h2 27.5 25.0 22.5 20.0 17.5 15.0 20.0 15.0 27.5 25.0 22.5 20.0 17.5 15.0g bulk (11/s2) 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.03 0.03 0.03 0.03 0.03 0.031.) (M2/S) 104 10-6 104 10-6 10-6 104 icrs 10-4 10-6 10-6 10-6 10-6 10-6 10-6Tn(theoretical) (s) 28.9 29.9 31.0 32.1 33.4 34.7 32.1 34.6 40.9 42.3 43.8 45.4 47.2 49.0Tnthumerican (s) 28.7 30.1 30.9 31.9 33.3 34.7 32.5 35.0 40.7 42.6 43.9 45.0 46.7 48.7Table (5-4). Forced internal seiche triggered by constant withdrawal, experiment 2;experimental parameters and results.74Parameter Exp 1 Exp 2 Exp 3 Exp 4L (m) 1.235 1.2 1.235 1.235h1 (m) 0.345 0.45 0.323 0.2258 (m) 0.044 0 0.075 0.200h2 (m) 0.136 0.15 0.127 0.100Pi (kg/m3) 1001.42 1000.00 1000.08 1000.08N (Hz) 1.21 — 0.95 0.58P2 (kg/m3) 1008.04 1004.00 1006.96 1006.96V (n%) 10-6 10-6 10-5 10-5AQ (m2 /s) 2.76*10-4 2.76*10-4 2.76*10-4 2.76*10-4F 0.029 0.067 0.034 0.020T11 (s) 32.3 34.9 32.2 34.2T12 (s) 150.8 — 118.2 83.4T13 (s) 281.5 — 220.3 143.4Table (5-5). Forced internal seiche triggered by unsteady withdrawal;experimental parameters.751.0epilimnionthermocline0.6 -0.4 -hypolimnion0.2 -0.0 10.8 -6.0^8.0^10.0^12.0^14.0^16.0Temperature (co)Figure (1-1). Temperature profile from Windermere northern basin showingthe typical temperature profile of a lake during summer (redrawn fromMortimer, 1952).76Qy=02h^Pihih2^P2x=0^ x=Ly=HSinkFigure (1-2). Defming sketch of a two-layered reservoir model770 10^30^SO^70TOCOLUMBIARIVERCORRA LINNDAMFigure (1-3). Location map for the Kootenay river drainage basin (fromCarmack & Gray, 1982).782HR^ 4HR 6HR^ 8HRLOCH EARNAug. 9,191110HR^ I2HR^ I4HRFigure (2-1). (a) Hourly positions of the thermocline bounded bythe 90 and 11°c isotherms on a longitudinal section of Loch Earnredrawn from Wedderburn (1912), (b) general flow pattern of freeinternal seiches in a two-layered lake, from Wedderburn &Williams (1911).79(c)P2h2L77=tFirst Horizontal Mode (n=1)^hi PiSecond Horizontal Mode (n=2)(b)Third Horizontal Mode (n=3)Figure (2-2). Interface profiles of internal seiches in a two-layered flukl with asharp interface; (a) first or fundamental horizontal mode (n=1), (b) secondhorizontal mode (n=2), (c) third horizontal mode (n=3).8002-kiNkTo= 2.5k/400= 5.0k/i/cro= 7.50.060.040.080.1or."0^0.4^f 0.6^0.80.54^tanh kh,Figure (2-3). The variation of tanhkht and tanhkh2 for which thevariation of natural frequency with wave amplitude changes. In regionA, the frequency decreases with increasing wave amplitude; In regionB, the frequency increases with increasing wave amplitude (fromThorpe, 1968a).0.0000^0.0002^0.0004^0.0006^0.0008^0.0010Viscosity (m2/s)Figure (2-4). Viscous effect on natural frequency of internal seiche in atwo-layered deep fluid.81yIH1yIHMode 2Mode 1(b) Figure (2-5). Schematic plot of a three-layered fluid system; (a) the density profile,(b) the associated two vertical internal wave modes, (c) the oscillation mechanism ofthe first vertical wave mode, (d) the oscillation mechanism of the second verticalwave mode (redrawn from Stevens & Imberger, 1993).82Y AH P1h2+8 -,-„.--„p=p1exp((N2g-1(h2+8-y)) --%%-,----..,h28p2=p1exp(N2g-18)lb-Figure (2-6). Density profile of a two-layered fluid system with aexponentially stratified interfacial layer.83(a)P = P0exP{-4P12P0 tanNY18))0.50.4 —0.3 —Cl.0.2 -First vertical wave modeSecond vertical wave modeThird vertical wave mode 0.1 —(b)2.00.0^ 4.0^6.0^8.0^10.0kSFigure (2-7). (a) Hyperbolic tangent density profile in deepfluid, (b) the associated dispersion relationship.841^1^t^I/DAMPING FACTOR (2)I^1,a ■%r 1r -\( a )7.4.01—20 0.2 04 06 0.8 1012 1.4 I 6 I 8 2 0 2.2 24 2 6 2.8 3 0RATIO OF FORCING AND NATURAL FREQUENCIES, COAT()g/IPIB_■__!_I_•_w..•...,_4F''?t:"-----.■- -;'"--^------ .- --iffir///..luioc, 5OMOM0.40CD*0.50,,,,, DAMPING FACTOR (213)I( b)___.,-,.—/.-'.'•PFRATIO OF FORCING AND NATURAL FREQUENCIES, ca/CrFigure(2-8). Linear damped vibration of mass-spring-dashpot system: (a)dynamic amplification factor /.2 and (b) phase angle a, as functions offrequency ratio co/cro and damping factor 2/3 (from Wilson, 1972).1807; 150CDCD-23 120U ; 90-J1.1.1 603000 0.2 04 06 0.8 1.0 I 2 I 4 1.6 I 8 2.0 2.2 2 4 2 6 2 8 3 0850-150-160-140-12Z:1\CN1^0-080-060-040-020-50^0-55^0-60^0-65^0-70^0-75a(rad I s)Figure (2-9). Response curves for internal seiches in deep fluids withdp/p2=4.7*10-3, where 2a/A is the wave steepness. The symbols representdifferent total plunger amplitudes (from Thorpe, 1968a).0-85^0-90^0-95 f^1 -00^105a(rad I s)Figure (2-10). Response curves for internal seiches in shallow lower fluid,kh2=0.48, and deep upper fluid, and with Ap/p2=4.7*10-3 (from Thorpe,1968a).0-05 -0-040-03C•10-020-0186Figure (3-1). Boundary conditions for stream function ty.87Axu(i+1,j+1)O v(i+1,j+1)c(i+ 1,j+ 1)u(i+ 1,j)v(i,j)Axu(i+2,j+ 1)v(i+2,j+1)c(i+2,j+ 1)u(i+2,j)v(i+2,j)^L.) c(i+2,j)v(i.;+i)u(i,j)A yu(i,j+ 1 )Figure (4-1). Staggered grid.I^ij+ I I i+jK th Pakicle•i-F.- Ir-1...--t=a; rgivirayrk N.. 4L.Ak 'I' 1 .ab,, \ I 1-43I- I- -I IIFigure (4-2). Bilinear interpolation.j+188PiP2^I-,1 250•ar•I3001— 11*41 grid- - 21*41 grid---. 41*41 gridi100^150^200Time (sec.)Figure (4-3a). A sketch showing locations of points A and B.Normalized Stream Function vs. Time (at point A)Figure (4-3b). Grid size dependence tests on dx; stream function versustime at point A .89(c)^ 21*21 grid21*41 grid 1.200 250^30041*41 grid41*61 grid^ -t^50^100^150(d)2015a 10Normalized Stream Function vs. Time (at point B)2520151050^100^150^200^250^300Time (sec.)Figure (4-3c). Grid size dependence tests on Ay, coarse mesh; stream functionversus time at point B.Time (cerFigure (4-3d). Grid size dependence tests on Ay, fine mesh; stream functionversus time at point B.901••••(a)120^14020^40^60^80— At = 0.0125 s^ At = 0.025 sAt– 0.05 s— At = 0.0125 s^ At = 0.025 s&=0.05s^t-i'—(b)4.Normalized Stream Function vs. Time (at point A)Time (sec.)Normalized Stream Function vs. Time (at point B)20^40^60^80^100^120^140Time (sec.)Figure (4-4). Time step At denpendence tests; (a) stream function versus timeat point A, (b) stream function versus time at point B.2.52.01 51.00.5910.170 16O *it 0.15El0.140.180.13100 300200Time (sec.)0.17--;^ .0 16= .2ttl4 . )› 0.15w0.140.180.13100 300200Time (sec.)hI , 'T I40.180.170.140.13v.100^200^300Time (sec.)Figure (4-5). Time history of the interface elevation at x =0 showing aliasingerrors increase with decreasing viscosity.92Ii400 — -500 --200600300 --• Numerical Results ve =10-4m2/s• Allan, 1993o Hyden, 1974A^Numerical Results v =10-6m2/s9080706050 ^40 ^30 ^20 ^20 30^40 50 60 70 80 90 200^300 400 500 600Tl(theoretical) (sec.)Figure (5-2). Free internal seiche in a two-layered fluid with a sharpinterface; comparison of observed fundamental natural periods with thetheoretical values based on the linear inviscid theory of equation (2-3).940.3150.310g 0.305S•-,•tcla0.3000.2950.290 .0I50^100^150Time (sec.)1 ^ 1•r_(a)Run 3MI = kh2 = 0.78I^I^I200 250 300Time (sec.)Figure (5-3). Free internal seiche in a two-layered fluid with a sharp interface,run 3 and run 9; time series of interface elevation at x=L showing the dampingprocess of free internal seiches.95Run 3Run 9Figure (5-4). Free internal seiche in a two-layered fluid with a sharpinterface, run 3 and run 9; surface plot of interface movement showingthe damping process of free internal seiches.96Elevation of Interface at X=0100 150 250200 300Horizontal Velocity at the Central Point of Upper Layer20-4-6z10.350100 200150 300250Horizontal Velocity at the Central Point of Lower Layer2-6x10-3500.3150.3100.3050.3000.2950.2900.28550 100 150 200 250 300Time (sec.)0.3150.3100.3050.3000.2950.2900.28550^100^150^200^250^300Time (sec.)Time (sec.)Time (sec.)Figure (5-5). Free internal seiche in a two-layered fluid with a sharp interface, run 3;(a) and (b) time series of density interface movement at and x=L, (c) and (d) timeseries of horizontal velocity at central point of upper and lower layers.971.0 .Horizontal Velociz Profile At X = 0.5L 0.6 —04 —0.004-0.002^0.000U (m/s)Hcrizontal Vclocityprofile At X = 0.5L0.0040.002-0.004 -0.002 0.000U (zo/4)0.8 —1.00.80.60.40.20.0-0.004 -0.002 0.000U (m/s)0.002 0.0040.004-0.002-0.004 0.0020.000U(m/4)Horizontal Velocityyrofile At X = 0.5LHorizontal Velocity.Profile At X = 0.5L1.00.80.6OA0.20.01.00.80.60.40.20.0Figure (5-6). Free internal seiche in a two-layered fluid with a sharp interface, run 3;horizontal velocity profiles at a time interval of 0.5T1 showing the large velocityshear at the density interface.98Figure (5-7). Free internal seiche in a two-layered fluid with a sharp interface, run 3;velocity vector plots over one complete cycle.99-1 -0.5 0 0.5 10.9 —0.8 —0.7 —0.6 —..,1 0.5 —0.4 —0.3 —0.2 —0.1Top LayerInterfa al Layer(a)IIIIII0.0 1.0 2.0 3.0 4.0 5.0 60Bottom LayerDensity Profile^Vertical Internal Wave Modesat ( kg/m3 ) Normalised UFigure (5-8). Forced internal seiche triggered by constant withdrawal,experiment 1, run 1; (a) piecewise linear density profile, (b) first four verticalinternal wave modes associated with the density profile of (a).1000.002^0.004^ -0.002^0.000^0.002^0.004U (m/s)-0.004 -0.002 0.000U (m/s)Horizontal Velocity Profile At T = 0.43 T11-0.004^-0.002^0.003^0.002^0.004U (m/s)Horizontal Velocity Profile At T =129 T111.0:^0.8 ^z 0.60A02 -0.0 -0.002^0.000^0.002^0.004U (m/s)Horizontal Velocity Profile At T = 2.15 T111.00.8z 0.60.40.20.01.00.8z 0.6)'• OA020.0-0.001 0.000U (m/s)Horizontal Velocity Profile At T = 1.72 T110.0040.0021.00.8z 0.6>" OA020.0-0.002^0.000^0.002^0.004UHorizontal Velocity Profile At T = 0.86 T111.00.8z 0.6• OA020.0Horizontal Velocity Profile At T = 2.58 T11 1.00.8 ^1z 0.6o4-02E.^0.0-0.004 -0.002 0.000(mis)0.002^0.0041.0 :0.8 ^0.6OA :-02 EOA^;S.1.00.8z 0.60.40.20.01.00.8z 0.6>* 0.40.20.0-0.004^-0.002^0.003^0.002^0.004U (0V.)Horizontal Velocity Profile At T = 4.73 T11Horizontal Velocity Profile At T = 5.59 T11^L^1.00.8z 0.6>" OA020.0Horizontal Velocity Profile At T = 3.01 T11Horizontal Velocity Profile At T = 3.87 T11Horizontal Velocity Profile At T = 3.44 T11 1.0z 0.60.80.40.2.0 -0.004^-0.002^0.000^0.002^0.004U (n)Horizontal Velocity Profile At T = 4.30 T110.8 ^ L____,........„._,....^ 40.6 iOA -:0.2^ 10.0 :Figure (5-9). Forced internal seiche triggered by constant withdrawal, experiment 1,run 1; horizontal velocity profiles at x = 0.5L.1010.000^0.072^0.004U (m/s)Horizontal Velocity Profile At T = 5.16 T11 -0002-^0.000^0.002^0.004U (to/s)Horizontal Velocity Profile At T = 6.02 T111.00.8z 0.60- OA020.01.0:0.8 0.6>"' 0A0.2 r0.0 -0.004^-0.002^0.000^0.002 0.004U (zo/s)Horizontal Velocity Profile At T = 6.45 T111.00.8z 0.60.40.20.0-0.002^0.000^0.032 0.004U (mu)0.8 ^0.6OA E.-020 .0 =-0.004 0.000^0.002^0.004U (m/s)-MOOS 0.000^0.032^0.004U(n)Horizontal Velocity Profile At T = 6.88 T11Frequency (Hz)Figure (5-10). Forced internal seiche triggered by constant withdrawal, experin -ent 1,run 1; (a) spectrum from the oscillation data of the center of interfacial layer in therange of 0.0-0.16 Hz, (b) spectra from the oscillation data of both upper interface andthe center of interfacial layer in the range of 0.0-0.015 Hz.1020.004-0.002^0.000^0.002U (mh)Horizontal Velocity Profile At T = 1.68 T11-0.004 -0.002 0.000^0.002^0.004U ('Ws)Horizontal Velocity Profile At T = 3.36 T110.0040.004 0.0020.002 -0.004-0.004 -0.002 -0.0020.000U (m/4)0.000U (m/s)Horizontal Velocity Profile At T = 0.42 T11 Horizontal Velocity Profile At T = 0.84 T111.00.80.6OA020.00.80.6OA020.0-0.004 -0.002 0.0040.0020.000U (rots)1.00.80.6 OA E-0.2^0.0 ^-0.004^-0.002^0.000^0.002^0.004U (m/s)Horizontal Velocity Profile At T = 1.26 T111.0Horizontal Velocity Profile At T = 2.10 T11 Horizontal Velocity Profile At T = 2.52 T111.00.80.6OA0.200orlonU (IWOHorizontal Velocity Profile At T = 2.94 T11-0.004 -0.002 0.0040.0021.00.80.6OA020.01.00.8 :-0.6 ^OA r^020.0^-0.004^ 0.000^0.002.^0.004U(u)-0.0041.00.80.6OA020.01.00.80.6OA020.0Figure (5-11). Forced internal seiche triggered by constant withdrawal, experiment 1,run 2; horizontal velocity profiles at X = 0.5 L.1031.61.5 — e1.4 —Numerical Results, ue=10-4m21sNumerical Results, .10-6m2/sFjeldstacl, 1933Phillips, 1966."."4^1.3-1.2 —•1.1(a)10.1 0.3 0.40.2OIL(.) 0.1 0.2 0.3 0.481,Numerical Results, .10-6m2/sFjeldstad, 19331.41.31.2 — .11.1 —(b)5•0•1.61.5Figure (5-12). Forced internal seiche triggered by constantwithdrawal, experiment 2; comparison of observed fundamentalfirst vertical mode natural periods with the theoretical values, (a)g 'buik=0 .06 m/s2, h1/L=h2/L=0.04-0.25, and (b) g 'buik=0 .03h1lL=h2/L=0.04-0.25.104Experimental Section Reservoir SectionFigure (5-13). Forced internal seiche triggered by unsteady withdrawal, experiment 1,comparison with physical experimnets; schematic of experimental setup (redrawn fromAllan, 1993).105Displacement (mm)Displacement (mm)8^8 8 .,' 000" 8o I _ §.0.§ t.., p - § g §1■,0Displacement (mm)Displacement (mm)Displacement (mm)oOs"1■•• o2010-2020 va^'Run 3, Tr= 40.0 s.I10-10( c )-20 ^0^50 100^150^200^250Time (sec.)20--20100 150 200 250Time (sec.)'Run 4, Tr= 42.0 s.I100^150^200^250Time (sec.)20101010-200^50^100 150 200 250Time (sec.)0^50^100^150 200 250Time (sec.)0^50^100 150 200 250Time (sec.)Figure (5-15). Forced internal seiche triggered by unsteady withdrawal, experiment 2;time series of interfacial oscillations at x = 0.108Horizontal Velocityprofile At T = 0.13 Ti^ Horizontal Veloci Profile At T = 0.26 Ti ^1.0^ 1.0^OA "z 0.6 z 0.61>'' OA^ )' OA02 02OA 1. ^OA-0.004^-0.002^0.000^0.002^0.001U (m/t)Horizontal Velocity Profile ALT = 0.39 Ti 1.0^ I0.8 Iz 0.6 zOA ^ I0.20.0 1 -0.001^-0.002^0.000^0.002^0.00$U (Mk)-0404^-0.002 0.003 0.002U (m/s)Horizontal Velocity Profile ALT = 0.52 Ti1.00.80.6OA0.20.0-0.004 -0.032 0.000U (iw.)aimHorizontal Velocity Profile At T = 0.65 Ti 1.0^ 1.00s- 0.8z o4- ^I ^0.6OA ^0.402^ 02Horizontal Velocity Profile At T = 0.78 T1oor q ^OA-0.004^-0.002^0.000^0.002^0.004 -0.004^-0.002^0.000^0.002U Ws) U (in's)0.032^0.004-0.002^0.000^0.002U (m/s)Horizontal Velocity Profile At T = 1.56 Ti1.00.8• 0.6OA0 0 : 02 1--0.004^-0.002^0.003^0.002C__0.00t^-0.001^-0.002^0.000^0.032^0.0040.20.0U (m/s)U (m/s)Horizontal Velocity Profile At T = L82 TiHorizontal Velocity Profile At T = 1.69 Ti0..01.0^ 10.0 ^2^02 ^I ^0.88--) ^z 0.6z 0.6>' 04>• OA-7^0.20.00.000^.002U (m/s)^0.002^ -0404^-0.002^U (m/s)0-0.004^-0.032^0.0001.0:0.8• 0.6OA1.00.8 -z 0.60A -0.20.0 ^-0.004 -0.002^0.000^0.002U (m/s)Horizontal Velocity Profile At T = 1.43 Ti1.00.8z 0.60.40.20.0-0.001^-0.002^0.000^0.002UHorizontal Velocity Profile At T = 0.91 Ti Horizontal Velocity Profile At T = L04 Ti1.00.80.6OA0.20.0-0.004 4.032 0.000U (m/s)Horizontal Velocity Profile At T ss 1.17 Ti^ Horizontal Veloci Profile At T = 1..30 Ti1.00.8z 0.6>" OA0.20.0-0.00470.004Horizontal Velocity Profile At T = 1.95 Ti1.0OAz 0.6>' OA02 ::-OA^-0.004 -0.032^0.003U (m/s)Horizontal Velocity Profile At T = 2.08 Ti 1.00.8• 0.40.2OA^-0.004 -0.002^0.003^0.032U (m/s)0.002Figure (5-16). Forced internal seiche triggered by unsteady withdrawal, experiment 2,run 3; horizontal velocity profiles at x = 0.5 L.109Ee^0.150oIo'cil 0.145^I I^1^100 150 200Time (sec.)1300250200 250 300central point of lower layer- • - - central point of upper layer,—1-6x10"50oTime (sec.)Figure (5-17). Forced internal seiche triggered by unsteady withdrawal, exexperiment 2,run 3; (a) time history of interface elevation at x,',1, (b) time history of horizontal velocityat central points of both upper and lower layers, (c) forcing discharge process.110Figure (5-18). Forced internal seiche triggered by unsteady withdrawal, experiment2, run 3; surface plot of density interface movement of resonant internal seiche.111200^.^300Time (sec.)20E 10-20200^300Time (sec.)400100Run 4= 31.4 s.3200^00 Time (sec.) 4005005004002010-10-20100^200^300^400^500Time (sec.)Figure (5-19). Forced internal seiche triggered by unsteady withdrawal, experiment 3,interfacial oscillation data at x=0.(continued)112100 200^300^400Time (sec.)100 200^300^400Time (sec.)20 ■^a 101I Run 5= 33.0 S.200300Time (see)20.1-20200300Time (sec.)20g 40-2020—g 10..,-20I^500100 400o 100 400 500o 500o 500Figure (5-19). Forced internal seiche triggered by unsteady withdrawal, experiment 3,interfacial oscillation data at x=0.113... .. .... .....:-Ifil0.0^01^02^03^0.4Frequency (Hz)Frequency (Hz)Figure (5-20). Forced internal seiche triggered by unsteady withdrawal,experiment 3; spectra of interfacial oscillation data at x=0.(continued)11410."1042(c)10-8... ... ... ... .. ... ..... ...I Run 8Tf = 54.0 s...... ...1.... . ^0.1^0.2Frequency (Hz)031043 — 1111111111^VI0.0^0.1 0.2^03^0.4Frequency (Hz)0.4Figure (5-20). Forced internal seiche triggered by unsteady withdrawal,experiment 3; spectra of inwrfacial oscillation data at x=0.1150.4 ^020.0:-0.004^-0.002^0.000^0.002^0.004U (m/s)0.004 0.0040.002 0.002-0.004 -0.004 -0.002-0.0O2 0.000U (Ws)0.000U (m/s)Horizontal Velocity Profile At t =0.77 T11-0.002UHorizontal Velocity Profile At t 1.08 T11 Horizontal Velocity Profile At t = 1.39 T111.0080.6OA0.20.00.002 0.004 -0.004 -0.002 0.000U(Ws)Horizontal Velocity Profile At t = 2.01 T11-0.004 0.000U (rawHorizontal Veloci Profile At t = 1.70 T11-0.002 0.002 0.0041.00.80.6OA020.01.0.OS0.60.40.20.0 0.002.^0.004Z.1.00.80.6OA020.01.0OS0.6OA0.200-0.002^0.000U (mh)-0.002 0.0040.002OA00U (ai/s)Horizontal Velocity Profile At t = 2.63 T11-0.0040.022^0.004-0.002 0.000U (m/s)0.002 0.004Horizontal Velocity Profile At t =2.32 T111.00.804 -0.40.20.00.032^0.004-0.002^0.000U (is)Figure (5-21). Forced internal seiche triggered by unsteady withdrawal, experiment 3,run 5; horizontal velocity profiles at x3.5L.Horizontal Velocity Profile At t = 2.94 T11 Horizontal Velocity Profile At t = 3.25 T111.0020.6OA020.0Horizontal Velocity Profile At t =0.46 T11Z.116^1500100 200 300 400Experiment 4, Run 1, Tr = 10.0 s.center of interfacial layerTime (sec.)Time (sec.)175170165160155150145 o 100 200^300 400 500Time (sec.)Figure (5-22). Forced internal seiche triggered by unsteady withdrawal, experiment 4,run 1; interfacial oscillation data of (a) the center of interfacial layer and (b) upper andlower interfaces at x,',I, (c) time series of interfacial layer thickness at x=0.(continued)11710-111 02center of interfacial layerRun 1, Tf= 10.0$ANA. 7(d) upper interface1117-111-1111ttArv■Ann.„,.., I34.6 s63 $ 11833 $ifA10.0 s20.0 s^1Vt,,^10.0^ 0.1^0.2^0.3Frequency (Hz)Figure (5-224). Forced internal seiche triggered by unsteady withdrawal,experirrent 4, run 1; spectra of interfacial oscillation data at x=0.118300200 400Time (sec.)lower interface ^upper interface200^300 400 500center of interfacial layerExperiment 4, Run 2, Tr = 31.4 s.Time (sec.)Time (sec.)Figure (5-23). Forced internal seiche triggered by unsteady withdrawal, experiment 4,run 2; interfacial oscillation data of (a) the center of interfacial layer and (b) upper andlower interfaces at x',31, (c) time series of interfacial layer thickness at x=0.(continued)119.^IIIII^IIIIIIIII^II^II^1^I^f^II^1^1^I0.1 0.2 03icr4 —31.4s10-5 Th.563 s if106_ 833s 1104 —104 —0110-9 I10-1° —1011 —(d)io-'2 —1....,0.0Run 2, Tf= 31.4 scenter of interfacial layerFrequency (Hz)Figure (5-23d). Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 2; spectra of interfacial oscillation data at x=0.12016501604.4^155center of interfacial layer100 200 300 400 500Time sec.)lower interface10-10100^200 300Time (sec.)200^300Time (sec.)400 500Experiment 4, Run 3, Tf = 33.0 s.Figure (5-24). Forced internal seiche triggered by unsteady withdrawal, experiment 4,run 3; interfacial oscillation data of (a) the center of interfacial layer and (b) upper andlower interfaces at x4), (c) time series of interfacial layer thickness at x=0.(continued)121center of interfacial layer0.1^ 0.2Frequency (Hz)03Figure (5-24d). Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 3; spectra of interfacial oscillation data at x=0.122Experiment 4, Run 4, Tf = 34.6 s.30 acenter of interfacial layer2010-10-20(a)-RAL,o 300100 200—Time (sec.)-upper interface200^300109^I500400„lower interface_Ei.4ig1.g00Time (sec.)Time (sec.)Figure (5-25). Forced internal seiche triggered by unsteady withdrawal, experiment 4,run 4; interfacial oscillation data of (a) the center of interfacial layer and (b) upper andlower interfaces at x4, (c) time series of interfacial layer thickness at x=0.(continued)123io4 — 34.6sRun 4, T1= 34.6 s—1563 $lo-6 — 3le —1 0-m —(d)io-" --1..,111.-..I.,..1, -I-0.0^ 0.1i^1I^II^I^I0.2 0.3Frequency (Hz)Figure (5-25d). Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 4; spectra of interfacial oscillation data at x=0.124Experiment 4, Run 5, Tf = 80.0 s.105-5-10200^300^400^500Time (sec.)10508-5-10100 200^300Time (sec.)400 5001751701650160E-■155150145100 200^300 400 500Time (sec.)Figure (5-26). Forced internal seiche triggered by unsteady withdrawal, experiment 4,run 5; interfacial oscillation data of (a) the center of interfacial layer and (b) upper andlower interfaces at x4), (c) time series of interfacial layer thickness at x=0.(continued)1001250.1^ 02^ 03Frequency (Hz)Figure (5-26d). Forced internal seiche triggered by unsteady withdrawal,experiment 4, run 5; spectra of interfacial oscillation data at x=0.126Horizontal Velocity Profile At t = 0.24 T12 Horizontal Velocity Profile At t = 0.48 T12o01.00.80.60.40.20.01.00.80.6OA0.20.01.00.80.6OA020.01.00.80.60.4020.01.00.80.6OA020.01.00.80.60.40.20.01.00.80.60.402001.00.8z 0.6020.01.0020.60.2 E -^^70.000^0.002^0.004U (m/s)-0.002^0.000^0.002^0.004U (m/s)Horizontal Velocity Profile At t =0.96 T12OA0.2 ^-0.0041.00.80.6OAF.-Horizontal Velocig Profile At t = 0.72 T121.00.80.60.40.20.00.000^0.002^MONU 03440Horizontal Velocity Profile At t =1.20 T12-0.004^-0.002^0.000^0.002^0.004U (m/s)-0.004^-0.002^0.000^0.002^0.004U (m/s)Horizontal Velocity Profile At t =2.16 T121.00.80.60.40.20.0-0.004^-0.002^0.000^0.002^0.004U (m/s)0.000^0.002^0.004U (m/s)Horizontal Velocity Profile At t = 1.44 T12 -0.002^0.000^0.002^0.004UHorizontal Velocity Profile At t = 1.92 T120.002^0.004Horizontal Velocity Profile At t =2.40 T120.002^0.004-0.002Horizontal Velocity Profile At t = 1.68 T12-0.004 -0.002 0.000U (m/s)-0.002 0.000U (m/s)1.00.80.60.40.20.01.00.80.6OA0.20.0Horizontal Veloc' Profile At t =2.88 T12 -0JO04^4002^0.000^0.002^0.004Horizontal Velocity Profile At t = 3.36 T12-0.004^-0.002^0.000^0.032^0.004U (WN)Figure (5-27). Forced internal seiche triggered by unsteady withdrawal, experiment 4,run 5; horizontal velocity profiles at x(3.5L.127Horizontal Velocity Profile At t =2.64 T12-0.004^-0.002^0.000U (m/s)Horizontal Velocity Profile At t = 3.12 T120.002 0.004-0.002 0.000U (m/s)0.002 0.004
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical modeling of resonant internal seiche triggered...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical modeling of resonant internal seiche triggered by unsteady withdrawal Gu, Li 1993
pdf
Page Metadata
Item Metadata
Title | Numerical modeling of resonant internal seiche triggered by unsteady withdrawal |
Creator |
Gu, Li |
Date Issued | 1993 |
Description | Most lakes and reservoirs are temperature stratified for some period during the year, at which time a typical fresh water body can be divided into a well-mixed, warm, surface layer called the epilimnion, which is separated from the colder and relatively stagnant bottom water, called the hypolimnion, by an intermediate region with strong temperature gradient known as the metalimnion, which in conditions of temperature stratification only, is sometimes called the thermocline. Internal standing waves or internal temperature seiches within a stratified lake or reservoir play an important role in the physical and biological dynamics. Most previous studies have considered wind-generated internal seiches. Little attention has been given to the problem of triggering internal seiche by unsteady withdrawal, which was first proposed by Imberger (1980) as a possible pump-back scheme to optimize the water quality in a lake or reservoir. The resonant internal seiche triggered by unsteady withdrawal is studied numerically in this thesis using the C/C (Cloud-In-Cell) method. The investigation is concerned with internal seiches in a viscous, incompressible and weakly stratified fluid which is confined within a rectangular tank with a horizontal line sink in an end wall. Two-layered fluid with a sharp interface and two-layered fluid with a diffuse interface are investigated. Two-layered fluid with a sharp interface approximates conditions of lakes and reservoirs in late summer or early autumn when the thermocline is quite thin, while two-layered fluid with a diffuse interface approximates conditions of lakes and reservoirs in spring and late autumn, when the thermocline can be quite thick. The CIC method used in present study combines the best features of both Lagrangian and Eulerian approaches, and is most advantageously applied to fluid dynamics problems with density stratification. The numerical model solves, in finite difference form, the governing equations in terms of stream function and vorticity, and therefore avoids the complexity of solving the pressure field. The free surface boundary conditions in this case are somewhat difficult to apply. However, for the case of small density variation, simplified free surface boundary conditions may be chosen without affecting the main features of the flow. A two-layered fluid with a sharp interface supports only one vertical internal wave mode. With the fluids in the upper and lower layers moving in the opposite directions, large velocity shear develops at the density interface. The introduction of a continuously stratified interfacial region between two homogeneous fluids enables more than one vertical internal wave mode to exist, although the energy is usually partitioned into the first few internal wave modes. It is found that the natural frequency of internal waves decreases with increasing interface thickness and viscosity. The natural periods of a two-layered fluid system with an exponentially stratified interfacial layer are well predicted by the numerical model, in comparison with the theoretical solutions of Fieldstad (1933). By discharging fluid from the tank, two distinct classes of internal waves are created, namely, forced and free internal waves. The forced internal wave follows the forcing discharge and sustains as long as the discharge is retained, while the free internal wave oscillates according to its natural frequency with decaying amplitude due to viscous damping. Periodic beating is observed as a result of the interaction between the forced and free internal waves. The resonant internal seiche is created when the frequency of unsteady forcing discharge corresponds to the natural frequency of internal waves. |
Extent | 6076626 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-09-10 |
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. |
IsShownAt | 10.14288/1.0050454 |
URI | http://hdl.handle.net/2429/1776 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1993-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1993_fall_gu_li.pdf [ 5.8MB ]
- Metadata
- JSON: 831-1.0050454.json
- JSON-LD: 831-1.0050454-ld.json
- RDF/XML (Pretty): 831-1.0050454-rdf.xml
- RDF/JSON: 831-1.0050454-rdf.json
- Turtle: 831-1.0050454-turtle.txt
- N-Triples: 831-1.0050454-rdf-ntriples.txt
- Original Record: 831-1.0050454-source.json
- Full Text
- 831-1.0050454-fulltext.txt
- Citation
- 831-1.0050454.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-0050454/manifest