A Finite Volume Implementation of the Shallow Water Equations for Boussinesq Gravity Currents by Thomas Meyer Hatcher Jr. A thesis submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Master of Science Auburn, Alabama August 4, 2012 Keywords: Gravity Currents, Lock-exchange, Finite Volume Method, Shallow Water Equations Copyright 2012 by Thomas Meyer Hatcher Jr. Approved by Jos?e Goes Vasconcelos Neto, Chair, Assistant Professor of Civil Engineering Prabakhar Clement, Professor of Civil Engineering Xing Fang, Associate Professor of Civil Engineering Abstract The shallow water equations (SWE) are a powerful tool for modeling the propagation of grav- ity currents (GC) because of their relative simplicity, computational efficiency and accuracy. Finite difference solutions, either based on the method of characteristics (MOC) or the implementation of numerical schemes such as Lax-Wendroff (LxW) have been traditionally used in such flow com- putations. On the other hand, the finite volume method (FVM) has been gaining popularity in several other hydraulic applications, being favored in cases when flow discontinuities are antic- ipated. This work is focused on an implementation of the finite volume method (FVM) to the solution of Boussinesq GC using the one and two-layer SWE models. The proposed two-layer mathematical model is a modification of the the work by [Rottman and Simpson, 1983], adapted to express such equations in a vectorial conservative format, amenable for FVM implementation. The traditional solution for the GC front boundary condition (BC), using a characteristic equation and a front condition, is compared to a new formulation that explicitly enforces local mass and momentum conservation. Linear numerical schemes (LxW and FORCE) and non-linear schemes based on the approximate solution of the Riemann problem (Roe and HLL) are implemented in this framework along with various front conditions. The proposed modeling framework is tested against experimental data collected by this investigation, and is also compared to previous inves- tigations. Results indicate that this proposed model has a comparably simple and robust imple- mentation, being flexible enough to be applied in a wide range of GC flow conditions and presents good accuracy. ii Acknowledgments I would like to thank my advisor, Professor Jose Goes Vasconcelos, for his experience, hard work and friendship throughout my time here at Auburn. He introduced me to CFD and experi- mental hydraulics and was instrumental in the development of this thesis. I also received a great deal of help from fellow students throughout the experimental portion of this work. Parker Ross assisted in the construction of the scale model of Mobile Bay. Kyle Moynihan and Kyle Strachan assisted in a large number of experiments along with Paul Simmons and Carmen Chosie. There contributions and friendship made the time I spent in the lab enjoyable. Bernardo Trindade and Gabriel Leite helped introduce me to numerical modeling and their friendship played an important role throughout this process. I would like to thank my family for all of their support these last couple of years and through- out my life. My best friend and wife, Jennifer Kelly Hatcher, has been the centerpiece in my life these last few years and has provided essential love and support. I would like to thank her for putting up with me throughout my time at Auburn. My parents, Tom and Sharon Hatcher, have provided incredible support throughout my life. I could not have asked for better role models to grow up around. I would like to congratulate my brother, Preston Hatcher, for posting a 1.73 ERA with 10 saves this past spring and to thank him for teaching me self control as a child. Finally, I would like to thank all of my friends in the civil engineering department who helped make my time at Auburn enjoyable. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Description of GC flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 Laboratory experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3 GC Flow stages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.4 Modeling approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.1 The description of the front . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.1.1 Theoretical investigations of GC front conditions . . . . . . . . . . . . . . 18 2.1.2 Empirical investigations of GC front conditions . . . . . . . . . . . . . . . 21 2.2 Experimental investigations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3 Analytical description of GC flows with integral models . . . . . . . . . . . . . . . 24 2.3.1 Deep ambient with constant Froude number . . . . . . . . . . . . . . . . . 24 2.3.2 Generalized solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.4 Numerical modeling of gravity currents . . . . . . . . . . . . . . . . . . . . . . . 27 2.4.1 One-Layer SWE model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.4.2 Two-Layer SWE model . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4.3 Navier-Stokes models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.5 Hyperbolic numerical schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.5.1 Linear schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 iv 2.5.2 Nonlinear schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3 Knowledge Gap and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.1 Experimental program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.1.1 Physical apparatus construction . . . . . . . . . . . . . . . . . . . . . . . 51 4.1.2 Initial salinity measurements . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.1.3 Front trajectory measurements . . . . . . . . . . . . . . . . . . . . . . . . 56 4.1.4 MicroADV measurements . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.2 One-layer SWE model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.3 Two-layer SWE model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.4 Aspects of the numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.4.1 Numerical schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.4.2 Boundary Conditions: DC and MOC . . . . . . . . . . . . . . . . . . . . 64 4.4.3 Stability and time step calculation . . . . . . . . . . . . . . . . . . . . . . 66 5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.1 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.2 Numerical model results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3 Comparison between experimental results and numerical predictions . . . . . . . . 76 5.4 Comparison with previous models . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6 Conclusions and future developments . . . . . . . . . . . . . . . . . . . . . . . . . . 97 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 v List of Figures 1.1 Aerial view of the front that forms between the Connecticut river and the Long Island Sound estuary. Source: environmentalheadlines.com . . . . . . . . . . . . . . . . . . 2 1.2 A particle-driven GC generated from the terrorist attacks of 9/11. Source: 911re- search.wtc7.net. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 A haboob or dust cloud smothering a residential area, my3monsters.com. . . . . . . . 5 1.4 High Reynolds number lock-exchange GC for Boussinesq fluids along a horizontal surface presented in [Slim, 2006]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.5 Schematic diagram of a partial depth lock-exchange experiment in a rectangular chanel. Adapted from [Rottman and Simpson, 1983]. . . . . . . . . . . . . . . . . . . . . . . 10 2.1 Schematic diagram of Benjamin?s formulation for the case of a heavier fluid propagat- ing underneath a lighter ambient fluid. Adapted from [Klemp et al., 1994]. . . . . . . . 20 2.2 Schematic diagram of a GC in an initially moving ambient system. Adapted from [Wright and Paez-Rivadeneira, 1996]. . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3 Schematic diagram for both fronts that is generated from a GC flowing into a shallow ambient (f > 0:5). Adapted from [Ungarish and Zemach, 2005]. . . . . . . . . . . . . 36 2.4 Two-wave structure for the HLL scheme presented in [Toro, 2001]. . . . . . . . . . . 45 4.1 The experimental apparatus with stationary ADV probes and a video camera moving with the GC front. The initial GC depth h0 is equal to the ambient depth H. . . . . . . 52 4.2 The photograph on the left is an aerial view of Mobile Bay, Alabama. Source: Landsat program. To the right is a snapshot of the depth-averaged velocities (arrows) and the water levels (contours) predicted by the estuarine and coastal ocean model. The vertical line in the middle of the bay represents the shipping canal. Presented in [Chen et al., 2005]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.3 Different components of the construction for the acrylic channel. . . . . . . . . . . . . 54 4.4 Photograph of the leading edge of a GC that illustrates the procedure in which the front trajectories were measured. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 vi 4.5 Schematic diagram of a GC that illustrates the cell structure for the DC BC at a given time step. The cell LE represents the last computational node that is completely filled by the front of the GC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.1 Dimensionless trajectory and velocity for the GC front. . . . . . . . . . . . . . . . . . 70 5.2 Velocity hydrograph comparison for 1%, 2%, and 3% salt in the 9.14 m tank. . . . . . 72 5.3 GC propagation for the one and two-layer SWE over a wide range of fractional depth: f = 0:01, 0:5 and 1:0. The DC BC is used with 400 discretization cells. . . . . . . . . 75 5.4 Front trajectory results for the one-layer SWE model using Roe scheme and DC BC with various front conditions and from experiments. . . . . . . . . . . . . . . . . . . . 77 5.5 Front trajectory results for the two-layer SWE model using Roe scheme and DC boundary conditions with various front conditions and from experiments. . . . . . . . 78 5.6 Slumping for the one and two-layer SWE models using the HLL and Roe scheme with two front conditions: Benjamin (B) and Huppert Simpson (HS). x0 = 1, h0 = 1 and e = 2% using 800 computational cells. . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.7 Velocity hydrographs with B and HS front conditions and 200 computational cells. The experimental data (e = 2:0%) is provided from MicroADV devices in the slump- ing stage: x = 7:68. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.8 Velocity hydrographs with the MRS front condition and 200 computational cells. The experimental data is provided from MicroADV devices at x = 7:68. . . . . . . . . . . 89 5.9 Comparison between measured ambient velocities and two-layer SWE predictions us- ing B and HS front conditions with 200 computational cells. The experimental data (e = 3:0%) is provided from MicroADV devices at x = 7:70. . . . . . . . . . . . . . 90 5.10 Comparison of GC profiles at various time steps between the proposed two-layer model and [Rottman and Simpson, 1983] model for the critical condition: f = 0:5. . . 92 5.11 Dimensionless slumping comparison between the proposed model and [Ungarish and Zemach, 2005] model using Roe scheme and modified FORCE scheme. . . . . . . . . 93 5.12 Dimensionless front velocity results for the two-layer model over a large range of fractional depth. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.13 Dimensional slumping comparison between the SWE and EFDC for a) t=1.8sec and b) t=39sec. The discretization size for the SWE models is 800 cells. . . . . . . . . . . 96 vii List of Tables 5.1 Experimental initial conditions and results for each experiment. Values for um were determined from a digital camera tracking the GC front. . . . . . . . . . . . . . . . . 69 5.2 Sensitivity of each modeling approach to the number of discretization nodes. The error for each model is compared to its computed propagation time at 1,000 nodes. The Roe scheme is used with B front condition and kDt = 0:1. . . . . . . . . . . . . . 74 5.3 Comparison between lock-exchange experiments and one-layer SWE model predic- tions using both BC approaches at the GC front, kDt = 0:1 and a 400-cell solution do- main. The time of propagation for the GC in the experiments was 93.67 s (e = 1:0%) and 67.63 s (e = 2:0%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.4 Comparison between lock-exchange experiments and two-layer SWE model predic- tions for different time steps and front conditions using the DC boundary condition. The GC time of propagation in the experiments was 93.67 s (e = 1:0%) and 67.63 s (e = 2:0%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.5 Two-layer SWE results using various values for b for the RS front condition with 400- cell solution domain. When b = 1:215, this front condition is referred to as MRS. The GC time of propagation in the experiments was 93.67 s (e = 1:0%) and 67.63 s (e = 2:0%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.6 Computational time comparison between the one-layer (1L) and two-layer (2L) SWE models and between the DC and MOC front BC strategies for various discretization sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 viii List of Notations The following symbols are used in this thesis: x0 = Initial length of denser fluid h0 = Initial depth of denser fluid h; h1 = Depth of GC h2 = Depth of ambient fluid hLE = Depth of leading edge or front of GC H = Initial ambient depth, total depth of system f = Initial depth ratio or fractional depth = h0=H fLE = Depth ratio of the GC leading edge u;u1 = Velocity of GC u2 = Velocity of the ambient fluid uLE = Velocity of the GC leading edge rc;r1 = Density of bottom GC ra;r2 = Density of ambient fluid Dr = Difference between the current and ambient densities =rc ra e = Relative density difference =(rc ra)=rc g0 = Reduced gravitational acceleration = ge um = Mean velocity at the GC front Re = Reynold?s number Rem = Mean Reynold?s number FrH = Total depth densimetric Froude number FrLE = Froude number at the GC leading edge b = Scaling factor for the front condition ix ~U = Vector of conserved variables ~F(U) = Vector of conserved variables fluxes ~S(U) = Vector of source terms A;B = Source term parameters Dx = Dimension of space discretization Dt = Dimension of time discretization kDt = Scaling factor for time step calculation Cr = Courant number MOC = Method of Characteristics DC = Dual-Cell boundary condition x Chapter 1 Introduction Gravity currents (GC) (i.e. buoyancy or density currents) are driven by a density difference between two or more fluids and generally travel in a quasi horizontal plane. In a large number of cases the difference in density is small so that the buoyancy forces are of equal or greater magnitude than the inertial forces. The primary discrepancies in the density are caused by two factors: the accumulation of dissolved solids or temperature differences between fluids. Examples of GC in which an accumulation of solids causes an increase in density include exchange flows in estuaries, sediment entrainment in rivers and the expansion of dense volcanic dust clouds. When a cold stream enters a warmer lake or vice versa, the density current is driven by temperature differences. In this thesis high Reynold?s number GC flows, which are formed from the discharge of saltwater into freshwater, are analyzed in the context of the shallow water equations (SWE). In water bodies such as rivers and estuaries, GC flows are a common occurrence and have a substantial impact on the flow field. When a river with a high discharge rate empties into an estuary with negligible tidal effects, a lighter GC may form at the free surface with a salt wedge (denser GC) propagating on the channel bottom into the river (Figure 1.1). In the Connecticut River during the spring months, there is a large discharge due to snow melt, which forms a salt-wedge regime with the Long Island Sound estuary (see Figure 1.1). Detailed measurements of the density profile have been recorded for the lighter GC that forms at the water surface [Simpson, 1997]. Sea-lochs are deep inlets from the ocean in which tidal and river effects are important. In systems with large amounts of river flow, well-mixed layers of brackish water develop on the river side of a sill or obstacle. If the tide rises high enough, a salt wedge may rise over this sill into the brackish region. This process has important effects on the pollution of the inlet and may determine the amount of dissolved oxygen in an unmixed zone that lies underneath the brackish 1 Figure 1.1: Aerial view of the front that forms between the Connecticut river and the Long Island Sound estuary. Source: environmentalheadlines.com region. Overturn due to dense water flowing over the sill may occur at rare intervals. In Loch Eil in Scotland, the overturn occurs at a frequency of about once a year, and this occurrence is caused by a GC flow [Simpson, 1997]. One of the first observed GC flows was related to volcanic eruptions in which hot gas and ash are expelled at high velocities. This ash cloud can travel hundreds of kilometers and is the most destructive component of the eruption [Baxter, 2000]. Once the buoyancy forces outweigh the initial momentum, the fluid rises vertically as a buoyant plume, which is outside the scope of quasi horizontal GC flows [Turner, 1973]. As the ash cloud moves higher and higher into the atmosphere, entrainment becomes more and more important [Slim, 2006]. Thus, the reduced gravity is decreased in magnitude until neutral buoyancy is achieved, which is inevitable because of the ambient stratification. At this depth, the plume begins to move horizontally as an intrusive GC. 2 Figure 1.2: A particle-driven GC generated from the terrorist attacks of 9/11. Source: 911re- search.wtc7.net. Similar examples occur at the ocean floor where 75% of the annual volcanic activity occurs. Small- scale emissions through hydrothermal vents and large-scale eruptions rise as buoyant plumes and have an important impact on ocean currents [Simpson, 1997]. Instead of forming a buoyant plume, a volcanic eruption can also spread as a pyroclastic GC when the density of the fluid is sufficiently large. Moreover, when the ash cloud is too dense to continue rising, the flow collapses and forms a ?ground-hugging? surge [Slim, 2006]. These pyroclastic flows can travel at speeds of up to 160 ms 1 with temperatures reaching 900 C [Blong, 2000]. Since the late 1700s, the primary cause of death in volcanic eruptions is due to pyroclastic flows (i.e. approximately 27% or 220,000 deaths) [Baxter, 2000]. A similar type of GC flow occured due to the terrorist attacks on 9/11/2001 (see Figure 1.2). A number of other types of GC flows can also occur in the context of volcanoes. When heavy rains re-mobilize the unconsolidated ash deposited from a volcanic eruption, a lahar or mud flow 3 is formed. This hazardous type of GC can travel at high velocities for several kilometers and can recur years after the initial eruption [Blong, 2000] and [Slim, 2006]. In 1985, the volcanic eruption of Nevado del Ruiz in Columbia claimed 25,000 lives. Four major lahars were formed, which traveled a total distance of 60 km at 10 ms 1. These mud flows over topped the banks a river channel in Armero and followed the bank into the town, which resulted in the majority of the casualties for this event [Simpson, 1997]. The GC examples presented above are higher velocity and lower viscosity flows. Thus, the Reynolds number is large (Re > 1;000). However, the flow of molten rock is an example of a low velocity and high viscosity flow. In addition, these low Reynolds number flows (Re < 10) are non-Newtonian meaning that there is a nonlinear relationship between stress and strain. However, the focus of this thesis is on the simpler Newtonian fluids, which are valid for gases and water flows. In addition, the buoyancy forces are balanced by inertial forces instead of by viscosity so that viscous forces are neglected. Another common atmospheric GC is a haboob or dust cloud, which is a cold atmospheric front several hundred meters in depth that can travel at over 20 ms 1. These types of GC flows are generated from thunderstorm outflows over sandy and arid landscapes. The turbulence from the outflow entrains dust, which has been measured at 40 mg=m3 [Simpson, 1997]. However, the difference in density due to the entrained dust particles is generally a small fraction of 1% compared to the effects of temperature. Therefore, it is common to model these haboobs only based on the temperature differences between the two fluids [Simpson, 1997]. GC flows are often categorized depending on the description of the density difference. When this difference is caused by temperature or dissolved solids, the GC is referred to as homogeneous. However, sometimes GC flows are driven by suspended particles that progressively settle out of the GC, and the density of the current is reduced. These flows are categorized as particle-driven whereas a combination of these two is known as particle-laden. A haboob is an example of a particle-laden GC; however, this type of flow is more homogeneous dependent. On the other 4 Figure 1.3: A haboob or dust cloud smothering a residential area, my3monsters.com. hand, a buoyant plume that occurs due to volcanic eruptions, ocean outfalls, etc. eventually move horizontally as a GC once neutral buoyancy is reached and is known as an intrusion [Slim, 2006]. In this thesis, the focus is on turbulent, homogeneous GC flows that travel along a horizontal surface in a Boussinesq system (i.e. ra=rc 1 where ra is the ambient density and ra is the current density). The thin layer assumption, which is widely used in the context of GC, is utilized in order to validate the SWE. Two thin-layer numerical models are proposed: one-layer SWE (the ambient fluid is neglected in the governing equations) and two-layer SWE (the governing equations are formulated based on the current and the ambient fluids). In addition, lock-exchange experiments are conducted to compare with the numerical models and to measure the GC velocity structure. 1.1 Description of GC flows GC flows are driven by a density difference in the vertical direction, which results in a pressure difference in the horizontal plane. This pressure difference is balanced by the dynamic velocity field in the horizontal direction. The gravitational effect associated with the density difference is 5 referred to as the reduced acceleration due to gravity, g0: [Ungarish, 2009]. g0= r c ra rc g =eg (1.1) in which ra is the ambient fluid density, rc is the density of the current and e is the relative density difference. The reduced gravity is a fundamental feature in the GC propagation and causes the decrease in velocity for GC compared to traditional dam break flows. This parameter replaces the gravity (g) in the SWE and plays a role in the front condition (described in 2.1). The first quantitative study of GC flows was done by von Karman, who evaluated the spread of poisonous gas [Huppert, 2006]. Subsequently, [von Karman, 1940] formulated a classical rela- tionship for the GC depth and velocity using Bernoulli?s theorem in which a heavier fluid advances into the lighter atmosphere: u LE (g0h LE)1=2 = Fr LE (1.2) in which hLE and uLE are the depth and velocity of the front or leading edge of the GC, respectively. FrLE is the Froude number of the GC front, which was evaluated by von Karman to equalp2 in the context of a relatively deep ambient. [Benjamin, 1968] argued that von Karman?s formulation was invalid because he applied Bernoulli?s theorem across a streamline characterized with head losses. Subsequently, Benjamin re-derived this expression (1.2) using the momentum integral, and ended up with the same results. This result surprised Benjamin; however, because both approaches used a different integral of the Euler equation, they could not have reached different results [Huppert, 2006]. The Froude number for a GC flowing into a deep ambient is actually closer to unity. More- over, [Huppert and Simpson, 1980] obtained FrLE from their experiments and [Bonnecase et al., 1992] stated that viscous drag and Reynold?s stresses are the reason for this discrepancy. [Shin et al., 2004] formulated a new expression for the GC front in which FrLE = 1 for the deep ambient scenario. Their result differs from Benjamin and von Karman?s value because [Shin et al., 2004] approach includes the energy transfer between the GC and the ambient fluid, which is important 6 for f!0. In the experiments conducted in this thesis f = 1; however, the numerical model is able to simulate any initial condition for the fractional depth (0 f 1). Two problems emerge from equation 1.2. For example, when a GC propagates into a shallow ambient fluid, the Froude number also becomes dependent on the ambient depth H or in effect, the depth ratio or fractional depth at the leading edge of the GC (i.e. fLE = hLE=H). In addition, (1.2) represents one equation for two unknowns (hLE and uLE) [Huppert, 2006] (further discussed in subsection 4.4.2). [Benjamin, 1968] solved the first potential problem by deriving an expression for FrLE that is a function of fLE by applying a flow force balance (section 2.1). For the lock-exchange flow that is simulated in this thesis, two approaches are utilized at the front BC in conjunction with Benjamin?s condition or an empirical alternative. The leading edge or front of a GC forms a distinct region that separates the advancing GC and ambient fluids. Mixing is prevalent at the GC front, but the entrainment between the two fluids is swept upstream in the form of shear instabilities. Towards the front of the GC a head or nose generally forms in which the flow depth is greater than that of the fluid upstream. This nose region is characterized by ??a zone of breaking waves and intense mixing?? [Simpson, 1997]. Initially, the nose of a GC moving along a horizontal surface is quasi-steady; however, in the case of an incline, the size of the GC nose increases with slope [Simpson, 1997]. For an inviscid GC flowing on a rigid and horizontal surface (Figure 1.4), the depth of the GC nose is dependent on the ambient depth, which varies from the deep ambient f = 0 to the full depth lock-exchange problem f = 1. In addition, the shape of the nose is dependent on the turbulence in the ambient fluid and the ambient velocity structure. This study is not concerned with viscous effects, but if viscosity becomes important, the rate of advance for the GC is decreased and the turbulence at the nose is lessened. For high velocity and low viscosity GC (i.e. Re > 1;000), the Reynolds number does not play a significant role in the propagation. For instance, in a thunderstorm the Reynold?s number is typically in the order of 1;000;000. Moreover, field observations suggest that the flow structure in these highly turbulent examples is similar for all Re > 1;000 GC flows [Simpson, 1997]. Viscosity can play an important 7 Figure 1.4: High Reynolds number lock-exchange GC for Boussinesq fluids along a horizontal surface presented in [Slim, 2006]. role for flows in which Re < 1;000; moreover, when Re < 10, viscosity tends to dominant over inertia in the GC flow [Simpson, 1997]. When a viscous oil slick spreads over top a water body, surface tension eventually exceeds viscosity in the level of importance [Fannelop and Waldman, 1971]. For high Reynold?s number GC flows, turbulent mixing is prevalent, and the instabilities travel upstream of the GC nose. This turbulent mixing feature is illustrated in Figure 1.4 in which most of the information is presented from [Britter and Simpson, 1978] and [Simpson and Britter, 1979]. The mixing in GC flows is a complex process that usually occurs in the following two forms: billows that roll up the nose of the GC and a pattern of lobes and clefts. The billows are formed from shear instabilities at the interface of the fluids because the two fluids flow in opposite directions. The lobes and clefts are formed by the influence of turbulence against the ground [Simpson, 1997]. Entrainment, the mixing between fluids, plays an important role in the propagation of the GC, and the importance of entrainment increases with slope. Moreover, the billows that entrain the ambient fluid play an important role in mass transfer, which decreases the front velocity. The velocity of the GC front is one of the key parameters and has received a great deal of attention from researchers. This value can be approximated by equating the mean potential energy 8 loss to the kinetic energy gain [Simpson, 1997]. For a lock exchange flow where the height of the lock equals h0, this relationship becomes: 1 2g 0h0 = 1 2u 2 (1.3) Isolating the velocity (u) in (1.3) yields: U = p g0h0 (1.4) in which U is an estimated potential velocity of the current and h0 is the initial GC depth. Although the velocity estimation in (1.4) is an unrealistic simplification, the results are typically of the same order of magnitude [Ungarish, 2009]. It is customary to normalize the results for the velocity with this velocity estimation. As discussed above, one of the key parameters in GC flows is the Reynolds number, which governs the importance of viscosity. An estimate of the magnitude for the Reynold?s number can be determined from equation 1.4 [Simpson, 1997]: Re =Uh0=n (1.5) in which n is the kinematic viscosity. For this study, we are concerned with high Reynolds number flows (i.e. Re >> 1;000), and the final viscous stage that may eventually develop is neglected [Simpson, 1997]. In addition, the present work is focused on thin layer systems (h0=x0 < 1) so that the SWE are more applicable because there is less interfacial mixing [Ungarish, 2009]. 1.2 Laboratory experiments Experiments are conducted to validate the proposed numerical model using a MicroADV device as well as digital cameras. One of the advantages of experimental research is that the GC flows are similar to large-scale environmental examples. Moreover, experiments are typically 9 easier to conduct than field measurements and it is easier to control parameters such as the effects of the ambient fluid velocity. Experimental investigations from [Britter and Simpson, 1978] and [Simpson and Britter, 1979] provided important information on the GC structure especially at the nose that forms near the leading edge. One of the first GC experiments (the lock-exchange problem) was conducted by removing a vertical barrier that separates freshwater from saltwater [Keulegan, 1957] and [Barr, 1967]. A similar approach was conducted in this study in which high definition digital cameras tracked the GC front while MicroADV probes measured the internal GC velocity distribution. Constant volume lock-exchange experiments are typically conducted by separating a denser fluid with a density r+Dr, length x0 and depth h0 from a lighter ambient fluid r of depth H in a long rectangular tank (see Figure 1.5). For full depth lock-exchange flows, h0 = H where as partial depth experiments are conducted when h0 < H. For partial depth releases, the ambient fluid of depth H h0 is added on top of the denser fluid so that the depths on both sides of the gate are equal (for a more comprehensive discussion, see [Shin, 2001]). Once this vertical gate is removed, a GC flow begins in which the denser fluid propagates underneath the lighter ambient. When h0 equals H, a GC flow also develops at the free surface moving in the opposite direction of the denser bottom current. Once the upstream moving depression wave (or hydraulic jump) that develops for the denser fluid reaches an upstream physical boundary, a reflection occurs, and a nose region is developed towards the front of the GC [Rottman and Simpson, 1983]. Experiments are often preferred to large-scale field measurements in the context of GC flows [Shin, 2001]. Laboratory experiments are easy to repeat in order to establish consistency with the data collection. They are less expensive and can be performed in a more controlled environment so that the data collection is of higher quality than for field trials. In addition, it is easier to isolate a certain aspect of the flow for a better understanding [Shin, 2001]. For field measurements, it is often difficult to eliminate the effects of wind, boating and tidal effects in estuaries. 10 Figure 1.5: Schematic diagram of a partial depth lock-exchange experiment in a rectangular chanel. Adapted from [Rottman and Simpson, 1983]. The differences between laboratory experiments and environmental flows can occur because of the difference in length scales. In order to achieve dynamic similarity between these two meth- ods, the relative importance of the controlling forces (e.g. inertial, viscous, etc.) should be main- tained [Shin, 2001]. The dimensionless Reynold?s number should be consistent for the laboratory experiments. UL n Prototype = ULn Experiment (1.6) where U is the velocity scale and n is the kinematic viscosity. The length scale L within the Reynold?s number is typically smaller in the laboratory; however, the problem is reduced if water is used instead of air for atmospheric flows [Shin, 2001]. Moreover, the kinematic viscosity of water is 0.01 cm2s 1 at 20 C and 0.15 cm2s 1 for air [Simpson, 1997]. Another important factor is that GC flows are weakly dependent on the Reynold?s number when the Reynold?s number is greater than 1,000 as previously described. The Peclet number is another relevant dimensionless parameter that is related to mass trans- port phenomena. Pe = ULk (1.7) 11 in which k is the diffusivity of the fluid. The Peclet number is the ratio of the advective and diffusive terms in the momentum equation [Shin, 2001]. In the majority of GC applications, the order of magnitude of the Peclet number is very large, sometimes greater than 1 million. Therefore, diffusion effects are negligible [Shin, 2001]. In laboratory experiments, which are presented in this thesis, the diffusivity of pure salt is 1.1 x 10 5cm2s 1 so that the Peclet number was also greater than 1 million [Shin, 2001]. 1.3 GC Flow stages The inviscid GC propagation in a rectangular channel is described by three stages [Ungarish, 2009]: slumping, transition, and self-similar. The transition stage is the least interesting of the three stages because it more or less serves as a bridge between the other two. However, the beginning of the transition stage (i.e. the end of the slumping stage) is an important feature that effects the GC depth and velocity (further discussed in 5.2). At the initial release that marks the beginning of the slumping stage, the depth decreases dramatically in comparison with the other stages. Moreover, the shape of the interface and the velocity of the GC undergo a transformation. Although the depth (h) drops to at least half of the initial depth (h0) after the initial release due to the energy conserving theory of [Benjamin, 1968], the velocity and the depth of the GC nose remain almost constant for the remainder of the slumping stage. [Ungarish, 2009] describes the slumping stage in two sub-categories starting with a motion analogous to a dam-break wave, which is characterized by the initial drop in depth and increase in velocity for the GC. This type of motion occurs at the beginning of the lock exchange flow when stationary fluid remains in the system (i.e. upstream of the release). Towards the front of the GC ?the interface has a negative slope, corresponding to a rarefaction (expansion) pressure wave, under which the fluid is accelerated? [Ungarish, 2009]. Next, the entire depth of the current begins to descend abruptly sending a disturbance towards the nose. A positive slope for this disturbance begins to form, which travels towards the front of the GC [Ungarish, 2009]. 12 The second stage is transitional in nature in which the GC ?is under the influence of the initial conditions and the front boundary conditions? [Ungarish, 2009]. Once the downstream propagating bore or wave (depending on f) that originated at a physical boundary reaches the GC front, the depth and velocity of the front begin to decrease. Moreover, the quasi steady motion that forms during the slumping stage is ended and there is an unsteady transition to the self-similar stage. The importance of the initial conditions begins to decrease as the GC propagates further downstream and enters this self-similar stage. In the final stage (i.e. self-similar), the inertial forces begin to balance with the buoyancy forces to retain an almost steady flow. For the inviscid assumption, the GC flow remains in the self- similar stage until there is contact with a physical boundary (wall) or the simulation is terminated after a certain minimum GC thickness (i.e. h!0) otherwise viscosity would have to be included in the computations. Because the depth and the velocity decrease with time, ??the importance of the viscous friction relative to the inertial terms increases monotonically during the self-similar propagation?? [Ungarish, 2009]. The following relationship formulated in [Huppert and Simpson, 1980] describes the time when the inviscid assumption becomes invalid: x =(x50h50gr=n2)17 (1.8) According to this relationship, the inviscid assumption is valid for all of the experiments and nu- merical simulations conducted in this study, 1.4 Modeling approaches One of the first applications in which GC were studied was oil spreading over top of a water body. Because of the relatively small thickness of oil compared to water in the context of oil spills within the ocean, the solution of the governing equations can be simplified to analytical solutions. First, it is customary to reduce the Navier-Stokes equations into the depth-averaged SWE, which are referenced throughout this paper. 13 [Fannelop and Waldman, 1971] and [Hoult, 1972] used the SWE in order to derive self-similar solutions. They focused on the three major phases in the propagation of oil on a water body: gravity-inertial, gravity-viscous, and surface tension-viscous. [Fannelop and Waldman, 1971] as- sumed that the SWE are time-independent, which is required in order to derive these analytical solutions, so that time derivatives vanish in the governing equations. In order to derive the SWE in the context of oil spreading using a single pair of governing equations (one-layer SWE), some simplifications are made. For example, the viscosity in the mo- mentum equation is neglected because the viscosity of crude oil is one to two orders of magnitude greater than water, so the horizontal spread of oil is considered independent of the vertical direc- tion. However, if the Reynold?s is greater than 1,000, viscosity is typically neglected for saltwater intrusions as well. The following derivation follows the formulation in [Fannelop and Waldman, 1971]; however, subtle changes are made in the notation for consistency with this thesis, which focuses on a heavier GC propagating underneath a lighter ambient fluid (e.g. saltwater into fresh- water). [Fannelop and Waldman, 1971] derives the 1-D one-layer SWE from the 2-D Navier-Stokes equations with the incompressible and Boussinesq assumptions: ?u ?x + ?v ?y + ju=x = 0 (1.9) ?u ?t +u ?u ?x +v ?u ?y = 1 r ?P ?x + m r ?2u ?y2 (1.10) ?v ?t +u ?v ?x +v ?v ?y = 1 r ?P ?y g (1.11) in which u and v represent the longitudinal and vertical velocity components, respectively. P is the pressure, r is the density and m is the shear stress. This set of equations (1.9-1.11) is valid for both 2-D and axisymmetric GC flows. For the planar case, j = 0, and for the radial case, j = 1. In order to derive the 1-D SWE from these 2-D Navier-Stokes equations, the thin layer assumption is utilized. Thus, equations (1.9-1.11) are integrated across the GC layer depth (h), 14 and v?u?y becomes zero in the x-momentum equation. In addition, the vertical accelerations are neglected (i.e. v ?v?y = 0). After integration, (1.11) is reduced to the hydrostatic relation between the pressure in the layer and the depth. Then, an expression is derived by differentiating this pressure relation with respect to x and assuming the weight of the GC is equal the weight of the displaced ambient fluid. ?P ?x =rg 0?h ?x (1.12) in which g0 = (r1 r2)=r1 where r1 is the density of the heavier current and r2 is the density of the ambient fluid. Equation 1.12 is substituted into equation 1.10 in order to obtain a single momentum equation. This thesis is concerned with inviscid GC, so the last term in the RHS of the x-momentum equation (1.10) is omitted and the resulting expression for momentum is: ?u ?t +u ?u ?x +g 0?h ?x = 0 (1.13) Following the formulation in [Fannelop and Waldman, 1971], equation (1.13) is not sufficient to solve for the two unknowns (h and u). Therefore, the continuity equation is required alongside of this momentum equation. In order to formulate the updated continuity equation from (1.9), the vertical accelerations are neglected. Then, the difference in the velocity (Dv) is equated to the continuity equation integrated along the depth: Dv = Z h 0 ?u ?x + ju=x dh = h?u?x juh=x (1.14) in which: Dv = ?h?t +u?h?x Therefore, the resulting continuity equation is written as: ?h ?t +u ?h ?x +h ?u ?x + ju=x = 0 (1.15) 15 Together, equations (1.13) and (1.15) form the 1-D SWE, which are written in primitive for- mat. This set of equations is amenable to the planar and axisymmetric scenarios, which are widely studied for both analytical models and SWE models. However, this thesis is not concerned with axisymmetric GC, so the radial term is omitted in the continuity equation. Moreover, the chain rule is utilized in order to write the resulting SWE in conservative format (i.e. in terms of h and uh), which is preferred near shocks: [LeVeque, 1992]: ?h ?t + ?uh ?x = 0 (1.16) ?uh ?t + ? ?x (uh)2 h + 1 2g 0h2 ! = 0 (1.17) 16 Chapter 2 Literature Review In this chapter, the theoretical and empirical GC front investigations are discussed. Some of the relavent experimental investigations into GC flows are presented. Subsequently, the use of analytical box models to simulate these GC flows are discussed before moving on to the more complex numerical models. The one and two-layer SWE models and the corresponding boundary conditions are analyzed in detail since this discussion leads to the proposed numerical models presented in the Methodology, Chapter 4. The advantages and disadvantages between the SWE models and the Navier-stokes models are presented in section 2.4.3, and this chapter is concluded with a discussion of hyperbolic numerical schemes that are used in SWE models. 2.1 The description of the front The description of the GC front (i.e. h and u / uh) has received considerable attention in recent decades. It is well known that the front velocity is very important to the GC propagation and interfacial mixing. In the widely used integral models (discussed in 2.3), the description of the front is one of the primary components along with volume continuity. Moreover, the front of the GC is typically utilized as a boundary condition in the SWE. The SWE are invalid at the GC front because of the importance of vertical accelerations, so the addition of a front condition is required to complete the model (described ahead in 4.4.2). One of the most traditional problems in unsteady, open channel hydraulics is the 1-D dam break problem, which is similar to the description of a GC front. In the dam break problem with a dry, frictionless bed scenario on one side, the front of the fluid (i.e. water) approaches a feather tip solution [Sturm, 2010]. Through characteristic analysis, Ritter was able to compute the velocity of the dam break front, which is described by the following relationship: 2u = c0 where c0 =pgh0 17 is the initial wave celerity [Klemp et al., 1994]. In addition, a similar expression was obtained for the velocity of the upstream moving depression wave: u = co. [Schoklitsch, 1917] compared these dam break results to the propagation of GC by means of laboratory experiments. It was apparent that the upstream moving wave of the GC compared well with the dam break analysis. However, the velocity of the GC front was about 1=2 of the dam break front velocity. This discrepancy lead to a new relationship (equation 1.2), which was formulated by [von Karman, 1940]. [Abbott, 1961] altered this relationship by replacing the Froude number with an empirical coefficient b that was determined from experiments: uLE =b p g0hLE (2.1) in which the subscript LE denotes the front or leading edge of the gravity current. One notices that this expression (2.1) equates to the dam break front velocity when b = 2 and the density of the ambient fluid is much smaller than the GC fluid density [Klemp et al., 1994]. In the description of GC propagation, the correct value for b or the Froude number has been the focus of a number of research investigations. Because of the importance and intrinsic complexity of interfacial mixing (specifically at the GC front), empirical relationships are often applied in SWE models. However, there have been many attempts to improve the theoretical description of the GC front. 2.1.1 Theoretical investigations of GC front conditions From experimental observations, it is known that the celerity of a GC front depends on the density difference (Dr) between the fluids and the depth ratio at the GC front region (i.e. fLE). In equation 2.1 the density difference is incorporated in the form of the reduced gravity. However, the depth ratio is not involved in this expression, so it is expected that b is a function of the depth ratio. One of the main theoretical breakthroughs for the description of the GC front was formulated by [Benjamin, 1968] in which the Froude-like number is a function of the fractional depth. 18 In order to formulate this expression at the leading edge of the GC, [Benjamin, 1968] equated the buoyancy forces to the drag forces between the two fluids at the front. For Ritter?s traditional dam break problem, the ambient fluid density is neglected as well as the shear with the ambient fluid, which leads to the feather tip front. In GC flows these simplifications are strictly invalid. Benjamin?s analysis was applied to a system in which an air cavity intruded over top of the water- filled pipe with a rectangular cross-section undergoing discharging at one end. The following development follows [Klemp et al., 1994], who applied Benjamin?s analysis in the context of a heavier fluid traveling underneath a lighter fluid (e.g. the propagation of a saltwater wedge into freshwater). In order to obtain the relationship for the GC front velocity, the flow is assumed steady and depth-averaged. The vertical momentum equation (1.11) is simply reduced to the hydrostatic relationship [Klemp et al., 1994]: ?P ?z = g 0 (2.2) The x-momentum equation is formulated from equation 1.10, and the effects of shear are neglected. This expression is integrated from a point far ahead r of the GC front to a point l far behind the GC across the channel of depth H, which yields the following flow force balance [Klemp et al., 1994]: Z H 0 P r +u2r dz =Z H 0 P l +u2l dz (2.3) in which ur is the ambient velocity at r and ul is the ambient velocity at l. Pl and Pr represent the pressures at l and r, respectively (see Figure 2.1). In order to compute this flow force balance, the continuity of volume is required: ul = urH=(H h). Moreover, Pl is determined by applying the vertical momentum equation (2.2) far behind the nose: Pl = 12u2r g0h. When these two expressions are applied to (2.3) along with the assumptions that Pr = 0 far ahead of the front and Ps = 12u2r at the stagnation point x f , the result is Benjamin?s formula: uLEp g0hLE = Fr(fLE)= (2 f LE)(1 fLE) 1+fLE 1=2 (2.4) 19 Figure 2.1: Schematic diagram of Benjamin?s formulation for the case of a heavier fluid propagat- ing underneath a lighter ambient fluid. Adapted from [Klemp et al., 1994]. This important expression represents one equation for two unknowns: uLE and hLE and is a solution for equation 2.1. From this relationship (2.4), it is clear that the velocity of the GC depends on the reduced acceleration due to gravity, the depth and the depth ratio. The implementation of this expression in numerical models is discussed in subsection 4.4.2. [Shin et al., 2004] extended Benjamin?s expression to include the effects of mass transfer between the front of the ambient fluid advance and the front of the heavier fluid. Their results are exactly the same as Benjamin?s for the lock-exchange problem because the interfacial waves (in which energy is transferred between the two fluids at the interface) generated from the ambient front are unable to reach the leading edge of the GC. However, the results from [Shin et al., 2004] expression improve in accuracy over Benjamin?s solution as the depth ratio decreases to a deep ambient problem [Shin et al., 2004]. For Boussinesq fluids, the expression formulated by [Shin et al., 2004] is: Up g0h = Fr(fLE)= 1 2 s h H 2 hH (2.5) This expression was reduced from the original non-Boussinesq relationship; furthermore, [Shin et al., 2004] expects the non-Boussinesq relationship to break down for large density dif- ferences. This reasoning is based on their assumption that the current and ambient fronts are connected by a horizontal interface with uniform flow conditions. This assumption is warranted 20 for Boussinesq systems; however, in non-Boussinesq fluids the interface is more complex and this assumption is invalid [Keller and Chyou, 1991]. 2.1.2 Empirical investigations of GC front conditions Huppert and Simpson (HS): This front condition was constructed empirically from the ex- periments of [Huppert and Simpson, 1980] and [Keulegan, 1957]. It performs better for full depth lock releases but is able to simulate all depth ratios [Ungarish, 2009]. Fr(fLE)= 8> >< >>: 1 2f 13 (0:075 f LE < 1) 1:19 (fLE 0:075) (2.6) Rottman and Simpson (RS): In this front condition, Benjamin?s theoretical condition is al- tered with an empirical front condition b. [Rottman and Simpson, 1983] adjusted b based on their partial depth experiments because their model was unable to simulate large initial values for f. On the contrary to the HS front condition, the RS front condition performs well for small fractional depths but not as well for full depth releases [Ungarish, 2009]. Fr(fLE)= bp2 (2 f LE)(1 fLE) 1+fLE 1=2 (2.7) Ungarish and Zemach (UZ): [Ungarish and Zemach, 2005] formulated a new front condition that bridges the HS and RS front conditions, so it is able to perform well for a large range of fractional depth. Fr(fLE)=(1+3fLE) 1=2 (2.8) Kranenburg (K): [Kranenburg, 1978] adjusted Benjamin?s theoretical front condition to al- low for an initially moving ambient. In addition, they utilized an empirical coefficient k that they estimated to equal 0.6. While the focus of this thesis is on counter-flows in which the ambient fluid has an initial condition of zero velocity, in many hydraulic applications the 21 Figure 2.2: Schematic diagram of a GC in an initially moving ambient system. Adapted from [Wright and Paez-Rivadeneira, 1996]. ambient fluid has non-zero velocity, either coflows or counter-current flows with respect to GC (e.g. cold discharges in river flows). The expression presented by [Kranenburg, 1978] is: C1 (grH)1=2 = Fr(h=H)= h=H(1 h=H)(2 h=H) 1+h=H +k(1 h=H) 1=2 (2.9) in which C1 = uLE q=H is equal to the velocity uLE of the GC for ambient fluids with initially zero velocity. Notice that the LHS of Equation 2.9 is a function of the ambient fluid depth instead of the GC depth. The schematic diagram for this system is illustrated in Figure 2.2. Wright (W): [Wright et al., 1990] developed a generalized relationship that is able to simu- late both co-flows and counter-flows. The depth ratio f in this relationship is based on the depth behind the nose (see Figure 2.2). Several experiments conducted for various ambient velocities demonstrated the effectiveness in this approach [Wright and Paez-Rivadeneira, 1996]. C1 (grH)1=2 = Fr(h=H)=[1 (1+qr)h=H] h=H(1 h=H)3 (1 h=H)3 q3r(h=H)3 1=2 (2.10) where qr is the velocity ratio between the ambient fluid and the GC (qr = q2=q1) in which q = uh. 22 2.2 Experimental investigations Along with analytical and numerical research, experimental investigations have provided im- portant insights into GC flows. [Simpson and Britter, 1979] studied a large range of dimensionless depths for GC flows and determined the resulting dimensionless velocities, rate of mixing, and the depth of the mixed layer behind the head. Their apparatus consisted of a moving floor and a downstream weir, which allowed them to halt the GC flow in order to obtain more accurate mea- surements. A similar procedure was utilized by [Wright and Paez-Rivadeneira, 1996] to analyze the effects of coflows and counterflows. For lock-exchange flows, [Rottman and Simpson, 1983] analyzed a large range of initial depth ratios (f) and compared the results their two-layer SWE model. For large initial depths, [Rottman and Simpson, 1983] concluded experimentally that the upstream moving jump occurs when f 0:7 and that the reflection off of the upstream boundary produces a hydraulic drop that eventually overtakes the GC front marking the beginning of the self similar stage. For smaller initial depths, the upstream moving jump is replaced by a depression wave [Rottman and Simpson, 1983]. In [Hacker et al., 1996] a digital image technique (DigImage, [Dalziel, 1993]) was used to de- termine the density structure for lock-release GC flows. They were able to observe the detrainment of dense fluid at the GC nose by breaking Kelvin-Helmholtz waves, which produced a stratified region behind the nose. The diluted fluid in this region was replaced by dense fluid behind the nose, which traveled toward the front of the nose. Thus, this process causes a re-circulation pro- cess and the entire GC becomes diluted. The work by [Hacker et al., 1996] utilized the following expressions to compute the total depth densimetric Froude number (FrH) and the mean Reynolds number (Rem) for the lock-exchange problem: FrH = umph 0 g0 Rem = 12 um h0n (2.11) where h0 is the initial depth of the denser fluid, g00 is the initial reduced gravity, um is the mean velocity of the GC front and nu is the kinematic viscosity. Rem is multiplied by 12h0 because this is 23 the energy conserving depth described by [Benjamin, 1968]. These dimensionless parameters are provided for the experiments conducted for this thesis in section 4.1. [Gerber et al., 2011] used a particle image velocimetry (PIV) to measure the Reynolds stress and shear production of turbulence of a stably stratified GC. Their experimental results compared well to their Reynold?s averaged Navier-stokes numerical model. [Firoozabadi et al., 2010] used MicroADV probes to measure the turbulence energy, Reynolds stresses and turbulence intensity of 3-D GC flows for various discharges, concentrations and slopes. As the discharge and the concentration increased, the turbulence kinetic energy also increased. Their results also indicate that the normalized turbulence intensity does not change with bed slope. A large number of experimental investigations have utilized sodium chloride (NaCl) in their experiments to generate Boussinesq fluids ([Rottman and Simpson, 1983], [Shin et al., 2004], etc.). The same approach was implemented in the experiments presented in this thesis. For non- Boussinesq lock-exchange flows, [Lowe et al., 2005] used sodium iodide (NAI) in addition to NACL in order to obtain density ratios (r2=r1) between 0.6 and 1.0. From their non-Boussinesq experiments, they determined that in most cases the lighter current retains its energy-conserving depth while the denser current is dissipative, so the depth decreases according to a decreasing density ratio. 2.3 Analytical description of GC flows with integral models For a quick and simple estimation of the GC propagation, integral models (also referred to as box models) can be a useful modeling option. In integral models only the resistance of the advance (the front condition) and the conservation of buoyancy are imposed. Thus, GC flows for an integral model are described by the spreading of an equal area rectangle, so horizontal variations in the flow are ignored [Slim, 2006]. While this is an unrealistic assumption; however, these modeling results can provide important insights into the average flow depth and velocity [Ungarish, 2009]. 24 2.3.1 Deep ambient with constant Froude number The first integral model is derived for a homogeneous GC propagating into a deep ambient f 0. However, similar models have also been proposed for particle-driven GC for both incom- pressible [Hogg et al., 1999] and compressible fluids [Timmermans et al., 2001]. In the context of the deep ambient, the Froude number, which is the ratio between the inertial and gravitational forces, is assumed constant (valid for the deep ambient problem). For the two-dimensional or pla- nar case, there are two key variables in the internal cell calculations: GC length xN(t) and height hN(t). As a result, two equations are required to complete the box model [Ungarish, 2009]. The first equation is the simple volume continuity equation: xN(t)hN(t)= V (2.12) in which V is the volume of the dense GC. For the fixed volume GC (e.g. the lock-exchange problem) with inertial-buoyancy balance (neglecting viscosity), the derivation begins with the front condition. This BC, which can be derived through dimensional analysis, with the Froude number (Fr) notation used in [Benjamin, 1968] is written as: uN = dxNdt = Fr(g0hN)1=2 (2.13) in which dxN=dt represents the GC front velocity; the Froude number is defined as Fr = u=pg0hN. This BC is the same relationship that is implemented for the SWE (described ahead in this paper). Unlike in Benjamin?s formulation, the Fr in (2.13) is assumed constant for this derivation. The combination of (2.12) and (2.13) yields the following relationship for xN(t) (or for hN(t)) [Ungar- ish, 2009]: Z t2 t1 = 1V 1=2 Z xN x 0 [Fr(f(x))] 1x 1=2dx (2.14) 25 where f(x)= V=(Hx) through algebraic manipulation. x N = x 0 1+ 32 Fr(g 0 V)1=2 x 0 3=2 t !2=3 (2.15) Once xN is determined from (2.15), h is computed from the continuity equation (2.12). The same strategy can be implemented to obtain an integral model (2.16) for the axisymmetric scenario: r = r0 1+2Fr(g 0 V)1=2 ro 2p1=2 ! (2.16) where r is the radius and r0 is the initial condition. 2.3.2 Generalized solution [Huppert and Simpson, 1980] developed an innovative integral model for both the two-dimensional and axisymmetric scenarios. Instead of using a constant Fr, they implemented the empirical front condition (HS) that was derived from their experiments and from [Keulegan, 1957]. As a result, they obtain a solution for both the slumping and the self-similar stage, which depends on the depth ratio (f). For the two-dimensional case, [Huppert and Simpson, 1980] substitutes the continuity equation (2.12) into the HS front condition. For the slumping stage, this relationship is written as: Fr = dxNdt =(g0 V=xN)1=2 = 12 V xNH 1=3 (0:075 hN=H < 1) (2.17) After integration and isolating xN, the slumping stage relationship is: xN = x0 1+ 712 g03 V H 2 x0 7 1=6 t !6=7 (x0 xN xs) (x0 x xs) (2.18) in which x0 is the initial lock length and xs is the slumping distance, which [Huppert and Simpson, 1980] define as the distance when (hN=H 0:075) 26 For the self-similar stage (i.e. hN=H 0:075), the authors propose that the Froude number is 1:19 according to the HS front condition. Following the same procedure as the slumping stage, the relationship for the self-similar stage is: xN = h x3=20 +1:78(g0 V)1=2(t ts) i2=3 (xs < x ) (2.19) in which x is the distance when viscous effects become larger than inertial effects (1.8) [Huppert and Simpson, 1980]: For the axisymmetric problem, which is derived in a similar manner to the two-dimensional case, the solution in the slumping stage is: [Huppert and Simpson, 1980]: r = r4=30 23p 1=6(g02 V H2)1=6t 3=4 (r0 r rs) (2.20) where rs and r in cylindrical coordinates are equivalent to xs and x for the two-dimensional problem. The relationship for the self-similar stage is: [Huppert and Simpson, 1980]: r =[r2s +2:37p 1=2(g0 V)1=2t]1=2 (rs r r ) (2.21) In summary, integral models are simple to implement and require minimal computational time as compared to the SWE and Navier-Stokes models since they lend themselves to analytical solutions. The development of such models does not require numerical methods or schemes, which can become extremely complex in other models (e.g. higher order Riemann solvers). The cost of this reduction in complexity is the use of assumptions: that GC spread in a quasi-steady manner with a constant depth in an exactly rectangular shape. Nevertheless, integral models can provide important insights into the flow of GC in a channel. The results for the deep ambient scenario are satisfactory depending on the context. However, when the initial depth ratio approaches the lock-exchange problem (f > 0:5), the integral model is unable to accurately capture the depth profile and other flow features although the Huppert and 27 Simpson formulation is a viable alternative. As a result, more complex models (e.g. SWE and Navier-Stokes) are required for the accurate simulation of these flows. 2.4 Numerical modeling of gravity currents Experimental studies (described in section 2.2) are an imperative tool for the understanding of GC flows. However, it is not practical to construct scale models for every type of scenario that is present in the field. It is very difficult to simulate experimentally several common flow features: tidal intrusions, complex geometries, non-Boussinesq systems, flows over porous substrates, etc. Analytical and numerical models are essential for a better understanding of GC propagation in which experimental data and/or field observations are used to calibrate these numerical models. Once calibrated, these numerical models can become powerful tools, yielding useful predic- tions of GC flows that are very important to the scientist or engineer in the field, who may need to determine the mixing mechanisms and/or flow field for a variety of issues such as the transport of contaminants. While analytical models are a viable option for these situations, they suffer from se- vere assumptions as described in the previous text. As a result, many GC flow applications require the use of numerical modeling. There are two primary methods for the numerical modeling of GC flows: SWE models and Navier-Stokes models, which are described below. 2.4.1 One-Layer SWE model The SWE models are generally separated into two groups: one-layer and two-layer. The two-layer SWE models differ in that they include the effects of the ambient fluid velocity in the internal cell calculations. Both of these models can provide an accurate estimation of the GC front and are computationally efficient compared to the more complex Navier-Stokes models. In the following text the details of the SWE models are analyzed and several of the important research investigations are included. The one-layer SWE models (4.2) are widely used for GC flowing into deep ambient fluids. [Ungarish and Huppert, 2004] used the one-layer SWE to simulate a GC traveling at the base of 28 a stratified ambient and determined that the stratification reduces the GC velocity. The key at- tractions in the single-layer SWE are the simplicity and the computational efficiency. Because the velocity of the ambient fluid is neglected, discontinuous solutions are generally not observed for flows in simplified rectangular geometries without obstacles. Details of this one-layer formula- tion are found in [Ungarish, 2009] for both Boussinesq and non-Boussinesq fluids. In this study, Boussinesq fluids are focused; however, it is relatively simple to formulate the one-layer model without the Boussinesq assumption. Therefore, the more generalized one-layer, non-Boussinesq SWE model is the focus of this discussion: ?h ?t + ?uh ?x = 0 ?uh ?t + ? ?x (uh)2 h + 1 2g 0h2 ! = 0 For this SWE model, the one-layer assumption (i.e. the ambient fluid velocity is neglected) poses inaccuracies for larger depth ratios. However, the single-layer SWE model provides fairly accurate results for the GC trajectory for depth ratios approaching unity f!1 even though the one layer assumption is least valid [Ungarish, 2007]. For the full depth lock-exchange problem, the ac- curacy is believed to be caused by a fortunate balance of errors from different simplifications (e.g. ambient velocity effects, entrainment, etc.) [Ungarish, 2007]. The ambient fluid velocity increases the GC velocity, but this increase is balanced by the reduction in velocity due to entrainment. The single-layer SWE model cannot simulate either one of these complex flow features, but the com- bined effect of neglecting these features equates to an accurate estimation for the GC front velocity. On the other hand, the flow depth in the initial slumping stage is inaccurately computed by this SWE model for larger initial depth ratios. Therefore, it is customary to move to the two-layer SWE model for this scenario. [Ungarish, 2007] formulated a single-layer SWE model for a large range of depth ratios and density differences. Furthermore, in this study Boussinesq GC flows for the lock-exchange prob- lem were analyzed as a special case of their non-Boussinesq model. Although there are several 29 important underlying assumptions in this approach, the generality, absence of empirical coeffi- cients and computational efficiency are attractive [Ungarish, 2007]. The characteristic equations for the one-layer model that are implemented in [Ungarish, 2007] are: g012 dhh1=2 du = 0 on dxdt = u (g0h)1=2 (2.22) in which the second expression represents the characteristic flow velocities. Although these equa- tions (2.22) can be used to compute the internal cell calculations, [Ungarish, 2007] implemented the Lax-Wendroff two step finite difference method. This numerical scheme is well documented and performs well in the absence of shocks [Chaudhry, 2008], which are generally not present in the one-layer model. However, the characteristic approach was still used in the computations by [Ungarish, 2007]. Since Benjamin?s theoretical front condition was used at the GC front BC, a characteristic equation was also used to solve for both unknowns (i.e. hLE and uLE). Therefore, the traditional solution to this problem has been to use a characteristic equation along with a front condition at this BC (described further in 4.4.2). For the other boundary condition that is located upstream of the release, either a reflective [Toro, 2001] or characteristic [Sturm, 2010] boundary condition may be used. One of the inconvenient components of SWE models for GC is the need for a boundary condi- tion at the GC front. Therefore, [Ruo and Chen, 2007] formulated the modified SWE for inviscid GC in order to eliminate the need for a front condition. Therefore, the only required boundary conditions are at the physical boundaries or end walls. In order to formulate their mathematical model, [Ruo and Chen, 2007] included the effects of the ambient resistance, which is normally accounted for by the front condition. The continuity equation for the modified SWE is the same as in the traditional SWE . How- ever, the momentum equation is altered in order to account for the ambient resistance. [Ruo and Chen, 2007] begins their derivation with equations (1.9) and (1.11). Then, the pressure term is 30 reduced from the following relationship, which is also used in the two-layer formulations: P(x;y;t)= Pi(x;t)+rg(H h) (2.23) where Pi(x;t) is the pressure at the fluid interface, and H and h are the ambient and current depth, respectively. In addition, equation (1.11) is implemented in order to develop a new x-momentum equation: ?u ?t +u ?u ?x +gr ?h ?x (2.24) However, a new problem emerges in this mathematical model: How to determine the pressure at the fluid interface. [Ruo and Chen, 2007] implements the unsteady Bernoulli equation along the instant streamline in order to overcome this issue. Their solution for Pi(x;t) is the following relationship: P i(x;t)= Pt ragh 12ra(eVi)2 ra?F?t y=h(x;t) (2.25) in which Vi is the absolute velocity of the ambient fluid at the interface, F is the velocity potential so that V = F, and Pt is the total pressure of the ambient fluid at a location far from the current [Ruo and Chen, 2007], which in this case is constant. When equation 2.25 is substituted into equation 2.24, it becomes clear that the new difficulties become computing eVi and f. Furthermore, [Ruo and Chen, 2007] breaks apart eVi into two components: tangential velocity wt and normal velocity wn. It is then assumed that ?(w2n)=?x >>?(w2t )=?x, which leads to the following relationships: Vi = wn = ?h=?tp1+(?h=?x)2 ??x ?F ?t = ??x (?h=?t)2 1+(?h=?x)2 (2.26) so that momentum equation for the modified SWE is [Ruo and Chen, 2007]: ?u ?t +u ?u ?x +gr ?h ?x = g 2 ? ?x (?h=?t)2 1+(?h=?x)2 (2.27) 31 in which g = ra=r. A look at this highly nonlinear momentum equation suggests that (2.27) is very similar to (1.17). However, there is an added source term on the RHS of (2.27), so [Ruo and Chen, 2007] solved this momentum equation using the perturbation method. Although they obtain accurate results that compare well to the more traditional SWE models for smaller depth ratios, [Ruo and Chen, 2007] does not believe that their approach can be extended to the lock-exchange problem. 2.4.2 Two-Layer SWE model The two-layer SWE have received more attention than the single-layer alternative for GC flow modeling. For instance, the lock-exchange problem is a relatively easy experiment to perform in which the two-layer SWE model is preferred. Numerous research investigations have focused on a wide range of GC flows and the use of the two-layer SWE. For example, [Bonnecaze et al., 1993] applied the two-layer SWE to particle-driven GC flows. However, this discussion focuses on homogeneous GC flows on a horizontal surface. As previously described, the two-layer SWE are preferred when f!1. When f!0, the results converge to the one-layer SWE model results. [Rottman and Simpson, 1983] developed the first two-layer SWE model for GC flows in the context of partial-depth releases. For larger fractional depth (i.e. f > 0:5), they encountered difficulties, which are described below. In the case of a two-layered model, the SWE represent a system with four PDEs to account for each phase?s mass and linear momentum conservation [Rottman and Simpson, 1983]. ?h1 ?t + ? ?x(u1h1)= 0 (2.28) ?h2 ?t + ? ?x(u2h2)= 0 (2.29) ?u1 ?t +u1 ?u1 ?x = 1 r1 ?P i ?x g ?h1 ?x (2.30) ?u2 ?t +u2 ?u2 ?x = 1 r2 ?P i ?x g ?h1 ?x (2.31) 32 The subscripts (1, 2) represent the lower (heavier) and upper (lighter) fluid layers, respectively. [Rottman and Simpson, 1983] implemented equation (2.28) as the sole continuity equation for their model. Then, the ambient variables (i.e. h2 and u2) were computed from the following rela- tionships: h2(x;t)= H h1(x;t) (Depth uniformity) (2.32a) u2 = u1h1h 2 (Local continuity) (2.32b) where equation 2.32a is deduced from the rigid lid assumption in which the total depth of the system does not change. The relationship for the ambient velocity (equation 2.32b) is formulated from the two continuity equations (2.28 and 2.29) along with 2.32a and is valid for constant volume gravity currents. In order to utilize one relationship for the equation of motion and to eliminate Pi (the pressure at the fluid interface), equations (2.30) and (2.31) are substituted together into a single momentum equation: (1+ra)?u1?t + 1 raH +h1H h 1 u1?u1?x + g0 r(1+a)3 u 21 H ?h 1 ?x = 0 (2.33) Furthermore, [Rottman and Simpson, 1983] implemented the Boussinesq simplification, which eliminates r (i.e. r2=r1), yielding: ?u1 ?t +(1 2a)u1 ?u1 ?x +(1 b)g 0?h1 ?x = 0 (2.34) The mathematical model that is formulated in [Rottman and Simpson, 1983] serves as the basis of the proposed model, which is presented in section 4.3. The two-layer mathematical model formulated in [Rottman and Simpson, 1983] is presented below in primitive format. By using the assumption of local flow continuity and overall depth uniformity, the governing equations are 33 reduced to 2 PDEs that represent the conservation of mass and momentum of the GC: ?h1 ?t + ? ?x(u1h1)= 0 ?u1 ?t +(1 2a)u1 ?u1 ?x +(1 b)g 0?h1 ?x = 0 (2.35) where a = h1H h 1 b = h1H + u 21 g0H 1 h1H 2 (2.36) [Rottman and Simpson, 1983] implemented the first expression in (2.35) as the sole continuity equation for their model along with the single momentum equation. However, one of the draw backs is clear from the parameters a and b in which the solution is undefined as the fractional depth approaches unity h!H. A solution to this problem is proposed in section 4.3. The context in which these two-layered models are most often implemented is in the release of dense fluids into relatively small ambient depths (e.g. the lock-exchange problem). The dense fluid is often originally contained by a barrier that is suddenly removed and a wall is located upstream of the release, which represents a physical boundary. Once the four PDEs are reduced to a single pair of PDEs, the numerical solution that is implemented in [Rottman and Simpson, 1983] involves characteristic analysis where these PDEs are transformed into the system of ODEs presented below: h1 du1dh 1 (1 2a)u1 +c = 0 c = u1(1 a) u21a2 +g0h1(1 b) 1=2 (2.37) One of the main benefits in this approach is computational efficiency and simplicity. However, the method of characteristics is not able to simulate shocks in which the characteristic lines of the same family intercept one another [Chaudhry, 2008]. In addition to the internal cell calculations, the c+ equation is implemented at the GC front boundary condition to provide closure to solve for the two local unknowns, hLE and uLE. The results in [Rottman and Simpson, 1983] compare well with experiments and exhibit minimal oscillations throughout the slumping stage; however, their front condition is only tested for smaller depth ratios 34 (f 0:5). Furthermore, their front condition (RS) decreases in accuracy for larger depth ratios, especially for the full depth lock-exchange problem [Ungarish and Zemach, 2005]. When f 0:5, the numerical results for [Rottman and Simpson, 1983] two-layer SWE model are non-physical (i.e. there is a multi-valued solution for the depth). Moreover, the upstream propagating shock that is generated for f > 0:5 becomes stronger as it travels at higher speeds, and the reflection off of a back wall triggers the formation of another sharp interface that eventually catches up with the GC front. While the numerical modeling suggests that the limit for this sharp interface to develop is f >0:5, experiments from [Rottman and Simpson, 1983] point that this limit is in fact closer to f > 0:7, albeit the transition is not exactly clear. The discrepancy, according to [Rottman and Simpson, 1983], is a result of the smooth profile in the velocity for the experiments and the mixing between the two fluids, which is neglected in the shallow water model. As a result, the two layer model developed by [Rottman and Simpson, 1983] is accurate quantitatively for f 0:5; the authors pointed that numerical simulations with f > 0:5 would require new developments in the theory of flow discontinuities. The work by [Klemp et al., 1994] extended the two-layer model developed by [Rottman and Simpson, 1983] in order to solve for the depth ratios greater than 0.5. Instead of using the method of characteristics, [Klemp et al., 1994] implemented a finite difference model using the second- order accurate Leapfrog scheme for the solution of the two-layered SWE. They enforce a bound- ary condition at the upstream propagating discontinuity based on an analytical relationship for hydraulic jumps. [Klemp et al., 1994] specified the pressure far downstream pr(d) according to the following expression: pr(d)= 12 V 2f u022 +g0D (2.38) in which Vf is the speed of the upstream moving jump, u02 is the velocity of the ambient fluid just downstream of the jump and D is the dissipation term (see Figure 2.3). Applying the flow force balance far upstream and downstream and using equation 2.38 and the continuity relations, [Klemp et al., 1994] obtained the following solution for the upstream moving 35 shock: V 2f g0H = f2 a2 2D0 [2f2 a=a]+(1 f)2(1 2a)=(1 a)2 (2.39) in which h f is the height of the upstream moving jump, f = h0=H, a =(h0 h f)=H and D0=D=H. [Klemp et al., 1994] set the dissipation term (D) equal to zero for consistency with Benjamin?s theory and based on the evaluation of weak hydraulic jumps. Equation 2.39 is not sufficient to solve for the moving jump BC, so [Klemp et al., 1994] implemented an additional condition for steady-state jumps based on the formulation in [Stommel and Farmer, 1952]: u 1 +Vf 2 g0h1 + u 2 +Vf 2 g0h2 = 1 (2.40) The schematic diagram for a GC flowing into a shallow ambient fluid is illustrated in Figure 2.3. This analytical relationship was required for their SWE model for f > 0:5 because the leading edge of the backward propagation wave travels slower than the trailing region until the front of the wave steepens into a discontinuity. [Klemp et al., 1994] describes this phenomenon as a hydraulic jump according to long wave theory; however, this jump differs from the classical hydraulic jump because it does not form an abrupt discontinuity. Instead, [Klemp et al., 1994] points that the slope of this jump remains finite due to non-hydrostatic effects, which are outside the scope of long wave theory. [Klemp et al., 1994] presented their results for a large range of depth ratios; however, they have focused their model on the initial stages of the current propagation, up to the instant when the upstream moving front contacts the end wall. The work by [Ungarish and Zemach, 2005] presents a model that extends the theory of [Rottman and Simpson, 1983] and [Klemp et al., 1994] in order to include the entire slumping stage in their formulation, extending the simulation beyond the instant when the backward prop- agating front contacts the upstream wall boundary. [Ungarish and Zemach, 2005] implemented a second-order accurate finite difference technique similar to [Klemp et al., 1994], but they used the Lax-Wendroff (LxW) two-step numerical scheme. The model performs an explicit tracking of the upstream propagating front before and after the reflection off the back wall, and these fronts are 36 Figure 2.3: Schematic diagram for both fronts that is generated from a GC flowing into a shallow ambient (f > 0:5). Adapted from [Ungarish and Zemach, 2005]. treated as discontinuities in the depth. In a certain way, the treatment of these fronts resembles a shock fitting approach used to track open channel bores presented in [Cunge et al., 1980]. When the shock contacts the end wall, the model by [Ungarish and Zemach, 2005] predicts the propaga- tion of a discontinuity towards the nose, which is also explicitly tracked. [Ungarish and Zemach, 2005] track this bore until the end of the dam break stage where the depth and the velocity of the current begin to decrease. One of the drawbacks in the approach by [Ungarish and Zemach, 2005] is the added model complexity in introducing two additional boundary conditions to track the flow discontinuities. In addition, the model requires two solution methods depending on the depth ratio. If f > 0:5, an explicit feature tracking approach is required; however, if f 0:5, the simpler approach is imple- mented in which no shock tracking methods are utilized except at the front of the GC [Ungarish, 2009]. In addition, experimental results indicate the discontinuities that form for two-layer GC flows are not represented by an exactly vertical interface between the two fluids. Therefore, it is an idealized approximation in which additional experimental confirmation is necessary. This motivated the development of an unified and more general shock-capturing approach that is com- putationally consistent for all depth ratios, which is formulated and discussed in section 4.3. 37 Although the BC at the GC front, upstream hydraulic jump and the reflected bore have been discussed, there is still another feature that needs to be addressed. For gravity currents flowing into shallow ambients (f 0:8 or larger depending on FrLE), a critical condition emerges at the GC nose [Ungarish, 2009]. This constraint occurs when b* < 1 (see equation 2.41) so that the GC front cannot travel faster than the c+ characteristic, which would cause cavitation [Klemp et al., 1994]. Therefore, the depth and velocity constraints are computed by setting b* = 1. This critical change of behavior is clearly seen from the following dimensionless, two-layer characteristic relationships [Ungarish, 2009]: 8 >>< >>: dh* du* = 1 1 b* h a*u* p(a*u*)2 +(1 b*)h* i ; on dx* dt* = c = u *(1 a*) p(a*u*)2 +(1 b*)h* (2.41) in which the asterisk denotes dimensionless variables. The parameters, a* and b*, are equivalent to the expressions in equation 2.36 in dimensionless form. In equation 2.41 and throughout Chapter 5, the dimensional variables are converted to dimensionless variables, denoted with an asterisk, according to the following procedure [Ungarish, 2009]: fx*; h*; t*; u*g=fx=x0; h=h0; t=T; u=Ug (2.42) where U = g0h0 1=2 T = x0U By using the critical nose approach described above, the constraints when implementing Ben- jamin?s condition are h = 0:347H and u* = 0:527. Thus, the depth and velocity cannot be greater than these constraints. More recently, [Adduce et al., 2011] formulated an alternative two-layer SWE model in which the conserved variables are: r1h1, r2h2, u1 and u2. Their primary objective was to include the effects of entrainment and shear forces, which decrease the GC front velocity. These features 38 enable [Adduce et al., 2011] to implement the theoretical front condition from [Shin et al., 2004]. However, their approach is more complex and is still largely dependent on empirical parameters. For the majority of SWE numerical models, the important effects of entrainment and shear forces require the use of an empirical front condition. However, [Adduce et al., 2011] avoided this approach by included these effects in the SWE formulation: ?(r1h1) ?t + ?(r1h1u1) ?x =r2ue ?(r2h2) ?t + ?(r2h2u2) ?x = r2ue ?u1 ?t + ? ?x u2 1 2 + r 2h2 +r1h1 r1 gcosq = gsinq t 1b +t 21r 1h1 ?u2 ?t + ? ?x u2 2 2 + r 1h1 r1 + r2h2 r2 gcosq = gsinq + t 2b +t 12r 2h2 (2.43) in which r1 and r2 are the densities of the denser GC and the ambient fluid, respectively. ue [L=T] is an entrainment parameter and q is the angle between the bottom of the tank and the horizontal [Adduce et al., 2011]. t1b and t2b represent the shear stress at the lower and upper boundaries, respectively. t12 and t21 are the interfacial shear stresses in each direction. The entrainment parameter (ue) was computed from an empirical relationship based on the formulation in [Ellison and Turner, 1959]. [Adduce et al., 2011] altered their expression in order to simulate the full depth lock-exchange system: ue u1 = k F2r F2r +5 (2.44) Then, [Adduce et al., 2011] computed the density of the GC r1(t) according to the following expression: r1(t)= M1 + Z t 0 dt Z xn(t) 0 r2ue dx V1 + Z t 0 dt Z xn(t) 0 ue dx (2.45) in which M1 is the mass and V1 is the volume of the GC. The advantage of this approach is that it takes into account the entrainment between the two fluids, which allows for a more accurate simulation of the GC front. However, the accuracy for this 39 approach is highly dependent on the empirical coefficient k, which was calibrated from experiments to be equal to 0.48. One questions whether other cases of GC flows would have the same k value. In addition, the model complexity and the computational time for this model should be substantially greater compared to the approach from [Ungarish and Zemach, 2005], or the model proposed in this thesis. 2.4.3 Navier-Stokes models For GC research, the majority of Navier-Stokes models are solved using Direct Numerical Simulation (DNS) or Large Eddy Simulation (LES). These models are able to simulate interfacial mixing and provide high-resolution results at the GC nose. A large number of insights into GC flows has been gained through high-resolution DNS models. [H ?artel et al., 2000] determined that the stagnation point is not at the GC front but at a point just upstream. [Necker et al., 2002] com- pared 2-D and 3-D models for particle-driven GC and concluded that 3-D models are required to accurately model the spanwise vortices or interfacial instabilities. These features survive over long times for 2-D simulations but are more rapidly broken down in the 3-D models. [Birman et al., 2005] investigated non-Boussinesq GC and confirmed that the denser GC front dissipates more for larger density differences whereas the lighter ambient front remains very close to the energy con- serving theory seen in Boussinesq systems. Moreover, the denser front is reduced to substantially less than half the channel height depending on the density difference while the light front remains near half the channel height. Although DNS models are the best choice theoretically, there are many drawbacks. The most important limitations are the intensive computational requirements and the increased complexity. Moreover, for 3-D simulations in present-day computing systems, the DNS method is limited to GC flows where Re = O(1;000) [Slim, 2006]. Another option for high resolution simulations of GC flows is Large Eddy Simulation (LES). [Ozgokmen et al., 2009] analyzed several LES models and concluded that a hybrid approach, which combined the dynamic eddy viscosity and approximate deconvolution approaches, com- pared best to the DNS results and was 1,000 times more computationally efficient. However, in 40 the first Navier-Stokes models that were used to simulate GC flows the details of the interfacial instabilities were not the primary objective. [Klemp et al., 1994] solved the two-dimensional non- hydrostatic incompressible Boussinesq equations in addition to the two-layer SWE to determine the validity of Benjamin?s front condition. The results from [Klemp et al., 1994] two-dimensional simulations indicate good agreement in comparison with the two-layer SWE. The main differences between their two models are the simulation of the upstream moving disturbance (for f > 0:5) and the ?hydrostatic circulation in the vicinity of the head? [Klemp et al., 1994]. [Ungarish and Hup- pert, 2004] also found good agreement between the one-layer SWE and Navier-Stokes simulations in the context of a GC propagating at the base of a stratified ambient fluid. [Klemp et al., 1994] also determined that Benjamin?s front condition performed well in com- parison with their model. Moreover, they concluded that the discrepancy between the SWE models and experiments is not caused by Benjamin?s theory but by other simplifications in the SWE mod- els the most important being the inability to simulate interfacial mixing. The validity of Benjamin?s theoretical expression is further discussed in section 5.3 with the use of MicroADV devices. 2.5 Hyperbolic numerical schemes Since the proposed numerical models for this thesis are implemented for the one and two-layer SWE, which is an example of a hyperbolic PDE, the focus is on hyperbolic numerical schemes. The majority of research and implementations for these numerical schemes is based on the more traditional water-air component of hydraulics. In the context of dam-breaks, an extensive compar- ison for both linear and nonlinear explicit numerical schemes is presented in [Zoppou and Roberts, 2003] for the SWE. Moreover, [Vasconcelos and Nunes, 2009] presented a detailed analysis for all of the numerical schemes that are implemented in this thesis: Lax-Wendroff (LxW), Modified FORCE (MFORCE), HLL, and Roe. They analyzed each of the numerical schemes based on error analysis (i.e. L1 and L?) and computational efficiency in the context of the dam break problem. Although the Roe scheme was generally the most accurate of the compared schemes, the HLL 41 scheme was considered the most appropriate choice based on the improved computational perfor- mance and similar accuracy. In the absence of shocks, the LxW scheme performed very well. The MFORCE scheme was accurate for each scenario, but the computational time rivaled the nonlin- ear schemes, which are considered superior [Vasconcelos and Nunes, 2009]. The implementation of these numerical schemes to inviscid GC using the finite volume method is further discussed in 4.4.1. Numerical schemes are applied on broad range of applications in hydraulics. More recently, nonlinear numerical schemes based on the solution of the initial value problem presented by Rie- mann (i.e. Riemann solvers) have been applied in a large range of hydraulic problems that include: one and two-dimensional dam-breaks [Zoppou and Roberts, 2003], the structural collapse of wa- ter supply reservoirs in urban areas [Zoppou and Roberts, 1999], the simulation of shocks in pipe filling events [Vasconcelos et al., 2009], etc. However, these schemes are rarely utilized in GC models. The main drawback for numerical schemes based on the Riemann problem is the in- crease in complexity without a clear increase in accuracy. A method to determine the adequacy of a numerical scheme is to assess computational efficiency, the possibility of shocks, etc. How- ever, whenever explicit schemes are used, it is essential that the CFL condition is satisfied. The maximum time step that is allowed for this condition is computed from the following expression: Dt =Cr Dx juj+c (2.46) in which Cr is the Courant number. In order to sustain model stability when using explicit schemes, Cr should be less than or equal to one. When the time step is well below this limit and shocks exist, linear schemes suffer either from numerical diffusion or oscillations depending on the type of scheme [Godunov, 1959]. On the other hand, nonlinear numerical schemes are not affected by this issue. For example, in the simulation of pipe filling events, it may be possible that the air phase requires small Courant numbers in certain locations along the pipe. In order to sustain a stable simulation, the Courant 42 condition (2.46) must be satisfied and nonlinear schemes are much more effective [Vasconcelos et al., 2009]. Implicit schemes can overcome this difficulty; however, if the Courant number is increased beyond unity, the accuracy of the simulation is compromised when shocks are present. In this section of the thesis, some of the most commonly used linear and nonlinear explicit schemes are presented and examples of their implementation are provided. This discussion begins with linear schemes, which are presented in the Methodology subsection 4.4.1, all in the context of the Finite Volume Methods (FVM). The linear schemes that are subsequently presented are all Centred methods. The advantage of Centred schemes with respect to alternatives (e.g. Upwind schemes) is that ?explicit local information on wave propagation for constructing approximations to the numerical flux? is not required [Toro, 2001]. In other words for more complex flow events in which is difficult to determine the location and the time at which the information is generated, Centered-cell schemes are preferred. Upwind schemes are superior for other problems (e.g. a contaminant is discharged at an upstream boundary) in which it is clear where the information generated. In a large number of recent models, the FVM is used because of its ability to implement nonlinear schemes, which are outside the scope of the FDM. While the FVM method is imple- mented in this study (presented in 4.4.1) for both linear and nonlinear schemes, the use of the FDM in previous GC investigations is more frequent. These numerical schemes are implemented for mathematical models of the following conservative form: ?~U ?t + ? ~F(U) ?x = ~S(U) (2.47) where ~U is the vector of the conserved variables, ~F(U) is the flux vector and ~S(U) represents the source terms. For the FVM, the conserved variables are updated according to the following expression: The various FVM schemes propose expressions for the interfacial fluxes ~F(U)n+1=2i+1=2 in order to solve this update equation for the conserved variables. ~Un+1i =~Uni + Dt Dx ~ F(U)n+1=2i 1=2 ~F(U)n+1=2i+1=2 +Dt ~S(U) (2.48) 43 2.5.1 Linear schemes For GC flows, the most widely used linear schemes are presented: Lax-Friedrich (LxF): This 1st order accurate, linear and monotone scheme is an extension of the unstable FTCS (forward-time centred-space) scheme. The general approach for GC problems is to move on to second-order accurate schemes due to the problems with numerical diffusion that are exhibited in first-order linear schemes. In addition, the second-order accu- rate schemes are more accurate for smaller discretization sizes. The Lax-Friedrich scheme is provided in [Toro, 2001] in FVM format: ~ULFi+1=2 = 1 2 h ~ F(U)ni + ~F(U)ni+1 i + 12 DxDt h~ Uni ~Uni+1 i (2.49) Lax-Wendroff (LxW): The LxW scheme is second-order accurate and is solved with predictor and corrector steps. This scheme is widely used in hydraulic problems and is one of the best options if shocks are not present [Vasconcelos and Nunes, 2009]. However, in the presence of flow discontinuities, the LxW scheme produces spurious non-physical oscillations in the results. Despite this drawback, this scheme is widely used in the context of GC with the aid of artificial viscosity [Ungarish and Zemach, 2005], which reduces the magnitude of the oscillations. One of the drawbacks of this approach is that the amount of diffusion that needs to be introduced is not easy to estimate a priori in order to avoid excessive smearing [LeVeque, 1992]. The Lax-Wendroff two-step scheme in FDM format is provided in [Sturm, 2010] for the St. Venant equations. In the following relationship, the Lax-Wendroff scheme is rewritten in a general format that is also applicable to the 1-D SWE. 8 >>< >>: ~Un+1=2i+1=2 = ~Uni +~Uni+12 + Dt2Dxh ~F(U)ni ~F(U)ni+1i ~Un+1i =~Uni + DtDxh ~F(U)n+1=2i 1=2 ~F(U)n+1=2i+1=2i+Dt ~F(U)ni (2.50) 44 MacCormick: This second-order accurate linear scheme yields results that are similar to the Lax-Wendroff scheme. However, [Zoppou and Roberts, 2003] determined that that the MacCormick scheme is slightly more accurate in the context of the dam break problem. In that study the error analysis was performed for the following parameters: normalized depth and normalized velocity. [Zoppou and Roberts, 2003] determined that the error associated with the combination of the MacCormick scheme and artificial viscosity is almost double the amount of the MacCormick scheme alone depending on the damping parameter. Unfor- tunately, the presence of spurious oscillations in the presence of shocks generally requires artificial viscosity. Despite this drawback, the MacCormick scheme has been utilized in a large number of GC models. [Adduce et al., 2011] implemented the MacCormack scheme in FDM format for their two-layer SWE model based on the implementation of [Garcia and Kahawita, 1986] for the St. Venant equations, which are closely related to the 1-D SWE. 8 >>> >>>< >>> >>> : ~Upi =~Uni DtDxh ~F(U)ni ~F(U)ni 1i+Dt ~S(U)ni predictor ~Uci =~Upi DtDxh ~F(U)pi+1 ~F(U)pi i+Dt ~S(U)pi corrector ~Un+1i = ~Uni +~Uci2 step n+1 (2.51) Leapfrog: This numerical scheme differs from the MacCormick and LxW schemes in that it is second-order accurate in time but not in space. In addition, this scheme is not as widely used as the previous two second-order accurate linear schemes. However, the Leapfrog scheme was implemented successfully in the work by [Klemp et al., 1994]. However, because they end their simulation when the upstream moving jump contacts the physical boundary, the performance of the Leapfrog scheme for the lock-exchange problem is not completely clear. It is expected that this scheme does not perform well after the reflection 45 Figure 2.4: Two-wave structure for the HLL scheme presented in [Toro, 2001]. from the upstream physical boundary, which is discussed in chapter 5. Details of this numer- ical scheme is described in [Sturm, 2010] in FDM format for the St. Venant equations: ~Un+1i =~Un 1i Dt Dx ~ F(U)ni+1 ~F(U)ni 1 (2.52) 2.5.2 Nonlinear schemes For the nonlinear schemes, the focus is on two of the most widely used implementations of the Riemann problem: HLL and Roe. The use of nonlinear schemes in GC models is not frequent, and for the lock-exchange problem, no examples were found. Therefore, the following analysis is based on more traditional open channels flows and follows the formulation in [Toro, 2001]. HLL: The HLL scheme is widely used in gas dynamics and hydraulics. However, it has rarely if ever been implemented for the SWE in the context of lock-exchange GC flows. The nonlinear HLL scheme assumes a two-wave structure that is correct for purely 1-D problems (see Figure 2.4). The HLL flux (FHLL) occurs in the intermediate or star region between the two wave speeds, SL and SR. If shear waves or species equations are introduced, the HLL scheme is inadequate because the middle wave is ignored [Toro, 2001]. 46 In order to utilize the FVM update expression (2.48), the interfacial fluxes ( ~F(U)i+1=2) must be computed. For the HLL scheme, ( ~F(U)i+1=2) is determined from the following relation- ship: ~F(U)i+1=2 = 8 >>> >>>< >>> >>> : FL FHLL = SRFL SLFR +SRSL(uR uL)S R SL FR (2.53) where SL and SR are the wave speeds, the subscript L denotes the ith cell and the subscript R denotes the cell i+1. These wave speeds can be computed from the following relationship that was suggested by [Toro, 2001]: SL = uL cLqL SR = uR +cRqR (2.54) and qK is given by: qK = 8 >>> < >>> : s 1 2 (h +hK)h h2K if h > hK; 1 if h hK (2.55) There are a multiple options available in which to determine the variables in the intermediate or star region (h and u ). [Toro, 2001] recommended the Two-Rarefaction Riemann Solver: 8 >>< >>: h = 1g0 12(aL +aR + 14(uL uR) 2; u = 12(uL +uR)+aL aR (2.56) Roe: For the SWE, the first-order accurate Roe scheme is presented in [Toro, 2001], which forms the basis of the following formulation. The accuracy of the Roe scheme compares well to the HLL scheme [Zoppou and Roberts, 2003]; however, Roe?s scheme is generally 47 more computationally intensive [Vasconcelos and Nunes, 2009]. The nonlinear Roe scheme was originally applied to the Euler equations, but was first applied to the SWE by [Glaister, 1988]. For the Roe scheme, the quasi linear form of (2.47) is approximated with a linear system: ?~U ?t +A ?~U ?x = 0 (2.57) where A is the Jacobian matrix. To compute the Roe averages, there are two common methods for SWE models: Roe-Pike [Roe and Pike, 1984] and Glaister [Glaister, 1988]. For the Roe-Pike approach presented in [Toro, 2001], the Roe-averages are computed from the following relationships: (Roe Pike) 8 >>>> >>< >>> >>> : ?u = uL ph L +uR ph Rp hL +phL ?h =phLhR ?c = q 1 2(c 2L +c2R) (2.58) where ?u, ?h and ?c are the Roe averages for the velocity, the depth and the celerity, respectively. For the Glaister approach the only difference is in the calculation of ?u, which is presented below: (Glaister) ?u = uL ph R +uR ph Lp hL +phR (2.59) Once the Roe averages are determined, the eigenvalues ( ?lj) and eigenvectors ( ?R(j)) are com- puted from the following relationships, which are presented in [Toro, 2001]: ?l1 = ?u ?a ?l2 = ?u+ ?a ?R(1) = h 1 ?u ?c i ?R(2) = h 1 ?u ?c i (2.60) 48 In addition, the wave strengths are computed from the Roe averages: wavestrengths 8 >>< >>: ?a1 = 12 h Dh ?h?cDu i ?a2 = 12 h Dh+ ?h?cDu i (2.61) The final step in the Roe scheme is to compute the inter-cell fluxes so that the conserved variables can be updated. ~F(U)i+1=2 = 1 2 ~ F(U)ni + ~F(U)ni+1 12 2 j=1 ?ajj?ljj?R(j) (2.62) 49 Chapter 3 Knowledge Gap and Objectives Although there have been significant insights into GC flows in recent years, there are still several unanswered questions: The link between the SWE models with theoretical front conditions from [Benjamin, 1968] and [Shin et al., 2004] and the experimental results is not clear. [Klemp et al., 1994] states that the error from the model is due to the SWE assumptions and not from the theoretical front conditions. Since the SWE models do not account for interfacial mixing, which slows down the GC advance, the velocity is overestimated. In experimental studies, it is difficult to determine the depth of the nose due to interfacial instabilities (e.g. Kelvin-Helmholtz). [Shin et al., 2004] presents a solution for this problem where h is a function of (g, r, H). However, it is difficult to determine the density variation from experiments. For state of the art two-layer SWE models, two modeling approaches are required depending on the depth ratio (f 0:5 and f > 0:5). For larger initial depth ratios, the model is much more complex and involves cumbersome, explicit tracking of flow features. A simpler and general model that is valid for all depth ratios is obviously desired. Nonlinear numerical schemes are becoming more and more popular in gas dynamics and hydraulics; however, there use in GC flows is still incipient. Therefore, the advantage of these nonlinear schemes on GC is not well known. The BC at the front of the GC has been traditionally solved with a front condition and a characteristic equation. However, for more complex mathematical models (e.g. the two- layer SWE), the characteristic equations also become more complex [Rottman and Simpson, 50 1983] and [Klemp et al., 1994]. In addition, the continuity errors are oftentimes unaccept- able (especially when depth constraints are implemented at the GC front). A relevant, robust solution for GC front calculations that improves continuity errors would be a clear improve- ment. This work is linked to the development of a new numerical model based on the FVM to simu- late the propagation of GC in the environment that can describe all relevant stages of the propaga- tion up to the self-similar stage. The proposed model should rely on a simple implementation and on conservative solutions of the SWE to simulate waves and flow discontinuities. The only explicit tracking that is enforced in the proposed model is at the GC front, and in this work several front condition alternatives were tested and compared. Details of the numerical model development as well as the experimental procedure are presented in the Methodology (Chapter 4). 51 Chapter 4 Methodology The procedure for the experimental investigations and the numerical modeling developments is discussed in the subsequent text. Both components are constructed in such a way to compare the front trajectory and to develop velocity hydrographs near the middle of the tank. The development is based on the full depth lock-exchange problem for Boussinesq fluids in a rectangular tank. Examples of the notation and a sketch of the experimental apparatus are provided in Figure 4.1. 4.1 Experimental program Lock-exchange experiments were conducted in this work for three different density differ- ences e = 1%;2%and3%) with the following initial conditions: x0 = 76:2cm and h0 = 40:6cm. There was no initial or constant ambient motion in the experiments conducted for this study. In- stead of analyzing a large range of initial parameters (e.g./ x0, h0, Dr, etc.), the experimental contributions are based on the analysis of the MicroADV probe results. However, the trajectory of the nose is also obtained from high-definition digital cameras for a comparison with the SWE models. Several experiments were conducted in smaller scale tanks in order to develop a consistent experimental program that involved the mixing procedure and the gate removal. Once consistent results for the GC trajectory were obtained in the smaller tanks, the experiments were performed in the 9.14 m tank that was used in this work. Two to three duplicates were performed for each run to ensure that the results were consistent and properly obtained. 4.1.1 Physical apparatus construction The 9.144 m long acrylic tank that is utilized in this thesis aims to represent a section of the navigation channel of Mobile Bay in Alabama, USA (see Figure 4.2. The full scale model was 52 Figure 4.1: The experimental apparatus with stationary ADV probes and a video camera moving with the GC front. The initial GC depth h0 is equal to the ambient depth H. used in a separate work to investigate the exchange flows between the Mobile and Tensaw rivers and the Gulf of Mexico that occur within the bay. The rectangular acrylic tank represents a channel 9.754 m that connects two reservoirs (see Figure 4.3). For this thesis, the channel was separated from the rest of the physical model with two wooden barriers approximately 2 cm thick with rubber attachments at each end. After these barriers were locked in place by two metal clamps, the length of the channel was reduced by 0.305 m at each end. In order to perform the experiments, a third gate was constructed to separate the two fluids, which initiated the GC flow upon removal. The acrylic channel was constructed in four segments because the acrylic sheets are only 2.44 m in length. In each component a 15 cm wide sheet that represents the channel bottom was welded to two side sheets (40.6 cm wide). Therefore, the maximum depth that could be achieved in the channel is 40.6 cm. After the side pieces were welded on top of the bottom sheet, the effective channel width was 12.7 cm. Before the the acrylic sheets were welded together, openings were drilled into the acrylic. Then the acrylic sheets were tightened together with metal screws in order to increase the welding 53 Figure 4.2: The photograph on the left is an aerial view of Mobile Bay, Alabama. Source: Landsat program. To the right is a snapshot of the depth-averaged velocities (arrows) and the water levels (contours) predicted by the estuarine and coastal ocean model. The vertical line in the middle of the bay represents the shipping canal. Presented in [Chen et al., 2005]. adhesion and to provide structural support. In the front of the physical model on the outside of the tank, small acrylic pieces were used to connect the adjacent acrylic sheets (Figure 4.3). On the inside of the tank, adhesive tape was placed between the acrylic pieces with silicon caulking applied along the edges in order to obtain a strong waterproof seal. As a result, the entire channel formed a cohesive and waterproof unit. 4.1.2 Initial salinity measurements The difference in density can be obtained by dissolving salt in one side of the tank CITE or by using a different temperature CITE. If the primary difference between the fluids is temperature, then the colder fluid will flow underneath the warmer fluid. In this study salt is dissolved in one side of the tank, so the saltwater will flow underneath the ambient freshwater. When dissolving the salt into the denser fluid, it is imperative that the resulting density is accurately computed. The methods that have been implemented to determine the density are: salinity tables, polynomial equations, hydrometers and refractive index measurements. Although 54 Figure 4.3: Different components of the construction for the acrylic channel. 55 reference tables are readily available, this method is not the most accurate because the volume of the fluid is known only to a few percent. In addition, the salt is never pure and may vary from one brand to another. Moreover, salt may contain certain amounts of water and may not completely dissolve in the fluid. An alternative and more accurate approach is to determine the refractive index of the solution, which depends on the amount of salt. [Shin, 2001] utilized an electronic refractometer, which is usually designed for medical use, with standard tables in which the density is determined from the value in Brix. According to [Shin, 2001], the error associated with this approach is 0.1 Brix or 0.005 gcm 3. The smallest value for the density difference in their experiments was 0.1 gcm 3, so the maximum error was less than 5%. The refractive index also depends on the temperature of the fluid, which can be calibrated by the instrument. In order to further reduce these errors due to temperature differences, [Shin, 2001] allowed the fluids to reach room temperature before the simulation. However, because of the simplicity and accurate results, the difference in density for this work was estimated from a polynomial expression in www.csgnetwork.com/h2odenscalc.html, which was provided by the University of Michigan and NOAA: r =r0 +Ac+Bc3=2 +0:00048314c2 (4.1) where A = 0:824493 0:0040899T +0:000076438T 2 0:00000082467T 3 +0:0000000053675T 4 B = 0:005724+0:00010227T 0:0000016546T 2 r0 = 1000 1 (T +288:9414)508929:2(T +68:12963) (T 3:9863)2 in which A, B and C are polynomial parameters, r is the density of the denser fluid in kgm 3, T is the temperature in C and c is the concentration in mg=l. The initial salt concentration was deter- mined from this relationship, and the resulting density difference was re-evaluated using precision 56 Figure 4.4: Photograph of the leading edge of a GC that illustrates the procedure in which the front trajectories were measured. hydrometers. These hydrometers determine the density with an error of approximately 5% or less. In addition to providing reliable data, experiments were performed relatively quickly. 4.1.3 Front trajectory measurements One of the most common data recording techniques that has been used in past decades is to record the trajectory of the GC front [Huppert and Simpson, 1980]. Although the position and the velocity of the front has been accurately obtained in previous studies, newer technology allows for a more accurate description. In the present study, digital cameras were used to record the position of the front at 30 frames per second. This, however, is not the first time that this type of quality of video camera has been used for GC front trajectories [Adduce et al., 2011]. Measurements of the trajectory of the front are very important in the evaluation of numerical models. In order to obtain values for the position of the GC nose at various positions along the channel, a grid (4 in by 1 in) was drawn on the front of the tank to record the longitudinal position of GC fronts. The video camera was operated in a way that followed the front of the GC. Because the 57 video camera also recorded the time, the velocity distribution of the GC front for each experiment was obtained. 4.1.4 MicroADV measurements For the velocity hydrograph measurements, MicroADV probes sampling at 20 Hz were placed near the middle of the tank in order to determine the longitudinal velocities for the denser GC and the ambient fluid. The lower ADV probe was placed low enough to measure the GC at 10.2 cm from the channel bottom. The upper probe was placed at 33.0 cm from the channel bottom in order to measure the ambient velocities. Although the ADV probes can measure the velocity in three-dimensions, the main focus of the data analysis was related to the longitudinal velocities. In this analysis, it was important to position the MicroADV probes as close to vertical as possible in order to accurately measure the velocity in the longitudinal direction. The probes were held in place with a wooden support, and adhesive tape was applied to hold the them in the vertical direction. The wooden support contains a cantilever section that hangs over the channel so that the ADV probe can be placed in the center of the tank. The data was recorded at 20 Hz for about 250 seconds or when the GC front reflected back and forth twice. 4.2 One-layer SWE model As mentioned, the objective was to develop a robust and generalized Boussinesq gravity cur- rent model for the lock-exchange problem. It was decided to minimize the use of explicit track- ing of flow features (i.e. discontinuities appearing when f > 0:5). Using the one and two-layer SWE written in physically conservative format, it was possible to implement linear and nonlinear shock-capturing schemes within a finite volume framework in an attempt to better describe the propagation of gravity currents, as presented below. The one-layer SWE models (equation 4.2) are widely used for GC flowing into a deep am- bient. In an attempt to establish clarity and simplicity, the one-layer SWE are written in vectorial 58 conservative format: ?~U ?t + ? ~F(U) ?x = ~S(U) ~U = 2 64 h uh 3 75 ~F(U)= 2 66 4 uh (uh)2 h + 1 2g 0h2 3 77 5 ~S(U)= 2 640 0 3 75 (4.2) This mathematical model is solved within the finite volume framework by integrating the previous PDE system, yielding the following expression that performs the updates on the conserved variables for each computational cell i at the time step n+1: ~Un+1i =~Uni + Dt Dx ~ F(U)n+1=2i 1=2 ~F(U)n+1=2i+1=2 where ~F(U)n+1=2i 1=2 are inter-cell fluxes (not to confuse with ~F(U)), which are calculated according to various numerical schemes formulations. Details of such formulations are presented in [Toro, 2001] and are omitted here for brevity. A major limitation of one-layer model is the inaccuracy when the ambient flow depth is not large in comparison with the depth of the current [Ungarish, 2009]. Moreover, the one-layer SWE cannot capture the slumping characteristics of gravity currents with large depth ratios where the effects of the ambient velocity are more important. In order to appropriately simulate the advance of gravity current fronts, SWE models (one and two-layered) require the use of a front condition. 4.3 Two-layer SWE model For the one-layer model, the velocity of the ambient fluid is neglected in the internal cell calculations. However, if the fractional depth becomes large enough (particularly when f > 0:5), 59 then the ambient velocity becomes an important parameter. The mathematical model for the two- layer SWE in a Boussinesq system was formulated in [Rottman and Simpson, 1983]: ?h ?t + ? ?x(uh)= 0 ?u ?t +(1 2a)u ?u ?x +(1 b)g 0?h ?x = 0 (4.3) where a = hH h b = hH + u 2 g0H 1 hH 2 Although equation 4.3 can produce robust results for certain GC applications, the previous equation has two limitations. First it is not expressed in terms of physically conserved variables (u instead of uh). Second, while amenable to be implemented in a finite difference framework, it is not in the conservative format used in finite volume models (i.e. equation 4.2). Thus, in cases of flow discontinuities, shocks will be simulated with incorrect velocities. In this work we propose to further manipulate the previous expression using the continuity equation and the chain rule (shown below) in order to obtain conservative variables. ?q ?t = ?uh ?t = ?u ?t + ?h ?t The terms that cannot be presented in the fluxes of the conservative variables ( ~F(U)) are isolated as source terms, which yields the following expression for the linear momentum conservation equation: ?uh ?t + ? ?x (uh)2 h + 1 2g 0h2 ! = 2auh?u?x +bg0h?h?x (4.4) 60 Now, combination of the equation (4.4) with the continuity equation for the dense layer yields a system of PDEs that is expressed in conservative, divergent format: ?~U ?t + ? ~F(U) ?x = ~S(U) ~U = 2 64 h uh 3 75 ~F(U)= 2 66 4 uh (uh)2 h + 1 2g 0h2 3 77 5 ~S(U)= 2 66 4 0 A?u?x +g0hB?h?x 3 77 5 (4.5) where A and B are written as: A = 2uh 2 H h B = (1 h=H)2 + u=pg0H 2 (1 h=H)2 (4.6) In the above equations, A and B are components of the source terms presented by [Rottman and Simpson, 1983] written in such a way to demonstrate that the numerator and denominator in both expressions approach 0 as f!1 and h!H, since the current velocity (u) approaches 0 in such conditions. For h > 0:95H, special handling of the source terms is required. A straight forward inspection of equation 4.6 indicates that the terms A and B appearing in the source term of equation 4.5 become undefined as h=H!1 where u!0, which in turn renders the numerical solution unstable. For the depth gradient term, this causes no problems except at f = 1; however,~S = 0 at f = 1. In this scenario B = 1 because u = 0 due to the horizontal free surface. For the velocity gradient component, the fix is not as simple. Moreover, instabilities are generated for h=H ranging from 0.95 to 1.0. In order to overcome this dilemma, L?Hopital?s rule is applied as h=H!1 so that a limiting value for A is formulated: lim h / H!1 A = 4uh (4.7) In the proposed two-layer SWE model, this use of L?Hopital?s rule was applied for f > 0:95. Moreover, the velocity gradient term is negligible for h=H > 0:9, so the possible side effects are negated. Using this approach, the two-layer SWE model is able to accurately simulate GC flows of 61 all fractional depth. In this implementation, the parameter A for the velocity gradient component of the source term is provided in the following expression. A = 8 >>< >>: 2uh2 H h (h=H < 0:95) 4uh (h=H 0:95) (4.8) The source terms are evaluated numerically using a finite difference calculation, and in this implementation a simple, central difference scheme scheme was used. The two-layer SWE may be solved by implementing a numerical scheme in a FVM update formulation (equation 2.48). In summary, the advantage of the proposed mathematical model is the simplicity since the same numerical scheme is applied everywhere in the solution domain at all time steps, except of course at the boundary conditions. Thus, the only boundary conditions that are used in the computations are the physical boundaries constraining the flow and the GC front. There is no explicit tracking of shocks except at the leading edge of the GC. 4.4 Aspects of the numerical solution A brief discussion on the strategy for the numerical solution implementation for these hyper- bolic PDE systems is presented, which include the selected numerical schemes; boundary condi- tions, including a new formulation to compute flow at the gravity current nose; and aspects related to the numerical stability and time steps that are required in the proposed model. 4.4.1 Numerical schemes As mentioned, to perform updates in the conserved variables using finite volume method, it is necessary to calculate the interface fluxes. While the proposed model has been tested with several linear and nonlinear schemes, four schemes are selected to perform the comparison in this work: the first-order accurate, linear FORCE scheme; second-order accurate, linear Lax-Wendroff 62 scheme; and two first-order accurate, nonlinear Roe and HLL schemes. All of these schemes are widely applied in CFD problems, and a detailed analysis is provided in [Toro, 2001]. Lax-Wendroff (LxW): The linear, two-step LxW scheme was presented in subsection 2.5.2 for the FDM. For the FVM, the expression for this second-order accurate scheme is [Toro, 2001]: 8 >>> >>> < >>> >>>: ~Un+1=2i+1=2 = 12(~Uni +~Uni+1)+ Dt2Dx( ~F(U)ni ~F(U)ni+1) ~F(U)n+1=2 i+1=2 = 2 66 4 uhn+1=2i+1=2 (uh)2 h + 1 2g 0h2 n+1=2 i+1=2 3 77 5 (4.9) Modified FORCE: The original FORCE scheme discussed in [Toro, 2001] is an average 50%-50% of the interfacial fluxes from the LxW and LxF schemes. This work follows the idea from [Vasconcelos and Nunes, 2009] and adds more weight to LxW scheme (i.e. 70%) so that diffusion is further reduced without a noticeable increase in oscillations at shocks. This results in the modified FORCE scheme adopted here: ~F(U)n+1=2 i+1=2 = 1 2 0:3 ~F(U)n+1=2i+1=2 LxF +0:7 ~F(U)n+1=2i+1=2 LxW (4.10) HLL: The nonlinear HLL scheme has performed well in the more common dam break prob- lem [Vasconcelos and Nunes, 2009] and should provide accurate and computationally effi- cient solutions for lock-exchange GC flows. The primary difference between this formula- tion and other SWE models is that the gravity g is replaced with the reduced gravitational acceleration g0 in the flux calculations. The procedure for the HLL scheme is provided in subsection 2.5.2. For the HLL scheme, the expression for the interfacial fluxes ( ~F(U)i+1=2) 63 that is used in the update expression (2.48) is: ~F(U)i+1=2 = 8 >>> >>>< >>>> >>: FL FHLL = SRFL SLFR +SRSL(UR UL)S R SL FR in which the subscript L denotes the ith cell and the subscript R denotes the cell i+1. The description of wave speeds (S) and the fluxes (F) is provided in subsection 2.5.2. In the proposed model, the intermediate or star variables are computed with the Two-Rarefaction Riemann Solver, which was recommended by [Toro, 2001]. However, similar results were obtained with the Riemann Solver based on exact depth positivity. Therefore, it is assumed that the choice of approximate-state Riemann solver for the intermediate variables is not crucial for these GC flows. Roe: The nonlinear Roe scheme was presented in subsection 2.5.2. The main difference in this implementation from other SWE models is that the gravity (g) is replaced with the reduced gravity (g0) due to Boussinesq gravity currents. The proposed numerical model was tested with the Glaister [Glaister, 1988] and Roe-Pike [Roe and Pike, 1984] approaches to compute the Roe averages, and the results were almost identical. For the Roe scheme, the expression for the interfacial fluxes (~F(~U)i+1=2) that is used in the update expression (2.48) is: ~F(U)i+1=2 = 1 2 ~ F(U)ni + ~F(U)ni+1 12 2 j=1 ?ajj?ljj?R(j) where the wave strengths ( ?aj), eigenvalues ( ?lj) and eigenvectors ( ?R(j)) are provided in subsection 2.5.2. 64 4.4.2 Boundary Conditions: DC and MOC At the boundaries of the solution domain (either between a moving gravity front and a physical boundary or between two physical boundaries) numerical schemes cannot be used to determine the two flow variables U=[h;uh]T and boundary condition calculations are required. At physical walls the alternatives to solve the problem are generally straightforward. One way is to either enforce zero flow (u = 0) and solve for the depth h using a relevant characteristic equation [Sturm, 2010]; an alternative would be the use of a reflective boundary condition as presented in [Toro, 2001]. The other boundary condition calculation is at the GC front, and two alternatives are im- plemented for this moving boundary in the case of the one-layer SWE model: the traditional characteristic approach (MOC) and a proposed approach based on local continuity and momentum conservation (DC) (described ahead in this subsection). Both alternatives enforce a constraint on the depth when the ambient depth H decreases and approaches the lock-exchange problem (e.g. 0:347 H for Benjamin?s front condition) [Klemp et al., 1994] and [Ungarish and Zemach, 2005]. The MOC BC that is only implemented here for the one-layer model involves the combination of a C+ characteristic equation and the enforcement of a front condition (as explained in section 2.1). This approach provides two equations that allow for the solution of the two unknowns [hLE;qLE]T at the gravity current cell LE, which denotes the front or leading edge. Assuming that the velocity of cell 1 (i.e. the wall boundary) will always be zero, the resulting set of equations are: ucell:1 +2 p g0hcell:1 = 2 p g0hcell:1 = uLE +2 p g0hLE uLE = p g0hLEFr LE (4.11) which are solved iteratively (for this study, the Bisection method was used). One notices that the characteristic equation transports the Riemann invariant value u+2c from the upstream wall (cell 1) until the GC nose according to the characteristic velocity: u+c. In addition, the solution for [hLE;uLE]T allows for the determination of uhLE. The boundary condition strategy for the two-layered system could also have followed the same strategy, but the more complex structure of the characteristic equation (as shown in [Rottman 65 and Simpson, 1983]), and consequently Riemann invariants, motivated the development of an alternative approach. This alternative method explicitly enforces continuity, linear momentum, the kinematic condition and a front condition at the leading edge of the GC. This front region is assumed to span over cell LE, a finite volume cell fully occupied by the dense fluid, and over cell LE+1, which undergoes the advance of the current and is not yet completely filled (see Figure 4.5). Because the front region exists between two computational cells, this method was named Dual-Cell (DC), and it requires five equations to solve for the unkown parameters: h LE, uh LE, h LE+1, uh LE+1 and Dx LE+1. At the upstream front BC cell (i.e. LE), continuity is enforced along with the x-momentum equation. The continuity equation and a front condition is implemented at the downstream cell (i.e. LE+1). The continuity equations that were used in the SWE models are provided below: Continuity 8 >>< >>: dALE dt =Dx dhLE dt = uh LE-1 uh LE+1 cell: LE dALE+1 dt = d(Dx LE+1 hLE+1) dt = uh LE cell: LE+1 (4.12) in which A is the cell area (A =DxhLE for cell LE and A =DxLE+1hLE+1 for cell LE+1). The continuity equation is slightly more complex for cell LE+1 because both the depth hLE+1 and the distance DxLE+1 change with time (see Figure 4.5). The x-momentum equation that was used in the DC boundary condition is provided from the following expression: Fx = cs ?mLE+1uLE+1 cs ?mLEuLE (4.13) in which ?m = rc uh is the mass flow rate. The implementation of the momentum and continuity equations for the SWE models is provided below along with the other two equations (i.e. front 66 condition and kinematic condition), which make up the DC boundary condition: DC 8 >>>> >>> >>>> >>>> < >>> >>>> >>> >>>> >: I r c uh n+1LE uh n+1 LE+1 h n+1LE uh n+1LE-1 h n+1LE +F 2 F1 = 0 II h n+1 LE = h n LE + Dt Dx(uh n+1 LE-1 uh n+1 LE+1) III uh n+1 LE+1 = h n+1 LE+1 pg0h n+1 LE+1 FrLE+1 IV Dx n+1 LE+1 =Dx n LE+1 +Dt uh n+1LE+1 h n+1LE+1 V uh n+1 LE = 1 Dt h Dx n+1LE+1(h n+1LE+1 h nLE+1)+h n+1LE Dx n+1LE+1 Dx nLE+1 i (4.14) where n is the time step index, F1 is the upstream hydrostatic force at the interface between cells LE-1 and LE, and F 2 is the downstream hydrostatic force at the interface between cells LE and LE+1. These five equations represent respectively: 4.14 I the linear momentum balance at cell LE; 4.14 II the local continuity at cell LE; 4.14 III the enforcement of a front condition relating depth and velocity at cell LE+1; 4.14 IV the kinematic condition, which updates the front position Dx n+1LE+1, at cell LE+1; and 4.14 V the local continuity at cell LE+1. These equations are solved iteratively for the pair of cells at the edge of the computational domain until the advance of the GC front at cell LE+1 (i.e x n+1LE+1, given by u LE+1Dt) exceeds the cell size Dx. When this happens, the front of the GC enters a new cell and Dx LE+1 is reset to zero. This alternative to compute the boundary conditions is shown to be accurate for both one-layer and two-layer SWE models, and have compared well with experimental data collected in this work as well with previous investigations. 4.4.3 Stability and time step calculation Due to the choice of applying explicit schemes for the numerical solution, there is an up- per limit for the computational time step Dt established by the Courant-Friedrichs-Lewy (CFL) condition: 67 Figure 4.5: Schematic diagram of a GC that illustrates the cell structure for the DC BC at a given time step. The cell LE represents the last computational node that is completely filled by the front of the GC. Dt Dxmax(u+pg0h) (4.15) The enforcement of this condition is not sufficient to render the model stable. What was determined is that in some cases, even when equation 4.15 is satisfied, the conditions become unstable at the front. When the time step was reduced and the number of time steps to traverse cell LE+1 increased, stability is obtained in the computations. Thus, a second limit for time step is also enforced in which the velocity scales with the celerity: Dt kDt Dxpg0h LE (4.16) in which kDt is a constant factor used in the computations ranging between 0.05 and 0.5, depending on the initial depth ratio and the density difference. However, stability was typically achieved with kDt = 0:1. Later in the model development, a change was made to the force calculations of the boundary conditions so that the Courant condition governed the stability. As a general rule, the 68 relationship that was implemented to compute the time steps is: Dt min Dx max(u+pg0h);kDt Dxp g0hLE (4.17) 69 Chapter 5 Results and Discussion 5.1 Experimental results Front trajectories It is convenient to display the results for lock-exchange GC flows in terms of dimensionless units (see section 2.4.2). Furthermore, the distance is normalized by the initial current length (x0), the depth by the initial GC depth (h0), the velocities by the initial celerity (pg0h0), etc. (see [Ungarish, 2009]). This procedure allows for a comparison for various initial conditions and verifies consistent experimental results. In addition, other researchers can easily compare there data to the dimensionless results. For the experiments conducted in this study, the initial conditions and some of the dimensionless results (see section 2.2) are presented in Table 5.1. It is clear from the results that the inviscid assumption is justified (Rem >> 1;000). The consistency for FrH demonstrate the consistency in which the experiments were performed. Table 5.1: Experimental initial conditions and results for each experiment. Values for um were determined from a digital camera tracking the GC front. Run x0 h0 r1 g0 um Rem FrH m m kgm -3 ms -2 ms -1 1 0.762 0.406 1010 0.107 0.098 19952 0.468 2 0.762 0.406 1020 0.202 0.135 27657 0.472 3 0.762 0.406 1028 0.277 0.158 32423 0.473 In Figure 5.1 dimensionless GC front trajectory and velocity results are provided for the 9.14 m tank. After an initial reduction in velocity, a near constant velocity region forms denoting the slumping stage (see section 1.3). The velocity results for the three experiments are almost identical 70 Figure 5.1: Dimensionless trajectory and velocity for the GC front. (i.e. steady) except at the beginning of the experiment where the removal of the gate is very impor- tant and could not be absolutely equal between repetitions. Likewise, the change in position with time reveals that the velocity is almost constant (actually very slowly decreasing) within the slump- ing stage. These results suggest that the experiments were conducted with sufficient accuracy to yield consistent results. MicroADV results MicroADV results were obtained in the 9.14 m segment of the channel at one location for all of the experiments: x = x=h0 = 7:68. The depth (h0) for the experiments was 40.6 cm and the initial lock length (x0) was 76.2 cm. Thus, the thin layer assumption for the shallow water theory is valid. For h0=x0 < 1, the mixing between the GC and the ambient fluids is not significant, and the majority of breaking waves are confined to the interface [Hacker et al., 1996]. The results for 71 the MicroADV probe are displayed in terms of dimensionless units in Figure 5.2. Thus, the data for all three experiments should fall along the same path. The results for the MicroADV probe provide a clear representation of the internal velocity structure for lock-exchange experiments. It is apparent from Figure 5.2 that pressure pulses are consequences of surface waves generated from the gate removal. An oscillatory pattern forms from this wave celerity; however, there is a clear interface caused by the arrival of the GC front. As the time increases, the pressure pulses are reduced. Furthermore, a parabolic decrease in velocity occurs after the GC front moves past the ADV probe until the end wall reflection generates negative velocities. The one and two-layer SWE models are able to simulate the GC propagation beyond the end wall reflection so that the velocity magnitudes are compared. The negative peaks that occur when t 40 represent the front of the reflected GC from the far wall so that the GC travels back towards its initial location. After the negative velocity front arrives at the MicroADV probe, one notices that the velocity quickly changes back to positive values before a negative curve develops. Thus a dense solitary wave was separated from the rest of the flow at the initial reflection with the far wall. This solitary wave was generated from a hydraulic jump, which was diminished due to the stably stratified layer that was left behind the GC nose from turbulent mixing [Simpson, 1997]. 72 Figure 5.2: Velocity hydrograph comparison for 1%, 2%, and 3% salt in the 9.14 m tank. 73 5.2 Numerical model results The results for the one and two-layer SWE models are illustrated in Figure 5.3 for a wide range of fractional depth. In Figure 5.3 f = 0:01, f = 0:5 and f = 1:0 for the deep ambient, critical condition and full depth scenario, respectively. The DC BC was implemented at the front of the GC with Benjamin?s front condition and 400 discretization nodes. Unless otherwise noted, kDt = 0:1 for the full depth and critical condition. For the deep ambient problem, kDt = 0:2 in order for the code to remain stable. It is clear from figure 5.3 that the two-layer SWE model approaches the one-layer alternative for the deep ambient scenario. The depth profile and the GC trajectory are almost identical at each point in time. In the deep ambient scenario, the velocity of the ambient fluid is much smaller in magnitude thus less momentum is transferred to the GC. However, as the depth ratio begins to increase, the two SWE models begin to diverge from one another. The resistance of the ambient fluid becomes more important, and the two-layer model is the better alternative [Ungarish, 2009]. The decrease in accuracy of the one-layer SWE model for smaller ambient depths is due to the increased importance of the ambient velocity. As the fractional depth increases from the deep ambient to the full depth lock-exchange problem, the velocity increases for both the one and two- layer models. Moreover, for the deep ambient scenario, a larger pressure from the ambient fluid is being exerted on the GC, which corresponds to a greater horizontal velocity. One of the interesting features in Figure 5.3 is that the GC in the two-layer model increases in velocity relative to the one-layer model as the depth ratio increases to the critical condition and then to the full depth lock-exchange problem. The reason for this discrepancy is due to the increased accuracy of the two-layer model in computing the length of the slumping stage, which is underestimated with the one-layer model. The two-layer model simulates the reflection with the upstream boundary with a sharp reflected bore, which eventually reaches the front. For the one-layer model, this reflection generates a wave, which reaches the GC front much faster, so the depth of the front and the resulting velocity are decreased quicker than in the two-layer model. 74 Although though the GC front in the two-layer model is constrained because of the more complex characteristics, the depth and velocity of the front are greater than in the one-layer model for the majority of the slumping stage. This constraint has been validated by high-resolution mod- els that suggest the GC front rarely moves faster than the c+ characteristic for Boussinesq fluids [Ungarish, 2009]. Because the two-layer model computes a longer slumping stage for the shallow ambient condition, the GC propagates at a faster. In table 5.2 the sensitivity of the DC boundary condition is analyzed in the context of the number of discretization nodes. For the one-layer model, the DC boundary condition is compared to the traditional MOC boundary condition. The results indicate that the DC BC is more accurate for simulations with smaller discretization nodes. However, both models converge in accuracy at approximately 400-800 nodes. In previous two-layer models using the shock-tracking approaches, a smaller number of dis- cretization nodes was often used to resolve the flow discontinuities. The discontinuous interfaces for these models are assumed vertical, so the numerical scheme does not need compute this feature even though post shock oscillations are a potential side effect [Ungarish, 2009]. By comparison, the shock-capturing approach proposed in this thesis requires more sophisticated numerical schemes and improved grid discretization in order to accurately reproduce the discontinuous interfaces. A comparison between these two approaches is presented in Figure 5.11. Table 5.2: Sensitivity of each modeling approach to the number of discretization nodes. The error for each model is compared to its computed propagation time at 1,000 nodes. The Roe scheme is used with B front condition and kDt = 0:1. Model Parameter Number of nodes 50 100 200 400 800 1000 GC propag. time (s) 77.6 71.1 68.1 68.0 67.2 67.31L: MOC Continuity error (%) +15.3 +5.5 +1.2 +1.0 -0.1 0 GC propag. time (s) 66.9 65.6 64.6 65.5 65.5 65.81L: DC Continuity error (%) +1.7 -0.3 -1.8 -0.4 -0.4 0 GC propag. time (s) 58.4 56.2 55.1 55.4 55.4 55.52L: DC Continuity error (%) +5.3 +1.4 -0.6 -0.1 -0.1 0 75 Figure 5.3: GC propagation for the one and two-layer SWE over a wide range of fractional depth: f = 0:01, 0:5 and 1:0. The DC BC is used with 400 discretization cells. 76 5.3 Comparison between experimental results and numerical predictions The assessment of the proposed numerical solutions for the one and two-layer SWE with regards to the experimental results is done in the following comparisons: 1. Experimental measurements of the GC front trajectory are compared to the results obtained from the one and two-layer models using the Roe scheme; 2. Experimental results for the overall GC propagation time are compared to numerical predic- tions with different numerical schemes, front conditions, front BC strategies (for the one- layer model) and time steps (for the two-layer model). Continuity errors in each alternative are also compared; 3. Velocity measurements obtained with MicroADV probes for both the denser GC and the upper ambient fluid are compared to the predicted velocity hydrographs obtained with the one and two-layer SWE models; Trajectory results Experimental results for the trajectory of the GC front in the 9.14 m tank were obtained from high definition digital cameras, recording at 30 frames per second. The advance of the front was measured against a grid placed on the acrylic walls, and typical results are presented in Figure 5.4. The GC flow was caused by the removal of a gate located at x0 = 0:762m, which separated saltwater from freshwater. The depth of both fluids prior to the gate removal was 0.406 m, and the relative density difference was approximately 2%. As described in [Simpson, 1997], the advance of the GC front is quasi-steady within the slumping stage, as one notices that the slope of the trajectory varies only slightly with time. Figure 5.4 presents numerical predictions for the GC front trajectory obtained with the one- layer SWE model using the Roe scheme to solve flow in internal cells associated with various front conditions. The DC front boundary condition equations were used in a 400-cell computa- tional domain. In the chart, the acronyms ?B?, ?HS? and ?UZ? represent for the front conditions 77 Figure 5.4: Front trajectory results for the one-layer SWE model using Roe scheme and DC BC with various front conditions and from experiments. proposed by Benjamin [Benjamin, 1968], Huppert and Simpson [Huppert and Simpson, 1980] and Ungarish and Zemach [Ungarish and Zemach, 2005], respectively. Results presented in this figure for the coordinate of the front are normalized by the original position of the gate x0, while the time was normalized by x0=pg0h0. One notices that the best prediction obtained with the one-layer model was the one obtained with Benjamin?s front condition followed by Huppert and Simpson. In general, the tested numerical alternatives were able to capture well the general trend of the experimental trajectory of the GC front. A similar comparison, this time obtained with the two layer model is presented in Figure 5.5; the options for the numerical solution are the same as for the one-layer model described above. Results are now reversed in that the predictions with Benjamin?s front condition were the worst, while results obtained with HS and UZ front conditions simulated the experimental measurements very well. 78 Figure 5.5: Front trajectory results for the two-layer SWE model using Roe scheme and DC bound- ary conditions with various front conditions and from experiments. Effects of the boundary conditions and numerical schemes The effect of the numerical scheme selection, boundary condition solution strategy (for one- layer model), time step size (for two-layer model) and front condition choice to the accuracy of the model predictions are compared and summarized in Tables 5.3 and 5.4 for the one and two-layer SWE models, respectively. The assessment of the model accuracy is measured using three criteria: 1) the difference between the experimental propagation time and the respective numerical predic- tion; and 2) the continuity (mass balance) errors; and 3) the computational efficiency involved in the calculations of the GC flow. Numerical results, all obtained with a 400-cell computational do- main, are compared to the experimental condition (x0 = 76:2cm; h0 = 40:6cm;Dr = 0:02gcm 3), but results for larger density differences follow the same general trend. One recalls that the application of the one-layer model is not recommended for the conditions used in experiments due to the large ratio between the GC depth and the ambient depth, particu- larly at the initial stages of the flow. This condition may be one explanation why the continuity errors obtained for this model associated with the traditional characteristic equation are too high, 79 regardless of the front condition, as it may be seen in Table 5.3. For this scenario, the constraint based on the energy conservation theory (f 0:5), causes an increase in the continuity error for the characteristic BC. The use of the alternative DC boundary condition equations, however, seems to address most of the issues with continuity errors. As far as the velocity of the GC front, the best results were obtained with Benjamin?s front condition for the one-layer model. In general, the worst numerical scheme with regards to the accuracy of the GC propagation time and continuity error is the second-order LxW scheme, possibly due to the artificial viscos- ity that must be introduced in the scheme so that oscillations are controlled. The best numerical scheme for this condition is the modified FORCE, but both non-linear schemes tested also pro- duced accurate results. The nonlinear Roe and HLL schemes performed well for all conditions that were tested. In addition, these nonlinear schemes are more stable for larger density differences and partial depth releases and perform better for smaller time steps than the linear schemes (Table 5.4). In Figure 5.6 the results with the Roe and HLL schemes are superimposed, and the depth profile and front trajectory are almost identical. The conditions tested in the experiments are more appropriately simulated with the use of a two-layer SWE formulation even though the front trajectory is also simulated accurately with the one-layer model. In general continuity errors are all under 1.4% for the two-layer model (Table 5.4), and the largest divergence between the measured and predicted propagation time is 18%. This error was expected because [Shin et al., 2004] stated that Benjamin?s theory over-predicts the propagation time by about 20%. As discussed in the Methodology chapter, a boundary condition for the GC front based on the characteristic equation was not implemented for the two-layer SWE model. Thus it was decided here to compare the effects of changing the time step size by using two values for the coefficient kDt presented in equation 4.17. The best agreement for the two-layer model was observed for the case when Huppert and Simpson?s front condition was used in terms of the GC front propagation time. With regards to the continuity error, both HS and UZ front conditions presented continuity errors consistently smaller than 1% while the results using Benjamin?s front condition were consistently the worst. There was 80 Figure 5.6: Slumping for the one and two-layer SWE models using the HLL and Roe scheme with two front conditions: Benjamin (B) and Huppert Simpson (HS). x0 = 1, h0 = 1 and e = 2% using 800 computational cells. 81 a slight increase in the continuity errors for the smaller kDt tested. All two-layer numerical schemes tested presented similarly accurate results, with a slight advantage to the non-linear schemes. The use of the LxW scheme resulted in instabilities that deemed this shock-capturing model unstable, and thus it was decided not to present the results here. Based on the error computed for the propagation time for the two-layer model, which was almost always greater in magnitude than 6%, an additional front condition was implemented. [Rottman and Simpson, 1983] proposed an approach based on Benjamin?s condition that allowed for a calibration constant. From their partial-depth experiments, they calibrated this parameter b to equal 1.0. However, this term changes for the full depth lock-exchange problem. Therefore, b is re-calibrated in order to determine a value that satisfied the experiments for this work (This condition is referred to as MRS). The results are summarized in Table 5.5, and b = 1:215 (see equation 2.7) provided the best results. For the full depth lock-exchange problem, it is impor- tant to note that the depth and velocity constraints at the GC front change depending on the front condition (see subsection 2.4.2). The front constraints for the MRS condition are: hLE 0:414H, uLE=pg0h0 0:448, respectively. Moreover, the error with experiments for the two-layer SWE model using the MRS front condition is constistently less than 1%. According to Table 5.3, the one-layer SWE are very accurate in conjunction with Benjamin?s theoretical front condition. It is speculated that the reason for these good results is that the errors associated with interfacial mixing/entrainment balance with the errors associated with neglecting the ambient velocity. In the actual GC propagation, the ambient velocity transfers momentum to the GC, which increases the GC velocity. On the other hand, turbulent entrainment, which is most prevalent at the GC front [Simpson, 1997], halts the advance of the GC. [Ungarish, 2007] similarly obtained accurate results for the full depth release using the one-layer SWE even though these equations are least valid for this depth ratio. When the two-layer SWE model is implemented, this balance of errors is changed because the ambient velocity is taken into account within the internal cell calculations. Therefore, because turbulent entrainment is neglected, the errors are greater compared to the one-layer model. 82 Table 5.3: Comparison between lock-exchange experiments and one-layer SWE model predictions using both BC approaches at the GC front, kDt = 0:1 and a 400-cell solution domain. The time of propagation for the GC in the experiments was 93.67 s (e = 1:0%) and 67.63 s (e = 2:0%). Experiment: e = 1:0% DC BC equations Characteristic BC equations B LxW FORCE HLL Roe LxW FORCE HLL Roe GC propag. time (sec) 91.54 89.90 90.21 89.87 92.43 91.74 92.02 93.71 Error w/ experim. (%) -2.3 -4.0 -3.7 -4.1 -1.3 -2.1 -1.8 +0.0 Continuity error (%) -5.0 -2.0 -2.1 -2.4 -8.3 -7.4 -7.7 -8.0 HS LxW FORCE HLL Roe LxW FORCE HLL Roe GC propag. time (sec) 108.39 107.42 107.49 107.49 109.94 109.66 109.76 109.63 Error w/ experim. (%) +15.7 +14.7 +14.8 +14.8 +17.4 +17.1 +17.2 +17.0 Continuity error (%) -7.4 -2.0 -2.3 -2.3 -12.9 -11.4 -11.8 -10.8 UZ LxW FORCE HLL Roe LxW FORCE HLL Roe GC propag. time (sec) 122.02 120.09 119.95 120.19 125.18 124.68 124.81 124.56 Error w/ experim. (%) +30.3 +28.2 +28.1 +28.3 +33.6 +33.1 +33.2 +33.0 Continuity error (%) -7.3 -2.0 -2.1 -2.3 -12.6 -11.0 -11.6 -10.7 Experiment: e = 2:0% DC BC equations Characteristic BC equations B LxW FORCE HLL Roe LxW FORCE HLL Roe GC propag. time (sec) 66.55 65.45 65.83 65.43 67.23 66.75 66.93 68.01 Error w/ experim. (%) -1.6 -3.2 -2.7 -3.3 -0.1 -1.3 -1.0 +0.6 Continuity error (%) -4.88 -2.03 -3.47 -2.23 -8.28 -7.43 -7.65 -8.18 HS LxW FORCE HLL Roe LxW FORCE HLL Roe GC propag. time (sec) 78.82 78.12 79.06 78.14 79.96 79.78 79.83 79.76 Error w/ experim. (%) +16.5 +15.5 +16.9 +15.5 +18.2 +18.0 +18.0 +17.9 Continuity error (%) -7.33 -2.04 -2.42 -2.22 -12.89 -11.51 -11.83 -11.02 UZ LxW FORCE HLL Roe LxW FORCE HLL Roe GC propag. time (sec) 88.73 87.30 87.56 87.32 91.05 90.68 90.77 90.62 Error w/ experim. (%) +31.2 +29.1 +29.5 +29.1 +34.6 +34.1 +34.2 +34.0 Continuity error (%) -7.27 -2.05 -2.10 -2.23 -12.63 -11.12 -11.58 -10.88 83 Table 5.4: Comparison between lock-exchange experiments and two-layer SWE model predictions for different time steps and front conditions using the DC boundary condition. The GC time of propagation in the experiments was 93.67 s (e = 1:0%) and 67.63 s (e = 2:0%). Experiment: e = 1:0% kDt = 0:1 kDt = 0:05 B FORCE HLL Roe FORCE HLL Roe GC propag. time (sec) 76.15 76.15 76.14 76.15 76.15 76.15 Error w/ experim. (%) -18.7 -18.7 -18.7 -18.7 -18.7 -18.7 Contintuity error (%) +1.0 +1.3 +1.3 +1.5 +1.3 +1.4 HS FORCE HLL Roe FORCE HLL Roe GC propag. time (sec) 100.26 99.65 99.63 99.56 99.68 99.65 Error w/ experim. (%) +7.0 +6.4 +6.4 +6.3 +6.4 +6.4 Contintuity error (%) +0.5 +0.6 +0.6 +0.8 +0.7 +0.7 UZ FORCE HLL Roe FORCE HLL Roe GC propag. time (sec) 107.48 106.36 106.28 106.28 106.39 106.32 Error w/ experim. (%) +14.7 +13.5 +13.5 +13.5 +13.6 +13.5 Contintuity error (%) +0.5 +0.6 +0.6 +0.8 +0.7 +0.7 Experiment: e = 2:0% kDt = 0:1 kDt = 0:05 B FORCE HLL Roe FORCE HLL Roe GC propag. time (sec) 55.38 55.38 55.38 55.39 55.39 55.38 Error w/ experim. (%) -18.1 -18.1 -18.1 -18.1 -18.1 -18.1 Contintuity error (%) +1.4 +1.3 +1.3 +1.6 +1.4 +1.4 HS FORCE HLL Roe FORCE HLL Roe GC propag. time (sec) 72.46 72.46 72.19 71.50 72.48 72.22 Error w/ experim. (%) +7.1 +7.1 +6.7 +0.57 +7.2 +6.8 Contintuity error (%) +0.7 +0.6 +0.7 +1.0 +0.7 +0.8 UZ FORCE HLL Roe FORCE HLL Roe GC propag. time (sec) 77.36 77.32 76.77 75.58 77.36 76.80 Error w/ experim. (%) +14.4 +14.3 +13.5 +11.8 +14.4 +13.6 Contintuity error (%) +0.7 +0.6 7 +1.0 +0.7 +0.8 84 Table 5.5: Two-layer SWE results using various values for b for the RS front condition with 400- cell solution domain. When b = 1:215, this front condition is referred to as MRS. The GC time of propagation in the experiments was 93.67 s (e = 1:0%) and 67.63 s (e = 2:0%). Experiment: e = 1% b based on equation 2.7 1.0 (RS) 1.1 1.2 1.215 1.3 p2 (B) GC propag. time (s) 131.35 112.81 95.73 93.33 83.64 76.16 Error w/ experim. (%) +40.2 +20.4 +2.2 +0.4 -10.7 -18.7 Continuity error (%) +0.51 +0.55 +0.70 +0.73 +0.91 +1.29 Experiment: e = 2% b based on equation 2.7 1.0 (RS) 1.1 1.2 1.215 1.3 p2 (B) GC propag. time (s) 95.53 82.08 69.45 67.88 60.83 55.40 Error w/ experim. (%) +41.3 +21.4 +3.0 -0.4 -10.1 -18.1 Continuity error (%) +0.51 +0.55 +0.70 +0.73 +0.91 +1.29 The computational time associated with the for SWE model alternatives is an important fea- ture that is compared in Table 5.6. This comparison is made between the one and two-layer models and between the DC and MOC boundary conditions for the one-layer model. The first section of Table 5.6 displays the computational time for the GC to travel the full length of the tank with the following initial conditions: x0 = 1m, Ltank = 10m, h0 = 1m and e = 2%. The second section of Table 5.6 compares the different modeling approaches. The results indicate that the DC BC is consistently more computationally efficient than the MOC BC. It is apparent from these results that the computational time is strongly dependent on the boundary condition solution strategy for the GC front. Also, the source terms that are present in the two-layer model have an important effect on the computational time (20%-50%, see Table 5.6). The two-layer model combined with the MOC BC and explicit tracking of the discontinuities is expected to provide even worse com- putational times as the one-layer model for comparable domain discretizations. Therefore, the two-layer SWE shock-capturing model proposed in this thesis is potentially attractive in terms of computational efficiency, fairly low continuity errors and above all conceptual simplicity and ease of implementation. It was found that the selection of the numerical scheme has a minor importance on the overall accuracy. It was determined that the nonlinear schemes are more stable for partial depth releases 85 Table 5.6: Computational time comparison between the one-layer (1L) and two-layer (2L) SWE models and between the DC and MOC front BC strategies for various discretization sizes. Num. of cells 100 200 400 800 2000 1L: MOC (s) 1.37 5.20 26.93 88.78 179.00 1L: DC (s) 1.02 2.07 2.62 9.58 36.32 2L: DC (s) 1.40 2.77 3.17 13.28 52.57 1L: MOC/DC 1.3 2.5 10.3 9.3 4.9 DC: 2L/1L 1.4 1.3 1.2 1.4 1.5 and for larger density differences. There was severe diffusion for the first-order Lax-Friedrich scheme. In addition, the second-order accurate linear schemes tested were unable to remain sta- ble for the two-layer SWE model in the context of the full depth lock-exchange problem due to excessive oscillations that form during the reflected bore. In previous models these second-order accurate schemes performed adequately because a shock-fitting procedure was used. In shock- capturing models, the choice of numerical scheme is more important since the discontinuity is not explicitly tracked. The most important factors in the simulation of GC propagation compared to lock-exchange experiments were: 1) selection of the appropriate mathematical model for the numerical solu- tion (best was two layer SWE); 2) selection of the appropriate front condition (modified RS for two-layer model); and 3) the selection of the nose BC solution strategy (most efficient was DC equations). Velocity measurements and predictions A comparison between the predicted and measured velocities at an intermediate point in the tank are presented in Figure 5.7. Velocity values are normalized by pg0h0 while time values are normalized by xo=pg0h0 (the dimensionless variables are denoted with asterisk). Velocity measurements were performed with the MicroADV probes sampling at a frequency of 20 Hz. The oscillations observed at the beginning of the velocity measurements were caused by the removal of the gate separating fresh and salt water. Suddenly, at about t = 15 the jump in the velocity 86 values corresponds to the arrival of the front at the ADV station, at x = x=h0 = 7:68. Velocity measurements indicate a gradual decline in the velocity after an initial peak that was caused by the arrival of the GC front. Measurements were continued after the gravity current advanced the entire length of the tank in order to observe the reflections. As the reflected current reached the MicroADV probes, the velocity becomes slightly negative, which was an expected outcome. Figure 5.7a) presents a comparison between the experimental measurements of the velocity with the predictions from the one-layer SWE model using the Roe scheme, DC boundary condi- tions and either Benjamin?s or Huppert and Simpson?s front conditions. These front conditions represent the theoretical implementation and the most accurate empirical condition that was tested for the full depth release. One notices good agreement between both simulations and the experi- mental data, with Benjamin?s front more accurately predicting both the arrival of the front and the reflected wave. The predicted velocity decay during the motion of the gravity current is milder than the measurements. Moreover, for the one-layer model the use an empirical front condition is not as important since the theoretical formulation in [Benjamin, 1968] is provides accurate results. A comparison between the predicted GC velocity using the two-layer SWE model and the experiments is presented in Figure 5.7b). As in the previous case, Roe scheme and DC boundary conditions were used along with Benjamin?s and Huppert and Simpson?s front conditions. From recent studies [Ungarish and Zemach, 2005], it is well known that the HS front condition performs well for the full depth lock-release, which is reinforced in 5.7b). On the other hand, Benjamin?s condition overestimates the GC front velocity by about 20% [Shin et al., 2004]. Therefore, one would expect that the HS front condition would provide a better estimate for the GC velocity structure. However, the magnitude (specifically at the velocity peaks) is better represented by Benjamin?s condition even though the front advance is better approximated by the HS condition. These results suggest that there is intense mixing and entrainment at the GC front (stated in [Simpson, 1997]), which decrease the velocities from their theoretical value. Although there have been many explanations regarding the inaccuracy of Benjamin?s theory [Klemp et al., 1994] and [Shin et al., 2004], Figure 5.7b) provides a clear understanding of the source of these problems in 87 the context of the velocity structure. The results in 5.7 indicate that the best agreement is observed with the HS front condition and the two-layer model. The results with the MRS front condition are presented in Figure 5.8 for the two-layer model to compare with two experiments: e = 1% and e = 2%. Although the HS front condition performed well in predicting the GC velocities, the results suggest that the MRS front condition performs best at predicting the arrival of both GC fronts: initial propagation and first reflection. It performs a little better than the HS front condition at predicting the peak velocities at the front but not as well as the B front condition. All things considered for the two-layer model, the MRS front condition is the best method tested for the Boussinesq, full depth lock-exchange problem. In Figure 5.9 the experimental results (Dr = 3:0%) for the upper MicroADV probe are com- pared to the two-layer SWE model predictions. The MRS and B front conditions are used with 200 computational cells. As seen in Figure 5.7b), the B front condition overestimates the arrival of both fronts: initial propagation and first reflection. The peak velocities compare well to both front conditions tested but are better represented by the MRS front condition. In addition, the MRS condition performs very well at predicting both GC front locations: t 15 and t 40. After the arrival of the GC front (t 20 30 in Figure 5.9), the ambient velocities computed with the two-layer SWE model underpredict the MicroADV results. 88 Figure 5.7: Velocity hydrographs with B and HS front conditions and 200 computational cells. The experimental data (e = 2:0%) is provided from MicroADV devices in the slumping stage: x = 7:68. 89 Figure 5.8: Velocity hydrographs with the MRS front condition and 200 computational cells. The experimental data is provided from MicroADV devices at x = 7:68. 90 Figure 5.9: Comparison between measured ambient velocities and two-layer SWE predictions using B and HS front conditions with 200 computational cells. The experimental data (e = 3:0%) is provided from MicroADV devices at x = 7:70. 5.4 Comparison with previous models This subsection compares the results of the proposed model to the ones obtained by previous models. Two models were selected for this comparison, both of these using two-layer formulations: Rottman and Simpson [Rottman and Simpson, 1983] and Ungarish and Zemach [Ungarish and Zemach, 2005]. The former model successfully simulated gravity current flows for partial depth releases (for f 0:5), but had limitations in representing the back-propagation of the wave during the slumping stage of gravity currents. The second model presented by [Ungarish and Zemach, 2005] is able to successfully simulate flows for cases when f > 0:5, but requires explicit tracking of the upstream moving hydraulic jump and the reflected bore. As it is shown below, the proposed model is able to simulate flows for any value of f without needing to perform explicit tracking of any flow feature except the GC front. 91 Figure 5.10 presents a comparison between the proposed model and the Rottman and Simpson model results for the case when f =0:5. The results presented are the profiles of the GC at different times along the simulation. All profiles are presented in a non-dimensional fashion, with depth normalized by h0, length by x0 and time by x0=pg0h0. [Rottman and Simpson, 1983] presented results using Benjamin?s front condition for a qualitative analysis; they obtained better agreement with experiments using their front condition. The proposed model results were obtained using Roe?s scheme and Benjamin?s front condition. One notices an excellent agreement between both models. A more challenging simulation is the case when f > 0:5 due to the formation of the backward moving feature resembling a discontinuity in the depth. A simulation for the case when f = 1 was performed with the proposed model in the conditions used in the work by [Ungarish and Zemach, 2005], and the results of both model?s GC profiles are presented in Figure 5.11. Two different modeling alternatives are presented for this comparison: Figure 5.11a) presents results obtained with the Roe scheme; and Figure 5.11b) presents results obtained with the Modified FORCE scheme. In both cases the HS front condition was used with the limiting value hLE = 0:427H at the GC front, following [Klemp et al., 1994]. Both modeling alternatives were successful in reproducing the model results presented by [Ungarish and Zemach, 2005]. Results obtained with the MFORCE scheme were slightly sharper (i.e better at simulating the vertical interface of shocks) than the results obtained with Roe scheme. For the Roe scheme, which performed similarly to the HLL scheme, Glaister?s approach [Glaister, 1988] was used to compute the Roe averages; however, the same results were obtained with the Roe-Pike approach (see [Toro, 2001]). One also notices that the few high frequency spurious oscillations presented in [Ungarish and Zemach, 2005] model (which uses LxW scheme) are not present in the proposed model results. In the previous discussion, the two-layer SWE model results were presented for three different initial depth ratios: f = 0, f = 0:5 and f = 1. In the following comparison, front conditions are tested for a larger number depth ratios in the context of the dimensionless front velocity. The results are presented in Figure 5.12 and were obtained using Roe?s scheme and three front conditions: B, 92 Figure 5.10: Comparison of GC profiles at various time steps between the proposed two-layer model and [Rottman and Simpson, 1983] model for the critical condition: f = 0:5. 93 Figure 5.11: Dimensionless slumping comparison between the proposed model and [Ungarish and Zemach, 2005] model using Roe scheme and modified FORCE scheme. 94 Figure 5.12: Dimensionless front velocity results for the two-layer model over a large range of fractional depth. HS and UZ. One notices that as the GC propagates into shallower ambient depths (thus larger f values) the velocity normalized by pg0h0 decreases. Results obtained with Benjamin?s front simulation were the largest, while results obtained with Huppert and Simpson?s front condition had an intermediate maximum at about f = 0:10. The results obtained with the UZ front yielded the smallest GC front propagation velocity. These results are compared with the numerical results presented by [Ungarish and Zemach, 2005] and with the experimental results by [Rottman and Simpson, 1983], and the best agreement with experimental results from the proposed model was observed for the UZ front condition. 95 A final comparison is made between the Environmental Fluid Dynamics code (EFDC) and the two SWE models (with the HS front condition), which is presented in [Hatcher et al., 2012]. If EFDC is used to compare the SWE models (Figure 5.13), the depth profile is much better ap- proximated by the two-layer model. However, both SWE models are able to accurately simulate the position of the front if the correct front condition is used. The positive results between the two-layer SWE models and EFDC indicate the possibilities of the thin layer models to GC flows. In addition, the position of the nose and the slope of the reflected bore are very similar for both approaches. It is important to note that the SWE model simulates a smaller magnitude for the length of the nose region and is unable to model the mixing between the two fluids. 96 Figure 5.13: Dimensional slumping comparison between the SWE and EFDC for a) t=1.8sec and b) t=39sec. The discretization size for the SWE models is 800 cells. 97 Chapter 6 Conclusions and future developments In this thesis Boussinesq lock-exchange flows were analyzed experimentally and numerically with the one and two-layer SWE models. These types of flows are common in canals, rivers, estu- aries, etc. and play an important role in mixing and contaminant transport. Thus, computationally efficient and accurate numerical models are important to determine the trajectory and velocities of the GC. Although the model was constructed for Boussinesq, 1-D lock-exchange flows, the pro- posed BC at the GC front can be applied to a larger range of flows (e.g./ non-Boussinesq). In order to test the numerical model, MicroADV probes were used to measure the velocities of the bottom current and the ambient fluid, and a high definition video camera was used to track the GC front. This study presents a modeling framework based on the FVM to simulate GC flows using the one and two-layer SWE. The motivation was the overall success in the application of the FVM in the context of free surface flow simulation. While the simulation of GC flows with numerical models has been performed for almost three decades, some conditions associated with GC flows are more challenging to be represented numerically. This is particularly true for the case when GC propagate into small ambient depths (f 0:5), as in such cases a backward propagating wave with a large depth discontinuity poses difficulties to SWE-based models. To date, such flows have been simulated with models that perform explicit tracking of this hydraulic jump and the reflected bore, which results in additional model complexity. The proposed model aims to present a simple, albeit robust, way to perform GC flow com- putation to any ambient depth. Previous two-layer SWE models required different approaches depending on the initial fractional depth. In the proposed two-layer SWE model, two main contri- butions are presented. The first is a numerical model in which a single shock-capturing framework 98 is used to simulate GC flows for any initial depth ratio. The second contribution is a new al- ternative to compute the GC front boundary condition enforcing explicitly mass and momentum conservation, referred to as Dual-Cell (DC). This method was generally able to increase model accuracy, particularly for smaller discretization sizes while being substantially more computation- ally efficient than the alternative based on the method of characteristics (MOC). For some of the discretization sizes that were tested, the DC BC was more than 10 times faster than the MOC. The characteristic backward moving flow discontinuity cannot be properly simulated with a single-layer model as is discussed in [Ungarish, 2009]. Results obtained with the proposed two-layer SWE model are much better qualitatively and were more successful in reproducing ex- perimental results with the MicroADV probes. Moreover, the two-layer SWE compared very well with previous numerical models both for small and large ambient depths. The development of this numerical solution involved the application of the SWE equation in conservative variables and the separation of some terms of the expression derived by Rottman and Simpson [Rottman and Simpson, 1983] to be handled as source terms. Such source terms required special treatment to overcome difficulties as f 1, and the use of L?Hopital rule in the two-layer model was successful for this purpose. The velocity gradient component of the two-layer source term vector required the use of L?Hopital?s rule for h=H > 0:95. This velocity gradient term was not important in the source term calculations for large fractional depth, so the effect on the numerical scheme was negligible. However, if L?Hopital?s rule is not used the numerical code was unstable. On the other hand, the depth gradient term played the most important role in establishing the discontinuity instead of the depression wave for large values of f. When this discontinuity (instead of the depression wave) reflects off of the upstream boundary, the reflection causes the GC nose to form in case of the two-layer model. Once the upstream reflection occurs, the source terms are only important at the interface of the reflected bore (where the depth gradient is large), which eventually reaches the front of the GC. 99 The resulting one and two-layer SWE models were tested with different front conditions, dif- ferent strategies to compute the nose boundary condition and for various fractional depth. The most important factors as far as accuracy is concerned were: 1) proper selection of the mathe- matical model; 2) selection of the front condition expression; and 3) selection of nose boundary condition calculation. Unlike the initial expectation, three out of four numerical models performed comparably well in the simulations. The Lax-Wendroff scheme (even with the use of artificial viscosity) produced the worst results due to the oscillations observed and the code crashing at the upstream reflection. The use of a modified FORCE scheme and the non-linear Roe and HLL schemes was gen- erally very successful in reproducing the experimental results collected in this investigation, par- ticularly when associated with the two-layer SWE formulation and the modified version of the RS front condition (MRS) that was presented in this thesis. The empirical coefficient b in this expres- sion was calibrated for the full depth lock-exchange experiments, and the best agreements with experiments was b = 1:215. This value is clearly larger than the value provided in [Rottman and Simpson, 1983] (b = 1:0) for partial depth releases. It is reiterated that an important feature of the proposed model is the ability to simulate the upstream moving discontinuity observed during the slumping stage without the need of explicit tracking. Moreover, the objective was to develop a shock-capturing two-layer SWE model that can simulate GC flows without the restriction of the ambient depth. Although there has been extensive research for GC flows in recent decades, further investi- gations are required for better understanding. The effects of entrainment and complex geometries are still poorly understood in the context of the SWE models. For laboratory experiments, it is difficult to accurately measure the GC depth. For this study, there was an attempt to utilize pres- sure transducers to overcome this problem; however, surface waves smear the results indicating that larger-scale GC flows are required if transducers are to be used. In future studies additional CFD models will be tested (e.g. OpenFOAM or FLOW3D). In addition, the proposed two-layer modeling approach should be extended to non-Boussinesq problems and axisymmetric GC flows. 100 Bibliography [Abbott, 1961] Abbott, M. B. (1961). On the spreading of one fluid over another, part II: the wave front. La Houille Blanche, 6:827?836. [Adduce et al., 2011] Adduce, C., Sciortino, G., and Proietti, S. (2011). Gravity currents produced by lock-exchange: experiments and simulations with a two layer shallow-water model with entrainment. J. Hydraulic Eng. [Barr, 1967] Barr, D. I. H. (1967). Densimetric exchange flow in rectangular channels. iii large scale experiments. Houille Blanche, 6:619?631. [Baxter, 2000] Baxter, P. J. (2000). Impacts of eruptions on human health. In Encyclopedia of Volcanoes, pages 1035?1043. Academic Press. [Benjamin, 1968] Benjamin, T. B. (1968). Gravity currents and related phenomena. J. Fluid Mech., 31:209?248. [Birman et al., 2005] Birman, V. K., Martin, J. E., and Meiburg, E. (2005). The non-boussinesq lock-exchange problem. part 2. high-resolution simulations. J. Fluid Mech., 537:125?144. [Blong, 2000] Blong, R. (2000). Volcanic hazards and risk management. In Encyclopedia of Vol- canoes (e. H. Sigurdsson), pages 1215?1227. Academic Press. [Bonnecase et al., 1992] Bonnecase, R. T., Huppert, H. E., and Lister, J. R. (1992). Particle-driven gravity currents. Journal of Fluid Mechanics, 250:339?369. [Bonnecaze et al., 1993] Bonnecaze, R. T., Huppert, H. E., and Lister, J. R. (1993). Particle-driven gravity currents. J. Fluid Mech. [Britter and Simpson, 1978] Britter, R. E. and Simpson, J. E. (1978). Experiments on the dynam- ics of a gravity current head. J. Fluid Mech., 88:233?240. [Chaudhry, 2008] Chaudhry, M. H. (2008). Open-Channel Flow. Springer, second edition. [Chen et al., 2005] Chen, Q., Zhao, H., Hu, K., and Douglass, S. (2005). Prediction of wind waves in a shallow estuary. J. Waterw., Port, Coastal, Ocean Eng., pages 137?148. [Cunge et al., 1980] Cunge, J. A., Holly, F. M. J., and Verwey, A. (1980). Practical Aspects of Computational River Hydraulics. Pitman Publishing Ltd., London, UK. [Dalziel, 1993] Dalziel, S. B. (1993). Rayleigh-taylor instability: experiments with image analy- sis. Dyn. Atmos. Oceans, 20:127?153. 101 [Ellison and Turner, 1959] Ellison, T. H. and Turner, J. S. (1959). Turbulent entrainment in strati- fied fluids. J. Fluid Mech., 6:259?274. [Fannelop and Waldman, 1971] Fannelop, T. K. and Waldman, G. D. (1971). Dynamics of oil slicks. AIAA, 10:506?510. [Firoozabadi et al., 2010] Firoozabadi, B., Afshin, H., and Bagerpour, A. (2010). Experimental investigation of turbulence specifications of turbidity currents. J. Applied Fluid Mech., 3(1):63? 73. [Garcia and Kahawita, 1986] Garcia, R. and Kahawita, R. A. (1986). Numerical solutions of the de saint venant equations with the mac-cormack finite difference scheme. International Journal for numerical methods in fluids, 6:259?274. [Gerber et al., 2011] Gerber, G., Diedericks, G., and Basson, G. R. (2011). Paricle image ve- locimetry measurements and numerical modeling of a saline density current. J. Hydraulic Eng., pages 333?342. [Glaister, 1988] Glaister, P. (1988). Approximate riemann solutions of the shallow water equa- tions. J. Hydraulic Res., 26. [Godunov, 1959] Godunov, S. K. (1959). Finite difference method for numerical computation of discontinous solution of the equations of fluid dynamics. Matematicheskii Sbornik, 47(3):271? 306. [Hacker et al., 1996] Hacker, J., Linden, P. F., and Dalziel, S. B. (1996). Mixing in lock-release gravity currents. Dynamics of Atmospheres and Oceans, 24:183?195. [Hatcher et al., 2012] Hatcher, T. M., KC, M., Vasconcelos, J. G., and Fang, X. (2012). A compar- ison between numerical modeling approaches for the simulation of gravity currents. Submitted to the EWRI conference 1/2012. [H ?artel et al., 2000] H ?artel, C., Meiburg, E., and Necker, F. (2000). Analysis and direct numerical simulation of the flow at a gravity-current head. part1. flow topology and front speed for slip and no-slip boundaries. J. Fluid Mech., 418:189?212. [Hogg et al., 1999] Hogg, A. J., Ungarish, M., and Huppert, H. E. (1999). A box model for non- entraining, suspension-driven gravity surges on horizontal surfaces. Eur. J. Mech. B - Fluids, 19:139?165. [Hoult, 1972] Hoult, D. P. (1972). Oil spreading on the sea. Anna. Rev. Fluid Mech., 4:341?368. [Huppert, 2006] Huppert, H. E. (2006). Gravity currents: a personal perspective. J. Fluid Mech., 554:299?322. [Huppert and Simpson, 1980] Huppert, H. E. and Simpson, J. E. (1980). The slumping of gravity currents. J. of Fluid Mech., 99:785?799. [Keller and Chyou, 1991] Keller, J. J. and Chyou, Y. P. (1991). On the hydraulic lock-exchange problem. Z. Angew. Math. Phys., 42. 102 [Keulegan, 1957] Keulegan, G. H. (1957). An experimental study of the motion of saline water from locks into fresh water channels. Nat. Bur. Stand. Rept. Technical Report 5168. [Klemp et al., 1994] Klemp, J. B., Rotunno, R., and Skamarock, W. C. (1994). On the dynamics of gravity currents in a channel. J. Fluid Mech., 269:169?198. [Kranenburg, 1978] Kranenburg, C. (1978). Internal fronts in two layer flow. J. Hydraul. Div. ASCE, 104(HY10):1449?1453. [LeVeque, 1992] LeVeque, R. J. (1992). Numerical Mehods for Conservation Laws. Birkhauser Verlag. [Lowe et al., 2005] Lowe, R. J., Rottman, J. W., and Linden, P. F. (2005). The non-boussinesq lock-exchange problem. part 1. theory and experiments. J. Fluid Mech., 537:101?124. [Necker et al., 2002] Necker, F., Hartel, C., Kleiser, L., and Meiburg, E. (2002). High-resolution simulations of particle-driven gravity currents. Intl J. Multiphase Flow, 28:279?300. [Ozgokmen et al., 2009] Ozgokmen, T. M., Lliescu, T., and Fischer, P. F. (2009). Large eddy simulation of stratified mixing in a three-dimensional lock-exchange system. Ocean Modell., 26:134?155. [Roe and Pike, 1984] Roe, P. L. and Pike, J. (1984). Efficient construction and utilisation of ap- proximate riemann solutions. Computing Methods in Applied Scence and Engineering. [Rottman and Simpson, 1983] Rottman, J. W. and Simpson, J. E. (1983). Gravity currents pro- duced by instantaneous releases of a heavy fluid in a rectangular channel. J. Fluid Mech., 135:95?110. [Ruo and Chen, 2007] Ruo, A. and Chen, F. (2007). Modified shallow water equations for inviscid gravity currents. Physical Review E, 75. [Schoklitsch, 1917] Schoklitsch, A. (1917). Uber dambruchwellen. Sitzungsber. Akad Wis- senschaf, 126(1489). [Shin, 2001] Shin, J. (2001). Colliding gravity currents. PhD thesis, University of Cambridge. [Shin et al., 2004] Shin, J., Dalziel, S., and Linden, P. F. (2004). Gravity currents produced by lock exchange. J. of Fluid Mech., 521:1?34. [Simpson, 1997] Simpson, J. E. (1997). Gravity Currents In the Environment and the Laboratory. Cambridge University Press. [Simpson and Britter, 1979] Simpson, J. E. and Britter, R. E. (1979). The dynamics of the head of a gravity current advancing over a horizontal surface. J. Fluid Mech., 94:477?495. [Slim, 2006] Slim, A. C. (2006). High Reynolds number gravity currents. PhD thesis, University of Cambridge. 103 [Stommel and Farmer, 1952] Stommel, H. and Farmer, H. G. (1952). Abrupt change in width in two-layer open channel flow. J. Mar. Res, 11:205?214. [Sturm, 2010] Sturm, T. W. (2010). Open Channel Hydraulics. McGraw-Hill, second edition. [Timmermans et al., 2001] Timmermans, M. E., Lister, J. R., and Huppert, H. E. (2001). Com- pressible particle-driven gravity currents. J. Fluid Mech., 445:305?325. [Toro, 2001] Toro, E. F. (2001). Shock-Capturing Methods for Free-Surface Shallow Flows. John Wiley and Sons. [Turner, 1973] Turner, S. J. (1973). Buoyancy Effects in Fluids. Cambridge Unveristy Press. [Ungarish, 2007] Ungarish, M. (2007). A shallow-water model for high-reynolds-number gravity currents for a wide range of density differences and fractional depths. J. Fluid Mech., 579:373? 382. [Ungarish, 2009] Ungarish, M. (2009). An Introduction to Gravity Currents and Intrusions. Chap- man and Hall/CRC, first edition. [Ungarish and Huppert, 2004] Ungarish, M. and Huppert, H. E. (2004). On gravity currents prop- agating at the base of a stratified ambient: effects of geometrical constraints and rotation. J. Fluid Mech., 521(69-104). [Ungarish and Zemach, 2005] Ungarish, M. and Zemach, T. (2005). On the slumping of high reynolds number gravity currents in two-dimensional and axisymmetric configurations. Euro- pean J. Mech. B/Fluids, 24:71?90. [Vasconcelos and Nunes, 2009] Vasconcelos, J. G. and Nunes, A. M. S. (2009). Evaluating nu- merical schemes for the shallow water equations applied to dam break flows. In ABRH-Brazil, editor, Proc. XVIII Brazilian Symposium in Water Resources, Campo Grande, Brazil. [Vasconcelos et al., 2009] Vasconcelos, J. G., Wright, S. J., and Roe, P. (2009). Numerical oscil- lations in pipe-filling bore predictions by shock-capturing models. J. Hydraulic Eng., 135(296). [von Karman, 1940] von Karman, T. (1940). The engineer grapples with nonlinear problems. Bull. M. Math Soc., 46:615?683. [Wright et al., 1990] Wright, S. J., Kim, Y., and Buhler, J. (1990). Density current propagation in flowing receiving fluid. Technical report, Int. Conf. on Physical Modeling of Transport and Dispersion. Boston, MA. ASCE, New York, pp. 12A19-12A24. [Wright and Paez-Rivadeneira, 1996] Wright, S. J. and Paez-Rivadeneira, D. (1996). Density in- trusions with large relative thickness. Dynamics of Atmospheres and Oceans, 24:129?137. [Zoppou and Roberts, 1999] Zoppou, C. and Roberts, S. (1999). Catastophic collapse of water supply reservoirs in urban areas. J. Hydraulic Eng., 125(7):686?695. [Zoppou and Roberts, 2003] Zoppou, C. and Roberts, S. (2003). Explicit schemes for dam-break simulations. J. Hydraulic Eng., 129:11?34. 104