University of Southern California Dissertations and Theses

Open channel flow instabilities: modeling the spatial evolution of roll waves

(USC Thesis Other)

Open channel flow instabilities: modeling the spatial evolution of roll waves

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content OPEN CHANNEL FLOW INSTABILITIES: MODELING THE SPATIAL EVOLUTION OF ROLL WAVES by Ziyi Huang A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (CIVIL ENGINEERING) December 2013 Copyright 2013 Ziyi Huang ii Dedication To my beloved family iii Acknowledgements The present work can never be accomplished by the author alone. Herein I would like to express my appreciation to the people who have made essential contributions to my Ph.D. study: Dr. Jiin-Jen Lee, my advisor, for providing me with this opportunity to pursue the academic research. During my ve years of USC study, he has continuously inspired me with creative ideas and guided me in the proper direction. Not only being a professor in academia, he always shares meaningful life philosophies with his students. My Ph.D. committee members, Dr. Carter Wellford, Dr. Sami Masri, Dr. Iraj Nasseri, Dr. Julian Domaradzki, and Dr. James Moore, for their invaluable sugges- tions on the mathematical equations and numerical methods. Dr. Qin Chen, Dr. Patrick Lynett, Dr. Jinhai Zheng for the brief discussions in conferences. Dr. Lucio Soilbelman for ocially hooding me; Dr. Erik Johnson, for arranging me the prerequisite courses and teaching assistantship; Jennifer Gerson and Jennie Craig for taking care of my student le; Dragana Davidovic and Janet Ng for their assistances on the nancial support aairs. My past and present research group mates: Dr. Mehrdad Bozorgnia, Dr. Chanin Chaun-Im, Amir Eftekharian, Maryam Ghadiri, Dr. Hyoung-Jin Kim, Dr. Zhiqing Kou, Shentong Lu, Dr. Yuan-Hung Tan, Dr. Ben Willardson, Dr. Xiuying Xing iv (alphabetic order of the last name) for the warm atmosphere they have provided in the oce. I have beneted from our interactive communications on a variety of research subjects. My friends who also study at USC: Heli Bao, Zhaohu Fan, Hao Feng, Yujing Hu, Nan Li, Weixuan Li, Haowang Wang, Zheng Yang, Dr. Xiaoming Zheng (alphabetic order of name) for the chats on recent progress and future plans. I wish you have a successful and wonderful life. Chinese-American Engineers and Scientists Association of Southern California (CESASC) for oering me 2011 Scholarship. USC CSSA (Chinese Students and Scholars Association) for organizing colorful activities that have relieved me of \permanent head damage". Especially, my parents, Mingfang Ding and Jianhui Huang, for their consistent encouragement which keeps me away from surrender. I sincerely wish you are accom- panied with health and happiness. Finally, Qian Ouyang, my anc ee, for your understanding on my endless research tasks, and love re ected by the beautiful smile. v Table of Contents Dedication ii Acknowledgements iii List of Figures vii List of Tables x Nomenclature xi Abstract xviii Chapter 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Objective and Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Chapter 2 Literature Review 10 2.1 Analytical Investigation . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.1 Theory of Stability . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.2 Mathematical Solution of Progressing Roll Waves . . . . . . . 16 2.1.3 Amplitude Evolution Equation . . . . . . . . . . . . . . . . . 17 2.2 Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3 Computational Dynamics . . . . . . . . . . . . . . . . . . . . . . . . 25 Chapter 3 Methodology of Model Construction 28 3.1 The Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Numerical Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2.1 Finite Volume Formulation . . . . . . . . . . . . . . . . . . . . 34 3.2.2 The Riemann Solver . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 48 3.2.4 Solution Procedure . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3 Model Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Chapter 4 Convergence and Validation Tests 59 4.1 Comparison with Transient Linear Theory . . . . . . . . . . . . . . . 59 4.1.1 Linearized Theory of Stability . . . . . . . . . . . . . . . . . . 62 4.1.2 The Transient Linear Theory . . . . . . . . . . . . . . . . . . 65 4.1.3 Numerical Experiment of Progressing Roll Waves . . . . . . . 67 4.1.4 Comparison Results . . . . . . . . . . . . . . . . . . . . . . . . 71 4.2 Comparison with Quasi-steady Nonlinear Theory . . . . . . . . . . . 76 4.2.1 Comparison Results . . . . . . . . . . . . . . . . . . . . . . . . 77 vi 4.2.2 Sensitivity Analyses . . . . . . . . . . . . . . . . . . . . . . . . 77 4.3 Comparison with Brock's Experiment . . . . . . . . . . . . . . . . . . 81 4.3.1 Time Series of Roll Waves . . . . . . . . . . . . . . . . . . . . 84 4.3.2 Time-averaged Wave Crest and Trough Depths . . . . . . . . 91 4.3.3 Sensitivity Analyses . . . . . . . . . . . . . . . . . . . . . . . . 92 4.3.4 Correction Length . . . . . . . . . . . . . . . . . . . . . . . . 102 4.3.5 Comparison Results . . . . . . . . . . . . . . . . . . . . . . . . 104 Chapter 5 Investigation of Roll-wave Properties 118 5.1 Demonstration of the Spatial Evolution of Roll Waves . . . . . . . . . 118 5.2 Wave-wave Interactions . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.3 Generality of the Spatial Evolution of Roll Waves . . . . . . . . . . . 131 5.4 Spectral Analysis of Roll Waves . . . . . . . . . . . . . . . . . . . . . 133 5.5 Eect of Bottom Roughness . . . . . . . . . . . . . . . . . . . . . . . 139 Chapter 6 Roll Waves in a Breaking-slope Open Channel 144 Chapter 7 Conclusions 163 7.1 On the Computational Model USC SWAT . . . . . . . . . . . . . . . 163 7.2 On Simulating the Spatial Evolution of Roll Waves . . . . . . . . . . 164 7.3 On Investigating Roll-wave Properties . . . . . . . . . . . . . . . . . . 165 7.4 On Mitigating the Spatial Evolution of Roll Waves . . . . . . . . . . 166 7.5 Future Research Considerations . . . . . . . . . . . . . . . . . . . . . 167 References 169 Appendix A Dressler's Mathematical Solution of Roll Waves 174 Appendix B Derivation of the Saint Venant Equations 181 B.1 Conservation of Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 B.2 Conservation of Momentum . . . . . . . . . . . . . . . . . . . . . . . 185 B.3 Disscusion on Approximation B.1 and B.2 . . . . . . . . . . . . . . . 193 vii List of Figures 1.1 Schematic of roll-wave evolution . . . . . . . . . . . . . . . . . . . . . 2 1.2 Photos of roll waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Interaction between water ow and channel bed . . . . . . . . . . . . 5 1.4 Sketch of roll-wave sections . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 Channel in Brock's experiment . . . . . . . . . . . . . . . . . . . . . 20 2.2 Inlet box in Brock's experiment . . . . . . . . . . . . . . . . . . . . . 21 3.1 Sketch illustration of open channel ow . . . . . . . . . . . . . . . . . 29 3.2 Schematic of computational domain discretization . . . . . . . . . . . 36 3.3 Schematic of formulated quantities . . . . . . . . . . . . . . . . . . . 38 3.4 Sweby's criterion for a second-order TVD scheme . . . . . . . . . . . 47 3.5 Schematic of applying rst-order extrapolation at in ow boundary . . 51 3.6 Structure of the model USC SWAT . . . . . . . . . . . . . . . . . . . 57 4.1 Boundary conditions for progressing roll waves . . . . . . . . . . . . . 69 4.2 Solution proles of progressing roll waves . . . . . . . . . . . . . . . . 73 4.3 Comparison with transient linear theory . . . . . . . . . . . . . . . . 74 4.4 Comparison with quasi-steady nonlinear theory . . . . . . . . . . . . 78 4.5 Sensitivity analyses in the numerical experiment of progressing roll waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.6 Simulated time series of roll waves: water depth . . . . . . . . . . . . 86 4.7 Contoured time series of roll waves: water depth . . . . . . . . . . . . 89 4.8 Time series of roll waves: ow velocity . . . . . . . . . . . . . . . . . 90 4.9 Sensitivity to the cell size in the simulation of real roll-wave ows . . 94 4.10 Sensitivity to the turbulent viscosity in the simulation of real roll-wave ows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 viii 4.11 Demonstration of corrected in ow boundary . . . . . . . . . . . . . . 103 4.12 Demonstration of time-averaged wave crest and trough depths before and after correcting upstream . . . . . . . . . . . . . . . . . . . . . . 104 4.13 Overview of comparison with Brock's experiment: Scenario (1) . . . . 107 4.14 Overview of comparison with Brock's experiment: Scenario (2) . . . . 108 4.15 Overview of comparison with Brock's experiment: Scenario (3) . . . . 109 4.16 Overview of comparison with Brock's experiment: Scenario (4) . . . . 110 4.17 Close view of comparison with Brock's experiment: Scenario (1) . . . 113 4.18 Close view of comparison with Brock's experiment: Scenario (2) . . . 114 4.19 Close view of comparison with Brock's experiment: Scenario (3) . . . 115 4.20 Close view of comparison with Brock's experiment: Scenario (4) . . . 116 5.1 Demonstration of the spatial evolution of roll waves . . . . . . . . . . 119 5.2 Demonstration of the spatial evolution of roll waves: separate channel sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.3 Demonstration of the wave overtaking process from the simulation of a real roll-wave ow . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.4 Demonstration of the wave absorption process from the simulation of a real roll-wave ow . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 5.5 Demonstration of the wave spawning process from the simulation of a real roll-wave ow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.6 Generality of the spatial evolution of roll waves . . . . . . . . . . . . 132 5.7 Simulated time series of roll waves for the spectral analysis . . . . . . 136 5.8 Amplitude of roll waves in the frequency domain from the present simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 5.9 Development of signicant wave height versus channel location for dif- ferent bottom roughness . . . . . . . . . . . . . . . . . . . . . . . . . 141 5.10 Development of normalized signicant wave height versus channel lo- cation for dierent bottom roughness . . . . . . . . . . . . . . . . . . 142 6.1 Prole geometry of the original open channel in the computational experiment of the breaking-slope open channel . . . . . . . . . . . . . 145 ix 6.2 Geometry of the modied open channel in the computational experi- ment of the breaking-slope open channel . . . . . . . . . . . . . . . . 146 6.3 Prole of normal ow conditions on the transition between Reach (1) and Reach (2) in the computational experiment of the breaking-slope open channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 6.4 Prole of normal ow conditions on the transition between Reach (2) and Reach (3) in the computational experiment of the breaking-slope open channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 6.5 Demonstration of time series of water depth on Reach (2) in the com- putational experiment of the breaking-slope open channel . . . . . . . 153 6.6 Comparison of time-averaged wave crest depth between constant-slope and breaking-slope channels . . . . . . . . . . . . . . . . . . . . . . . 156 6.7 Comparison of time-averaged wave trough depth between constant- slope and breaking-slope channels . . . . . . . . . . . . . . . . . . . . 157 6.8 Comparison of signicant wave height between constant-slope and breaking-slope channels . . . . . . . . . . . . . . . . . . . . . . . . . . 160 A.1 Dressler's mathematical solution of roll waves . . . . . . . . . . . . . 178 B.1 Schematic for conservation of mass . . . . . . . . . . . . . . . . . . . 183 B.2 Schematic for conservation of momentum . . . . . . . . . . . . . . . . 187 x List of Tables 4.1 Comparison of computational time among dierent cell sizes . . . . . 96 4.2 Percentages of increase of wave crest depths due to decreasing the value of turbulent viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.3 Percentages of decrease of wave trough depths due to decreasing the value of turbulent viscosity . . . . . . . . . . . . . . . . . . . . . . . . 101 4.4 Hydraulic parameters applied in Brock's experiment and the present simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.5 Computational parameters applied in the present simulation for com- parison with Brock's experiment . . . . . . . . . . . . . . . . . . . . . 106 5.1 Hydraulic parameters for investigating the eect of bottom roughness 140 6.1 Hydraulic parameters for three reaches in the modied, breaking-slope open channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 xi Nomenclature Following symbols are used in this study: A surface area A 0 initial maximum wave amplitude a acceleration a j lower coecients in the tridiagonal matrix B system some property of a system B transitional matrix for the vector of ow variables b some property in a control volume b j middle coecients in the tridiagonal matrix C arbitrary constant C c Ch ezy coecient C d drag coecient C transitional matrix for the vector of convective uxes Cr Courant number c wave celerity c j upper coecients in the tridiagonal matrix c 0 j updated upper coecients in the tridiagonal matrix c w wave speed of progressing roll waves c 0 kinematic wave speed c slowest kinematic wave speed c + fastest kinematic wave speed D dissipation vector in the TVD Lax-Wendro scheme xii d components in the dissipation vector in the TVD Lax-Wendro scheme e mathematical constant 2:71828 e j coecients in the vector associated with the tridiagonal system e 0 j updated coecients in the vector associated with the tridiagonal system e 1 rst right eigenvector of the Jacobian matrix of the vector of convective uxes e 2 second right eigenvector of the Jacobian matrix of the vector of convective uxes F force F p force due to water pressure F b drag force due to bottom shear stress F xx viscous force due to normal stress F 1 c rst component of the vector of convective uxes F 2 c second component of the vector of convective uxes F Froude number F 0 undisturbed Froude number F vector of uxes F c vector of convective uxes F d vector of diusive ux f Darcy-Weisbach friction factor f n current value of the friction factor in the the Newton-Raphson's method f n+1 future value of the friction factor in the the Newton-Raphson's method f 0 initial value of the friction factor in the the Newton-Raphson's method g gravitational acceleration H water depth in the progressing coordinate xiii H A transitional water depth in the progressing coordinate H B transitional water depth in the progressing coordinate H b water depth at wave back in the progressing coordinate H c critical water depth in the progressing coordinate H f water depth at wave front in the progressing coordinate H n some progressing roll wave H n+1 progressing roll wave next to H n H 1/3 signicant wave height h water depth h c time-averaged wave crest depth h t time-averaged wave trough depth h 0 normal water depth I identity matrix i square root of -1 J Jacobian matrix of the vector of convective uxes j index of control cells in a computational domain K progressing discharge per unit width K c critical progressing discharge per unit width K 1 coecient in the Colebrook-White equation K 2 coecient in the Colebrook-White equation K 3 coecient in the Colebrook-White equation k wavenumber k s bottom roughness height L wavelength ` index of time level M momentum M system momentum of a system xiv m mass m system mass of a system N total number of control cells in a computational domain NR function associated with the Newton-Raphson's method P eigenvector matrix in the TVD Lax-Wendro scheme p water pressure Q vector of ow variables Q mirror vector of ow variables on the mirror cell q ow rate per unit width q 0 normal ow rate per unit width R h hydraulic radius < Reynolds number R + all positive real numbers r ratio for the ux limiter r p modulus of the complex expression for the initial ow velocity S external sources S vector of external sources t time U depth-averaged ow velocity in the progressing coordinate U c critical depth-averaged ow velocity in the progressing coordinate u ow velocity depth-averaged throughout cross-section u 0 normal ow velocity depth-averaged throughout cross-section V volume W weight w vertical ow velocity w 1 rst component in the parameter vector w 2 second component in the parameter vector xv w parameter vector x spatial coordinate along channel bed x interface channel location on a cell interface x max maximum channel location Z + all positive integers z spatial coordinate vertical to channel bed vector of wave strength 1 rst component of the vector of wave strength 2 second component of the vector of wave strength coecient in the quadratic equation for the complex angular frequency coecient in the quadratic equation for the complex angular frequency Q value dierence of the vector of ow variables t time step u value dierence of the ow velocity w value dierence of the parameter vector x size of control cells in a computational domain discriminant of the quadratic equation for the complex angular frequency small tolerance " small nonnegative quantity progressing coordinate s shock location in the progressing coordinate small perturbation of the water depth angle of channel bed inclination p argument of the complex expression for the initial ow velocity eigenvalues of the Jacobian matrix of the vector of convective uxes xvi 1 rst eigenvalue of the Jacobian matrix of the vector of convective uxes 2 second eigenvalue of the Jacobian matrix of the vector of convec- tive uxes kinematic viscosity t turbulent viscosity small perturbation of the ow velocity depth-averaged throughout cross-section O order of truncation error water density h standard deviation of the time series of the water depth b shear stress between water boundary layer and channel bottom xx normal stress entropy violation prevention function ux limiter ! complex angular frequency ! I imaginary part of the complex angular frequency ! R real part of the complex angular frequency xvii General top signs not specied above: e approximate value arithmetic mean value General superscripts not specied above: v index for the component in the dissipation vector in the TVD Lax-Wendro scheme General subscripts not specied above: i.b. quantity on the in ow boundary left quantity on the left hand side of a cell interface right quantity on the right hand side of a cell interface quantity on cell interfaces xviii Abstract This study is concerned with \roll waves", which are a series of intermittent shock- like waves occurring in turbulent water ows down a wide, rectangular open channel due to the natural instability. Most literature on the subject of roll waves have concentrated on theoretical investigations. Only two of them made eorts on compu- tationally simulating real roll-wave ows, but no validation tests were implemented. A computational model capable of simulating the spatial evolution of roll waves is constructed in the present work. The equations of motion are the one-dimensional viscous Saint Venant equations, incorporating the turbulently viscous force to describe the momentum diusion resulting from tremendous gradient and curvature of roll waves. The equations are solved by a high-resolution, shock-capturing numerical scheme. The scheme is based on the conservative nite volume formulation and adopts the second-order TVD Lax-Wendro Riemann solver. The computational model, written in FORTRAN 95, is named as \Shallow Water Analyzed by a TVD Scheme at University of Southern California" with an acronym \USC SWAT". It can be applied to analyze a wide range of shallow water ows. The numerical scheme is tested through comparisons with both transient linear and quasi-steady nonlinear theories. The numerical solutions of maximum amplitude of progressing roll waves match the linear theory at time levels in which the wave amplitude is small. The numerical solution of converged progressing roll waves match a mathematical solution of roll waves deduced from the governing equations in a xix progressing coordinate. The constructed computational model is utilized to simulate real roll-wave ows in a constant-slope channel. In the validation test, the spatial development curves of time-averaged wave crest and trough depths achieved from the simulation are com- pared with those from experimental data. It is shown that the value of the turbulent viscosity needs to be appropriately selected so as to agree with the experimental data. The comparison results manifest the simulation satisfactorily predicts the spatial evo- lution of wave crest and trough depths. Some of important roll-wave properties are investigated. The simulation discovers three types of wave-wave interactions, which impact characteristics of the roll-wave evolution. The spatial evolution of roll waves obeys a generality despite dierent hydraulic conditions. The nonperiodicity of roll waves increases from upstream to downstream channel locations. The extent of the ow instability is weakened in channels with rough bottom. A concept of the breaking-slope open channel is tested through a computational experiment. This concept is found to be useful in suppressing growth of the roll-wave amplitude. This computational study covering the comparison with experimental data en- riches the knowledge of computational mechanics regarding free-surface water ows, and provides a complete insight on simulating the spatial evolution of roll waves using the Saint Venant equations. Findings from the computational experiments on miti- gating the roll-wave evolution could help overcome diculties caused by open channel ow instabilities. 1 Chapter 1 Introduction The uniform, free surface, turbulent water ows down an inclined, wide, rectan- gular open channel may become unstable spontaneously, if some sucient conditions are satised. As a result, steady hydraulic conditions are interrupted, including the water depth (or pressure) and ow velocity. When the channel is suciently long, the small instabilities on the water surface will develop to a series of intermittent shock-like waves separated by sections of the gradually varying ow. This natural phenomenon is known as \roll wave". The schematic of the spatial evolution of roll waves is presented in Figure 1.1. 1.1 Background Roll waves have been widely discovered and recorded in literatures. Mostly they are mentioned to occur in turbulent ows down free surface inclined channels such as spillways from dams, open aqueducts, run-o channels, and creeks. Cornish (1934) presented detailed observations of roll waves in 1910 and published photographs of them in a long rectangular open conduit. Rouse (1938) described the fashion of pro- ceeding slugs and mentioned that the process is readily observed in nature, although he did not perform analyses of the roll-wave motion. Brock (1969) published a photo of progressing roll-wave trains down the Santa Anita Wash in Acadia, California. More recently it is found that water roll waves can be occurrences in other natural 2 Water flow 0 h u c Horizontal θ Figure 1.1. Schematic prole of the spatial evolution of roll waves in an inclined open channel. h 0 is the normal water depth,u is the depth-averaged ow velocity through- out cross-section, c is the wave celerity, and is the angle of channel inclination. 3 environments, for example supraglacial meltwater channel ows (Carver et al. 1994), dense current in an incline (Alavian 1986; Cenedese & Whitehead 2004), propagating periodic bores in the ocean (Swaters 2003), and o-shore cascading slugs in the lake (Fer et al. 2002). Additionally, various other uid ow instabilities analogous to roll waves have been detected and studied, such as the bentonite suspension (Engelund & Wan 1984), granular ow (Forterre & Pouliquen 2003), gravity-driven ow in col- lapsible tubes (Brook et al. 1999), viscoplastic mud ow (Ng & Mei 1994; Liu & Mei 1994; Coussot 1997; Balmforth & Liu 2004), multiphase and power-law ow (Woods et al. 2000; Pascal 2006; Pascal & D'Alessio 2007). Similar disturbances are also visible in laminar sheet ows (Juliean & Hartley 1986), which is considered as the counterpart of roll waves in turbulent ows. Lee took some photographs of roll-wave trains down a rectangular free surface creek in 1993. The creek located in city of Sierra Madre, California conveys water from the San Gabriel mountain. The photos are presented in Figure 1.2. The left three photos show multiple intermittent roll waves passing a portion of the creek while the right two show closer views of a single wave. The photos show the water appears to be white at the vicinity of the wave front, demonstrating the existence of turbulent eddies in the travelling shock-like waves. The formation of the instabilities depends on the interaction between the ow sheet and the rough boundary wall. It is analytically proved that roll waves do not establish if the channel resistance is zero and that the size of the waves will approach zero when the resistance is suciently reduced (Dressler 1949). Rouse (1938) men- tioned that roll waves do not occur if the channel bed is suciently irregular, which suggests too much resistance also prevents the formation of perturbations. Although 4 Figure 1.2. Photographs of roll waves in a long rectangular creek. The creek is located in the city of Sierra Madre, California conveying water from the San Gabriel mountain. Left three photos show multiple roll waves passing a portion of the creek and right two photos show closer views of a single roll wave. Photographs taken by Jiin-Jen Lee. roll waves have some analogies to bores in the ocean environment, they are domi- nated by the resistance whilst bores are not. In addition, bores will break whereas roll waves propagate without distortion. Figure 1.3 shows the interaction between the water ow and the rough channel bed. The roughness is illustrated by the uneven bottom. When the resistance leads to the instability, the ow kinematic-wave velocity becomes greater than the propagation velocity of shallow water waves. Consequently, the small perturbations expand themselves as owing downward, and eventually grow to large, shock-like waves. The Froude number,F, is known as the most important dimensionless parameter 5 h Water flow θ u Horizontal Figure 1.3. Micro-view of interaction between the water ow and rough channel bed. h is the water depth above the bottom, u is the depth-averaged ow velocity throughout cross-section, and is the inclination angle of the channel bed. in free surface, gravity driven ows. It is dened as F = u p g cosh ; (1.1) whereu is the mean ow velocity,g is the gravitational acceleration, is the angle of channel inclination, andh is the water depth. The denominator of the Froude number is considered as the phase speed of shallow water waves. According to the dimensional analysis, the Froude number is the ratio of the inertial force to the gravitational force. The Froude number classies open channel ows into three categories: 6 Horizontal Water flow Figure 1.4. Sketch illustration of roll-wave transition between subcritical and super- critical ow in a continuous section. Two continuous sections are separated by the discontinuous hydraulic jump. subcritical ow, ifF< 1; supercritical ow, ifF> 1; critical ow, ifF = 1. Thomas (1937) observed that when roll waves propagate downward, the water ow continuously transits from subcritical to supercritical ow, and transits back to sub- critical ow again through a moving hydraulic jump (or a shock), as shown in Fig- ure 1.4. However, the absolute mean ow velocity is greater at the subcritical part than the supercritical region, and is everywhere less than the wave speed. Hence, what Thomas observed is the ow velocity with respect to progressing roll waves. 7 1.2 Objective and Scope When a steep conduit is built for conveying water, the roll-wave problem always comes to the attention of open channel hydraulic engineers. Such disturbances are an undesirable occurrence. The original channel eligibility of accommodating the designed water discharge may be violated if the uniform ow becomes the roll-wave ow. The water will over ow the channel sides with disturbances, and thus results in ood threats to adjacent environments. Besides, these instabilities lead to sud- den variations of the pressure and momentum, which can impact the robustness of hydraulic structures. Similar problems also occur in run-o channels of mud ows and residential gutter ows. In recent decades, clinical engineers have found unstable blood ows in vessels of humans and animals. These phenomena greatly enhance the probability of vessel burst due to abrupt changes of the blood pressure. The objective of the present study is to computationally explore major properties of roll-wave trains down a wide, rectangular open channel. When the water is expected to ow in a natural or articial channel, the initial concern is whether the small perturbations, which result from turbulences and the resistance, will grow to large roll waves. If the channel is such that instabilities self-grow, the maximum water depth behind some certain section will denitely be greater than the original normal depth. Therefore, it is indispensable to estimate the maximum water depth along the whole channel so as to protect environments near the water ow from ooding. For a real water ow that is highly nonlinear and subject to eddy viscous eects, there does not exist a fully exact solution. A physical hydraulic model has to be installed for each open channel case if the ow is studied experimentally, and thus 8 impact the eciency. The computational model should be the most proper approach to simulate a complete water ow characterized by instabilities and roll waves. Once the model is established and validated, it can be easily applied for the interpretation of roll-wave properties. A thorough ow process, specifying the spatial evolution of roll waves from the normal ow condition, destabilization to small perturbations, to developed large, shock-type waves, can be achieved in a computer based upon equations of motion and an appropriate numerical scheme. For simulating channels having dierent dimensions and designed discharges, the computing task can be easily accomplished by adjusting some necessary parameters. Besides the maximum water depth, further interests may also include the mini- mum water depth, ow velocity, wave period, and wave-wave interactions. Moreover, available approaches for attenuating the roll-wave evolution are proposed in this study. The rst approach is to increase the bottom roughness so that the ow resistance is suciently large to prevent instabilities. The second approach is to modify the chan- nel slope. The original channel with a constant slope is modied to a breaking-slope channel in which no large roll waves are generated. Generally, the principal goal of attenuating the roll-wave evolution is to reduce the wave amplitude and restore the stability of open channel ows. The center of the present study is therefore located at the construction of a com- putational model capable of simulating the spatial evolution of roll waves in a wide, rectangular open channel. The equations of motion describe the conservation of phys- ical properties of water ows and conform to the shallow water theory. The equations form a system of nonlinear, nonhomogeneous, parabolic-hyperbolic type partial dif- 9 ferential equations. We need to develop a conservative high-resolution numerical scheme capable of capturing discontinuous shocks and retaining physical properties in the computational domain. In Chapter 2, a comprehensive literature survey will be conducted including ana- lytical, experimental, and computational roll-wave studies. Chapter 3 will discuss the method of the model construction, including the equations of motion and the numer- ical scheme. In Chapter 4, the accuracy of the numerical scheme is veried through convergence tests. Then the simulation results on the spatial evolution of roll waves are compared with experimental data. Chapter 5 will investigate some of important roll-wave properties using the constructed computational model. A concept of the breaking-slope open channel is introduced and tested by a computational experiment in Chapter 6. Finally, major conclusions from the present study and proposed future works will be stated in Chapter 7. 10 Chapter 2 Literature Review In this chapter, we review important studies pertinent to roll waves. It should be emphasized that the present study concentrates on roll waves occurring in turbulent water ows down a inclined open channel. Those literatures studying roll waves or similar physical processes of other uids associated with other mediums, as mentioned in Chapter 1, are not considered in this chapter. 2.1 Analytical Investigation Since Cornish (1934) described the detail roll-wave phenomenon, most researches in this eld have been analytical. Those theories are not able to predict real wa- ter ows in an open channel, because the equations of motion have nonlinearity and are not analytically solvable. However, theories approximated from the equations of motion adequately identify important necessary conditions for the formation of roll waves. The theoretical solution of periodic steady-state progressing roll waves have been used for several numerical scheme tests, even though they are not related to natural ows. Amplitude evolution equations deduced from the asymptotic approxi- mation had been used to study the open channel ow instabilities. 11 2.1.1 Theory of Stability Jereys (1925) rst obtained a stability criterion by imposing small sinusoidal perturbations on a smooth water surface. He started with the uni-dimensional shallow water theory, augmented by the Ch ezy formula that addressed the gravitational force and bottom resistance: @h @t + @(hu) @t = 0; (2.1) @u @t +u @u @x +g @h @x =g tan f 8 u 2 h ; (2.2) where h was the water depth, u was the depth-averaged ow velocity, g was the gravitational acceleration, was the angle of channel inclination, andf was the Darcy- Weisbach friction factor, x and t were arguments denoting the spatial coordinate and time, respectively. The dimensionless form of the above equations were derived and linearized by neglecting any products of dependent variables. Hence, a single linearized equation was achieved: (1 1 F 2 0 ) @ 2 @(x 0 ) 2 + 2 @ 2 @x 0 @t 0 + @ 2 @(t 0 ) 2 + ( tan F 2 0 )( h 0 )(3 @ @x 0 + 2 @ @t 0 ) = 0; (2.3) by substituting f = 8 tan F 2 0 : (2.4) was the dimensionless perturbation from the unity normal depth, was the di- mensional wavelength,h 0 was the dimensional normal depth,F 0 was the undisturbed Froude number, and the prime sign represented dimensionless groups. A sinusoidal 12 perturbation could be expressed by either the real or imaginary part of = 0 exp[i2(x 0 Ct 0 )]; (2.5) where C was the complex wave speed, and could be written as C =C R +iC I ; (2.6) where C R denoted the phase speed, and C I was the imaginary speed relating the change of the wave amplitude. Hence, equation 2.5 became = 0 exp(2C I t 0 ) exp[i2(x 0 C R t 0 )]: (2.7) The stability occurred when the amplitude of the prescribed small perturbation did not magnify as time increases. Obviously, this happened whenC I was equal to or less than zero. In order to nd out the stability criterion, equation 2.5 was substituted into 2.3 to give C R = 3Y + 4C I 2Y + 4C I ; (2.8) and F 0 = 1 q 1 (2C R +C 2 I +C 2 R +C I Y / ) ; (2.9) where Y = tan F 2 0 h 0 : (2.10) 13 Therefore, Jereys derived the linearized stability criterion: C R = 3 2 ; (2.11) F 0 = 2: (2.12) Jereys's study discovered an important linearized theory of stability. It indicates that water ows having a undisturbed Froude number of less than or equal to 2 are stable, otherwise disturbances will develop to large-amplitude roll waves. In addition to linearization approximation, the wave amplitude 0 exp(2C I t 0 ) must be small so that the exponential sinusoidal perturbation is consistent with the original nonlinear equations. This theory has proved to be accurate in most cases through many other similar approaches, both analytically and numerically. Dressler (1949) also discovered two necessary conditions for roll waves to form. He applied the same equations as Jeerys's study but constructed a progressing coor- dinate system that incorporated arguments of space and time. Using the progressing coordinate he was able to obtain an open channel ow prole equation in a dierential form. At the critical condition in the progressing coordinate, he showed that no roll waves formed if the bottom resistance was equal to zero. Furthermore, after substi- tution of the critical values into the prole equation, the other criterion of existence of roll waves was devised as 4r 2 <m (2.13) wherer 2 denoted the bottom resistance coecient andm was the channel slope. This expression indicated if the resistance was too large compared to the channel slope, no 14 roll waves would form. The two necessary conditions for the growth of the perturbation obtained by Dressler are also valid in most cases, as have been demonstrated by many other researches. The second instability condition was also derived by Thomas (1937) in a somewhat dierent fashion. Needham and Merkin (1984) used an approach similar to Jereys's work to derive a linearized theory of stability, except that they established an inequality by requiring the real part of the wave frequency to be non-negative for the stability, while Jereys required the imaginary part of the wave speed to be non-positive. However, they still derived an identical sucient and necessary condition to Jereys for the ow to be stable: F6 2. Yu and Kevorkian (1992) also used the small perturbation method to derive a linearized system. The system was solved analytically to show that the initial per- turbations decay exponentially if 0 6F < 2 whereas if F > 2 the solution presents an exponentially growing disturbance. They emphasized it is crucial to consider the cumulative eect of nonlinearities to obtain a correct description of roll waves over long time. The Jereys's method was duplicated by Balmforth and Mandre (2004) except that the innitesimal perturbations were represented by a truncated Fourier series with a Floquet multiplier. With this approach, they were able to nd out the critical Froude number of both turbulent and laminar ows. Moreover, they explored that the critical value of Froude number for the stability was in uenced by the wavenumber of the uneven bottom topography and strength of the viscous eect. 15 In the work of Noble (2007), the asymptotic linear stability of roll waves specied the viscous Saint Venant equations was considered. The viscous eect was described by the second derivative of ow variables multiplied by a real viscosity. The aim of the study was to prove a strong spectral stability implying the linear stability of roll waves. This was accomplished by computing the pointwise estimate on Green functions. This study also deduced the necessary conditions for the spectral stability of roll waves in the low frequency regime via computing an expansion of the Evans function. Extending Noble's work, Johnson et al. (2011) showed the spectral stability that implied nonlinear stability of spatially periodic viscous roll-wave solutions to the one- dimensional Saint Venant equations. They treated the equations in Lagrangian coor- dinates from which large-amplitude nonlinear damping estimates could be obtained. One of diculties to reach the purpose of this study was the nonconservative form of the equations which could not be dealt using the estimates of the conservative case. This was resolved by observing the decaying amount of the undierentiated quadratic source terms. The linearized theory of stability for roll waves in open channel ows derived by Jeereys and Dressler can be obtained by another approach, which is not covered by any of above literatures. Similarly, this approach starts from the one-dimensional hyperbolic Saint Venant equations and applies the small perturbation method. After having the linearized equation, the low order approximation needs to be performed. The adjustment of the roll-wave stability is dependent on the sign of the momentum diusivity. The specic derivation process of this approach will be given in Chapter 4 16 of this report. 2.1.2 Mathematical Solution of Progressing Roll Waves In Dressler's work (1949), besides stability theories, a progressing roll-wave pro- le was also deduced. This result was the special mathematical solution to the one- dimensional Saint Venant equations in Dressler's progressing coordinate. Those iden- tical sections of a continuous roll-wave prole were joined together to construct dis- continuous periodic solutions. The discontinuities connecting roll-wave prole were determined by hydraulic shock conditions. The method contributed by Dressler enables one to establish an exact solution of converged progressing roll waves, if the channel slope, bottom drag, wavelength, and the wave speed are prescribed. Although this method cannot predict natural water ows, it is commonly used to test the convergence of a numerical scheme solving for the shallow water theory. Dressler's solution is also used for testing the convergence of the numerical scheme of the present study, and the comparison results will be given in Chapter 4. Both the detail derivation of the continuous roll-wave prole and how the discontinuous periodic solutions are constructed are reproduced in Appendix A. Needham and Merkin (1984) applied the equations similar to Dressler, but further augmented by an turbulently viscous dissipation term 0 @ 2 v/@x 2 shearing normal to the ow. The constant 0 is indeed known as the turbulent viscosity, but it was not dened in this work. Hence they were able to construct continuous, periodic solutions by including the turbulent viscosity term to dissipate energy near the shock. Hwang and Chang (1987) only analyzed their model equations and discovered a new family 17 of roll-wave solutions using the dynamic singularity theory. Analytical solutions of roll waves were found in a inclined channel of an arbitrary cross-section (Boudlal & Liapiadevskii 2004). The starting equations and the progressing coordinate is the same as Dressler's study, except the assumption of the nonconvex pressure function. There have been other analytical roll-wave solutions governed by amplitude evo- lution equations. These equations derived from the well-known shallow water con- servation laws not only specify the quasi-steady solutions as t!1, but they are also capable of predicting transient solutions as disturbances develop with respect to time. Hence, the literatures related to amplitude evolution equations are discussed separately in next subsection. 2.1.3 Amplitude Evolution Equation The amplitude evolution equation is a model equation. It provides transient solutions from small disturbances on a uniform open-channel ow to large roll waves. It is derived in such a way that dependent variables of original equations of motion are approximated using asymptotic expansions. Due to the complicacy, the detail derivation procedures for amplitude evolution equations will not be given herein. Yu and Kevorkian (1992) obtained a set of amplitude evolution equations by concentrating on weakly unstable problems, i.e. 0<F 2 1. As typical as others, they started from the dimensionless shallow water theory and assumed a multiple scale expansion form of dependent variables. The expansion was substituted into the shallow water theory to obtain the evolution expression for ow variables. The results of the evolution expression compared quite well with those of the original equations. 18 The analytical solutions were weakly nonlinear due to the truncation approximation applied from the asymptotic expansion. Similar to Needham and Merkin (1984), Kranenburg (1992) added a turbulent diusion term into the shallow water theory so as to account for the non-uniform velocity prole near shocks. He scaled the dimensionless dependent variables by the asymptotic expansions with a second-order approximation. These expansions were substituted into the full equations of motion and the second-order amplitude evolu- tion equation was obtained by requiring all terms to be of the same order of magni- tude. The evolution equation derived by Kranenburg was shown to be able to predict roll-wave evolution with prescribed small, spatially periodic sinusoidal waves. Nev- ertheless, ill conditions of the equation was found when a subharmonic disturbance, which has smaller wavenumber and amplitude, was added to the sinusoidal initial condition. The solution with a subharmonic disturbance presented displacements of the original waves from their equilibrium positions. Consequently he concluded that the asymptotic method was not suitable for describing the continuing spatial wave evolution. Balmforth and Mandre (2004) considered eects of uneven bottom topography besides sources of the gravitational force, frictional resistance, and eddy viscous dissi- pation. The ow variables were expanded asymptotically with considering the topo- graphic eect. The nal evolution equation governed a new dependent variable which represented the ow velocity. The evolution equations was in the canonical form and the turbulent viscosity was related to the rst and third derivative of the dependent variable. The solutions of the evolution equations was shown to be consistent with 19 those of the full Saint Venant equations. The eect of the uneven bottom topography on the ow stability was explored utilizing the derived equation. Jin and Kim (2001) considered a model roll-wave equation that arose as the weakly nonlinear asymptotic approximation of the shallow water theory. Since this work mainly coped with the diculty of the numerical approximation to the model, it will be summarized in next section. Asymptotic expansions allow the derivation of amplitude evolution equations, which reduce the eort of directly solving the full shallow water theory for roll-wave solutions. They have proved to be able to converge to the solution of the original equations of motion. This method is, though, limited in analyzing the amplitude growth from small disturbances as a function of time. In addition to stability theories, mathematical solutions of progressing roll waves, and amplitude evolution equations, there have been several other analytical studies related to roll waves. Cristo et al. (2010) investigated the eect of the variability of the Darcy-Weisbach friction factor on the development of roll waves. They suggested that the assumption of constant friction factor led to underestimate of the roll-wave growth, according to models of time-asymptotic solutions of the linearized shallow water theory. Julien et al. (2010) analytically examined the growth of disturbances in a steep channel using a dynamic model of the shallow water theory. Their conclusion was that waves with smaller wavelength had the greatest amplication over a xed channel length, but the results had not been veried by experimental or eld data. 20 2.2 Experiment Brock (1969) conducted an experiment using a steel rectangular laboratory chan- nel. Figure 2.1 shows the overview of the experimental open channel. The upstream end of the channel was connected to a inlet box (Figure 2.2) feeding the steady water discharge. Three dierent channel slopes were designed in the experiment. Two or three dierent water ow rates were carried out under each scenario of the channel slope. Both naturally evolving processes of water surface perturbations and large quasi-periodic roll waves were considered in the study. Figure 2.1. General view of the laboratory open channel with water ows in Brock's experiment. Roll waves naturally evolved from the upstream to the end of the channel. 21 Figure 2.2. View of the inlet box and laboratory open channel near the inlet in Brock's experiment. The inlet box fed steady water discharge to the channel. The normal depth was measured by a wire gauge that indicated the turbulent water surface was higher or lower than a prescribed height. In other words, when the wire was set at a know distance from the channel bottom, the wire gauge could only tell whether the wire was in the water or not. If the wire was in the water half of the time series recorded by the gauge, the normal depth was determined as that wire elevation. Major roll-wave properties including the wave crest depth, wave trough depth, wave period, and wave celerity were measured and analyzed in the experiment. The wave trough depth, wave period, and wave celerity along the channel were measured from water pressure records as a function of time at the channel bottom. Those pressure records were taken by a pressure transducer and shown by an oscillograph 22 system. Pressure records were taken concurrently at two separated stations in order to measure the wave celerity. The time-averaged values of those properties were obtained at various channel locations by calculating the arithmetic means of the time series. The time-averaged wave crest depth of small perturbations were also obtained from pressure records. However, the crest depth of large shock-type waves were determined by a wire gauge and pressure records simultaneously. The reason was that the pressure distribution at a non-breaking jump was not hydrostatic. Therefore the wave crest depth could not be determined even if the bottom pressure was accurately measured. Both pressure records and wire gauge readings were taken at the middle of channel sections. 200 values were used to determine the arithmetic mean and frequency distribution of the wave trough depth and wave celerity while the time-averaged wave period at each station was based upon 1,500-2,000 values. For the wave crest depth, the frequency distribution was rst constructed and the arithmetic mean was then computed from it. Development curves for the major properties of roll waves were obtained through approaches introduced above. They were described as a function of the channel location. Those properties were made dimensionless by the transformation h c 0 = h c h 0 ; (2.14) h t 0 = h t h 0 ; (2.15) T 0 = sinT r g h 0 ; (2.16) c 0 = c p gh 0 : (2.17) 23 h denoted the water depth, subscripts \c" and \t" represented the wave crest and trough depth respectively,h 0 was the normal depth,T was the wave period, was the channel inclination,g was the gravitational acceleration, andc was the wave celerity. The upper bar sign denoted the time-averaged quantity while the prime sign implied the dimensionless value. Development curves were constructed to specify the relation- ship between those dimensionless roll-wave properties and the dimensionless channel location,x 0 =x/h 0 , wherex was the spatial coordinate along the channel. The char- acteristics of development curves of four roll-wave properties were summarized herein on the basis of the Brock's experimental outcome. Time-averaged wave crest depth h c 0 : The h c 0 curve showed that the growth rate with respect to the channel location was low near the inlet, and increased as a quasi- linear behaviour at more downstream. Finally the curve became concave downward presenting low growth rate. It was found that the location where shock waves began to form almost coincided with where the wave overtaking started. The dimensionless coordinate value x 0 for the starting point of the overtaking was the same for ows of dierent discharges at a xed channel slope. Time-averaged wave trough depth h t 0 : The development curve for h t 0 showed a rapid decay initially, and slid smoothly to a constant level. The measurement of h t 0 was not suciently complete, as it was missing at either upstream or downstream for some cases. Time-averaged wave period T 0 : The T 0 curve showed a constant value at the beginning, and grew linearly right after the wave overtaking began. This indicated the increase of the wave period resulted from coalescences of multiple perturbations. 24 Time-averaged wave celerity c 0 : The development feature of c 0 was shown to be consistent with the growth ofh c 0 , which was expected for shallow water waves. More- over, c 0 at a given x 0 decreased when the normal depth of the ow was reduced. The developments of h c 0 , h t 0 , and T 0 versus x 0 all presented a unique fashion for varying Froude number at a xed channel slope. The growth rate ofh c 0 as a function of x 0 was greater at steeper channels as the Froude number was also larger. In all of the runs within this experiment, it was found that the time-averaged wave crest depth was still increasing at the end of the channel. This implied roll waves had not fully evolved because the experiment was limited by the space of the laboratory. Therefore, it was not determined whether the roll-wave amplitude would stop growing and where was the stopping point if it existed. The experimental data, particularly the development of h c 0 and h t 0 versus x 0 , obtained by Brock will be used for validation tests of the present computational study. The comparison results will be given and discussed in Chapter 4. Plumerault et al. (2010) performed an experiment to study shallow water ows over an inclined sinusoidal-shape bottom. A parametric analysis was implemented based upon experiments with varying Reynolds number and channel slope. They distinguished two types of travelling waves: roll wave and pulse wave. They also established ranges of parameters of Re and slope I in which one type of wave was dominant over the other or two types coexisted. 25 2.3 Computational Dynamics Since the roll-wave phenomenon had been discovered, there were only a few com- putational studies on this subject. Here we brie y summarize their approaches and outcomes. Detail information including the governing equations, numerical schemes, and discussions of results are provided by referred works. Brook et al. (1999) developed a code to explore the possibility of roll waves occur- ring in collapsible tubes. The conservation laws in a tube were applied as governing equations of this study. The second-order Godunov scheme specied by Falle (1991) was used to solve the equations. The nonlinear Riemann solver in this scheme was devised by van Leer (1977). Their code was tested by exact solutions of the dam break problem and a nonlinear articial tube law and the comparison with MacCor- mack schemes. They performed a numerical experiment of progressing roll waves in an inclined open channel. The numerical solutions were compared with the linearized theory of amplitude growth and decaying rates and Dressler's mathematical solution (1949). At the end, they studied instabilities in inclined collapsible tubes. In order for the convergence test of the numerical scheme, the present study utilizes the pro- cedure of the numerical experiment of progressing roll waves as well as comparisons between theories and numerical solutions described in their paper. The detail of the convergence test will be given in Chapter 4. As mentioned in Subsection 2.1.3, Jin and Kim (2001) analyzed the behaviour of numerical solutions to a model roll-wave equation which arose as a weakly nonlinear asymptotic approximation of the shallow water theory. A so called round-o error would build up and destroy the long time behaviour of the roll-wave solution, if the 26 cell-averaged source was used. In order to overcome this diculty, they utilized a numerical scheme developed by Jin (2001) in which the interface method was applied for approximating the source term. Hence they achieved steady-state mass preserving due to the consistency between the ux and source term. The interface method proved to be more accurate than the cell-averaged source method, especially for quasi-steady numerical solutions. Que and Xu (2006) used a kinetic scheme, originally developed for the Boltzmann equation, to solve the Saint Venant equations. The BGK kinetic model (Bhatnagar et al. 1954) was utilized to approximate the ux at cell interfaces after the integration based on the nite volume method. The identical procedure of progressing roll waves used by Brook et al. (1999) was conducted to test their numerical scheme. The schemes was also applied to simulate the run-up of solitary waves and the solutions were compared with Synolakis's experimental data (1987). The above three computational studies merely concentrated on the evolution from small amplitude perturbations to progressing roll waves over a one-dimensional sim- ulation domain. Those approaches were eligible to verify the accuracy of a numerical scheme, but it was not related to natural roll-wave ows which develop spatially in an open channel. Next we review two literatures numerically simulating real roll-wave ows. Chang et al. (2000) numerically solve the viscous Saint Venant equations used by Needham & Merkin (1984). They found the time-averaged wave properties at every channel location is shown to be self-similar. This self-similar coarsening dynamics was caused by coalescences among propagating perturbations. They also suggested 27 a generic cascade coalescence model for wave separations to capture the self-similar family. Zanuttigh and Lamberti (2002) numerically analyzed the spatial development of roll waves by modeling the equations for a non-prismatic open channel. The weighted average ux method was adopted to establish the numerical scheme. In the validation tests, the code was applied to nd the steady ow conditions in a open channel with variable width and slope. The code was utilized to simulate roll-wave ows in Brock's experiment. However the computed roll-wave properties were not compared with the experimental results. They carried out sensitivity analyses of their numerical model to the mesh size. Both constant and variable friction factor were used in their study. 28 Chapter 3 Methodology of Model Construction As mentioned in Chapter 1, the center of this study is to establish a computational model. The computational model should be capable of simulating entire processes of the generation of instabilities and the spatial evolution of roll waves in watercourses down a wide, rectangular open channel. The model basically comprises the equations of motion and a numerical scheme. The equations of motion mathematically spec- ify water ows in a free surface channel and the numerical scheme is committed to solving the equations. In this chapter, we discuss in detail the approach of the model construction. 3.1 The Equations of Motion Water ows in a wide, rectangular open channel are illustrated by Figure 3.1. The problem is two-dimensional because the transverse component of the ow velocity is small compared to the longitudinal component. The hydraulic conditions including the water depth and ow velocity do not alter at dierent proles. Therefore the appropriate system to model this problem is the one-dimensional viscous Saint Venant equations: @h @t + @ @x (hu) = 0; (3.1) @ @t (hu) + @ @x hu 2 + 1 2 g cosh 2 =g sinh b + 1 @F xx @x : (3.2) 29 Water flow x z g g ( , ) h x t θ ( , ) u x t Horizontal Figure 3.1. Sketch illustration of the open channel ow specied by the Saint Venant equations. The uneven channel bed indicates the bottom roughness. h (x;t) is the water depth, u (x;t) is the depth-averaged ow velocity, represents the angle of channel inclination, and g denotes the gravitational acceleration. x and z are the spatial coordinates. t denotes time,x is the spatial coordinate along the channel bed,g is the gravitational acceleration, represents the angle of channel inclination with a slope tan, and is the water density. The dependent variables of the system are the water depth h (x;t) and depth-averaged ow velocity throughout cross-section u (x;t). b is the stress pertinent to the bottom drag between the water boundary layer and channel bed. F xx is the turbulently viscous force shearing normal to cross-section. The water depth and ow velocity in the equations are time-averaged such that the contributions from the turbulent uctuations are not taken into account. The conservation of mass (equation 3.1) can be deduced from the Reynolds trans- port theorem balancing the rate of change and net eux of mass. The conservation of momentum (equation 3.2) is derived from Newton's second law balancing the total 30 force and the summation of the rate of change and net eux of momentum. Two approximations are adopted in the derivation: the ow velocity is depth-averaged and uniform throughout cross-section; the water pressure obeys a hydrostatic distribution and this is how the term @ 1 2 g cosh 2 @x comes from. The second approximation is inferred based on the assumption that the vertical com- ponent (i.e. z direction) of the particle acceleration has a negligible eect on the water pressure. Those approximations are commonly known as the shallow water theory. It is so called because the approximations are valid when the horizontal dimension is much larger than the vertical scale. Since the derivation of the Saint Venant equa- tions requires integrations over the water depth, the system can be also obtained from depth-integrating the incompressible Navier-Stokes equations. Not only can the system model open channel ows, but it is able to describe other shallow water ows, for instance, tidal bore propagation and wave run-up on a sloping beach. The conven- tional derivations of the Saint Venant equations are given by Stoker (1958), Whitham (1999), or Billingham (2000). External sources consist of the gravitational force and bottom drag. The water particles in an open channel are subject to the gravitational acceleration due to the inclination. The gravitational acceleration substantially leads to the uid motion along the channel bed. The eect of the gravity is modeled by g sinh. To model the bottom drag, we utilize the relationship among the mean ow veloc- ity, hydraulic radius (water depth in wide, rectangular channels), and channel slope 31 specied by the Ch ezy formula for turbulent ows: b =C d ujuj; (3.3) whereC d is the dimensionless drag coecient and can be evaluated using the Darcy- Weisbach friction factor C d = f 8 (3.4) which is an empirical function of the ow variables and bottom roughness height. In many previous roll-wave literatures, the drag coecientC d was treated as a constant. Nevertheless this hypothesis is only valid in the steady-uniform ow where ow vari- ables do not change. In real roll-wave ows, ow variables have sudden variations at the wave front region. It is essential to use variable C d in the computational simu- lation of roll-wave ows. The detail calculation of C d will be provided in Chapter 4. The bottom drag as the ow resistance always acts in the opposite direction of the ow velocity, and this is the reason for adding the absolute value sign on one of u in equation 3.3. In most studies related to open channel ows, the viscous eect is neglected from the governing equations. The equations are hyperbolic without the diusive viscosity and produce discontinuous shock-type solutions. Disregarding the viscous eect is reasonable when the ow variables are near uniform along the channel. Roll-wave ows are very non-uniform with large slope of water depth and ow velocity. Hence the viscous force should be added to the equations. It is observed that considerable turbulent eddies occur on the crest of roll waves. The diusive transport and dissi- 32 pation of energy due to eddies are signicant and can be specied by the turbulent viscosity term. Jeereys (1925) did not add the viscous eect as deriving the linearized theory of stability. Dressler (1949) omitted the viscosity to construct piecewise continuous roll-wave solutions which are connected by discontinuities. Yu & Kevokian (1992) also neglected the viscous eect to obtain the weakly nonlinear roll-wave solutions of an amplitude evolution equation. In the numerical simulation of real roll-wave ows, Zanuttigh & Lamberti (2002) did not consider the viscosity in their governing equations. The rst work emphasizing the role of the viscous eect in analyzing roll waves were performed by Needham & Merkin (1984). They represented the eect of energy dissipation by the diusion of ow velocity. Consequently the work of Dressler was extended and they were able to acquire continuous periodic roll-wave solutions. Following their study, several other works included the turbulent viscosity, for example, Hwang & Chang (1987), Kranenburg (1992), Chang et al. (2000), and Balmforth & Mandre (2004). Noble (2007) and Johnson et al. (2011) analyzed the theoretical stability of viscous roll waves. We model the turbulently viscous force using the common approach for dening the normal stress of uids of incompressible ow but replacing the kinematic viscosity by a turbulent viscosity with the same dimension F xx = Z h 0 t @u @x dz = t h @u @x (3.5) where t is the constant turbulent viscosity andz is the vertical coordinate. This ex- pression is consistent with what had been used by Kranenburg (1992) and Balmforth 33 & Mandre (2004) who provided more details about the theoretical background of its rationality. Developed open channel ow instabilities are shock-type waves with tremendous adverse gradient of water depth (near discontinuity) and curvature of ow velocity at the vicinity of the wave front region. Thus the normal stress specied by equation 3.5 must have a signicant contribution on the momentum diusion, and the magnitude of diusivity is t weighted by the water depth h. The problem of using the turbu- lent viscosity is the absence of empirical information for quantifying t . Although Needham & Merkin (1984), Kranenburg (1992), and Chang et al. (2000) suggested variations in t did not aect nonlinear solutions in a more qualitative sense, we shall show in Chapter 4 simulation results of the roll-wave spatial evolution are substan- tially in uenced by the choice of the turbulent viscosity. We shall also nd the value of t that makes model results agree with Brock's experimental data of the wave amplitude satisfactorily. After substitution from equations 3.1, 3.3, and 3.5, momentum conservation equa- tion 3.2 can be written as @u @t +u @u @x +g cos @h @x 1 h @ @x t h @u @x =g sinC d ujuj h : (3.6) The complete derivation of the Saint Venant equations 3.1 and 3.6 for water ows in a wide, rectangular open channel is given in Appendix B. The derivation is based on the Reynolds transport theorem and Newton's second law and is somewhat dierent from approaches in the literatures. The equations 3.1 and 3.6 form a system of nonlinear, nonhomogeneous, parabolic- 34 hyperbolic partial dierential equations. The full analytical solution to the system is intractable unless certain terms are truncated in special cases. We next develop a conservative high-resolution numerical scheme for the system capable of capturing discontinuous shocks and retaining physical properties in the computational domain. 3.2 Numerical Scheme 3.2.1 Finite Volume Formulation The purpose of this section is to develop a numerical scheme to nd the conser- vative ow variables described by the equations of motion 3.1 and 3.6. We consider the conservation laws @Q (x;t) @t + @F c (Q) @x + 1 h @F d (Q; Q 0 ) @x = S (Q) (3.7) where Q (x;t) is the vector of ow variables, F c (Q) is the vector of the convective ux, F d (Q; Q 0 ) is the vector of the diusive ux, and S (Q) represents the source vector. According to equations 3.1 and 3.6, we can dene: Q = 0 B B @ h u 1 C C A ; (3.8) F c = 0 B B @ hu 1 2 u 2 +g cosh 1 C C A ; (3.9) 35 F d = 0 B B @ 0 t h @u @x 1 C C A ; (3.10) and S = 0 B B @ 0 g sinC d ujuj h 1 C C A : (3.11) The equations of motion governing water ows in a wide, rectangular open channel is a one-dimensional system. The only argument of space isx. The problem specied by equation 3.7 is supposed to be mathematically well-posed. Depending on proper initial and boundary conditions, the solution to Q (x;t) can be uniquely obtained inside the domain x2 [0;x max ] and for time t> 0, where x max 2R + . Figure 3.2 presents the schematic of how the computational domain is discretized. The discretization is based on the nite volume method. The computational domain is divided into N control cells with equal spacing x, where N2Z + . Each control cell has two interfaces. The cell interfaces separate adjacent two control cells. For instance, the jth (j = 1; 2; ;N) cell has the space x2 x j1/2 ;x j+1/2 where x j1/2 = (j 1) x (3.12) is the location of the interface between the j 1th and jth cell. The conception of the nite volume method is attributed to Godunov's scheme (Godunov 1959). The ow variables are assumed to be uniform inside each control cell. The primary problem is to determine adequate numerical ux function that approximates the exact ux on the cell interface, and the only data available are the 36 Cell interface Outflow boundary j 1 2 j x − Control cell Inflow boundary 1 2 j 1 − 1 2 j x + Figure 3.2. Sketch illustration of discretization of the one-dimensional computational domain. The discretization is based on the nite volume method. The domain is divided into many control cells with equal spacing. Two adjacent control cells are separated by a cell interface. cell-averaged ow variables (Leveque 2002). To derive the discretized equation using the nite volume method, equation 3.7 is integrated over the jth control cell to give dQ j dt + 1 x h (F c ) j+1/2 (F c ) j1/2 i + 1 x 1 h j h (F d ) j+1/2 (F d ) j1/2 i = S j (3.13) where Q j = 1 x Z x j+1/2 x j1/2 Qdx (3.14) is the cell-averaged vector of ow variables on the jth cell, h j = 1 x Z x j+1/2 x j1/2 hdx (3.15) 37 is the cell-averaged water depth throughout the jth cell, and S j = 1 x Z x j+1/2 x j1/2 Sdx (3.16) is the cell-averaged source vector on the jth cell. (F c ) j1/2 and (F d ) j1/2 denote the convective and diusive ux on the cell interface x j1/2 . Figure 3.3 shows the formulated quantities specied by equation 3.13. It should be emphasized that the source term is only dependent on the ow variables. Therefore the source term is also assumed uniform inside control cells. The nite volume formulation also requires integrating equation 3.13 over a time step from t ` to t `+1 , which results in Q `+1 j = Q ` j 1 x Z t `+1 t ` h (F c ) j+1/2 (F c ) j1/2 i dt 1 x Z t `+1 t ` 1 h j h (F d ) j+1/2 (F d ) j1/2 i dt + Z t `+1 t ` S j dt (3.17) where Q `+1 j represents the vector of ow variables on the jth cell at the future time level t `+1 . The merit of using the nite volume formulation in the present study is the con- servation of physical properties specied by equations 3.1 and 3.6. The uxes of those physical properties through the in ow and out ow boundaries are retained, as long as the uxes leaving each control cell and entering the adjacent cell are evaluated in a consistent manner (Eymard et al. 2000; Versteeg & Malalasekera 2007). The cell-averaging treatment of the nite volume formulation results in piecewise constant ow variables over control cells with a single discontinuity on interfaces. Our goal is to update Q `+1 j provided ow variables at the current time level t ` are known 38 x Δ ( ) 1 2 c j− F ( ) 1 2 d j− F 1 2 j x − 3 2 j x − , j j Q S 1 1 , j j − − Q S x Δ ( ) 1 2 c j+ F ( ) 1 2 d j+ F 1 2 j x + 3 2 j x + , j j Q S 1 1 , j j + + Q S x Figure 3.3. Sketch illustration of the formulated quantities described by equation 3.13, based on the nite volume method. Q j1 , Q j and Q j+1 are the cell-averaged values of the vector of ow variables on each control cell, while S j1 , S j and S j+1 denote the cell-averaged values of the source vector. (F c ) j1/2 , (F c ) j+1/2 , (F d ) j1/2 , and (F d ) j+1/2 represent the convective and diusive uxes on cell interfaces. on all control cells. We need to calculate the convective ux associated with the hyperbolic part by solving the Riemann problem locally at cell interfaces (Leveque 2002). 3.2.2 The Riemann Solver A typical hyperbolic system of conservation laws without source is written as @Q @t + @F @x = 0 (3.18) where Q is the dependent vector and F is some vector-valued function of Q. The Riemann problem for the system equation 3.18 is dened to be the initial value 39 problem with special initial data. The initial data, Q (x; 0), are discontinuous with a single jump as following: Q (x; 0) = 8 > > < > > : Q left ; x<x interface Q right ; x>x interface : (3.19) Seeking the dependent vector at subsequent time levels is actually studying the break-up of the discontinuity. The initial data decomposes into a shock and a rar- efaction which propagate at opposite directions. In order to obtain the information of the break-up, we require a Riemann solver to compute the state of F at the location of x interface . More detail descriptions about the Riemann problem are given by Toro (2009). The convective ux F c represents the hyperbolic part of the conservation laws equation 3.7. The knowledge of wave speeds and jump conditions is usually essential to solve the Riemann problem for a nonlinear system. Let J (Q) = @F c @Q = 0 B B @ u h g cos u 1 C C A (3.20) be the Jacobian matrix of the convective ux, based on equations 3.8 and 3.9. Then the wave speeds are the eigenvalues of J and the projections of Q = Q left Q right onto the eigenvectors of J are the jumps associated with the state on the cell interface. Q left and Q right represent the ow variables on control cells at the left and right of the interface. However the ow variables are unknown on the interface, and so is J. Roe (1981) considered the way of exact solutions to an approximate problem. The 40 attempt was made to nd e J (Q left ; Q right ) denoting the approximate Jacobian matrix on the cell interface. e J was acceptable only if the following properties were satised: (i) It constitutes a linear mapping from the vector space Q to the vector space F c ; (ii) As Q left ! Q right ! Q, e J (Q left ; Q right )! J (Q); (iii) For any Q left , Q right , e J (Q left ; Q right ) (Q left Q right ) = (F c ) left (F c ) right ; (iv) The eigenvectors of e J are linearly independent. Those properties were called property U, meaning the intention to ensure uniform validity across discontinuities. Hence the analytical expressions of the eigenvalues and eigenvectors of e J could be obtained. Yet the direct construction of e J tended to be dicult. Even though the formulae of e J was derived, the algebra became almost impossible to carry through without error. Therefore, a technique of the parameter vector was inspired. We introduce the parameter vector w = 0 B B @ p h u 1 C C A = 0 B B @ w 1 w 2 1 C C A (3.21) such that any jump in the space of Q and F c can be specied in terms of their images in the space of w. The vector of ow variables is described by the image w as Q left Q right = e B (w left w right ) (3.22) where e B = 0 B B @ 2 w 1 0 0 1 1 C C A : (3.23) 41 The upper bar stands for the arithmetic mean of values at the left and right of the interface, i.e. w 1 = (w 1 ) left + (w 1 ) right 2 = p h left + p h right 2 : (3.24) Likewise for F c we can write (F c ) left (F c ) right = e C (w left w right ) (3.25) where e C = 0 B B @ 2 w 1 w 2 w 2 1 2g cos w 1 w 2 1 C C A : (3.26) It is easy to detect e J = e C e B 1 . In order to nd eigenvalues e we can formulate e C e B 1 e I Q = 0 (3.27) which is equivalent to e C e e B w = 0 (3.28) and consequently the characteristic equation e C e e B = 0: (3.29) I represents the identity matrix. Then the eigenvalues of e J are dened as: e 1 = u c; e 2 = u + c (3.30) 42 where u = w 2 = u left +u right 2 (3.31) and c = q g cosw 2 1 = r g cos h left +h right 2 : (3.32) The corresponding right eigenvectors are obtained by solving equation 3.28 with found e and mapping into Q-space. The results are e e 1 = 0 B B @ 1 g cos c 1 C C A ; e e 2 = 0 B B @ 1 g cos c 1 C C A : (3.33) There are numerous Riemann solvers that provide the numerical function of con- vective ux on the cell interface. They can be categorized into rst- and higher-order methods in terms of the truncation error. Higher-order schemes usually produce more accurate solutions since rst-order methods induce a high level of false diu- sion. Nevertheless many numerical studies have corroborated regular higher-order schemes for hyperbolic conservation laws have stability problems, e.g. (Burguete & Garc a-Navarro 2001; Gasc on & Corber an 2001). Their solutions are oscillatory, espe- cially when they are used to solve turbulent problems. Those non-physical oscillations tend to occur near discontinuous solutions in most cases. Therefore neither rst-order schemes nor regular higher-order schemes are shock-capturing. Harten (1983) suggested the monotonicity property of a Riemann solver was spec- ied as: (i) No new local extrema in the computational domain may be created; 43 (ii) The value of a local minimum is nondecreasing, and the value of a local maxi- mum is nonincreasing. Those schemes short of the monotonicity property generate spurious oscillations (wig- gles) on the numerical solutions. Harten proved the numerical scheme for hyperbolic conservation laws must be TVD (total variation diminishing) so as to guarantee the monotonicity property. The total variation of a numerical scheme is dened as TV (Q) = N X j=2 jQ j Q j1 j; N2Z + : (3.34) A scheme is TVD if the total variation does not increase with time: TV Q `+1 6TV Q ` ; `2 0 &Z + : (3.35) We require a second- or higher-order TVD Riemann solver to cope with shock- type roll waves and discontinuous initial data existing on cell interfaces. We utilize an explicit TVD Lax-Wendro scheme that has been successfully tested for solving the Euler's equations in aerodynamics (Yee 1987) and the Saint Venant equations in open channel ows (Delis & Skeels 1998). The predecessor of this scheme was devised by Sweby (1984) and Roe (1984) to solve the linear scalar convection equation. It was extended to hyperbolic systems of conservation laws by Yee (1987) and can be viewed as the usual Lax-Wendro scheme supplemented by the addition of a conser- vative numerical dissipation. The numerical dissipation counteracts those spurious oscillations. The scheme is spatially symmetric due to the central-dierencing algo- rithm yet holds the transportiveness property because of the ux limiter built into 44 the numerical dissipation term. The scheme has the same accuracy as the original Lax-Wendro scheme which is second-order in space and time. Therefore the time integral of the convective ux on cell interfaces in equation 3.17 is written as 1 t Z t `+1 t ` (F c ) dt = (F c ) ` (3.36) where t is the computational time step. According to the TVD Lax-Wendro scheme, (F c ) = 1 2 h (F c ) left + (F c ) right e P D i (3.37) where (F c ) left and (F c ) right are the convective ux on control cells at the left and right of the interface. e P = 0 B B @ 1 1 g cos c g cos c 1 C C A (3.38) is a 2 2 matrix composed of the right eigenvectors of e J ,e e 1 ande e 2 . D is a 2 1 dissipation vector whose elements are expressed as d v = t x e v 2 v + e v (1 v ) v (3.39) with v = 1; 2. e is a continuously dierentiable function to prevent entropy violation (Harten & Hyman 1983) and has a form e = 8 > > < > > : e ; e >" e 2 +" 2 2"; e <" : (3.40) " is a small nonnegative quantity which measures the violation of the entropy con- 45 dition. Motivated by Oleinik's entropy condition, " varies at dierent locations and time levels and is determined by " = max 0; e left ; right e : (3.41) left and right are the eigenvalues of the exact Jacobian matrix of the convective ux at the left and right of the interface. Based on equation 3.20, the eigenvalues of J are written as 1 =uc; 2 =u +c (3.42) where c = p g cosh. is the wave strength of the local characteristic variables dened as = 0 B B @ 1 2 1 C C A = e P 1 (Q r Q l ) (3.43) and e P 1 = 0 B B @ 1 2 1 2 c g cos 1 2 1 2 c g cos 1 C C A : (3.44) In equation 3.39, represents the ux limiter. The current TVD Lax-Wendro scheme is a one-parameter family so that a wide group of ux limiters can be ac- cepted. There have been many studies dealing with the development of ux limiters since 1970s. Some of them include van Leer (1974), Chakravarthy & Osher (1983), Sweby (1984), Roe (1986), and Lien & Leschziner (1994). All of the developed ux limiters are described as functions of a ratio r , i.e. = (r ). For nonlinear sys- tems of governing equations, r should be the backward ratio of the product of the approximate eigenvalue and wave strength on consecutive interfaces. For example, 46 on the cell interface x j+1/2 r j+1/2 = e j1/2 j1/2 e j+1/2 j+1/2 : (3.45) Since the original purpose of adding the ux limiter is to ensure TVD and coun- teract the wiggles, the criterion of ux limiters proposed by Sweby (1984) must be satised for producing a second-order TVD scheme. The criterion constrains the relationship between the ux limiter (r ) and the ratio r as following: r 6 (r )6 2r ; if 0<r < 1 2 ; r 6 (r )6 1; if 1 2 6r < 1; 16 (r )6r ; if 16r < 2; 16 (r )6 2; if r > 2; and (r ) = 0; if r 6 0: This criterion is presented by Figure 3.4, which is a graph specifying the relationship between (r ) and r . We choose the ux limiter invented by van Leer (1974): = (r ) = r +jr j 1 +jr j : (3.46) 47 3 2 * * 2r ψ = * * r ψ = * ψ * r 0 1 2 4 6 0 1 * ψ Figure 3.4. The graph of the criterion for a second-order TVD scheme proposed by Sweby (1984). The criterion describes the relationship between the ux limiter (r ) and the ratio r . The ux limiter must lie within the shadow region to ensure that a scheme is second-order and TVD. This limiter is a smooth, monotone increasing function of r and has a symmetric property (r ) r = 1 r : (3.47) Its r curve settles inside the shadow region of Figure 3.4, so it complies with the second-order TVD criterion. Note that r goes to innity when the denominator of equation 3.45 is zero and (r ) lies on the upper boundary of Sweby's second-order TVD region, which is 2. Also (r ) remains at zero when r is nonpositive. 48 Owing to the CFL stability condition for the explicit scheme max v j t x 6Cr< 1; (3.48) the time step alters and is determined before the computation at each time level. The selection of the Courant number Cr will be discussed under the convergence sensitivity analysis in Chapter 4. 3.2.3 Boundary Conditions The computational model is constructed for simulating real roll-wave ows in a wide, rectangular open channel. Large roll waves evolve from small perturbations on the water surface of the steady-uniform ow. The small perturbations exist be- cause turbulent ows are characterized by uctuations of ow variables. The ow variables in the present equations of motion are time-averaged without considering time-dependent uctuations. Nevertheless small perturbations must be imposed at the in ow boundary so as to simulate roll-wave ows. Hence the water depth at the in ow boundary is expressed as the normal water depth supplemented by a small perturbation: h i.b. =h 0 + i.b. (t): (3.49) h i.b. denotes the water depth at the in ow boundary, h 0 is the normal water depth, and i.b. (t) stands for the time-dependent small perturbation to the normal depth at the in ow boundary. The small perturbation is set to be random with a normal 49 probability distribution function. The mean of i.b. (t) is zero, resulting in the normal depth at the in ow boundary after the averaging treatment over a sucient time period. The ow rate at the in ow boundary is constrained to be q 0 =h 0 u 0 whereq 0 andu 0 are the ow rate and ow velocity in the steady-uniform condition, respectively. Therefore the ow velocity at the in ow boundary is calculated as u i.b. = q 0 /h i.b. . Since the ow variables are prescribed at the in ow boundary, the convective ux F c at the in ow boundary can be easily obtained without using the approximate Riemann solver discussed in Subsection 3.2.2. As specied in Subsection 3.2.2, the application of the TVD Lax-Wendro scheme requires the knowledge of the ux limiter on cell interfaces. The ux limiter on the cell interfacex 3/2 , which is the downstream interface of rst control cell, is computed as a function of r 3/2 and r 3/2 = 1/2 1/2 e 3/2 3/2 : (3.50) To nd the ratio r 3/2 , the quantities of the eigenvalue and wave strength must be determined on interfaces x 1/2 and x 3/2 in which x 1/2 is the location of the in ow boundary. The denominator of equation 3.50 is determined through the method introduced in Subsection 3.2.2. In the numerator, 1/2 is evaluated exactly based on equation 3.42 because the ow variables are prescribed at the in ow boundary. The wave strength at x 1/2 is computed as 1/2 = P 1 1/2 (Q 1 Q mirror ) (3.51) where P 1 1/2 is the inverse of a matrix comprising the right eigenvectors of the exact 50 Jacobian matrix of the convective ux. The exact right eigenvectors can be obtained according to equation 3.20, so P 1 1/2 = 0 B B @ 1 2 1 2 c g cos 1 2 1 2 c g cos 1 C C A 1/2 : (3.52) Q mirror represents the vector of ow variables at the imagined mirror cell. The mirror cell is located at the upstream of the in ow boundary. The ow variables at the mirror cell are evaluated using the rst-order extrapolation. That is Q mirror = 2Q 1/2 Q 1 : (3.53) The rst-order extrapolation treatment at the in ow boundary is shown by Figure 3.5. We can notice that the convective ux at x 3/2 is calculated based on Q mirror , Q 1/2 , Q 1 , and Q 2 . The most downstream of the computational domain is the out ow boundary. The out ow boundary is articially imposed because the equations can be only solved in a bounded domain. The out ow boundary should absorb any signal reaching the boundary. In other words, outgoing waves should leave the domain cleanly so that no spurious re ections are generated from the out ow boundary into the computational domain. Such out ow boundary condition is known as Non-Re ecting Boundary Condition (NRBC). The NRBC is achieved through applying the zeroth-order ex- trapolation at the out ow boundary. We imagine the control cell at the downstream of the out ow boundary is the N + 1th cell. The ow variables are equal at the Nth andN +1th cell: Q N = Q N+1 . The zeroth-order extrapolation is recommended since 51 Inflow boundary ( ) 3 2 mirror 1 2 1 2 , , , F Q Q Q Q mirror Q 1 2 Q 1 Q x 1 2 x 1 3 2 x 2 ( ) 3 2 mirror 1 2 1 2 , , , F Q Q Q Q Mirror cell Figure 3.5. Schematic of applying the rst-order extrapolation on the ow variables at the in ow boundary of the computational domain. The mirror cell at the upstream of the in ow boundary is imagined in order to compute the convective ux at the location x 3/2 . The ow variables at the mirror cell is calculated using the rst-order extrapolation technique. it does not have stability problems (Leveque 2002). Hence the convective ux at the out ow boundary, which is the interface x N+1/2 , is equal to the ux at the Nth cell: (F c ) N+1/2 = (F c ) N . 3.2.4 Solution Procedure The water depth on all control cells at the future time level,h `+1 j , can be computed once the mass and momentum convection (F c ) is obtained. This is accomplished by explicitly solving the rst part of the vector equation 3.17 and the remaining problem is the second part. 52 The time integrals of the diusive ux and source vector in equation 3.17 are approximated using the trapezoidal method, that is 1 t Z t `+1 t ` 1 h j (F d ) dt = 1 2 " 1 h `+1 j (F d ) `+1 + 1 h ` j (F d ) ` # (3.54) and 1 t Z t `+1 t ` S j dt = 1 2 S `+1 j + S ` j : (3.55) In other words, the arithmetic mean values of 1/h j (F d ) and S j at `th and ` + 1th time levels are evaluated. The velocity gradient@u/@x inside the diusive ux on the cell interface is treated by the central dierencing method. For example, @u @x j1/2 = u j u j1 x : (3.56) Hence the semi-implicit Crank-Nicholson scheme (1947) is established for the diusive ux. The accuracy of this scheme is second-order in space and time in terms of the Taylor series truncation error. A von Neumann stability analysis manifests the Crank-Nicholson scheme is unconditionally stable (Flethcer 1991). Though we still comply with the CFL stability condition as this scheme may disobey the boundedness property if the time step is too large (Versteeg & Malalasekera 2007). The operation of the trapezoidal approximation on the time integral of the source term in the momentum equation gives rise to S `+1 j where S =g sinC d ujuj h : (3.57) 53 S `+1 j contains the nonlinear ow variable as well as drag coecient at the future time level. We perform the second-order conservative linearization using a local Taylor expansion about u ` j to express S `+1 j : S `+1 j =S ` j + @S @u ` j u `+1 j +O t 2 (3.58) in which @S/@u =2C d u=h and u `+1 j =u `+1 j u ` j . Finally, the water depth on the interface h , as a component of the momentum diusivity among the diusive ux, is obtained by taking advantage of the convective ux vector (F c ) calculated beforehand through equation 3.37. In detail, let F c = 0 B B @ hu 1 2 u 2 +g cosh 1 C C A = 0 B B @ F 1 c F 2 c 1 C C A ; (3.59) so a cubic equation is formulated for u as following u 3 2F 2 c u + 2g cosF 1 c = 0: (3.60) The above equation is analytically solved by the Cardano's method and the water depth h is calculated via h =F 1 c /u. Therefore a system of linear discretized equations, governing the ow velocity at the future time level u `+1 j , is constructed. It comprises the second part of the vector equation 3.17 on each control cell. Thus the number of equations is N. In order to solve the discretized equations, they are arranged into a matrix form which is written 54 as 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 b 1 c 1 a 2 b 2 c 2 : : : a j b j c j : : : a N1 b N1 c N1 a N b N 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 u `+1 1 u `+1 2 : u `+1 j : u `+1 N1 u `+1 N 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 e 1 e 2 : e j : e N1 e N 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ; (3.61) where a j , b j , c j , and e j are scalar coecients and blank elements denote zero values. The size of the coecient matrix is NN. This type of matrix system is known as the tridiagonal system that can be solved using the Thomas algorithm (Fletcher 1991) composed of two steps: forward sweeping and back substitution. The objective of the forward sweeping is to generate an upper-triangular coecient matrix from the original system. This operation eliminates all elements of a j and replaces b j with unity on the diagonal. It also updates elements of c j and e j . The forward sweeping produces 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1 c 0 1 1 c 0 2 : : 1 c 0 j : : 1 c 0 N1 1 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 u `+1 1 u `+1 2 : u `+1 j : u `+1 N1 u `+1 N 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 e 0 1 e 0 2 : e 0 j : e 0 N1 e 0 N 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 : (3.62) 55 c 0 1 = c 1 /b 1 and e 0 1 = e 1 /b 1 . For other equations, c 0 j = c j b j a j c 0 j1 and e 0 j = e j a j e 0 j1 b j a j c 0 j1 . Having the above modied system,u `+1 j is solved by the back substitution: u `+1 N = e 0 N and u `+1 j =e 0 j u `+1 j+1 c 0 j . The Thomas algorithm is a direct matrix solver which is more computationally ecient than iterative methods. It can be utilized due to the special structure of the tridiagonal system. The operation count for this algorithm is only 5N 4, which is clearly preferable to N 3 for the full Gaussian elimination. Nevertheless, there are unknowns inside the coecient matrix of the tridiagonal system, which are the water depth on cell interfaces at the future time level, h `+1 . This is because the time integral approximation equation 3.54 generates the implicit diusive ux, and (F c ) is calculated only at the current time level. To overcome this problem, we carry out the following predictor-corrector iterative approach: 1. Assume the water depth on the interface is equal at the current and future time level, i.e. h `+1 =h ` ; 2. Substituteh `+1 (orh ` ) into the coecient matrix of the tridiagonal system and solve for u `+1 j ; 3. Use the predicted u `+1 j along with the known h `+1 j to nd the convective ux (F c ) `+1 by the Riemann solver, and thus the water depth h `+1 on the interface; 4. Substitute the updatedh `+1 into the coecient matrix of the tridiagonal system and correct u `+1 j ; 5. Compare the corrected u `+1 j with the predicted u `+1 j . If they are equal, then the correctedu `+1 j has converged; if not, let the correctedu `+1 j be the predicted one, and repeat the step 3, 4 and 5 until u `+1 j converges. 56 3.3 Model Structure The content of this section provides an overview of the outcome of the present model construction, in the hope of beneting other researchers who may adopt this model and the author's future work on the model extension to other problems. A computer code, written in the language of FORTRAN 95, is developed to imple- ment the numerical scheme solving the equations of motion specied in this chapter. This computer program is utilized for analyzing the spatial evolution of roll waves in the current study. It is capable of simulating roll-wave ows in a wide, rectangular open channel given the ow rate, channel slope, and bottom roughness. However the equations can also describe other shallow water ows such as the dam-break wave, tidal bore propagation, and wave run-up on a sloping beach. The model can be scaled for problems of two-dimensional shallow water ows without much diculty. The Rie- mann problem is solved using a TVD Lax-Wendro scheme. Hence this program is named as \Shallow Water Analyzed by a TVD Scheme at University of Southern California" with an acronym \USC SWAT". The structure of USC SWAT is clearly arranged such that this model consists of a main program, six function subprograms, and ten subroutines, as shown by Figure 3.6. The obligation of the main program is to perform iterations with time until the prescribed maximum time level is reached. The main program also calculates the water depth and conducts the predictor-corrector iterative approach for the ow ve- locity at the future time level. Within each time step, the main program calls the function subprograms and subroutines to obtain required values. The function \INBH" and \INBU" calculates the random water depth and ow 57 MAIN Function Subprograms INBH INBU INBFC2 Subroutines NORMALCONDITION VELOCELE WAVESTRENGTH FLUXLIMITER1 MAIN INBFC2 INBC INB_ALPHA1 INB_ALPHA2 FLUXLIMITER2 CELLFACEFLUXC CELLFACEH DRAGCOEF FACTORIZATION SOLVE Figure 3.6. Illustration of the model structure of USC SWAT. The \MAIN" program is responsible for performing time-iterations and calling function subprograms and subroutines. Those function subprograms and subroutines execute respective assign- ments when they are called. Function subprograms and subroutines return to the \MAIN" program with computed values after execution. velocity at the in ow boundary. The function \INBFC2" calculates the second com- ponent of the convective ux at the in ow boundary. The function \INBC" cal- culates the wave speed at the in ow boundary. The function \INB ALPHA1" and \INB ALPHA2" calculates the wave strengths at the in ow boundary. The subroutine \NORMALCONDITION" computes the drag coecient, normal water depth and normal ow velocity at the steady-uniform condition. The sub- routine \VELOCELE" computes the ow velocity and wave speed on cell interfaces. The subroutine \WAVESTRENGTH" computes the wave strengths on cell interfaces. The subroutine \FLUXLIMITER1" computes van Leer's ux limiter (1974) on cell 58 interfaces and \FLUXLIMITER2" computes Sweby's ux limiter (1984) on cell in- terfaces. The subroutine \CELLFACEFLUXC" computes the convective ux on cell interfaces. The subroutine \CELLFACEH" computes the water depth on cell inter- faces using Cardano's method. The subroutine \DRAGCOEF" computes the drag coecient based on the roughness height and ow variables. The subroutine \FAC- TORIZATION" performs the forward sweeping for the tridiagonal matrix system and \SOLVE" executes the back substitution to compute the ow velocity on each control cell at the future time level. 59 Chapter 4 Convergence and Validation Tests Once the equations of motion are formulated and the numerical scheme is es- tablished, we need to verify the accuracy of the numerical scheme. In this chapter, tractable theoretical solutions associated with the current equations of motion are obtained and compared with numerical results. The solutions obtained through ana- lytical approaches consist of transient linear theory and quasi-steady nonlinear theory. After the convergence tests, the numerical results of real roll-wave ow simulations are compared with Brock's experimental data. 4.1 Comparison with Transient Linear Theory The comparison procedure between the theoretical and numerical solutions used by Brook et al. (1999) are adopted in this study. In order to derive the theory, we begin with the Saint Venant equations: @h @t +h @u @x +u @h @x = 0; (4.1) @u @t +u @u @x +g cos @h @x =g sinC d u 2 h : (4.2) This theory is only available when the normal stress term is eliminated from the momentum equation. We assume the ow is always in non-negative x direction such that u 2 =ujuj. The angle of channel inclination and drag coecient C d are taken 60 to be constants. In reality the drag coecient depends on ow characteristics and bottom roughness and it cannot be constant. Yet here we are comparing the theoret- ical and numerical solutions to a same set of equations. This enables us to make the assumption. The derivation of the transient liner theory requires us to impose small perturba- tions to the steady-unifomr state. The ow variables are written as the summation of quantities at the steady-uniform state and small perturbations: h (x;t) =h 0 + (x;t); (4.3) u (x;t) =u 0 + (x;t) (4.4) where constants h 0 and u 0 denote the water depth and ow velocity in the steady- uniform state. Variables (x;t) and (x;t) stand for small perturbations to the normal water depth and ow velocity. We substitute the above expressions of ow variables into equations 4.1 and 4.2 and notice h 0 and u 0 are constants. Therefore the equations specifying the small perturbations are obtained as @ @t +h 0 @ @x + @ @x +u 0 @ @x + @ @x = 0; (4.5) @ @t +u 0 @ @x + @ @x +g cos @ @x +C d u 2 0 (h 0 +) +C d 2 (h 0 +) +C d 2u 0 (h 0 +) g sin = 0: (4.6) In equations 4.5 and 4.6, any term including multiplications of small perturbations are considered to be nonlinear. We neglect all nonlinear terms to get the linearized 61 equations @ @t +h 0 @ @x +u 0 @ @x = 0; (4.7) @ @t +u 0 @ @x +g cos @ @x +C d u 2 0 (h 0 +) +C d 2u 0 (h 0 +) g sin = 0: (4.8) This step of linearization implies the reason why the theory we are looking for is linear. The elimination of nonlinear terms of small perturbations is valid when the magnitude of small perturbations is tiny. Therefore the linear theory cannot stand as the wave amplitude magnies, which will be corroborated in the following comparison results. When the ow is at the normal condition, the gravitational force on the water is balanced by the drag force exerted by the channel bottom. Hence we have C d =g sin h 0 u 2 0 (4.9) which is substituted into equation 4.8 to give @ @t +u 0 @ @x +g cos @ @x +g sin 2 u 0 h 0 = 0: (4.10) The conservation of mass equation 4.7 can be rewritten as @ @x = 1 h 0 @ @t +u 0 @ @x (4.11) which can be integrated to have = 1 h 0 Z @ @t +u 0 @ @x dx + (t) +C (4.12) 62 where (t) represents a possible function of t and C is an arbitrary constant. Then we substitute equation 4.12 into equation 4.10 and take rst derivative with x to yield @ @t +c @ @x @ @t +c + @ @x + 2g sin u 0 @ @t +c 0 @ @x = 0 (4.13) where c =u 0 p g cosh 0 ; (4.14) c + =u 0 + p g cosh 0 ; (4.15) and c 0 = 3 2 u 0 : (4.16) c 0 represents the kinematic wave speed. c and c + represent the slowest and fastest traveling speeds. Equation 4.13 is the linearized equation governing the small per- turbation of water depth . 4.1.1 Linearized Theory of Stability Since we are exploring the open channel ow instabilities, it would be interesting to theoretically analyze the stability criterion from equation 4.1 and 4.2. Equation 4.13 is rewritten as @ @t +c 0 @ @x = u 0 2g sin @ @t +c @ @x @ @t +c + @ @x : (4.17) 63 Applying the lower order approximation (Whitham 1974) @ @t 'c 0 @ @x (4.18) on the right hand side of equation 4.17, we are able to have @ @t +c 0 @ @x = u 0 2g sin (c 0 +c ) (c 0 +c + ) @ 2 @x 2 (4.19) which is equivalent to @ @t +c 0 @ @x = h 0 u 0 2 tan u 0 2g sin (c 0 u 0 ) 2 @ 2 @x 2 (4.20) when equation 4.14 and 4.15 are substituted. Equation 4.20 is a typical convection- diusion equation with a wave speed of c 0 and a diusivity of h 0 u 0 2 tan u 0 2g sin (c 0 u 0 ) 2 : Obviously the ow described by equation 4.20 would be stable if there is positive diusion. In other words, small perturbations of the water ow do not amplify if the diusivity is positive. Therefore the stability criterion becomes h 0 u 0 2 tan u 0 2g sin (c 0 u 0 ) 2 > 0 (4.21) or u 0 p g cosh 0 < 3 2 u 0 <u 0 + p g cosh 0 : (4.22) 64 u 0 p g cosh 0 < 3/2u 0 is automatically satised. From 3/2u 0 < u 0 + p g cosh 0 , we have u 0 < 2 p g cosh 0 (4.23) which is equivalent to u 0 p g cosh 0 =F 0 < 2 (4.24) where F 0 is the undisturbed Froude number at the steady-uniform condition of the water ow. Equation 4.24 manifests that water ows down a wide, rectangular open channel are stable when the Froude number associated with the normal condition is less than two. At the equilibrium state, the gravitational force is balanced by the bottom drag, so we have g cosh 0 =C d u 2 0 tan (4.25) which is substituted into equation 4.24 to give tan< 4C d : (4.26) Equation 4.26 is an equivalent stability criterion to equation 4.24, and it implies the channel slope must be small compared to the bottom drag force in order to achieve the stability. The open channel ow is supercritical whenF> 1. When open channel ow instabilities and roll waves occur, the ow is highly supercritical, i.e. F 0 > 2, and the highly supercritical ow results from very steep slope and smooth channel bottom. 65 4.1.2 The Transient Linear Theory The linearized equation governing the small perturbation of water depth has been inferred from the nonlinear Saint Venant equations, as shown by equation 4.13. The solution of the small perturbation of water depth can be taken to be a complex trigonometric function of wavenumber and angular frequency: =A 0 e i(kx!t) ; (4.27) whereA 0 denotes the initial maximum wave amplitude,k is the wavenumber and! = ! R +i! I is the complex angular frequency. Substituting the complex trigonometric function into equation 4.13, we obtain the following dispersion relation ! 2 +! + = 0 (4.28) in which =2u 0 k +i 2g sin u 0 (4.29) and = u 2 0 g cosh 0 k 2 i3g sink: (4.30) Above is the quadratic equation for !. This theory invented by Brook et al. (1999) species the maximum wave amplitude with prescribed initial perturbations, wavenum- ber and hydraulic characteristics at the normal condition. Hence the angular fre- quency must be quantitatively determined. In order to solve for !, we utilize the 66 general quadratic formula ! = + p 2 (4.31) where the expression for the discriminant is = 2 4 = 4g cosh 0 k 2 4(g sin) 2 u 2 0 +i4g sink: (4.32) The complex discriminant can be also written as =jjfcos [arg ()] +i sin [arg ()]g; (4.33) and so p =jj 1/2 cos 1 2 arg () +i sin 1 2 arg () : (4.34) Hence the solution for ! is expressed as ! =u 0 k + 1 2 jj 1/2 cos 1 2 arg () +i g sin u 0 + 1 2 jj 1/2 sin 1 2 arg () : (4.35) Since the angular frequency! is a complex, we can rewrite the small perturbation of water depth as =A 0 e ! I t e i(kx! R t) : (4.36) Above representation of manifests the wave amplitude changes with respect to time. Also the wave amplitude is dependent on the imaginary part of !. The wavenumber k as well as the real part of ! only in uence the phase shift of wave propagation. For future comparison between the theory and numerical results, we take the natural 67 logarithm of : ln () =! I t + ln (A 0 ) +i (kx! R t): (4.37) The maximum wave amplitude is then expressed as following by taking the real part: ln [max (jj)] =! I t + ln (A 0 ): (4.38) The natural logarithm of maximum magnitude of the small water depth perturbation is the transient linear theory. It is analytically deduced from the nonlinear Saint Venant equations through linearization. The theory denes the natural logarithm of maximum wave amplitude is a linear function of time with a slope of the imaginary part of the angular frequency. Next we perform a numerical experiment to compute the time-dependent amplitude of progressing roll waves, so that the comparison can be executed. 4.1.3 Numerical Experiment of Progressing Roll Waves The numerical experiment simulates how initial perturbations grow to progressing roll waves down an inclined channel resulting from the instability. It is implemented by solving equation 4.1 and 4.2 that the linear theory originates from and using the nite volume formulation and TVD Lax-Wendro scheme described in section 3.2. The equations are solved inside a one-dimensional computational domain. The undisturbed Froude number F 0 , drag coecient C d , and the ow rate q 0 at the normal condition will be prescribed such that the channel slope tan, water depthh 0 and ow velocityu 0 at the steady-uniform state can be computed due to the 68 balance between the gravitational force and bottom drag. As previously explained, the drag coecient is assumed to be constant because the purpose of this numerical experiment is to verify the accuracy of the numerical scheme. At the steady-uniform state, we have g sin =C d u 2 0 h 0 (4.39) that leads to tan =C d u 2 0 g cosh 0 =C d F 2 0 (4.40) which is used to calculate the channel slope. The normal water depth h 0 and ow velocity u 0 are then obtained by substituting u 0 = q 0 /h 0 into equation 4.39. The initial maximum wave amplitudeA 0 and the wavenumber k are also prescribed. The boundary conditions associated with the numerical experiment of progressing roll waves are dierent from those specied in Subsection 3.2.3. The boundary con- ditions should make initial perturbations and roll waves propagate innitely without any interference. Hence we choose that the 1st control cell is located at the down- stream of the Nth cell next to the out ow boundary while the Nth control cell is located at the upstream of the 1st cell next to the in ow boundary. Consequently the solution at the out ow boundary is imposed as the solution at the in ow boundary. In other words, the out ow boundary coincides with the in ow boundary. These boundary conditions are illustrated by Figure 4.1. The initial condition of water depth is composed of the normal water depth and a series of monochromatic sinusoidal perturbations: h (x; 0) =h 0 + (x; 0) (4.41) 69 Inflow boundary Outflow boundary 1 − Inflow boundary 1 2 x Outflow boundary Figure 4.1. Illustration of boundary conditions for the numerical experiment of pro- gressing roll waves. The out ow boundary coincides with the in ow boundary. The 1st control cell is located at the adjacent downstream of the Nth cell while the Nth control cell is located at the adjacent upstream of the 1st cell. and (x; 0) =A 0 sin (kx): (4.42) Brook et al. (1999) suggested the initial ow velocity must be corresponding to the prescribed initial water depth through the linearized equation of mass conservation for small perturbations. Therefore the ow velocity should be obtained by substituting =A 0 [sin (kx!t)i cos (kx!t)] (4.43) into equation 4.7 to deduce = A 0 h 0 ! k u 0 [sin (kx!t)i cos (kx!t)]: (4.44) 70 Then the initial perturbation of the ow velocity turns out to be (x; 0) = A 0 h 0 ! k u 0 [sin (kx)i cos (kx)]: (4.45) Since ! =! R +i! I is complex, we can write !/ku 0 in the polar form as ! k u 0 =r p e ip (4.46) where r p = ! k u 0 (4.47) and p = arg ! k u 0 : (4.48) After substitution from equation 4.46 and taking the real part, the initial perturbation of the ow velocity becomes (x; 0) = A 0 h 0 r p sin (kx + p ) (4.49) and the initial ow velocity is u (x; 0) =u 0 + (x; 0): (4.50) We should notice, from equation 4.49, the initial ow velocity is related to the angular frequency which is computed through equation 4.35. The solution proles of progressing roll waves at dierent instantaneous time levels 71 are presented by Figure 4.2. The horizontal axis indicates the channel locationx and the vertical axis represents the water depth h. The parameter values selected in this numerical experiment are: F 0 = 3:0,C d = 0:006,q 0 = 0:001 m 2 /s,A 0 = 0:005h 0 ,k = 10. The total channel distance is 1 m which is discretized into 1000 control cells. The Courant numberCr is chosen to be 0.75. The ow direction is from left to right. It is seen from the proles that the initial sinusoidal perturbations spontaneously magnify as they propagate and begin sharpening after they reach some certain amplitude. Eventually they transform to shock-like roll waves with extremely steep fronts. The growth of the wave amplitude reduces abruptly after the development of roll waves. We found developed roll waves travel without any change of wave amplitude and shape after t = 12:0 s approximately. 4.1.4 Comparison Results The water depth throughout the computational domain are computed at various time levels in numerical experiments specied above. To execute the comparison, we take the peak wave amplitude as a function of time from the numerical solutions and calculate the natural logarithm. The time-dependent linear theory is computed by equation 4.38 with a same set of parameters applied in numerical experiments. The comparison results between the theoretical and numerical solutions for the scenarios ofF 0 = 2:25,F 0 = 2:5, andF 0 = 3:0 are shown by Figure 4.3(a). Solutions with dierent undisturbed Froude number share a same set of parameters: C d = 0:006, q 0 = 0:001 m 2 /s, A 0 = 0:005h 0 , and k = 10. The vertical axis stands for the natural logarithm of maximum wave amplitude while the horizontal axis is time. 72 h(m) 0 0.2 0.4 0.6 0.8 1 0.0018 0.0022 0.0026 0.003 t=0.0 s h(m) 0 0.2 0.4 0.6 0.8 1 0.0018 0.0022 0.0026 0.003 t=2.0 s h(m) 0 0.2 0.4 0.6 0.8 1 0.0018 0.0022 0.0026 0.003 t=3.8 s x(m) h(m) 0 0.2 0.4 0.6 0.8 1 0.0018 0.0022 0.0026 0.003 t=5.0 s See caption on the next page. 73 h(m) 0 0.2 0.4 0.6 0.8 1 0.0018 0.0022 0.0026 0.003 t=5.8 s h(m) 0 0.2 0.4 0.6 0.8 1 0.0018 0.0022 0.0026 0.003 t=7.0 s h(m) 0 0.2 0.4 0.6 0.8 1 0.0018 0.0022 0.0026 0.003 t=8.4 s x(m) h(m) 0 0.2 0.4 0.6 0.8 1 0.0018 0.0022 0.0026 0.003 t = 12.0 s Figure 4.2. Presentation of numerical solutions of progressing sinusoidal pertur- bations and roll waves. The plots describe the water depth h against dier- ent channel locations within x 2 [0; 1] m, at instantaneous time levels: t = 0:0; 2:0; 3:8; 5:0; 5:8; 7:0; 8:4; 12:0 s. Associated computational parameters are F 0 = 3:0, C d = 0:006, q 0 = 0:001 m 2 /s,A 0 = 0:005h 0 , and k = 10. 74 t(s) ln [max (⏐η⏐ )] (ln (m)) 0 5 10 15 20 25 30 -11 -10 -9 -8 -7 =3.0 =2.5 =2.25 (a) t(s) ln [max (⏐η⏐ )] (ln (m)) 0 5 10 15 20 25 30 -14 -13 -12 -11 -10 =2.0 =1.8 =1.5 (b) Figure 4.3. Comparison results between transient linear theory and numerical solu- tions. The vertical axis describes the natural logarithm of maximum wave amplitude while the horizontal axis is time. The dash lines represent the linear theory and the solid lines represent numerical solutions. The solutions are for the undisturbed Froude numbers of F 0 = 2:25, F 0 = 2:5, F 0 = 3:0 in (a) and F 0 = 1:5, F 0 = 1:8, F 0 = 2:0 in (b). 75 The dash lines represent the linear theory and the solid lines represent nonlinear numerical solutions. It is seen that the theoretical and numerical solutions have good agreement when both of them develop exponentially. When the perturbation magnitude is small, the numerical results also present linear functions with a slope of the imaginary part of the angular frequency. Nevertheless nonlinear eects start dominating and the exponential growth of numerical solutions declines as the wave amplitude increases. Eventually the wave amplitude remains at a constant value. If we compare Figure 4.2 and 4.3(a), it is not dicult to discover the time level when waves begin to sharpen (e.g. t = 5:8 s) coincides with the time level when the linear behaviour of numerical solutions stops. In other words, the linear function dominates when waves are sinusoidal whereas nonlinear eecs play the more signicant role when waves become hydraulic jumps. Figure 4.3(b) shows comparison results between the linear theory and numerical solutions for the scenarios ofF 0 = 1:5,F 0 = 1:8, andF 0 = 2:0. Again the numerical results are capable of matching the linear slope closely. Since the stability criterion is satised, ! I is calculated to be zero whenF 0 = 2:0 and negative whenF 0 = 1:5 and F 0 = 1:8. Analogously, the computed wave amplitude nearly holds when F 0 = 2:0 and attenuates when F 0 = 1:5 and F 0 = 1:8. The nonlinear eects of numerical solutions do not present as the wave amplitude keeps small and no roll waves form. It is also interesting to mention the growth slope of wave amplitude increases when the undisturbed Froude number F 0 turns greater. For scenarios of F 0 > 2:0, the nonliear eects begin dominating and the growth rate begins decaying after a shorter time period when F 0 is greater. For example, the linear slope ceases at 76 t 6 s for F 0 = 3:0 whereas at t 18 s for F 0 = 2:25. Also the nal wave amplitude is larger with greater F 0 when a constant value is reached. For example, the nal wave amplitude hits ln [max (jj)]7:33 ln (m) for F 0 = 3:0 but only ln [max (jj)]9:12 ln (m) forF 0 = 2:25. 4.2 Comparison with Quasi-steady Nonlinear Theory In addition to the case of the time-dependent amplitude, the comparison of quasi- steady solutions after full evolution of progressing roll waves must be considered. Dressler (1949) constructed a mathematical solution of progressing roll waves by piecing together sections of a continuous wave prole through a series of discontinuous shocks. The wave prole was analytically deduced from the shallow water theory as equation 4.1 and 4.2 with prescribed channel slope, bottom drag coecient, wave speed and wavelength. This discovery was based on a progressing coordinate system that incorporates the eect of space and time. The specic derivation of Dressler's mathematical solution is given in Appendix A. As mentioned in Subsection 2.1.2, the method used by Dressler is unable to predict natural water ows, but it is commonly applied for testing the convergence property of a numerical scheme solving the Saint Venant equations. We also apply this quasi-steady nonlinear theory to test the present numerical scheme. The large roll waves have been shown by Figure 4.2 to emerge from initial per- turbations. As mentioned previously, the waves will evolve to a quasi-steady state in which the size and shape do not change. The numerical solution of the time level t = 19:1 s at which roll waves have fully converged is picked for comparison. On 77 the other hand, we use the same parameter values in the numerical experiment of progressing roll waves to calculate Dressler's mathematical solution. The wave speed is measured to be 0.6 m/s according to the computed time series from the numer- ical experiment. Thus several sections of the calculated Dressler's wave prole are connected through vertical jumps. 4.2.1 Comparison Results Figure 4.4 presents how the numerical solution of quasi-steady converged roll waves is compared with Dressler's theory. The agreement is satisfactorily good. Both the wave height and wavelength match very closely. Noting these two solutions are derived from a thoroughly identical mathematical model, we can see the numerical scheme is capable of accurately catching the wave prole. So far we have veried the accuracy of the present numerical scheme for com- puting roll waves in both transient and quasi-steady cases. Therefore the scheme is dependable to be used for solving the Saint Venant equations and simulating real roll-wave ows. 4.2.2 Sensitivity Analyses In the eld of computational uid dynamics, the eciency is always considered. The eciency of computation implies decreasing computational time while keeping satisfactory accuracy. Once a numerical scheme has been selected, the most usual ways of enhancing the computational eciency are to increase the cell size or time step. Hence before the real roll-wave ows are simulated, the most ecient compu- 78 x(m) h(m) 0 0.2 0.4 0.6 0.8 1 0.0018 0.0022 0.0026 0.003 Figure 4.4. Comparison of quasi-steady converged roll-wave solution between nu- merical and Dressler's analytical approaches. The vertical axis is the water depth while the horizontal axis is the channel location. The solid line stands for the present numerical solution while Dressler's theory is represented by the circle symbol. The numerical solution is picked at the instantaneous time level t = 19:1 s. tational approach must be determined. The present numerical scheme is an explicit one in computing the convective ux, so the CFL stability condition described by equation 3.48 must be satised. Since the maximum eigenvalue is variable and can be quite large, the Courant number rather than the time step must be a constant to ensure the stability condition. Hence the values of the cell size and Courant number are prescribed before each computation. We can then increase the cell size, x or Courant number, Cr so as to seek for a possible way which is more ecient. The procedure of the sensitivity analysis is then implemented by changing the cell size while keeping a constant Courant number or by changing the Courant number while keeping a constant cell size. 79 First the sensitivity to the cell size is tested. The Courant number is xed to be 0.75. Three numerical experiments are conducted using the same parameter values as what are used in Subsection 4.1.3 but dierent cell sizes: 0.001 m, 0.002 m, and 0.005 m. The results of a single, fully converged roll wave along with Dressler's theory are compared in Figure 4.5(a). The numerical solutions are all picked at the instantaneous time level t = 19:1 s. It is obvious that the solution associated with the cell size of 0.001 m has the best convergence property. When the cell size is increased to 0.002 m or 0.005 m, the wave shape inclines upstream and deviates from the original position. Also the solutions with coarse cell sizes show deformed shocks that are o from the vertical direction. Therefore the present numerical scheme is sensitive to the cell size, and the ne size of 0.001 m should be applied in roll-wave simulations despite its more computational time. A similar analysis is carried out for the sensitivity to the Courant number. The cell size is xed to be 0.001 m and three dierent values of Courant number are tested: 0.5, 0.75, and 0.95. The results are presented by Figure 4.5(b). The numerical solutions are also picked at the time levelt = 19:1 s. The numerical scheme seems not sensitive to the Courant number as the solutions associated with dierent Courant numbers all have satisfactory accuracy. The solutions with the Courant number of 0.5 and 0.75 are almost identical. However the solution with the Courant number of 0.95 has a sharp tip at the crest resulting in a slightly overestimated wave height. The computed crest depth will be compared with experimental data in the next section. Thus the Courant number of 0.95 is not considered. The Courant number of 0.75 will be used in roll-wave simulations because it has less computational time than the 80 x(m) h(m) 0.1 0.2 0.3 0.0018 0.0022 0.0026 0.003 Dressler’s theory Δx = 0.001 m Δx = 0.002 m Δx = 0.005 m (a) x(m) h(m) 0.1 0.2 0.3 0.0018 0.0022 0.0026 0.003 Dressler’s theory Cr =0.5 Cr =0.75 Cr =0.95 (b) Figure 4.5. Sensitivity analyses of the present scheme applied in the numerical experi- ment of progressing roll waves. The vertical axis is the water depth and the horizontal axis is the channel location. The numerical solutions are picked at the instantaneous time level t = 19:1 s. (a) shows the sensitivity to the cell size with a xed Courant number of 0.75. (b) shows the sensitivity to the Courant number with a xed cell size of 0.001 m. 81 Courant number of 0.5. 4.3 Comparison with Brock's Experiment In this section we adopt the computational model USC SWAT constructed in Chapter 3 to simulate roll-wave ows observed in the real life. The free surface ows down an inclined, wide, rectangular channel are simulated with prescribed normal ow rate per unit width q 0 , channel slope tan, and absolute bottom roughness k s . The water depth h 0 and ow velocity u 0 at the steady-uniform condition are calculated based on the discharge equation q 0 =h 0 u 0 ; (4.51) the Ch ezy formula u 0 =C c p R h tan (4.52) where C c denotes the Ch ezy coecient and R h is the hydraulic radius, and the Colebrook-White equation 1 p f =K 1 log k s K 2 R h + K 3 4< p f (4.53) where f is the Darcy-Weisbach friction factor. < is the Reynolds number evaluated by < = R h u (4.54) 82 in which stands for the kinematic viscosity. K 1 ,K 2 , andK 3 are constant coecients applied in the equation. Using the Darcy-Weisbach friction factorf, we can evaluate the Ch ezy coecient by C c = s 8g cos f (4.55) (Sturm 2001). Since the present study concentrates on ows in a wide, rectangu- lar channel in which the channel width is much greater than the water depth, the hydraulic radius is approximately equal to the water depth. Therefore the Ch ezy formula is rewritten as u 0 = s 8g cos f h 0 tan: (4.56) Some values of the constant coecients are listed by Yen (2002). In this study we utilize the set referred by Henderson (1966): K 1 = 2:0, K 2 = 12:0, and K 3 = 2:5. Thus the Colebrook-White equation is rewritten as 1 p f =2 log k s 12h + 2:5 4< p f : (4.57) The drag coecient, C d , within the term describing the bottom resistance needs to be calculated on each control cell at each time level. C d is calculated byC d =f/8, as mentioned in Section 3.1. The Darcy-Weisbach friction factorf is computed based on equation 4.57. This equation states that f is dependent on the relative roughness and Reynolds number, which alter among each control cell. The equation is implicit onf. Hence the Newton-Raphson's iterative method is applied to obtain the friction 83 factor. A function is established as NR (f) = 1 2 p f + log k s 12h + 2:5 4< p f : (4.58) The rst derivative with respect to f of the above function is taken to be NR 0 (f) = 1 4 f 3/2 + log (e) ks 12h + 2:5 4 Re p f 2:5 4< 1 2 f 3/2 (4.59) wheree denotes the mathematical constant with a value of 2:71828 . The iterative procedure is executed as: 1. Guess an initial value of f 0 ; 2. Compute NR (f n ) and NR 0 (f n ), n = 0; 1; 2; ; 3. Compute f n+1 =f n NR (f n ) NR 0 (f n ) ; 4. Ifjf n+1 f n j6jf n+1 j, then f n+1 is the nal answer; Ifjf n+1 f n j>jf n+1 j, then letf n =f n+1 and repeat step 2, 3, and 4 until the solution converges. is the small tolerance. The above algorithm for the computation of the drag coecient is incorporated into USC SWAT as a subroutine named \DRAGCOEF". The boundary conditions have been specied in Subsection 3.2.3. Normal ow conditions can be set to be the initial conditions. We can also set no- ow condition initially. The cell size in the computation is selected to be 0.001 m and the Courant number is selected to be 0.75, as specied in Subsection 4.2.2. The Compaq Visual 84 Fortran Standard Edition 6.6.0 is adopted to compile the computer code. The com- putations run in a Hewlett-Packard desktop computer with the Windows Vista 64-bit operating system. The processor is Intel Core 2 Quad CPU Q8300 at the speed of 2.5GHz and the memory is 6.0 GB. 4.3.1 Time Series of Roll Waves The complete outcome of the roll-wave simulation using USC SWAT is time series of ow variables at all channel locations. Figure 4.6 exhibits nite sections of computed time series of the water depth at selected channel locations. The associated parameters used in the simulation are q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. The undisturbed Froude number is calculated as F 0 2:60. The water depth is recorded at eight dierent channel locations from upstream to downstream: 10, 20, 30, 40, 50, 60, 80, and 100 m away from the channel entrance (the in ow boundary of the computational domain). The interesting time period is t2 [90; 110] s for all locations. The time series of water depth represent the ow conditions at each channel location. In general, the perturbations are small near the channel entrance whereas they become larger at downstream locations. This corroborates small perturbations magnify spontaneously as they propagate downstream because of the instability. The wave behaviors are quite uniform at one specic channel location. For instance, waves are always small and smooth at upstream locations such as x6 30 m whereas all of them become large and shock-like at downstream locations, particularly at x > 80 m. In other words, no developed roll waves can be found at upstream locations and 85 h(m) 90 94 98 102 106 110 0.01 0.02 0.03 x = 10 m h(m) 90 94 98 102 106 110 0.01 0.02 0.03 x = 20 m h(m) 90 94 98 102 106 110 0.01 0.02 0.03 x = 30 m t(s) h(m) 90 94 98 102 106 110 0.01 0.02 0.03 x = 40 m See caption on the next page. 86 h(m) 90 94 98 102 106 110 0.01 0.02 0.03 x = 50 m h(m) 90 94 98 102 106 110 0.01 0.02 0.03 x = 60 m h(m) 90 94 98 102 106 110 0.01 0.02 0.03 x = 80 m t(s) h(m) 90 94 98 102 106 110 0.01 0.02 0.03 x = 100 m Figure 4.6. Presentation of time series of the water depth in the simulation of a real roll-wave ow. The vertical axis is the water depth while the horizontal axis is time. The time series are recorded within t2 [90; 110] s at channel locations of 10, 20, 30, 40, 50, 60, 80, and 100 m. The solutions are for the parameters of q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. 87 no small-amplitude perturbations can be found at very downstream locations. The small perturbations begin to sharpen and have steeply sloping fronts at some certain channel section. In this case, shock-like roll waves start to form at x 40 m. The location at which a small perturbation become a roll wave is not stationary, as the perturbations with larger amplitude develop to roll waves more upstream than those with smaller amplitude. It is found that no new shock-like waves are formed behindx 60 m in that all perturbations have completed transformations. The well- evolved roll waves possess steep fronts while retaining smooth backs. Theoretically, the formation of shock-like roll waves in the simulation results is attributed to the nonlinearity of the equations of motion. From upstream to downstream locations, the wave period has a tendency to in- crease in that the wavenumber over a certain time period decreases. This tendency is caused by wave coalescences that appear right after the formation of shock-like waves. The occurrences of wave coalescences decrease the wavenumber or enlarge the wavelength. Hence, within a certain time period, the wavenumber presented at downstream locations is less than that at upstream locations. Also the combination of multiple waves results in a novel wave with increased size. Thus the wave coalescence is the other reason besides the instability for the increase of the wave amplitude. The present simulation discovered several types of wave-wave interactions including the wave coalescence and they will be specically discussed in Chapter 5. The simulated naturally growing roll waves are nonperiodic at all channel lo- cations. The wave height and wave period alter among dierent waves at a given location. The nonperiodicity of the natural roll-wave ow had been demonstrated 88 by Brock's experiment. He showed that the periodic ow could only be achieved by periodic motions of the wave maker at the channel entrance. The time series of water depth at channel locations of x = 10, 20, 30, 40, 60, and 80 m are re-displayed by Figure 4.7 in a contour format for a better comparison. The ow velocity was not analyzed in Brock's experiment. Using the present computations, the ow velocity prole of roll waves is compared with the water depth and the results are presented in Figure 4.8. The results include time series of the water depth, ow velocity, and momentum ux per unit width per unit density over a particular time period at a given channel location. The computational parameters are set to be what was applied for the time series of water depth presented in Figure 4.6. The interesting time period is t2 [80; 100] s and they are recorded at the location x = 80 m. It is seen that the behavior of the ow velocity prole is entirely consistent with that of the concurrent water depth. When the water depth is at the crests, the ow velocity reaches the peaks; when the water depth drops to the troughs, the ow velocity also decreases to the low points. Even the transitions on wave backs and fronts have a consistent pattern. When the water depth is gradually varied at the backs, the change of ow velocity is smooth; when the water depth runs into a shock, the ow velocity also presents a jump that is nearly discontinuous. Therefore the momentum ux prole, calculated as hu 2 , is analogous to the fashion of the water depth and ow velocity as well. Thomas (1937) observed that, in sections of the continuous wave prole, the ow were supercritical near the wave troughs and subcritical near the wave crests. This 89 t (s) 90 95 100 105 110 x(m) 10 20 30 h (m) 0.01 0.02 t (s) 90 95 100 105 110 x(m) 40 60 80 h (m) 0.01 0.02 Figure 4.7. Contour of time series of the water depth in the simulation of a real roll- wave ow. The vertical axis is the water depth while two horizontal axes represent time and channel locations. The time series are recorded within t2 [90; 110] s at channel locations of 10, 20, 30, 40, 60, and 80 m. The solutions are for the parameters of q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. 90 h(m) 80 85 90 95 100 0.01 0.02 0.03 u (m/s) 80 85 90 95 100 0.6 1 1.4 t(s) hu 2 (m 3 /s 2 ) 80 85 90 95 100 0.01 0.03 Figure 4.8. Presentation of comparison between the time series of ow velocity and water depth in the simulation of a real roll-wave ow. The results show the water depth, ow velocity, and momentum ux per unit width per unit density within t2 [80; 100] s at x = 80 m. The solutions are for the parameters of q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. 91 is questionable since all of the roll-wave ows are unstable and highly supercritical. Dressler (1949) mentioned the ow velocity was everywhere less than the wave speed c w . Hence, in the progressing coordinate, the ow velocity computed as u c w is negative denite toward the opposite direction of travelling waves. As such the observation by Thomas regarding the transition between supercritical and subcritical ows is based on the progressing coordinate, and the ow velocity is that with respect to the wave speed. According to the simulation results in the present study,juc w j is greater near the wave troughs than the wave crests. Meanwhile the water depth is smaller near the wave troughs than the wave crests. Thus Thomas declared his observation. 4.3.2 Time-averaged Wave Crest and Trough Depths One of the major objectives of the present study is to computationally analyze the development of the time-averaged wave crest and trough depths versus the channel location, which were researched experimentally by Brock (1969). The experimental procedure and principal discoveries are introduced in Section 2.2. In order to preserve the consistency for the comparisons, the time-averaged wave crest and trough depths are computed through analogous approaches to Brock's experiment. The time series of water depth at various channel locations introduced in Subsec- tion 4.3.1 are utilized to obtain the time-averaged wave crest and trough depths. The computational procedure is as following: 1. At an interesting channel location, the time series of water depth characterizing either small perturbations or large roll waves is simulated using model USC 92 SWAT. The time period should be suciently long such that more than 200 waves are included, because Brock's experiment also picks 200 waves as a sample for the calculation of the arithmetic average value; 2. Find the crest and trough depths for each single wave out of 200 waves in the time series, and then obtain the time-averaged wave crest and trough depths at this location by calculating the arithmetic mean value; 3. Repeat steps 1 and 2 at various interesting channel locations to form devel- opment curves of the time-averaged wave crest and trough depths versus the channel location. It should be emphasized that the report of Brock's experiment does not provide the absolute roughness height of the steel laboratory channel, which is an essential parameter for calculating the drag coecient in the model USC SWAT. Hence, be- fore the computation is executed, we need to nd out the roughness height of the laboratory channel. For each run with dierent channel slope or discharge in the experiment, the angle of inclination, normal water depth, and normal ow velocity are measured and provided by the report. Therefore the friction factor f at the steady-uniform condition can be calculated by equation 4.56. Then the calculated friction factor is substituted into equation 4.57 to obtain the roughness height of the laboratory channel. 4.3.3 Sensitivity Analyses In Subsection 4.2.2, the TVD Lax-Wendro numerical scheme has been corrob- orated to be sensitive to the cell size. The coarse cell sizes lead to deformed shapes 93 of progressing roll waves. The sensitivity of the model USC SWAT to the cell size in the simulation of real roll-wave ows is analyzed herein. The results are shown in Figure 4.9. The time-averaged wave crest and trough depths are plotted as func- tions of the channel location. All quantities are normalized by the water depth at the steady-uniform condition. The wave crest and trough depths are obtained with a spatial interval of 10 m and the most upstream interesting location is 10 m. The total channel length is set to be 151 m. The maximum time level of the time series reaches 1300 s so that at least 200 waves are included. The hydraulic parameters are = 0:05013,h 0 = 0:00523 m,u 0 = 0:783 m/s, and t = 0:26 m 2 /s. Four dierent cell sizes are tested: x = 0:001, 0.002, 0.005, and 0.01 m. The simulated results with coarse cell sizes tend to overestimate the amplication rate of wave height comparing with the result of x = 0:001 m. When the cell size is increased, the spatial growth rate of wave crest depth increases. There is a similar manner for the spatial decaying rate of wave trough depth. Nevertheless, at the very downstream stations where the amplication rate of wave height almost reduces to zero, the amplitude of fully evolved roll waves is quite close for dierent cell sizes. The sensitivity analyses in Subsection 4.2.2 indicate the coarse cell sizes result in deformed wave shapes, but the wave amplitude is not aected. Yet in the simulation of real roll-wave ows coarse cell sizes induce overestimated amplication rate. For a computational domain with a particular length, the ner cell size leads to more control cells, and thus leads to more computational time. The computational times to generate the time series of water depth, which are required for the time- averaged wave crest and trough depths presented in Figure 4.9, are recorded for all 94 x/h 0 ⎯ h c /h 0 0 5000 10000 15000 20000 25000 30000 1 1.3 1.6 1.9 2.2 2.5 2.8 Δx = 0.001 m Δx = 0.002 m Δx = 0.005 m Δx =0.01 m (a) × × × × × × × × × × × × × × × x/h 0 ⎯ h t /h 0 0 5000 10000 15000 20000 25000 30000 0.5 0.6 0.7 0.8 0.9 1 Δx = 0.001 m Δx = 0.002 m Δx = 0.005 m Δx =0.01 m × (b) Figure 4.9. Sensitivity of the computational model USC SWAT to the cell size in the simulation of real roll-wave ows. The horizontal axis is the normalized channel location while the vertical axis is the normalized time-averaged wave crest depth in (a) and normalized time-averaged wave trough depth in (b). The associated hydraulic parameters are = 0:05013, h 0 = 0:00523 m, u 0 = 0:783 m/s, and t = 0:26 m 2 /s. 95 cases of the dierent cell sizes. They are shown in Table 4.1. When the cell size is decreased as retaining other computational parameters, such as the Courant number, total channel length, and maximum time level of the time series, the computational time dramatically increases. For instance, the computational time associated with x = 0:01 m is 0 days 3 hours and 0.20 minutes. The computational time associated with x = 0:001 m reaches 9 days 10 hours and 2.77 minutes, which is approximately seventy-ve times longer than the case of x = 0:01 m. The eciency of any computational model comprises two aspects: speed and accuracy. Enhancing the speed with losing the accuracy is never acceptable. The comparison between the theoretical and numerical solutions has manifested the coarse cell sizes result in bad convergence. The computed time-averaged wave crest and trough depths are also lack of agreement among dierent cell sizes. Hence, although the coarse cell sizes greatly cut down the computational cost, we still adopt the nest cell size x = 0:001 m because of its better accuracy. As discussed in Section 3.1, we use the common approach for dening the normal stress of uids of incompressible ow to model the momentum diusion. However the kinematic viscosity must be replaced by a turbulent viscosity t because the roll- wave ows specied in this study are turbulent. In the eld of turbulence modeling, particularly Reynolds-averaged Navier-Stokes (RANS) equations modeling, the tur- bulent viscosity (or eddy viscosity) is commonly used to close the Reynolds stresses. Turbulence models have been devised to compute the turbulent viscosity, for example mixing length model, k" model, and Reynolds stress model, etc (Fletcher 1991; Wilcox 2006; Versteeg & Malalasekera 2007). Of those turbulence models, the mixing 96 Cell size (m) Courant number Channel length (m) 0.01 0.75 151 0.005 0.75 151 0.005 0.75 151 0.002 0.75 151 0.001 0.75 151 Channel length (m) Maximum time level (s) Computational time 151 1300 0 days + 3 hours + 0.20 minutes 151 1300 0 days + 10 hours 151 1300 + 10 hours + 55.46 minutes 151 1300 2 days + 15 hours + 5.54 minutes 151 1300 9 days + 10 hours + 2.77 minutes Table 4.1. Presentation of comparison of computational time among dierent cell sizes in the simulation of real roll-wave ows. The computational times to generate the time series are indicated for dierent cell sizes of x = 0:01, 0.005, 0.002, and 0.001 m, but identical Courant number 0.75, channel length 151 m, and maximum time level 1300 s. length andk" models are by far the most widely used, since they are well validated by controlled experiments. Nevertheless they do not provide satisfactory results in some circumstances of computing complex turbulent ows. It is demonstrated that none of those models is suitable for all kinds of ow conditions. The turbulent viscos- ity applied in the Saint Venant equations can not be computed by those turbulence models associated with RANS modeling due to their distinguishing mathematical contexts. Also they have not been approved to be used in the simulations of real roll-wave ows. Moreover no empirical information regarding the computation of turbulent viscosity in the Saint Venant equations can be found in literatures. Thus there is no proposed model that has been validated by experiments for quantifying 97 the turbulent viscosity in the present study, and the sensitivity analysis must be implemented. Figure 4.10 shows the results of the sensitivity of the computed time-averaged wave crest and trough depths to the turbulent viscosity. The plots are produced in an identical manner as the sensitivity analysis about the cell size. The computational and hydraulic parameters are also what are used in the sensitivity analysis about the cell size except that ve dierent values of the turbulent viscosity are tested: t = 2:6e-01, 1:0e-01, 1:0e-02, 1:0e-03, and 1:0e-04 m 2 /s. In general, greater values of the turbulent viscosity result in smaller amplitude of wave crest and trough. At any channel location, the wave crest depth decreases when the turbulent viscosity is increased. With greater values of turbulent viscosity, the amplication of the wave crest depth and transformation from small perturbations to shock-like roll waves due to the instability occur at more downstream stations. In other words, with greater values of turbulent viscosity, the channel section where small perturbations dominate is longer starting from the in ow boundary. For instance, with t = 2:6e-01 m 2 /s, the amplication due to the instability occurs at x/h 0 2 (10000; 22000) roughly while with t = 1:0e-01 m 2 /s the amplication due to the instability occurs at x/h 0 2 (5000; 12000) roughly. The amplication with t = 1:0e-02, 1:0e-03, and 1:0e-04 m 2 /s occurs at even more upstream locations. Also the spatial growth rate due to the instability is smaller as the turbulent viscosity is increased. Consequently for the cases of t = 1:0e-03 and 1:0e-04 m 2 /s, the waves are already not small perturbations and the wave crest depths grow to high levels even at the front of the most upstream interesting location of x = 10 m. Therefore these 98 x/h 0 ⎯ h c /h 0 0 5000 10000 15000 20000 25000 30000 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 ν t = 2.6e-01 m 2 /s ν t = 1.0e-01 m 2 /s ν t = 1.0e-02 m 2 /s ν t = 1.0e-03 m 2 /s ν t = 1.0e-04 m 2 /s (a) × × × × × × × × × × × × × × × x/h 0 ⎯ h t /h 0 0 5000 10000 15000 20000 25000 30000 0.4 0.5 0.6 0.7 0.8 0.9 1 ν t = 2.6e-01 m 2 /s ν t = 1.0e-01 m 2 /s ν t = 1.0e-02 m 2 /s ν t = 1.0e-03 m 2 /s ν t = 1.0e-04 m 2 /s × (b) Figure 4.10. Sensitivity of the time-averaged wave crest and trough depths to the turbulent viscosity in the simulation of real roll-wave ows. The horizontal axis is the normalized channel location while the vertical axis is the normalized time-averaged wave crest depth in (a) and normalized time-averaged wave trough depth in (b). The associated hydraulic parameters are = 0:05013, h 0 = 0:00523 m, u 0 = 0:783 m/s. 99 distinctive features regarding the amplication due to the instability lead to the gaps of the wave crest depths among dierent values of turbulent viscosity. Behind the amplication due to the instability, the spatial growth rate of wave crest depth due to wave coalescences is alike for dierent values of turbulent viscosity. The spatial decaying rate of wave trough depth has an analogous fashion to the growth rate of wave crest depth, except the downstream wave trough depths with t = 1:0e-02, 1:0e-03, and 1:0e-04 m 2 /s become larger than those with smaller values of turbulent viscosity. Since the variation of the value of turbulent viscosity has a signicant impact on the results of time-averaged wave crest and trough depths, we should numerically summarize how much dierence is caused by decreasing the value of turbulent vis- cosity from t = 2:6e-01 m 2 /s. The percentages of increase of wave crest depths with t = 1:0e-01, 1:0e-02, 1:0e-03, and 1:0e-04 m 2 /s from the wave crest depth with t = 2:6e-01 m 2 /s at each interesting channel location are tabulated in Table 4.2. The percentages of decrease of wave trough depths with t = 1:0e-01, 1:0e-02, 1:0e-03, and 1:0e-04 m 2 /s from the wave trough depth with t = 2:6e-01 m 2 /s are tabulated in Table 4.3. For the wave crest depths, the percentage of increase is greater when the value of turbulent viscosity is smaller at any channel location, like what is shown by Figure 4.10(a). Hence Table 4.2 also corroborates the wave crest depths are increased as the value of turbulent viscosity is decreased. For each value of turbulent viscosity, the greatest percentage occurs at x = 50 or 60 m where waves with smaller value of turbulent viscosity have already magnied due to the instability and transformed to shock-like roll waves, whereas waves with t = 2:6e-01 m 2 /s are still small per- 100 ν t (m 2 /s) x (m) 1.0e−01 1.0e 10 0.34 20 1.66 30 7.52 40 37.58 50 95.75 60 113.99 70 100.09 70 100.09 80 66.59 90 46.17 100 30.66 110 21.39 120 16.55 130 15.10 140 15.95 150 17.51 1.0e−02 1.0e−03 1.0e−04 12.15 76.20 97.36 102.40 144.71 196.45 142.28 172.12 252.39 166.52 199.90 276.70 190.95 207.05 280.22 182.74 206.44 266.78 160.94 180.37 239.04 160.94 180.37 239.04 118.25 134.45 182.44 94.55 108.09 153.91 66.69 79.02 122.98 53.57 63.77 107.13 47.04 55.80 99.02 44.60 52.88 96.36 42.65 50.73 94.95 42.41 51.06 97.12 Table 4.2. Presentation of percentages of increase of wave crest depths with t = 1:0e-01, 1:0e-02, 1:0e-03, and 1:0e-04 m 2 /s from the wave crest depth with t = 2:6e-01 m 2 /s. The total channel length is 151 m. The rst interesting channel location is 10 m and the interval is 10 m. turbations. When t = 1:0e-04 m 2 /s the wave crest depth can increase by nearly three times at x = 50 m. For the wave trough depths, the percentage of decrease is greater when the value of turbulent viscosity is smaller at x < 50 m. At x > 50 m, the wave trough depth with t = 1:0e-01 m 2 /s has the greatest percentage of de- crease. The percentages of decrease become smaller as the channel location increases for t = 1:0e-02, 1:0e-03, and 1:0e-04 m 2 /s. The reason for these is that the wave 101 ν t (m 2 /s) x (m) 1.0e−01 10 0.23 20 1.21 30 6.40 40 18.68 50 44.09 60 48.08 70 43.80 70 43.80 80 38.74 90 31.62 100 25.96 110 20.55 120 17.59 130 16.30 140 17.93 150 17.91 1.0e−02 1.0e−03 1.0e−04 9.85 41.19 47.22 47.86 47.91 55.73 51.28 48.28 56.75 51.39 46.14 53.97 48.47 45.28 52.60 42.95 41.92 48.89 35.62 35.87 43.53 35.62 35.87 43.53 25.39 26.55 34.60 15.47 16.54 23.53 7.00 10.24 14.93 0.37 2.72 7.29 −1.40 −0.32 5.30 −4.52 −3.88 2.01 −1.84 −2.02 2.76 −1.46 −1.85 2.73 Table 4.3. Presentation of percentages of decrease of wave trough depths with t = 1:0e-01, 1:0e-02, 1:0e-03, and 1:0e-04 m 2 /s from the wave trough depth with t = 2:6e-01 m 2 /s. The total channel length is 151 m. The rst interesting channel location is 10 m and the interval is 10 m. trough depths for t = 1:0e-02, 1:0e-03, and 1:0e-04 m 2 /s start to increase at x = 50 m. The sensitivity analysis about the turbulent viscosity indicates when the value of turbulent viscosity is increased, the spatial growth rate of wave amplitude is decreased. The decease of the growth rate of wave amplitude implies existences of wave energy dissipation. The wave energy dissipation can be observed in real roll-wave ows 102 that show turbulent eddies at the vicinity of wave fronts. The energy dissipation is mathematically modelled by the application of the momentum diusion which adopts the turbulent viscosity in the present equations of motion. 4.3.4 Correction Length The value of turbulent viscosity has to be selected so that the spatial growth and decaying rates of computed time-averaged wave crest and trough depths are able to match Brock's experimental data. However, when the turbulent viscosity is opti- mized for the comparison, the length of the upstream section starting from the in ow boundary in which small perturbations dominate is usually too large. Consequently the locations where waves magnify due to the instability are excessively downstream and the simulation results lose agreements with the experimental data. Hence we need to neglect a channel section at the very upstream locations and dene a new in ow boundary more downstream. The demonstration of the corrected in ow boundary is presented by Figure 4.11. The associated hydraulic parameters are = 0:05013, h 0 = 0:00523 m,u 0 = 0:783 m/s, and t = 0:26 m 2 /s. The total channel length is 180 m. The prole solution is picked at the instantaneous time level t = 230 s. In this case, the corrected in ow boundary is 50 m more downstream than the original in ow boundary. So after correction all channel locations are 50 m less than the original channel location. The demonstration of time-averaged wave crest and trough depths before and after correcting upstream are shown in Figure 4.12. The normalized wave crest and trough depths are correlated with the normalized channel location. The total channel length 103 x(m) h(m) 0 20 40 60 80 100 120 140 160 180 0.004 0.008 0.012 0.016 x = 50 m: corrected inflow boundary Figure 4.11. Demonstration of corrected in ow boundary in the simulation of a real roll-wave ow. The vertical axis is the water depth and the horizontal axis is the channel location. The prole solution is picked at the instantaneous time level t = 230 s. The corrected in ow boundary is indicated as x = 50 m. The hydraulic parameters are = 0:05013, h 0 = 0:00523 m, u 0 = 0:783 m/s, and t = 0:26 m 2 /s. is 151 m before correction. The rst interesting location is 10 m before correction and the interval is 10 m. The maximum time level of the time series of water depth reaches 1300 s in order to include 200 waves. The hydraulic parameters are those used in Figure 4.11. Compared with crest and trough depths before correction, the total channel length becomes 101 m and the rst interesting location changes to 0 m after correction. The total number of interesting locations reduces to 11, because other four upstream channel locations fall to negative. After correction the excessive channel section with small perturbations is eliminated and the amplication of wave amplitude due to the instability occurs at more upstream locations. As such, with appropriate selections of the value of turbulent viscosity and the correction length, the comparisons between the simulation results and experimental data can be conducted. 104 × × × × × × × × × × × × × × × × × × × × × × × × × × x/h 0 ⎯ h/h 0 0 6000 12000 18000 24000 30000 0.5 1 1.5 2 2.5 3 Crest depth Trough depth Corrected crest Corrected trough × × Before correction After correction Figure 4.12. Demonstration of time-averaged wave crest and trough depths before and after correcting upstream in the simulation of a real roll-wave ow. The vertical axis is the normalized time-averaged wave crest and trough depths and the horizontal axis is the normalized channel location. The crest and trough depths before and after correction are indicated. The hydraulic parameters are = 0:05013,h 0 = 0:00523 m, u 0 = 0:783 m/s, and t = 0:26 m 2 /s. 4.3.5 Comparison Results In Brock's experiment, three dierent channel slopes are used and two or three dierent ow rates are tested for each channel slope. Thus there are totally eight dierent scenarios in the experiment. One of major outcomes from the experiment is the spatial development of time-averaged wave crest and trough depths. The ap- proaches on how to obtain this experimental outcome are specied in Section 2.2. Four scenarios out of eight in the experiment are selected to perform the comparison. 105 Scenario Angle of inclination Discharge (m 2 /s) Normal depth (m) Normal velocity (m/s) Froude number Turbulent viscosity (m 2 /s) Correction length (m) (1) 0.05013 0.004098 0.00523 0.783 3.46 0.26 50 (2) 0.05013 0.008264 0.00798 1.04 3.71 0.26 46 (3) 0.08439 0.01108 0.00798 1.39 4.97 0.36 38 (4) 0.1201 0.006827 0.00533 1.28 5.60 0.52 50 Table 4.4. Presentation of hydraulic parameters applied in Brock's experiment and the present simulation. Four scenarios are analyzed in the present study, including three dierent channel slopes and four dierent discharges. The selections of the value of turbulent viscosity and correction length are also indicated. The hydraulic parameters for the selected four scenarios are summarized in Table 4.4. Due to the designed angles of inclination and discharges, the resulting Froude num- bers are fairly large and the ows are highly supercritical. The experimental data of time-averaged wave crest and trough depths versus channel locations associated with these four scenarios can be acquired from Brock's Ph.D. dissertation report (1967). For those scenarios which are not selected, the experimental data of either wave crest or trough depth are incomplete. The wave crest or trough depth was not recorded at a number of channel locations in those scenarios. Hence they are inadequate to be selected for the purpose of comparisons. The hydraulic parameters in the four scenarios are applied in the present simula- tion of real roll-wave ows, which achieves computed time-averaged wave crest and trough depths at various channel locations. The selections of the value of turbulent viscosity and correction length are shown in Table 4.4. The computational param- 106 Scenario Channel length (m) Maximum time level (s) (1) 151 1300 (2) 191 1300 (2) 191 1300 (3) 191 800 (4) 151 800 Courant number Cell size (m) Computational time 0.75 0.001 9 days + 10 hours + 2.77 minutes 0.75 0.001 12 days + 6 hours 0.75 0.001 + 6 hours + 39.52 minutes 0.75 0.001 18 days + 8 hours + 55.10 minutes 0.75 0.001 13 days + 22 hours + 12.97 minutes Table 4.5. Presentation of computational parameters applied in the present simulation for comparison with Brock's experiment. Computational parameters, including the computational time, for all of the four scenarios are indicated. eters, including the computational time to generate sucient time series of water depth, are shown in Table 4.5 for all of the four scenarios. The computation for each scenario takes extremely long time for the production of more than 200 waves. The interval between adjacent interesting channel locations is 10 m. The overviews of comparison results are presented by Figure 4.13 Figure 4.16 for the Scenario (1) Scenario (4). The vertical axes are the normalized time-averaged wave crest and trough depths and the horizontal axes are the normalized channel lo- cation. The normalization is executed through a division by the normal water depth. It can be seen that the simulated wave crest and trough depths are capable of fol- lowing the development trend of those obtained by Brock' experiment for all of the four scenarios. Both the wave crest and trough depths match closely with Brock's experimental data at various channel locations. For the wave crest depth, the wave 107 × × × × × × × × × × × x/h 0 ⎯ h/h 0 0 2500 5000 7500 10000 12500 15000 17500 20000 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 Crest depth by present simulation Crest depth by Brock’s experiment Trough depth by present simulation Trough depth by Brock’s experiment × Figure 4.13. Overview of comparison with Brock's experiment: Scenario (1). The vertical axis is the normalized time-averaged wave crest and trough depths and the horizontal axis is the normalized channel location. 108 × × × × × × × × × × × × × × × x/h 0 ⎯ h/h 0 0 2500 5000 7500 10000 12500 15000 17500 20000 0.4 0.7 1 1.3 1.6 1.9 2.2 2.5 2.8 3.1 Crest depth by present simulation Crest depth by Brock’s experiment Trough depth by present simulation Trough depth by Brock’s experiment × Figure 4.14. Overview of comparison with Brock's experiment: Scenario (2). The vertical axis is the normalized time-averaged wave crest and trough depths and the horizontal axis is the normalized channel location. 109 × × × × × × × × × × × × × × × × x/h 0 ⎯ h/h 0 0 2500 5000 7500 10000 12500 15000 17500 20000 0 0.5 1 1.5 2 2.5 3 3.5 4 Crest depth by present simulation Crest depth by Brock’s experiment Trough depth by present simulation Trough depth by Brock’s experiment × Figure 4.15. Overview of comparison with Brock's experiment: Scenario (3). The vertical axis is the normalized time-averaged wave crest and trough depths and the horizontal axis is the normalized channel location. 110 × × × × × × × × × × × x/h 0 ⎯ h/h 0 0 2500 5000 7500 10000 12500 15000 17500 20000 0 0.5 1 1.5 2 2.5 3 3.5 4 Crest depth by present simulation Crest depth by Brock’s experiment Trough depth by present simulation Trough depth by Brock’s experiment × Figure 4.16. Overview of comparison with Brock's experiment: Scenario (4). The vertical axis is the normalized time-averaged wave crest and trough depths and the horizontal axis is the normalized channel location. 111 amplitude as well as the growth rate of wave amplitude are quite small near the chan- nel inlet. The growth rate magnies behind a certain channel location or a certain level of crest depth. The large amplication rate is caused by the open channel ow instability. The wave crest depth becomes large as small perturbations transform to large shock-like roll waves. Since the space of the laboratory where the experiment was carried out was limited, the experimental data for ows after the amplication due to the instability were not recorded. Therefore the further downstream ow char- acteristics were not studied by Brock. Nevertheless the present simulation manifests the wave crest depth continues to grow behind the end of experimental data, but the growth rate decays because of the balance between the instability and nonlinear eect. The continuous but smaller growth of crest depth is induced by wave coales- cences. The growth of wave crest depth seems not to stop even at very downstream locations, because there are still wave coalescences that result in increase of the wave size. The wave coalescence will be discussed in Chapter 5. The simulated wave crest depth can magnify to large levels at the downstream locations, for instance, close to 4 times of the normal water depth in Scenario (4). Yet the experimental wave crest depth can only increase to roughly 2.5 times of the normal depth due to the limited laboratory space. The development of wave trough depth is somewhat analogous to that of wave crest depth. The available experimental data only show the phase in which the trough depth decays, but they are insucient to completely display how the trough depth evolves spatially. The simulation results show the initial wave am- plitude and decaying rate of trough depth are small. The decaying rate magnies due to the instability and eventually becomes very small. The decaying rate of trough 112 depth is smaller than the growth rate of crest depth. Therefore the amplitude of trough of roll waves is smaller than that of crest. Moreover, at the location where the crest depth still grows due to the instability, the trough depth stops decaying and almost remains at a xed level. In order to make more detailed comparisons of the spatial amplication rate, we extract the channel sections where experimental data are available and magnify the characteristics of development curves from Brock's experiment and the present simu- lation. The close views of comparison results are shown in Figure 4.17 Figure 4.20 for the Scenario (1) Scenario (4). The interesting channel lengths are shorter than those in the overviews so that the characteristics of development curves can be inter- preted more clearly. The interval among interesting channel locations is decreased to 5 m. The vertical axes are the normalized time-averaged wave crest and trough depths and the horizontal axes are the normalized channel location. Although the simulation results agree well with experimental data from the overviews, the close views reveal slight dierences, particularly for the crest depth. The simulated spatial evolutions of wave trough depth still match closely with experimental data in the close views. For the crest depth, the agreements at very upstream locations where small pertur- bations dominate are quite good for all of the four scenarios. Yet the growth rates of simulation results become smaller than those of experimental data, when small per- turbations begin magnifying and transforming to roll waves. Hence the simulation results slightly underestimate the crest depth at the intermediate phases of the exper- imental data. Then the growth rates of experimental data start to decrease whereas the growth rates of simulation results keep increasing or do not change. Therefore 113 × × × × × × × × × × × × × x/h 0 ⎯ h/h 0 0 2000 4000 6000 8000 10000 12000 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 Crest depth by present simulation Crest depth by Brock’s experiment Trough depth by present simulation Trough depth by Brock’s experiment × Figure 4.17. Close view of comparison with Brock's experiment: Scenario (1). The vertical axis is the normalized time-averaged wave crest and trough depths and the horizontal axis is the normalized channel location. 114 × × × × × × × × × × × x/h 0 ⎯ h/h 0 0 1000 2000 3000 4000 5000 6000 7000 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 Crest depth by present simulation Crest depth by Brock’s experiment Trough depth by present simulation Trough depth by Brock’s experiment × Figure 4.18. Close view of comparison with Brock's experiment: Scenario (2). The vertical axis is the normalized time-averaged wave crest and trough depths and the horizontal axis is the normalized channel location. 115 × × × × × × × × × x/h 0 ⎯ h/h 0 0 1000 2000 3000 4000 5000 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 Crest depth by present simulation Crest depth by Brock’s experiment Trough depth by present simulation Trough depth by Brock’s experiment × Figure 4.19. Close view of comparison with Brock's experiment: Scenario (3). The vertical axis is the normalized time-averaged wave crest and trough depths and the horizontal axis is the normalized channel location. 116 × × × × × × × × x/h 0 ⎯ h/h 0 0 1000 2000 3000 4000 5000 6000 7000 0.4 0.7 1 1.3 1.6 1.9 2.2 2.5 2.8 Crest depth by present simulation Crest depth by Brock’s experiment Trough depth by present simulation Trough depth by Brock’s experiment × Figure 4.20. Close view of comparison with Brock's experiment: Scenario (4). The vertical axis is the normalized time-averaged wave crest and trough depths and the horizontal axis is the normalized channel location. 117 the simulated crest depth eventually rises over the experimental data. If we look back at the overviews, the locations of simulation results where the growth rates of wave crest depth begin decreasing are more downstream than those of experimental data. Generally, the comparison results between Brock's experiment and the present simulation are satisfactory, in spite of slight dierences about the growth rates of the wave crest depth. The validation tests using the spatial evolution curves of wave crest and trough depths from Brock's experiment manifest the present computational model can be utilized for simulating real roll-wave ows. 118 Chapter 5 Investigation of Roll-wave Properties In this chapter the tested computational model USC SWAT is utilized to explore some of important properties of roll waves. The spatial evolution of roll waves is demonstrated by simulation results throughout the whole channel at an instantaneous time level. The generality of the spatial evolution of roll waves is described following the demonstration. Several phenomena of wave-wave interactions in roll-wave ows are discovered by the present computational model. The spectral analysis is also implemented to study the distribution of wave energy over the frequency domain at various channel locations. At the end how the magnitude of roughness height on the channel bottom in uences the signicant wave height of roll waves is analyzed. 5.1 Demonstration of the Spatial Evolution of Roll Waves The outcome of the computational simulation of a real roll-wave ow has been shown by time series at selected channel locations in Subsection 4.3.1. The outcome can also be viewed by picking the simulation results throughout the whole channel at an instantaneous time level. As such the spatial evolution of roll waves can be dis- played in a more straightforward manner, though the time level is xed. Figure 5.1 illustrates the demonstration of the spatial evolution of roll waves down an inclined, wide, rectangular open channel. The vertical axis is the water depth and the hori- zontal axis is the channel location. Associated hydraulic parameters are q 0 = 0:008 119 x(m) h(m) 0 10 20 30 40 50 60 70 80 0.005 0.01 0.015 0.02 Figure 5.1. Demonstration of the spatial evolution of roll waves. The vertical axis is the water depth and the horizontal axis is the channel location. The prole solution is picked at the instantaneous time level t = 900 s. The hydraulic parameters are q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. The undisturbed Froude number is calculated as F 0 2:60. The total channel length is 80 m. The prole solution is picked at the instantaneous time levelt = 900 s. This instantaneous prole is capable of representing simulation results throughout the whole channel at other time levels. The small-amplitude perturbations dominate at the upstream locations within x2 [0; 20] m. Owing to the ow instability those small perturbations begin to magnify at x = 20 m approximately and keep growing considerably until x = 40 m. This dramatic amplication process results in that the wave amplitude at the location of x = 40 m is much greater than the upstream locations. Behind x = 40 m the increase of wave amplitude becomes less remarkable. Figure 5.1 can only inform the variation of wave amplitude at dierent channel locations. The whole channel can be separated into several channel sections so that the features of wave shape and wavelength can also be claried. Hence the demon- 120 stration of the spatial evolution of roll waves with separate channel sections is shown by Figure 5.2. The whole channel is separated into four sections with equal length of 20 m. It still illustrates the increase of wave amplitude from upstream to downstream locations. In addition, it shows the wave shape begins to sharpen within x2 [36; 40] m. Consequently the waves are smooth at the upstream locations of x = 36 m but have steep fronts and relatively mild backs at the downstream locations of x = 40 m. This process accompanied with the increase of wave amplitude leads to the transfor- mation from small perturbations to large, shock-like roll waves. The demonstration indicates the wavelength starts to increase at the downstream locations of x = 50 m. Brock (1969) found that the location where wave coalescences began coincided with the location where the increases of both wavelength and wave period began. There- fore the increase of the wavelength results from existences of the wave coalescence. The wave coalescences also cause slight increases of the wave amplitude. 5.2 Wave-wave Interactions Brock (1969) observed occurrences of the wave overtaking process right after the formation of shock-like roll waves. The wave overtaking that led to growth of the wave size was prevalent after the eect of the spontaneous instability terminates. The present simulation discovers three types of wave-wave interactions in real roll-wave ows: wave overtaking, wave absorption, and wave spawning. The denition of the wave overtaking process is that a wave with greater size, including amplitude and wavelength, chases down another wave with smaller size, and they merge into a new wave. Figure 5.3 illustrates a demonstration of the wave 121 h(m) 0 4 8 12 16 20 0.005 0.01 0.015 0.02 h(m) 20 24 28 32 36 40 0.005 0.01 0.015 0.02 h(m) 40 44 48 52 56 60 0.005 0.01 0.015 0.02 x(m) h(m) 60 64 68 72 76 80 0.005 0.01 0.015 0.02 Figure 5.2. Demonstration of the spatial evolution of roll waves: separate channel sections. The vertical axis is the water depth and the horizontal axis is the channel location. The prole solution is picked at the instantaneous time levelt = 900 s. The hydraulic parameters are q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. 122 overtaking process from the simulation results. The vertical axis is the water depth and the horizontal axis is the channel location. The interesting channel section is x2 [55; 67] m. The time level in which a wave overtaking process is about to occur is picked as the initial time level. The proles are recorded at the instantaneous time level 0.0, 2.2, 3.8, 5.1, 5.4, 5.6, 5.7, and 6.2 s. The hydraulic parameters are q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, t = 0:04 m 2 /s, and F 0 2:60. At the initial time level, a roll wave at the most upstream location is bigger than the wave at right downstream position. At the subsequent time levels, the bigger wave overtakes and combines with the smaller one and gradually merge into a single wave. In the whole process the combination point is always located at the front of the bigger wave and the back of the smaller wave. The process generates a new wave that is even bigger than both of the previous two waves. We nd in the simulation that a roll wave can overtake not only one but many smaller waves before it is overtaken by another even larger wave. The wave overtaking process is induced by the dierence of the ow velocity between two waves. We have deduced in Subsection 4.3.1 that the magnitude of the ow velocity is consistent with that of the water depth. Hence the wave with bigger size must possess faster ow velocity. Furthermore the wave with bigger size has faster propagation celerity according to the shallow water wave theory. The wave absorption has not been mentioned by any previous study. Yet it is found in the present simulation of a real roll-wave ow. Figure 5.4 illustrates a demonstration of the wave absorption process from the present simulation. The ver- tical axis is the water depth and the horizontal axis is the channel location. The 123 h(m) 55 57 59 61 63 65 67 0.01 0.02 0.03 t=0.0 s h(m) 55 57 59 61 63 65 67 0.01 0.02 0.03 t=2.2 s h(m) 55 57 59 61 63 65 67 0.01 0.02 0.03 t=3.8 s x(m) h(m) 55 57 59 61 63 65 67 0.01 0.02 0.03 t=5.1 s See caption on the next page. 124 h(m) 55 57 59 61 63 65 67 0.01 0.02 0.03 t=5.4 s h(m) 55 57 59 61 63 65 67 0.01 0.02 0.03 t=5.6 s h(m) 55 57 59 61 63 65 67 0.01 0.02 0.03 t=5.7 s x(m) h(m) 55 57 59 61 63 65 67 0.01 0.02 0.03 t=6.2 s Figure 5.3. Demonstration of the wave overtaking process from the simulation of a real roll-wave ow. The vertical axis is the water depth and the horizontal axis is the channel location. The proles are recorded at the instantaneous time level 0.0, 2.2, 3.8, 5.1, 5.4, 5.6, 5.7, and 6.2 s. The solutions are for the parameters of q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. 125 h(m) 50 53 56 59 62 65 0.01 0.02 0.03 t=0.0 s h(m) 50 53 56 59 62 65 0.01 0.02 0.03 t=3.3 s h(m) 50 53 56 59 62 65 0.01 0.02 0.03 t=5.0 s x(m) h(m) 50 53 56 59 62 65 0.01 0.02 0.03 t=5.8 s See caption on the next page. 126 h(m) 50 53 56 59 62 65 0.01 0.02 0.03 t=6.3 s h(m) 50 53 56 59 62 65 0.01 0.02 0.03 t=6.7 s h(m) 50 53 56 59 62 65 0.01 0.02 0.03 t=7.1 s x(m) h(m) 50 53 56 59 62 65 0.01 0.02 0.03 t=8.3 s Figure 5.4. Demonstration of the wave absorption process from the simulation of a real roll-wave ow. The vertical axis is the water depth and the horizontal axis is the channel location. The proles are recorded at the instantaneous time level 0.0, 3.3, 5.0, 5.8, 6.3, 6.7, 7.1, and 8.3 s. The solutions are for the parameters of q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. 127 interesting channel section is x2 [50; 65] m. The time level in which a wave ab- sorption process is about to occur is picked as the initial time level. The proles are recorded at the instantaneous time level 0.0, 3.3, 5.0, 5.8, 6.3, 6.7, 7.1, and 8.3 s. The hydraulic parameters are the same as those used in the demonstration of the wave overtaking process. There are two roll waves quite adjacent to each other to the most upstream locations at the initial time level. The upstream one is smaller than the downstream one. Then the entire body of the smaller wave is gradually absorbed by the back of the bigger wave. Eventually the smaller wave submerges into the back of the bigger one and a new wave is generated that is larger than both of them. The wave absorption also begins after the formation of shock-like roll waves. The wave absorption process is dierent from the wave overtaking process in that the size of the upstream wave is smaller than the size of the downstream wave. Indeed the wave absorption process is also induced by the dierence of the ow velocity. Although the downstream wave is bigger, the water depth and thus the ow velocity of the crest of the smaller wave is greater than the back of the bigger wave. Therefore the front of the smaller wave is able to catch up with the back of the bigger one. However the entire process seems that the smaller wave is absorbed by the bigger wave, so we term it the wave absorption. The wave absorption can only be found when the distance between two waves is very short. Both the wave overtaking and absorption processes result in the wave coalescence. Wave coalescences start to appear right after the formation of roll waves and exist prevalently in roll-wave ows. The combination of multiple waves generates a novel wave with increased size. Hence the wave coalescence is the other reason besides the 128 natural instability for the increase of the wave amplitude from upstream to down- stream locations. The occurrences of the wave coalescence also enlarge the wavelength and wave period. The wave spawning process has not been observed in real roll-wave ows. Balm- forth & Mandre (2004) numerically explored the existence of the wave-spawning in- stability by solving the Saint Venant equations and their amplitude equation. A demonstration of the wave spawning process from the present simulation is illus- trated by Figure 5.5. The vertical axis is the water depth and the horizontal axis is the channel location. The interesting channel location is x2 [10; 20] m. The time level in which a wave spawning process is about to occur is picked as the initial time level. The proles are recorded at the instantaneous time level 0.0, 1.2, 2.2, 3.5, 4.9, 5.7, 6.3, and 7.4 s. The hydraulic parameters are the same as those used in the demonstration of the wave overtaking process, except that the value of the turbulent viscosity is set to be 0.001 m 2 /s. There is a small perturbation, which is articially imposed at the in ow boundary, to the very upstream locations at the initial time level. The perturbation self-magnies as propagating forward because of the natural instability. Meanwhile the upstream end of the perturbation begins to subside below the normal water depth and a trough depression is generated. A swell perturbation subsequently emerges at the upstream end of the depression. All of the perturbations continue to magnify as the swell grows to a new wave. As such a wave spawning process has accomplished. We nd one of the necessary conditions for the occurrence of the wave spawning is that the distance between two articially imposed waves must be very large. In other words, the original wave should have an extremely long tail at 129 h(m) 10 12 14 16 18 20 0.009 0.01 0.011 t=0.0 s h(m) 10 12 14 16 18 20 0.009 0.01 0.011 t=1.2 s h(m) 10 12 14 16 18 20 0.009 0.01 0.011 t=2.2 s x(m) h(m) 10 12 14 16 18 20 0.009 0.01 0.011 t=3.5 s See caption on the next page. 130 h(m) 10 12 14 16 18 20 0.009 0.01 0.011 t=4.9 s h(m) 10 12 14 16 18 20 0.009 0.01 0.011 t=5.7 s h(m) 10 12 14 16 18 20 0.009 0.01 0.011 t=6.3 s x(m) h(m) 10 12 14 16 18 20 0.009 0.01 0.011 t=7.4 s Figure 5.5. Demonstration of the wave spawning process from the simulation of a real roll-wave ow. The vertical axis is the water depth and the horizontal axis is the channel location. The proles are recorded at the instantaneous time level 0.0, 1.2, 2.2, 3.5, 4.9, 5.7, 6.3, and 7.4 s. The solutions are for the parameters of q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:001 m 2 /s. 131 the upstream side. Hence, successive wave spawning processes will occur to produce a series of new waves at the upstream side of the original wave until the long tail van- ishes. Instead of adding perturbations at each time level, we impose low-frequency perturbations with relatively large time interval (e.g. 1 s) at the in ow boundary so as to show the wave spawning process. However, the water surface of turbulent ows always has perturbations in reality, so the wave spawning process has not been observed. Based on the demonstration the wave spawning denotes the birth of a new wave at the upstream side of an existing wave. In the simulation, the wave generated by the wave spawning process comes out spontaneously due to the nonlinear wave dynamics. It does not belong to the perturbations articially imposed at the in ow boundary. The present demonstration shows the wave spawning process at upstream locations where roll waves have not formed. The wave spawning can also occur for shock-like roll waves in a numerical experiment like what is shown by Balmforth & Mandre. In the present simulation of real roll-wave ows, the wave-wave interactions de- scribed above characterize the spatial evolution of roll waves. Spontaneous processes of the wave coalescence, including the wave overtaking and absorption, under the case of shorter distance and the wave spawning under the case of longer distance help determine the scale of the wave amplitude and wavelength. 5.3 Generality of the Spatial Evolution of Roll Waves There is a generality for the spatial evolution of roll waves, which is indepen- dent on the hydraulic parameters such as the ow rate, channel slope, and bottom 132 h Initial phase Transition phase h h t h 0 h Final phase Transition phase c h x c h t h Figure 5.6. Sketch illustration for the generality of the spatial evolution of roll waves. Based on the development curves of wave crest and trough depths versus the channel location, the spatial evolution of roll waves is characterized by the initial phase, transition phase, and nal phase. roughness height. The generality, presented by Figure 5.6, is on the basis of the development curves of wave crest and trough depths versus the channel location. The generality contends that the spatial evolution of roll waves with any hydraulic parameters consists of three phases: initial phase, transition phase, and nal phase. The initial phase starts from the channel inlet. Both the amplitude and growth rate of perturbations on the normal depth are small in the initial phase. The small perturbations maintain their smooth shape without sharpening. The initial phase is free of the wave coalescence, so the wave period and wavelength nearly do not change. 133 The inceptive boundary of the transition phase is the location where the growth rate of the wave amplitude suddenly becomes large. In the transition phase the wave amplitude magnies signicantly due to the spontaneous instability. All of the small perturbations from the initial phase transform into large, shock-like roll waves with steep fronts. The wave coalescence begins to occur right after the formation of roll waves. The end of the transition phase is the location where the natural growth of the wave amplitude due to the instability terminates. In the nal phase, the evolved roll waves propagate without breaking, but there are plenty of existences of the wave coalescence. The wave crest depth continues to grow with small rates whereas the wave trough depth almost remains at a xed level. The wave coalescence in the transition and nal phases is the other factor that causes the increase of the wave amplitude besides the natural instability. The wave period and wavelength are also increased due to occurrences of the wave coalescence. 5.4 Spectral Analysis of Roll Waves The simulated time series in Subsection 4.3.1 corroborate the natural roll waves are nonperiodic at any channel location. Herein we take samples of simulated time series of the wave amplitude of roll waves at selected channel locations, and perform the spectral analysis on those time series. As such the characteristics of variations of the frequency or period of roll waves along the whole channel can be discovered. The samples of time series are presented in Figure 5.7. The vertical axis is the wave amplitude, obtained by subtracting the normal depth from the computed water depth, and the horizontal axis is time. For the convenience of performing the spectral 134 analysis, the time step of this computation is xed to be 0.0004 s that is suciently small to guarantee the CFL stability condition Cr < 1. The representative channel locations are x = 10, 30, 50, 80, 100, and 120 m and the interesting time period is t2 [150; 200] s. Associated hydraulic parameters are q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, t = 0:04 m 2 /s, andF 0 2:60. The fast Fourier transform is executed on the time series presented in Figure 5.7 so as to achieve the features of the wave amplitude distribution on the frequency domain. The spectral analysis results corresponding to the time series in Figure 5.7 are shown in Figure 5.8. The wave amplitude is described as a function of frequency. The dominant wave period(s) at each location is(are) indicated. The spectral analysis results show the wave amplitude increases when the channel location increases at the front of x = 50 m. This is consistent with the time series showing the natural growth of the wave amplitude due to the instability. The wave period has a tendency to increase from upstream to downstream locations. For in- stance, the dominant wave period is around 1.52 s at x = 10 m while shifts to 4.55 s atx = 120 m. We nd the wave period nearly does not change at locations ofx< 30 m but begins to increase behind. Also x = 30 m is roughly the location where small perturbations transform to shock-like roll waves and waves begin to coalesce after the formation of roll waves. The shock-like roll waves have larger wavelength than small perturbations and the wave coalescence results in a novel wave with greater wave- length. Thus the increase of the wave period should be attributed to the formation of shock-like roll waves and existences of the wave coalescence. Brock (1969) also manifested the location where the wave overtaking began coincided with the location 135 150 160 170 180 190 200 −5 0 10 20 x 10 −3 Amplitude (m) x = 10 m 150 160 170 180 190 200 −5 0 10 20 x 10 −3 Amplitude (m) x = 30 m 150 160 170 180 190 200 −5 0 10 20 x 10 −3 Time (s) Amplitude (m) x = 50 m See caption on next page. 136 150 160 170 180 190 200 −5 0 10 20 x 10 −3 Amplitude (m) x = 80 m 150 160 170 180 190 200 −5 0 10 20 x 10 −3 Amplitude (m) x = 100 m 150 160 170 180 190 200 −5 0 10 20 x 10 −3 Time (s) Amplitude (m) x = 120 m Figure 5.7. Presentation of simulated time series of the water depth of roll waves for the spectral analysis. The vertical axis is the wave amplitude and the horizontal axis is time. The time series are recorded within t2 [150; 200] s at channel locations of 10, 30, 50, 80, 100, and 120 m. The solutions are for the parameters of q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. 137 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 x 10 −5 Amplitude (m) 1.52 sec x = 10 m 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 2 4 6 x 10 −4 Amplitde (m) 1.72 sec x = 30 m 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 x 10 −3 Frequency (Hz) Amplitude (m) 2.5 sec 2.27 sec 1.92 sec x = 50 m See caption on next page. 138 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 x 10 −3 Amplitude (m) 2.38 sec x = 80 m 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 x 10 −3 Amplitude (m) 3.13 sec 2.63 sec x = 100 m 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 x 10 −3 Frequency (Hz) Amplitude (m) 4.55 sec 3.33 sec x = 120 m Figure 5.8. Presentation of the amplitude of roll waves in the frequency domain corresponding to the time series shown in Figure 5.7. The vertical axis is the wave amplitude and the horizontal axis is frequency. The dominant wave period(s) is(are) indicated. The spectral results are analyzed at locations of 10, 30, 50, 80, 100, and 120 m. The solutions are for the parameters of q 0 = 0:008 m 2 /s, tan = 0:05, k s = 0:001 m, and t = 0:04 m 2 /s. 139 where the wave period began to increase. The spectral analysis results illustrate the wave energy is distributed on wider frequency range at more downstream locations. Except for the locations of x = 10 and 30 m, multiple spikes exist on the whole frequency range and the amplitude of those spikes can be considerably high, espe- cially for the locations of x = 100 and 120 m. The wide frequency distribution of the wave amplitude manifests roll-wave ows are nonperiodic and the diversity of the wave period becomes more signicant at more downstream locations. Brock (1967) suggested the standard deviation of the wave period rose up when the channel loca- tion increased. Hence the experiment also deduces more various wave periods exist in addition to the average wave period at more downstream locations. 5.5 Eect of Bottom Roughness The absolute roughness height at the channel bottom is a signicant factor that impacts the characteristics of open channel ows. In this section, we investigate how the spatial evolution of roll waves is in uenced by the bottom roughness. We perform a numerical experiment using the present computational model with constant discharge per unit width, channel slope, and turbulent viscosity, but variable bottom roughness height. Five scenarios are considered and the hydraulic parameters are presented by Table 5.1. The normal depth and velocity are computed according to equation 4.51, 4.56, and 4.57. When the bottom roughness increases from 0.2 mm to 10 mm, the resulting normal depth increases whereas the normal velocity decreases. Consequently the undisturbed Froude number associated with normal ows decreases and the ow becomes less supercritical. 140 Bottom roughness (mm) Discharge (m 2 /s) Channel slope Turbulent viscosity (m 2 /s) Normal depth (m) Normal velocity (m/s) Froude number 0.2 0.008 0.05 0.04 0.00841 0.952 3.32 0.5 0.008 0.05 0.04 0.00915 0.874 2.92 0.5 0.008 0.05 0.04 0.00915 0.874 2.92 1 0.008 0.05 0.04 0.00988 0.810 2.60 2 0.008 0.05 0.04 0.0108 0.742 2.28 10 0.008 0.05 0.04 0.0140 0.572 1.55 Table 5.1. Presentation of hydraulic parameters applied in the numerical experiment for investigating the eect of bottom roughness on the spatial evolution of roll waves. Five scenarios are considered with constant discharge per unit width, channel slope, and turbulent viscosity, but variable bottom roughness. Resulting normal depth, normal velocity, and undisturbed Froude number are also indicated. In order to investigate the eect of the bottom roughness, we analyze the signif- icant wave height at interesting channel locations for each scenario. The signicant wave height H 1/3 , measured as the average of the highest one-third of waves (trough to crest) in a given time period, represents the average of largest waves. Most waves are smaller than the signicant wave height because of the approach of its measure- ment. Hence time series of the water depth, including more than 200 waves, are generated by the computational simulation at interesting channel locations. Then the signicant wave height is calculated as H 1/3 = 4 h (5.1) where h stands for the standard deviation of a sample time series of the water surface elevation. Figure 5.9 shows the development of the signicant wave height 141 x(m) H 1/3 (m) 10 20 30 40 50 60 70 80 0 0.004 0.008 0.012 0.016 0.02 k s =0.2 mm k s =0.5 mm k s =1 mm k s =2 mm k s =10 mm Figure 5.9. Presentation of development of the signicant wave height versus the channel location for dierent values of the bottom roughness. Five scenarios of the bottom roughness are analyzed: k s = 0:2, 0.5, 1, 2, and 10 mm. Associated hydraulic parameters are q 0 = 0:008 m 2 /s, tan = 0:05, and t = 0:04 m 2 /s. versus the channel location associated with ve scenarios of dierent values of the bottom roughness. The minimum interesting channel location is x = 10 m while the maximum isx = 80 m and the interval is 10 m. The maximum time level reaches 800 s to include 200 waves. The development of the signicant wave height is substantially analogous to the development of the wave crest depth. The signicant wave height for the scenarios of k s = 0:2, 0.5, 1, 2 mm grows spatially due to the spontaneous instability as well as wave coalescences. On the other hand, the signicant wave height becomes smaller when the bottom roughness is increased. We have shown 142 x(m) H 1/3 /h 0 10 20 30 40 50 60 70 80 0 0.4 0.8 1.2 1.6 2 2.4 k s =0.2 mm k s =0.5 mm k s =1 mm k s =2 mm k s =10 mm Figure 5.10. Presentation of development of the normalized signicant wave height versus the channel location for dierent values of the bottom roughness. Five scenarios of the bottom roughness are analyzed: k s = 0:2, 0.5, 1, 2, and 10 mm. Associated hydraulic parameters are q 0 = 0:008 m 2 /s, tan = 0:05, and t = 0:04 m 2 /s. the resulting undisturbed Froude number decreases when the bottom roughness is increased. Therefore the ow becomes less unstable and the growth rate of the wave height also decreases. The signicant wave height for the scenario of k s = 10 mm lies on zero at all channel locations. This implies there are no waves on the water surface. The resulting undisturbed Froude number associated with k s = 10 mm is less than two, manifesting the bottom roughness is suciently large to prevent the destabilization while attenuate the wave amplitude and stabilize the water ow. The relationship between the normal depth and signicant wave height can be 143 interpreted by analyzing the normalized signicant wave height. The normalized signicant wave height is dened as the ratio of the signicant wave height to the normal depth. Figure 5.10 presents the normalized signicant wave height versus the channel location. The results are quite similar to the absolute signicant wave height, except that the dierences among each scenario are more remarkable. When the bottom roughness is increased, the signicant wave height decreases but the normal depth increases. Therefore the relatively small signicant wave height becomes even smaller after it is normalized. The wave height for k s = 10 mm is still zero because of the stabilization. For k s = 0:2 mm, the signicant wave height is already 2 times of the normal depth at x = 50 m and nearly reaches 2.3 times at x = 80 m. Yet the normalized signicant wave height is only around 1.5 for k s = 2 mm at the most downstream location. Hence the reduction of the ow stability attributed to the increase of the bottom roughness is clearly perceived. According to the above analyses regarding the signicant wave height we empha- size that the increase of the bottom roughness is a possible way of mitigating the roll-wave evolution in unstable open channel ows, when other hydraulic parameters, such as the discharge and channel slope, can not be modied. 144 Chapter 6 Roll Waves in a Breaking-slope Open Channel When a wide, rectangular open channel is suciently steep and smooth, the water ows could be spontaneously unstable. The natural instability will lead to large, shock-like roll waves on the water ow surface. Roll-wave ows are undesirable because they could bring about ood threats to adjacent environments and disrupt the original purpose of designing the channel. Therefore we should avoid the occurrence of roll-wave ows utilizing the natural mechanism of the spatial evolution of roll waves. In Section 5.5 we have explored the roll-wave evolution can be mitigated by increasing the magnitude of the bottom roughness. In this chapter, a concept of the breaking-slope open channel to prevent the formation of roll waves is specied and tested through a computational experiment. Assume a wide, rectangular open channel is designed to accommodate a water ow with a discharge per unit width 0.008 m 2 /s. The prole geometry of the channel is illustrated by Figure 6.1. The horizontal distance in the longitudinal direction is 100 m and the elevation drop is 5 m, so the channel slope is 0.05 and the total channel length is 100.125 m. The roughness height on the channel bottom is xed to be 1 mm. Hence, using equation 4.51, 4.56, and 4.57, the normal conditions can be obtained and the resulting undisturbed Froude number is F 0 2:6. The ow is therefore unstable and large roll waves will occur on the water surface. When the discharge and bottom roughness are invariable, the only approach to prevent the formation of roll waves is to make the channel slope less steep. Herein 145 5 m 100 m 100 m Figure 6.1. Sketch illustration of the prole geometry of the original open channel. The elevation drop and horizontal distance in the longitudinal direction are indicated. we replace the original channel of a constant slope with a modied, breaking-slope channel with two dierent slopes. The modied channel is illustrated by Figure 6.2. The hydraulic conditions along with the channel length and elevation drop for each reach are listed in Table 6.1. This breaking-slope channel is composed of three con- terminous reaches. Reach (1) and Reach (3) at two ends of the channel retain the slope of the original channel, except that the horizontal distance in the longitudinal direction is shortened to 30 m. As a result the length and elevation drop of them are also shortened to about 30.038 m and 1.5 m, respectively. The undisturbed Froude number in Reach (1) and Reach (3) is therefore still 2.60, which indicates the existence of the natural instability. However Reach (2) in the middle is the portion making the modied channel distinctive from regular constant-slope channels. The slope of Reach (2) is adjusted to be 0.015, which is less steep than Reach (1) and Reach (3). For the whole channel, the total elevation drop is xed to be 5 m and the total horizontal distance in the longitudinal direction is xed to be 100 m. Hence the elevation drop of Reach (2) 146 Reach (1) Reach 5 m 30 m (a) (b) Reach (2) Reach (3) 40 m 30 m 63.596 m Figure 6.2. Sketch illustration of the geometry of the modied, breaking-slope open channel. The whole channel is divided into three reaches. (a) shows the prole view in which the total elevation drop and the horizontal distance in the longitudinal direction for each reach are indicated. (b) shows the plan view of the channel path in which the horizontal scale of Reach (2) in the transverse direction is indicated. 147 Reach Channel slope Channel length (m) Elevation drop (m) (1) 0.05 30.038 1.5 (2) 0.015 133.348 2 (3) 0.05 30.038 1.5 Elevation drop (m) Normal depth (m) Normal velocity (m/s) Froude number 1.5 0.00988 0.810 2.60 2 0.0141 0.568 1.53 1.5 0.00988 0.810 2.60 Table 6.1. Presentation of hydraulic parameters for three reaches in the modied, breaking-slope open channel. The channel length and elevation drop for each reach are also indicated. must be 2 m; the horizontal distance of Reach (2) in the longitudinal direction must be 40 m. Because of its slope and elevation drop, the channel length of Reach (2) should be around 133.348 m. If Reach (2) is straight, the elevation could not drop 2 m within the horizontal distance of 40 m, and also the channel length could not reach 133.348 m. Therefore we set up three bending points by utilizing the horizontal space extendible in the transverse direction. As illustrated by Figure 6.2(b), the bends are located at the junction between Reach (1) and Reach (2), the junction between Reach (2) and Reach (3), and the middle point of Reach (2). The consequent path of the whole channel appears to be axisymmetric. The paths of Reach (1) and Reach (3) only stretch in the longitudinal direction, whereas Reach (2) covers distances in both longitudinal and transverse directions. The horizontal distance extending in the transverse direction is 63.596 m approximately. As such the length of Reach (2) reaches 133.348 m and the elevation drop is able to reach 2 m. Table 6.1 indicates the normal depth, normal velocity, and undisturbed Froude number associated with Reach (2) have changed due to the variation of the channel slope. The normal depth 148 has increased whereas the normal velocity has decreased. Thus the Froude number becomes 1.53 which implies the stabilization of the water ow. The present computational model, USC SWAT, is utilized to simulate water ows in the constant-slope channel and breaking-slope channel. The discharge, bottom roughness, and turbulent viscosity are q 0 = 0:008 m 2 /s, k s = 0:001 m, and t = 0:04 m 2 /s in both channels. The cell size is selected as 0.001 m and the Courant number is 0.75. In the computation of the breaking-slope channel ow, the total channel length is 193.424 m comprising all of the three reaches. The channel slope is assigned to the computational domain of each reach according to the design of the breaking-slope channel. Any eect on the ow conditions that might be caused by the bends in the transverse horizontal direction is disregarded. Before simulating the roll-wave ow in the breaking-slope channel, we are inter- ested in investigating the gradually varied ows resulting from the slope changes. Hence the normal ow conditions are simulated in the breaking-slope channel with no perturbations on the water surface. Figure 6.3 shows the transitions of the normal water depth in (a) and normal ow velocity in (b) between Reach (1) and Reach (2). The location where the channel slope breaks is indicated. The Froude number in both reaches are greater than 1 and the slope in Reach (2) is smaller than the slope in Reach (1). Therefore the transition from a steep slope to a less steep slope appears. When water passes through the breaking point, the lower depth from Reach (1) has to climb up the higher depth in Reach (2). The transition water depth between two normal depths presents a convex prole with negative curvatures. Thus this is a S- 3 type gradually varied ow (Sturm 2001). The water depth on Reach (1) does not 149 change. However the upstream ow controls and exerts its in uence on the water ow on Reach (2). Note the water surface slope approximates to zero at both the upstream and downstream ends of the transition. The gradually varied ow is nonuniform but steady ow. The ow velocity should corresponds to the water depth so as to main- tain the normal ow rate. Hence the ow velocity in the transition presents a concave prole before turning to the normal condition on Reach (2). Although the Reach (2) is less steep than Reach (1), the ow on it is still supercritical. There does not exist a hydraulic jump between Reach (1) and Reach (2). Similarly the transitions of the normal water depth and normal ow velocity between Reach (2) and Reach (3) are shown in Figure 6.4. Since the slope in Reach (2) is smaller than the slope in Reach (3), it is a transition from a less steep slope to a steep slope without a hydraulic jump. When water passes through the breaking point, the higher depth from Reach (2) has to drop to the lower depth in Reach (3). This is a S-2 type gradually varied ow in which the transition depth presents a concave prole with positive curvatures. The upstream ow in Reach (2) nearly keeps unchanged while controlling and exerting its in uence on the water ow in Reach (3). The corresponding ow velocity presents a convex prole in Reach (3) so that the ow rate is maintained. The small perturbations are imposed on the normal ow at the in ow bound- ary so that we can simulate the roll-wave evolution in both the constant-slope and breaking-slope channels. In the constant-slope channel, we nd small perturbations magnify and become large, shock-like roll waves due to the natural instability. In the breaking-slope channel, small perturbations also magnify in Reach (1) and Reach (3). Nevertheless the ow characteristics are quite distinctive in Reach (2). 150 x(m) h 0 (m) 29 30 31 32 33 34 0.01 0.011 0.012 0.013 0.014 0.015 Reach (2) Reach (1) (a) x(m) u 0 (m/s) 29 30 31 32 33 34 0.55 0.6 0.65 0.7 0.75 0.8 0.85 Reach (2) Reach (1) (b) Figure 6.3. Prole of normal ow conditions on the transition between Reach (1) and Reach (2). The vertical axis is the normal water depth in (a) and normal ow velocity in (b) while the horizontal axis is the channel location. The junction between Reach (1) and Reach (2) where the channel slope breaks is indicated. The hydraulic parameters are q 0 = 0:008 m 2 /s, k s = 0:001 m, t = 0:04 m 2 /s, and tan = 0:05 in Reach (1) and 0.015 in Reach (2). 151 x(m) h 0 (m) 162 163 164 165 166 0.01 0.011 0.012 0.013 0.014 0.015 Reach (3) Reach (2) (a) x(m) u 0 (m/s) 162 163 164 165 166 0.55 0.6 0.65 0.7 0.75 0.8 0.85 Reach (3) Reach (2) (b) Figure 6.4. Prole of normal ow conditions on the transition between Reach (2) and Reach (3). The vertical axis is the normal water depth in (a) and normal ow velocity in (b) while the horizontal axis is the channel location. The junction between Reach (2) and Reach (3) where the channel slope breaks is indicated. The hydraulic parameters are q 0 = 0:008 m 2 /s, k s = 0:001 m, t = 0:04 m 2 /s, and tan = 0:015 in Reach (2) and 0.05 in Reach (3). 152 Figure 6.5 presents a sample for the simulated time series of the water depth on Reach (2) of the breaking-slope channel. The vertical axis is the water depth and the horizontal axis is time. The proles are recorded within t2 [150; 200] s at channel locations x = 38, 48, 58, and 68 m. The junction between Reach (1) and Reach (2) is located at x = 30:038 m. Hence the prole at x = 38 m shows the amplied perturbations from Reach (1) where the water ow is destabilized. Those perturbations have not transformed to shock-like roll waves since Reach (1) is not suciently long. The water depth proles at x = 48, 58, and 68 m in Reach (2) manifest the wave amplitude decreases from upstream to downstream locations. The attenuation of the wave amplitude is attributed to the stabilization of the water ow in Reach (2). Reach (2) has a less steep channel slope that results in the less supercritical ow. The perturbations in Reach (2) maintain their smooth shapes without sharpening. The wave amplitude in Reach (2) keeps attenuating spatially until the junction between Reach (2) and Reach (3) at x = 163:386 m. However the wave amplitude never returns to zero and there are always perturbations on the water surface. After the perturbations enter Reach (3), they start to magnify but do not transform to shock-like roll waves because of the limited channel length. The spatial development curves for the time-averaged wave crest and trough depths are obtained from the simulation results in the constant-slope and breaking- slope channels. In order to include more than 200 waves for the arithmetic mean at each channel location, the maximum time level in the simulations reaches 800 s for the constant-slope channel and 1000 s for the breaking-slope channel. The comparison of the time-averaged wave crest depth between the constant- 153 h(m) 150 160 170 180 190 200 0.013 0.014 0.015 x = 38 m h(m) 150 160 170 180 190 200 0.013 0.014 0.015 x = 48 m h(m) 150 160 170 180 190 200 0.013 0.014 0.015 x = 58 m t(s) h(m) 150 160 170 180 190 200 0.013 0.014 0.015 x = 68 m Figure 6.5. Demonstration of time series of the water depth on Reach (2) where the ow is stabilized. The vertical axis is the water depth and the horizontal axis is time. The proles are recorded at the channel locationsx = 38, 48, 58, 68 m. The hydraulic parameters are q 0 = 0:008 m 2 /s, k s = 0:001 m, t = 0:04 m 2 /s, and tan = 0:015 in Reach (2). 154 slope and breaking-slope channels is shown in Figure 6.6(a). The vertical axis is the time-averaged wave crest depth while the horizontal axis is the channel location. In the constant-slope channel, the most upstream interesting location for recording the wave crest depth is at x = 3 m and the most downstream interesting location is at x = 98 m; the interval among adjacent interesting locations is 5 m. In Reach (1) of the breaking-slope channel, the most upstream interesting location is atx = 3 m while the most downstream interesting location is at x = 28 m; the interval among adjacent interesting locations is 5 m. In Reach (2) of the breaking-slope channel, the most upstream interesting location is at x = 38 m while the most downstream interesting location is at x = 158 m; the interval among adjacent interesting locations is 10 m. In Reach (3) of the breaking-slope channel, the most upstream interesting location is at x = 168 m while the most downstream interesting location is at x = 188 m; the interval among adjacent interesting locations is 5 m. The breaking-slope channel is longer than the constant-slope channel, although they have identical elevation drop and horizontal distance in the longitudinal direction. The junctions between adjacent reaches in the breaking-slope channel are indicated. In the constant-slope channel, the initial growth of the wave crest depth is quite small. Behindx = 20 m roughly, the wave crest depth increases signicantly to a high level. The spatial growth returns to small rates near the end of the channel. Hence the spatial evolution of roll waves in the constant-slope channel conforms to the generality discussed in Section 5.3 because of the natural instability and wave coalescence processes. In the breaking-slope channel, the evolution curve in Reach (1) is identical to the constant-slope channel curve, since they have a common slope. In Reach (2), the normal water depth is higher than that 155 of Reach (1), so the wave crest depth climbs to a level that is slightly higher than the normal depth in Reach (2). The wave crest depth begins to decrease in Reach (2) due to the stabilization of the water ow. Then it almost remains at a constant level quite close to the normal depth until the end of the reach. In Reach (3), the wave crest depth drops back to the lower normal depth. Because the wave amplitude attenuates in Reach (2), the wave crest depth is very close to the normal depth at the most upstream location of Reach (3). The wave crest depth starts to increase in Reach (3) due to the natural instability. Even though the normal depth in Reach (2) is relatively high, the maximum water depth represented by the time-averaged wave crest depth is much higher in the constant-slope channel, resulting from the occurrence of large roll waves. Figure 6.6(b) shows the comparison of the normalized time-averaged wave crest depth. The time-averaged wave crest depth is divided by the normal water depth in the corresponding channel or reach. The wave crest depth in the constant-slope channel almost increases to 2.2 times of the normal depth at the end of the channel. However the maximum value of normalized wave crest depth is constrained beneath 1.2 in the breaking-slope channel. Even though the water ow is unstable in Reach (1) and Reach (3), the wave crest depth does not magnify remarkably since these reaches have limited lengths. In Reach (2) of the breaking- slope channel, the wave crest depth after the normalization drops to the level just above 1.0 and keeps constant until the end of the reach. Therefore the suppressing eect on the wave crest depth in Reach (2) ensures the normalized wave crest depth in Reach (3) is extremely close to 1.0 initially and does not increase to a high level at the end of the breaking-slope channel. 156 × × × × × × × × × × × × × × × × × × × × Reach (2) Reach (1) Reach (2) Reach (3) x(m) ⎯ h c (m) 0 20 40 60 80 100 120 140 160 180 200 0.01 0.012 0.014 0.016 0.018 0.02 0.022 Constant-slope channel Breaking-slope channel × (a) × × × × × × × × × × × × × × × × × × × × Reach (2) Reach (1) Reach (2) Reach (3) x(m) ⎯ h c /h 0 0 20 40 60 80 100 120 140 160 180 200 1 1.2 1.4 1.6 1.8 2 2.2 Constant-slope channel Breaking-slope channel × (b) Figure 6.6. Comparison of the time-averaged wave crest depth between the constant- slope and breaking-slope channels. The vertical axis is the time-averaged wave crest depth in (a) and normalized time-averaged wave crest depth in (b) while the horizontal axis is the channel location. The junctions between adjacent Reaches in the breaking- slope channel are indicated. The common hydraulic parameters are q 0 = 0:008 m 2 /s, k s = 0:001 m, and t = 0:04 m 2 /s. 157 × × × × × × × × × × × × × × × × × × × × Reach (2) Reach (1) Reach (2) Reach (3) x(m) ⎯ h t (m) 0 20 40 60 80 100 120 140 160 180 200 0.006 0.008 0.01 0.012 0.014 Constant-slope channel Breaking-slope channel × (a) × × × × × × × × × × × × × × × × × × × × Reach (2) Reach (1) Reach (2) Reach (3) x(m) ⎯ h t /h 0 0 20 40 60 80 100 120 140 160 180 200 0.6 0.7 0.8 0.9 1 Constant-slope channel Breaking-slope channel × (b) Figure 6.7. Comparison of the time-averaged wave trough depth between the constant-slope and breaking-slope channels. The vertical axis is the time-averaged wave trough depth in (a) and normalized time-averaged wave trough depth in (b) while the horizontal axis is the channel location. The junctions between adjacent Reaches in the breaking-slope channel are indicated. The common hydraulic param- eters are q 0 = 0:008 m 2 /s, k s = 0:001 m, and t = 0:04 m 2 /s. 158 The comparison of the time-averaged wave trough depth between the constant- slope and breaking-slope channels is shown in Figure 6.7(a). The vertical axis is the time-averaged wave trough depth while the horizontal axis is the channel location. The selections of the interesting locations for recording the wave trough depth are the same as those for the wave crest depth, in both constant-slope and breaking- slope channels. In the constant-slope channel, the decaying rate of the wave trough depth is quite small near the channel inlet, and becomes fairly large behind x = 20 m. Eventually the spatial decaying returns to small rates near the end of the channel. Hence the development curve of the wave trough depth in the constant- slope channel indicates the spatial evolution of large roll waves occurs owing to the natural instability. In Reach (1) of the breaking-slope channel, the time-averaged wave trough depth coincides with that in the constant-slope channel, since they have the same slope. Entering Reach (2), the wave trough depth rises upward because the normal water depth in Reach (2) is higher than Reach (1). The curve increases to a level just lower than the normal depth, manifesting the existence of perturbations amplied in Reach (1). Then the wave trough depth begins increasing slightly and almost holds at a level quite close to the normal depth. In Reach (3), the curve drops back to the lower normal depth, and the wave trough depth starts to decrease due to the destabilization of the water ow. Figure 6.7(b) shows the comparison of the normalized time-averaged wave trough depth. The normalized wave trough depth in the constant-slope channel drops to a level that is lower than 0.6. Nevertheless the wave trough depth in the breaking-slope channel only decreases to a level which is just below 0.9 times of the normal depth. In Reach (1) of the breaking-slope channel, 159 the normalized wave trough depth decays due to the natural instability. However it does not fall signicantly since the length of Reach (1) is restricted. In Reach (2) the water ow is stabilized, so the normalized wave trough depth rises up and remains at a level close to 1.0 until the end of the reach. The attenuation of the wave amplitude in Reach (2) ensures the normalized wave trough depth in Reach (3) is extremely close to 1.0 initially and does not decay to a low level at the end of the breaking-slope channel. The comparison of the signicant wave height development curves between the constant-slope and breaking-slope channels is shown in Figure 6.8(a). The signicant wave height, obtained by the means described in Section 5.5, is plotted as a function of the channel location. The selections of the interesting locations for recording the signicant wave height are the same as those for the wave crest depth, in both constant-slope and breaking-slope channels. In the constant-slope channel, the spatial growth of the signicant wave height is quite analogous to that of the time-averaged wave crest depth, manifesting the existence of large, evolved roll waves. In Reach (1) of the breaking-slope channel, the signicant wave height is identical to the curve in the constant-slope channel, because the water ow is destabilized. Entering Reach (2), the signicant wave height begins to decrease and then remains at a level close to zero. In Reach (3), it starts growing again indicating the amplication of the wave height. It can be clearly seen the absolute signicant wave height in the breaking- slope channel is much smaller than that in the constant-slope channel. Figure 6.8(b) shows the comparison of the normalized signicant wave height. The signicant wave height is divided by the normal water depth in the corresponding channel or reach. 160 ×× × × × × × × × × × × × × × × × ××× Reach (2) Reach (1) Reach (2) Reach (3) x(m) H 1/3 (m) 0 20 40 60 80 100 120 140 160 180 200 0 0.003 0.006 0.009 0.012 0.015 0.018 Constant-slope channel Breaking-slope channel × (a) ×× × × × × × × × × × × × × × × × ××× Reach (2) Reach (1) Reach (2) Reach (3) x(m) H 1/3 /h 0 0 20 40 60 80 100 120 140 160 180 200 0 0.3 0.6 0.9 1.2 1.5 1.8 Constant-slope channel Breaking-slope channel × (b) Figure 6.8. Comparison of the signicant wave height between the constant-slope and breaking-slope channels. The vertical axis is the signicant wave height in (a) and normalized signicant wave height in (b) while the horizontal axis is the channel location. The junctions between adjacent Reaches in the breaking-slope channel are indicated. The common hydraulic parameters are q 0 = 0:008 m 2 /s, k s = 0:001 m, and t = 0:04 m 2 /s. 161 The signicant wave height magnify to 1.9 times of the normal depth at the end of the constant-slope channel. Nevertheless the normalized signicant wave height is constrained beneath 0.3 in the whole breaking-slope channel. The wave height grows in Reach (1) and Reach (3) due to the natural instability. Yet the amplications are not considerable since the lengths of Reach (1) and Reach (3) are limited. In Reach (2), the normalized signicant wave height attenuates and holds at a low level owing to the stabilization of the water ow. Hence the perturbation amplitude is quite small at the entrance of Reach (3). Through this computational experiment depicted above, we prove growth of the roll-wave amplitude can be signicantly suppressed in the breaking-slope channel. With the same set of hydraulic conditions including the discharge and bottom rough- ness, but dierent channel slopes, the wave amplitude in the breaking-slope channel is much smaller than in the constant-slope channel. Indeed shock-like roll waves do not occur in the whole breaking-slope channel. Hence the robustness of the channel is protected from abrupt variations of the water pressure induced by roll waves. Also the maximum water depth in the breaking-slope channel is smaller, even if the normal water depth is relatively high in Reach (2). In this computational experiment, the slope of Reach (1) and Reach (3) in the breaking-slope channel is chosen to be the same as that of the constant-slope chan- nel. As such the merit of the breaking-slope channel over the constant-slope channel can be revealed in an evident fashion, through the comparisons of those properties (i.e. time-averaged wave crest and trough depths as well as signicant wave height). However the horizontal space in the transverse direction has to be utilized for design- 162 ing the breaking-slope channel, otherwise the elevation could not drop suciently in Reach (2). The breaking-slope channel is not necessarily extended in the transverse direction, if Reach (1) and Reach (3) are steeper than the constant-slope channel. Meanwhile the breaking-slope channel is still able to retain the function of suppress- ing growth of the roll-wave amplitude. 163 Chapter 7 Conclusions The subject of this study is to explore the availability of computationally simulat- ing the spatial evolution of roll waves down a wide, rectangular open channel. Major conclusions from the present study are drawn in this chapter. They are divided into four categories and each category is discussed in a single section. Future research considerations are also stated. 7.1 On the Computational Model USC SWAT A computational model is constructed with the purpose of simulating the roll- wave evolution. The equations of motion in this model are the viscous Saint Venant equations which describe the conservation of mass and momentum in open channel water ows. The equations adopt the shallow water theory and seek for the depth- averaged ow velocity throughout cross-section. The external sources in the momen- tum equation are modeled by the Ch ezy formula along with the Colebrook-White equation. The turbulently viscous force is also incorporated into the momentum equation. The equations of motion are solved by a high-resolution, shock-capturing numerical scheme. The scheme is based on the nite volume formulation. The convec- tive uxes on cell interfaces are evaluated using the TVD Lax-Wendro scheme (Yee 1987) in which the eigenvalues and eigenvectors are obtained from Roe's parameter vector approach (1981). The Crank-Nicholson method (1947) is used to cope with the 164 diusive ux. This model is named as \Shallow Water Analyzed by a TVD Scheme at University of Southern California" with an acronym \USC SWAT". The numerical scheme in the computational model has been tested through con- vergence tests concerning progressing roll waves. The time-dependent numerical so- lutions on the maximum amplitude of progressing roll waves are found to match the transient linear theory proposed by Brook et al. (1999) when the wave amplitude is small. The numerical solution of converged progressing roll waves are found to match Dressler's quasi-steady nonlinear theory. The results of the convergence tests corroborate the numerical scheme is capable of providing accurate solutions to the equations of motion. 7.2 On Simulating the Spatial Evolution of Roll Waves The computational model is utilized to carry out simulations of natural roll- wave ows in the real world. The small perturbations on the normal ow variables are imposed at the in ow boundary and the out ow boundary is made thoroughly non-re ecting. The simulation results are compared with Brock's experimental data (1967). Findings in this category are summarized as following: 1. The computational model along with the present boundary conditions can gen- erate simulated natural roll waves that are observed in the laboratory channel and eld channels. 2. The development curves of time-averaged wave crest and trough depths are obtained from the simulated time series of the water depth at various chan- nel locations. The simulation results agree satisfactorily with the experimental 165 data by Brock, for the spatial development curves of both wave crest and trough depths. The computational simulation is able to provide results at more down- stream locations which are not studied in the experiment due to the space limit. 3. The simulation results are sensitive to the size of control cells in a computa- tional domain. The cell size should be diminished to 0.001 m so as to gain de- cent agreements between the simulation results and experimental data. Coarser cell sizes result in deformed wave shapes as well as overestimated wave ampli- tude throughout the whole channel. However the computational time greatly increases when the ne cell size is used. 4. The simulation results are sensitive to the value of the turbulent viscosity. Greater values induce smaller spatial growth of the wave amplitude and smaller amplitude of evolved roll waves. In the comparison between the simulation re- sults and experimental data, the values of the turbulent viscosity reside in the range of [0:26; 0:52] m 2 /s. 7.3 On Investigating Roll-wave Properties The computational simulations are performed to investigate some of important properties of natural roll waves. Findings in this category are summarized as follow- ing: 1. The magnitude of the ow velocity is consistent with that of the water depth. 2. The spatial evolution of natural roll waves obeys a generality in spite of dierent hydraulic conditions. This generality is composed of three phases: initial phase, transition phase, and nal phase. 166 3. The simulated natural roll-wave ows are characterized by three types of wave- wave interactions: wave overtaking, wave absorption, and wave spawning. Only the wave overtaking process is mentioned to be observed in the real world ac- cording to previous literatures. The wave overtaking and absorption processes result in the coalescence of multiple waves. The wave spawning process implies natural roll-wave ows must occur as wave trains without the occurrence of an individual wave. 4. Natural roll waves are nonperiodic. The dominant wave periods tend to increase from upstream to downstream channel locations. Also the nonperiodicity of natural roll waves increases from upstream to downstream channel locations. 7.4 On Mitigating the Spatial Evolution of Roll Waves Two approaches that could help overcome diculties caused by roll waves are considered in this study. How the variation of the bottom roughness height in uences the signicant wave height is investigated. The water ow in a breaking-slope open channel is analyzed, and the comparisons of important properties related to the wave amplitude are conducted between the constant-slope and breaking-slope channels. Findings in this category are summarized as following: 1. The water ow instability can be weakened by using channels with rough bot- tom. 2. Large values of the bottom roughness result in small wave height throughout the channel. Hence, in the engineering application, the uncertainty regarding the water depth prediction is minimized by increasing the value of the bot- 167 tom roughness. Even though the normal water depth is relatively high with large values of the bottom roughness, the water depth can be readily predicted without signicant disturbances on the water surface. 3. There is no hydraulic jump on the transition from the steep reach to the less steep reach in the breaking-slope channel, because the undisturbed Froude num- ber in the less steep reach is 1.53 and the water ow is still supercritical. 4. Shock-like roll waves do not form in the breaking-slope channel. 5. The wave amplitude attenuates spatially in the less steep reach where the water ow is stabilized. 6. When the ow rate and bottom roughness are invariable, the utilization of the breaking-slope channel is able to eectively suppress growth of the wave amplitude compared with the constant-slope channel. Due to the existence of large roll waves in the constant-slope channel, the maximum water depth is greatly decreased by utilizing the breaking-slope channel, even if the normal water depth in the less steep reach is relatively high. 7.5 Future Research Considerations According to the ndings and progress from the present work, following prospec- tive research directions are proposed: 1. Experiments can be carried out to produce and validate a mathematical model, which evaluates the turbulent viscosity in the Saint Venant equations for shallow water ows. 168 2. The Manning's equation can be applied to model the external sources in the Saint Venant equations. The simulation results concerning the spatial evolution of roll waves associated with the Manning's equation can be compared with those from the current equations using the Ch ezy formula. 3. The eect of the channel width change on the spatial evolution of roll waves can be investigated both computationally and experimentally. In the computational study about the channel width change problem, two-dimensional shallow water equations may have to be used, since the water ow in the transverse direction may not be negligible. 169 References [1] Alavian, V. 1986. Behaviour of density currents on an incline. ASCE J. Hydraul. Eng. 112(1):27-42. [2] Balmforth, N.J., Liu, J.J. 2004. Roll waves in mud. J. Fluid Mech. 519:33-54. [3] Balmforth, N.J., Mandre, S. 2004. Dynamics of roll waves. J. Fluid Mech. 514:1- 33. [4] Bhatnagar, P.L., Gross, E.P., Krook, M. 1954. A model for collision processes in gases. I: small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94:511-525. [5] Billingham, J., King, A.C. 2000. Wave Motion. Cambridge University Press, UK. [6] Bogoliubov, N.N., Mitropolski, Y.A. 1961. Asymptotic Methods in the Theory of Nonlinear Oscillations. Hindustan Publishing Corp., India. [7] Boudlal, A., Liapidevskii, V.Y. 2004. Roll waves in a non-regular inclined chan- nel. Euro. Jnl of Applied Mathematics 15:257-271. [8] Brock, R.R. 1967. Development of Roll Waves in Open Channels. Technical re- port No. KH-R-16, California Institute of Technology. [9] Brock, R.R. 1969. Development of roll-wave trains in open channels. ASCE J. Hydraul. Div. 95:1401-1428. [10] Brook, B.S., Falle, S.A., Pedley, T.J. 1999. Numerical solutions for unsteady gravity-driven ows in collapsible tubes: evolution and roll-wave instability of a steady state. J. Fluid Mech. 396:223-256. [11] Burguete, J., Carc a-Navarro, P. 2001. Ecient construction of high-resolution TVD conservative schemes for equations with source terms: application to shal- low water ows. Int. J. Numer. Methods Fluids 37:209-248. [12] Carver, S., Sear, D., Valentine, E. 1994. An observation of roll waves in a supraglacial meltwater channel, Harlech Gletsher, East Greenland. J. Glaciol- ogy 40(134):75-78. [13] Cenedese, C., Whitehead, J.A. 2004. A dense current owing down a sloping bottom in a rotating uid. J. Phys. Oceanogr. 34(1):188-203. [14] Chakravarthy, S.R., Osher, S. 1983. High resolution applications of the Osher upwind scheme for the Euler equations. Proc. AIAA Computational Fluid Dy- namics Conference, Danvers, MA 363-372. 170 [15] Chang, H.C., Demekhin, E.A., Kalaidin, E. 2000. Coherent structures, self- similarity, and universal roll wave coarsening dynamics. Phys. Fluids 12:2268- 2278. [16] Cornish, V. 1934. Ocean Waves and Kindred Geophysical Phenomena. Cam- bridge University Press, London. [17] Coussot, P. 1997. Mud ow Rheology and Dynamics. Balkema, Netherlands. [18] Crank, J., Nicholson, P. 1947. A practical method for numerical evaluation of so- lutions of partial dierential equations of the heat-conduction type. Proc. Camb. Phil. Soc. 43:50-67. [19] Cristo, C.D., Iervolino, M., Vacca, A., Barbar, Z. 2010. In uence of relative roughness and Reynolds number on the roll-waves spatial evolution. ASCE J. Hydraul. Eng. 136:24-33. [20] Delis, A.I., Skeels, C.P. 1998. TVD schemes for open channel ow. Int. J. Numer. Meth. Fluids 26:791-809. [21] Dressler, R.F. 1949. Mathematical solution of the problem of roll-waves in in- clined open channels. Commun. Pure Appl. Math. 2:149-194. [22] Engelund, F., Wan, Z. 1984. Instability of hyperconcentrated ow. ASCE J. Hydraul. Eng. 110(3):219-233. [23] Eymard R, Gallou et T, Herbin R. 2000. Finite volume methods. Handbook of Numerical Analysis 7:713-1020. [24] Falle, S.A.E.G. 1991. Self-similar jets. Mon. Not. R. Astron. Soc. 250:581-596. [25] Fer, I., Lemmin, U., Thorpe, S.A. 2002. Winter cascading of cold water in Lake Geneva. J. Geophys. Res. 107(C6):3060. [26] Fletcher, C.A.J. 1991. Computational Techniques for Fluid Dynamics. 1: Fun- damental and General Techniques. Springer, Germany. [27] Forterre, Y., Pouliquen, O. 2003. Long-surface-wave instability in dense granular ows. J. Fluid Mech. 486:21-50. [28] Gasc on, L., Corber an, J.M. 2001. Construction of second-order TVD schemes for nonhomogeneous hyperbolic conservation laws. J. Comput. Phys. 172:261-297. [29] Godunov, S.K. 1959. A dierence method for numerical calculation of discontin- uous solutions of the equations of hydrodynamics. Mat. Sb. (N.S.) 47:271-306. [30] Harten, A. 1983. High resolution schemes for hyperbolic conservation laws. J. Comput. Phys. 49:357-393. [31] Harten, A., Hyman, J.M. 1983. Self adjusting grid methods for one-dimensional hyperbolic conservation laws. J. Comput. Phys. 50:235-269. 171 [32] Henderson, F.M. 1989. Open Channel Flow. Macmillan, New York. [33] Hwang, S.H., Chang, H.C. 1987. Turbulent and inertial roll waves in inclined lm ow. Phys. Fluids 30:1259-1268. [34] Jereys, H. 1925. The ow of water in an inclined channel of rectangular section. Philosphical Magazine Series 6 49:793-807. [35] Jin, S. 2001. A steady-state capturing method for hyperbolic systems with source terms. Math. Model. Numer. Anal. 35(4):631-645. [36] Jin, S., Kim, Y.J. 2001. On the computation of roll waves. Math. Model. Numer. Anal. 35(3):463-480. [37] Johnson, M., Zumbrun, K., Noble, P. 2011. Nonlinear stability of viscous roll waves. SIAM. J. Math. Anal. 43:577-611. [38] Julien, P.Y., Friesen, N., Duan, J.G., Eykholt, R. 2010. Celerity and amplica- tion of supercritical surface waves. ASCE J. Hydraul. Eng. 136:656-661. [39] Julien, P.Y., Hartley, D.M. 1986. Formation of roll waves in laminar sheet ow. J. Hydraul. Res. 24:5-17. [40] Kranenburg, C. 1992. On the evolution of roll waves. J. Fluid Mech. 245:249-261. [41] Lax, P., Wendro, B. 1960. System of conservation laws. Commun. Pure Appl. Math. 13:217-237. [42] Leonard, B.P. 1979. A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comput. Methods Appl. Mech. Eng. 19:59- 98. [43] Leveque, R.J. 2002. Finite-Volume Methods for Hyperbolic Problems. Cambridge University Press, UK. [44] Lien, F.S., Leschziner, M.A. 1994. Upstream monotonic interpolation for scalar transport with application to complex turbulent ows. Int. J. Numer. Methods Fluids 19:527-548. [45] Liu, K.F., Mei, C.C. 1994. Roll waves on a layer of a muddy uid down a gentle slope a Bingham model. Phys. Fluids 6:2577-2590. [46] Needham, D.J., Merkin, J.H. 1984. On roll waves down an open inclined channel. Proc. R. Soc. Lond. A 394:259-278. [47] Ng, C., Mei, C.C. 1994. Roll waves on a shallow layer of mud modelled as a power-law uid. J. Fluid Mech. 263:151-183. [48] Noble, P. 2007. Linear stability of viscous roll waves. Communications in Partial Dierential Equations 32:1681-1713. [49] Pantan, R.L. 2005. Incompressible Flow. John Wiley & Sons, New Jersey. 172 [50] Pascal, J.P. 2006. Instability of power-law uid ow down a porous incline. J. Non-Newtonian Fluid Mech. 133:109-120. [51] Pascal, J.P., D'Alessio, S.J.D. 2007. Instability of power-law uid ows down an incline subjected to wind stress. Appl. Math. Model. 31(7):1229-1248. [52] Plumerault, L.R., Astruc, D., Thual, O. 2010. High-Reynolds shallow ow over an inclined sinusoidal bottom. Phys. Fluids 22(5):054110. [53] Que, Y.T., Xu, K. 2006. The numerical study of roll-waves in inclined open channels and solitary wave run-up. Int. J. Numer. Meth. Fluids 50:1003-1027. [54] Reynolds, O. 1903. Papers on Mechanical and Physical Subjects. Cambridge Uni- versity Press, London. [55] Roe, P.L. 1981. Approximate Riemann solvers, parameter vectors, and dierence scheme. J. Comput. Phys. 43:357-372. [56] Roe, P.L. 1984. Generalized formulation of TVD Lax-Wendro schemes. Tech- nical report ICASE 84-53, NASA Langley Research Center, Virginia, USA. [57] Roe, P.L. 1986. Characteristic-based schemes for the Euler's equations. Ann. Rev. Fluid Mech. 18:337-365. [58] Rouse, H. 1938. Fluid Mechanics for Hydraulic Engineers. McGraw-Hill, New York. [59] Segal, L.A. 1980. Mathematical Models in Molecular and Cellular Biology. Cam- bridge University Press, London. [60] Sewby, P.K. 1984. High resolution schemes using ux limiters for hyperbolic conservation laws. SIAM J. Numer. Anal. 21:995-1011. [61] Stoker, J.J. 1958. Water Waves: The Mathematical Theory with Applications. John Wiley & Sons, New York. [62] Sturm, T.W. 2001. Open Channel Hydraulics. McGraw-Hill, New York. [63] Swaters, G.E. 2003. Baroclinic characteristics of frictionally destablized abyssal over ows. J. Fluid Mech. 498:349-379. [64] Synolakis, C.E. 1987. The runup of solitary waves. J. Fluid Mech. 185:523-545. [65] Thomas, H.A. 1937. The propagation of stable wave congurations in steep chan- nels. Technical report, Carnegie Institute of Technology, Pittsburgh, PA. [66] Toro, E.F. 2009. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer, Germany. [67] van Leer, B. 1974. Towards the ultimate conservative dierence scheme II. Mono- tinicity and conservation combined in a second-order scheme. J. Comput. Phys. 14:361-370. 173 [68] van Leer, B. 1977. Towards the ultimate conservative dierence scheme. IV: a new approach to numerical convection. J. Comput. Phys. 23:276-299. [69] Versteeg H.K., Malalasekera W. 2007. An Introduction to Computational Fluid Dynamics. Pearson Education Limited, England. [70] Whitham, G.B. 1974. Linear and Nonlinear Waves. Wiley-interscience, New York. [71] Wilcox, D.C. 2006. Turbulence Modeling for CFD. DCW Industries, California. [72] Woods, B.D., Hurlburt, E.T., Hanratty, T.J. 2000. Mechanism of slug formation in downwardly inclined pipes. Intl. J. Multiphase Flow 26(6):977-998. [73] Yee, H.C. 1987. Construction of explicit and implicit symmetric TVD schemes and their applications. J. Comput. Phys. 68:151-179. [74] Yen, B.C. 2002. Open channel ow resistance. ASCE J. Hydraul. Eng. 128:20-39. [75] Yu, J., Keorkian, J. 1992. Nonlinear evolution of small disturbances into roll waves in an inclined open channel. J. Fluid Mech. 243:575-594. [76] Zanuttigh, B., Lamberti, A. 2002. Roll waves simulation using shallow water equations and weighted average ux method. J. Hydraul. Res. 40:610-622. 174 Appendix A Dressler's Mathematical Solution of Roll Waves Dressler (1949) constructed a mathematical solution of progressing roll waves which has been commonly cited and used for testing the accuracy of a numerical scheme solving Saint Venant equations. The comprehensive derivation of the quasi- steady nonlinear theory is specied by Dressler (1949), in this appendix we brie y summarize how the solution is achieved. The fundamental theory from which the mathematical roll-wave solution origi- nates is the Saint Venant equations: @h @t +h @u @x +u @h @x = 0; (A.1) @u @t +u @u @x +g cos @h @x =g sinC d u 2 h : (A.2) When progressing roll waves travel down an open channel with a constant speed of c w , a single progressing coordinate can be established to incorporate the space and time, which is expressed as =xc w t: (A.3) So we can deduce the following relation u (x;t) =U (xc w t) =U (); (A.4) h (x;t) =H (xc w t) =H () (A.5) 175 in which H is the water depth and U is the ow velocity in the progressing coor- dinate. Having the progressing coordinate and using the substitution of @h/@t = c w (@H/@t), @h/@x =@H/@x and likewise for u, equation 4.1 and 4.2 become H dU d + (Uc w ) dH d = 0; (A.6) (Uc w ) dU d +g cos dH d =g sinC d U 2 H : (A.7) Solving these two equations that govern ow variables in the moving system gives dU d = (Uc w ) (g sinC d U 2 /H) (Uc w ) 2 g cosH ; (A.8) dH d = H (g sinC d U 2 /H) (Uc w ) 2 g cosH : (A.9) Dividing above equations leads to dH dU = H c w U (A.10) which integrates to (c w U)HK: (A.11) The constant, K, is treated as the progressing discharge rate when viewed from the moving coordinate. From equation A.11, we know U = (c w HK)/H that is substituted into equation A.9 to generate dH d = H g sinC d (c w HK) 2 H 3 (Uc w ) 2 g cosH : (A.12) 176 This equation is the dierential prole equation valid for progressing waves. It can be integrated to determine what possible water depth proles are obtainable, but none of them are periodic and capable of giving roll-wave solutions. Nevertheless Dressler devised a special solution that is constructed by joining together sections of a continuous wave prole through a series of discontinuous shocks. In order to do this, two steps must be implemented: rst, a continuous water depth prole must be found; second, the locations of the shocks or the maximum and minimum water depth at the shocks must be found. When U is the ow velocity and c w is the wave speed, we can consider Uc w as the ow velocity relative to the travelling waves. Since c w is everywhere greater than U, the relative ow velocity Uc w is negative denite. When an open channel ow is critical, we have (U c c w ) 2 =g cosH c (A.13) attributed to F = 1 where the subscript \c" represents quantities evaluated when the relative ow is critical. It can be noted that the denominator on the right hand side of equation A.12 vanishes if the relative ow is critical. Dressler corroborated that roll-wave solutions can be formed only from the exceptional case in which the numerator is also required to be zero. Therefore, by setting both denominator and numerator to be zero we are able to quantify the critical parameters: U c = c w 1 + p C d /sin ; (A.14) H c = 1 g cos c w 1 + p sin/C d ! 2 ; (A.15) 177 and K c = 1 g cos c w 1 + p sin/C d ! 3 : (A.16) With the aid of obtained critical values, equation A.12 can be written as dH d = sin H 2 + [H 0 (c 2 w C d )/(g sin)]H +C d H 2 0 /sin H 2 +H 0 H +H 2 0 (A.17) which is a modied dierential form of the wave prole. Through integration on the above equation, we can obtain the continuous roll-wave prole associated with H (0) =H c (i.e. the relative ow is critical at = 0): = 1 sin (HH c ) + H 2 A +H A H c +H 2 c H A H B ln HH A H c H A H 2 B +H B H c +H 2 c H A H B ln HH B H c H B (A.18) where H A;B = p C d 2 sin ( p C d + 2 p sin p C d + 2 p sin 2 4 sin 1/2 ) H c (A.19) are designated to be two roots of the numerator on the right hand side of equa- tion A.17. If we name the wave described by equation A.18 as H n , then the next wave at the downstream H n+1 is expressed as =L + 1 sin (HH c ) + H 2 A +H A H c +H 2 c H A H B ln HH A H c H A H 2 B +H B H c +H 2 c H A H B ln HH B H c H B (A.20) whereL denotes the wavelength. Wave proles represented by equation A.18 and A.20 178 H n H 1 n H + b H H c H 0 ζ = S ζ ζ Back Front L f H B H A H c H Figure A.1. Sketch illustration of two consecutive progressing roll-wave proles de- vised by Dressler (1949). The vertical axis is the water depth while the horizontal axis is the progressing coordinate. Two roll waves H n and H n+1 that have identical prole and wavelength are separated by a discontinuous shock at the location s . The back of the shock H b represents the maximum water depth at the crest while the shock front H f represents the minimum water depth at the trough. are illustrated by Figure A.1. It should be emphasized that the wave prole given by equation A.18 is innitely long, if the wave trough and crest depths are not determined. According to conserva- tion laws applied at hydraulic shocks, the shock front and back possess the relation of H f = 1 2 " H 2 b + 8K 2 c g cosH b 1/2 H b # (A.21) in which H b is the shock back and H f is the shock front. The shock front and back correspond to the trough and crest point respectively. Moreover the back of the wave H n and the front of the next wave H n+1 share the same shock location s . Hence we 179 can utilize this characteristic to derive the other relationship between H b and H f as s =H n (H b ) =H n+1 (H f ): (A.22) After substitution of the prole equation A.18, the above relation becomes H 2 A +H A H c +H 2 c H A H B ln H b H A H f H A H 2 B +H B H c +H 2 c H A H B ln H b H B H f H B sinL +H b H f = 0: (A.23) Therefore the values of H b and H f are obtained by solving equation A.21 and A.23 iteratively. In conclusion, Dressler's mathematical solution of progressing roll waves can be es- tablished if four parameters are prescribed: the angle of channel inclination, bottom drag coecientC d , progressing wave speed c w and wavelength L. The quantication of these parameters should be correlated to demands of the comparison between an- alytical and numerical solutions. The values used in Dressler's theory should be consistent with what are used in the numerical experiment of progressing roll waves. Nevertheless the wave speed is not specied in the numerical experiment, so it pro- duces some more complication. To measure the wave speed, we record a computed time series of progressing roll waves at any location inside the computational domain and acquire its wave period T using the fast Fourier transform. As such the wave speed is calculated by c w = L/T . Physically the mathematical solution represents fully developed roll waves that do not alter in size and shape but merely progress forward. If viewed from the progressing coordinate , those waves even remain at 180 the same location. Thus Dressler's solution is known as the quasi-steady nonlinear theory of progressing roll waves. 181 Appendix B Derivation of the Saint Venant Equations In order to construct a computational model for simulating the spatial evolu- tion of roll waves, we need to have a set of equations which governs water ows in a wide, rectangular open channel. In this appendix, we derive the equations of mo- tion providing the quantication of ow variables including the water depth and ow velocity as functions of space and time. The equations are a system of partial dier- ential equations describing the conservation of mass and momentum inside a specic control volume. Before deriving the equations of motion, we rst recall the Reynolds transport theorem (Reynolds 1903): the rate of change of an extensive property of a system is equal to the summation of the rate of change of the property in a control volume and the net out ow rate of the property through the surface of the control volume. The theorem can be written into a formula as dB system dt = d dt Z c.v. bdV + Z c.s. budA (B.1) where B system stands for the property of a system and b is the property in a control volume. t is the time,V is the volume, andA denotes the area of the control surface. is the density and u represents the velocity. c.v. and c.s. are the acronyms for control volume and control surface, respectively. The Reynolds transport theorem connects the Lagrangian and Eulerian approaches for analyzing uid problems. It is 182 the fundamental concept for using the control volume approach. Since we consider water ows in a wide, rectangular open channel, the problem is two-dimensional. The ow conditions almost do not change at dierent proles. The ow velocity in the meridional direction is small compared to the zonal direction. Therefore we assume the selected control volume has a unity width in the meridional direction. Turbulent ows are characterized by uctuations of ow variables. The ow is unsteady even under the steady-uniform condition because of the uctuations. The ow variables, water depthh and ow velocityu, mentioned below are treated to be time-averaged values such that the contributions from the turbulent uctuations are disregarded. B.1 Conservation of Mass The conservation of mass for a control volume is illustrated by Figure B.1. The space enclosed by the dash line is the selected control volume with a length of x and height of h, letting x be the spatial coordinate along the channel bed and h be the water depth. The locations of the upstream and downstream control surfaces are x 1 and x 2 . Note the lateral ow is not considered in the current study. The extensive property is the mass in this case. We assume the mass of the entire system is steady, so dm system /dt = 0 where m system is the mass of the system. The property in the control volume b on the right hand side of equation B.1 is evaluated as the mass per unit mass. That is b = m system m system = 1: (B.2) 183 Water flow h 1 q 2 q z x Δ 1 x 2 x x Figure B.1. Schematic for the conservation of mass in the selected control volume for water ows down a wide, rectangular open channel. The control volume is indicated by the dash line. x andz are the spatial coordinates. x is the length of the control volume and the water depthh is the height of the control volume. The control volume has a unity width. x 1 andx 2 represent the locations of the upstream and downstream control surfaces. q 1 and q 2 stand for the volume uxes through the control surfaces. The density of water is constant because the open channel ow is incompressible. Therefore the Reynolds transport theorem equation B.1 reduces to d dt Z c.v. dV + Z c.s. udA = 0: (B.3) Since the control volume has a unity width, the total volume of the control volume is evaluated as Z c.v. dV = Z x 2 x 1 hdx: (B.4) 184 The net volume eux through the control surfaces can be written as Z c.s. udA = Z h 2 0 u 2 dz Z h 1 0 u 1 dz = Z h 0 udz: (B.5) h 1 andh 2 are the water depth on the control surfaces. u 1 andu 2 are the ow velocity on the control surfaces. z is the vertical coordinate. Substituting equations B.4 and B.5 into equation B.3 yields d dt Z x 2 x 1 hdx + Z h 0 udz = 0: (B.6) We take the length of the control volume to be innitely small: lim x!0 x =dx, such that the water depth is constant inside the control volume. To derive the equation of the mass conservation, we need to adopt an approximation on the ow velocity. The approximation will be discussed in Section B.3. Approximation B.1 The ow velocity is depth-averaged and uniform throughout cross-section. Hence equation B.6 becomes dh dt dx +d (hu) = 0: (B.7) Dividing the above equation by the length of the control volume and using the partial dierentiation operator, we achieve the conservation of mass in a dierential form @h @t + @ @x (hu) = 0: (B.8) 185 B.2 Conservation of Momentum The Newton's second law states that the particle acceleration is proportional to the summation of all forces on that particle in a specic direction. The formula of the law is written as X F =ma (B.9) where F stands for the force, m is the mass, and a denotes the acceleration. Letting the mass be constant, the right hand side of the formula can be interpreted as the time rate of change of momentum. Hence the formula can be rewritten as X F =m du dt = d dt (mu) = dM dt (B.10) where M is the momentum. The law can also be formulated for an entire system to give X F = dM system dt (B.11) where M system is the momentum of the system. In order to apply the Reynolds transport theorem, we let the extensive property of the system be the momentum. The corresponding intensive property b becomes the momentum per unit mass, that is b = M m =u: (B.12) 186 Therefore the Reynolds transport theorem for the momentum is dM system dt = d dt Z c:v: udV + Z c:s: u 2 dA: (B.13) Combining equation B.11 and B.13, we are able to deduce the integral form of the momentum equation as X F = d dt Z c:v: udV + Z c:s: u 2 dA: (B.14) This equation dictates that, for a specic control volume, the summation of all forces is equal to the summation of the rate of change of momentum and the net momentum eux through the control surfaces. The conservation of momentum for a control volume specied by equation B.14 is shown in Figure B.2. The control volume is chosen for water ows in a wide, rectangular open channel. The dimension of the control volume is similar to that of mass conservation. The momentum inside the control volume is evaluated as Z c:v: udV = Z x 2 x 1 Z h 0 udzdx: (B.15) The net momentum eux through the control surfaces is Z c:s: u 2 dA = Z h 2 0 u 2 2 dz Z h 1 0 u 2 1 dz = Z h 0 u 2 dz: (B.16) In open channel ows, forces for the system as well as acting on the control volume along the channel bed direction have four components: the force due to the 187 Water flow h x z ( ) 1 p F ( ) 2 p F ( ) 1 xx F τ sin W θ ( ) 1 qu ρ ( ) 2 qu ρ x Δ 1 x 2 x x b F τ ( ) 2 xx F τ W Figure B.2. Schematic for the conservation of momentum in the selected control volume for water ows down a wide, rectangular open channel. The control volume is indicated by the dash line. x andz are the spatial coordinates. x is the length of the control volume and the water depth h is the height of the control volume. The control volume has a unity width. x 1 and x 2 represent the locations of the upstream and downstream control surfaces. (qu) 1 and (qu) 2 stand for the momentum uxes through the control surfaces. F p is the force due to the water pressure. F xx is the viscous force pertinent to the normal stress. F b is the frictional force due to the bottom drag. W denotes the weight of the control volume. is the angle of channel inclination. water pressureF p , the viscous force pertinent to the normal stressF xx , the frictional force due to the bottom shear stressF b , and the gravitational forceW sin. W is the weight of the control volume and is the angle of channel inclination. The shear stress among water particles is neglected, because the strain rate @u/@z is zero resulting from Approximation B.1. Also the rate of strain of the vertical velocity is negligible. The pressure force and viscous force act on the upstream and downstream control surfaces. The frictional force acts on the bottom of the control volume between the 188 water boundary layer and channel bed. The gravitational force is the body force. In an equation, the total force on the control volume is the summation of those four components: X F = X F p + X F xx F b +W sin: (B.17) The minus sign at the front of the frictional forceF b indicates it acts on the negative x direction. The aggregate pressure force is expressed as X F p = Z h 1 0 p 1 dz Z h 2 0 p 2 dz = Z h 0 pdz: (B.18) p denotes the water pressure. p 1 and p 2 are the water pressure on the upstream and downstream surfaces. Here we need to adopt the other approximation on the water pressure. The rantionality of the approximation is discussed in Section B.3. Approximation B.2 The water pressure obeys a hydrostatic distribution, namely @p @z =g cos (B.19) where g is the gravitational acceleration. Multiplying the above hydrostatic law by @z on both sides and integrating from z to h generates Z h z @p = Z h z g cos@z: (B.20) 189 After evaluating the integrals, we have (p) h (p) z =g cos (hz): (B.21) (p) h is the water pressure on the top surface, which is zero. (p) z is the pressure we are seeking at some vertical locationz. Hence the water pressure is calculated by the following: p =g cos (hz): (B.22) By substituting equation B.22 into equation B.18, we obtain the total pressure force as X F p = Z h 0 g cos (hz)dz = 1 2 g cosh 2 : (B.23) The aggregate viscous force is X F xx = Z h 2 0 ( xx ) 2 dz Z h 1 0 ( xx ) 1 dz = Z h 0 xx dz (B.24) where xx represents the normal stress on the vertical control surfaces. ( xx ) 1 and ( xx ) 2 are the normal stress on the upstream and downstream surfaces. We utilize the common approach for dening the normal stress of uids of incompressible ow: two times the product of the dynamic viscosity and strain rate. The water ows in an open channel are turbulent ows. Thus the kinematic viscosity is replaced by a turbulent viscosity. The normal stress is written as xx = 2 t 0 @u @x = t @u @x : (B.25) 190 t is the constant turbulent viscosity and t = 2 t 0 . After substitution from equa- tion B.25, the aggregate turbulently viscous force equation B.24 becomes X F xx = Z h 0 t @u @x dz: (B.26) The gravitational force leading to the uid motion along the channel bed is written as W sin = Z x 2 x 1 Z h 0 g sindzdx: (B.27) The frictional force due to the bottom drag is expressed as F b = Z x 2 x 1 b dx (B.28) where b stands for the shear stress between the water boundary layer and channel bottom. The four components forming the total force on the control volume are specied by equation B.23, B.26, B.27 and B.28. Substituting those equations into equation B.17 produces X F = 1 2 g cosh 2 + Z h 0 t @u @x dz Z x 2 x 1 b dx + Z x 2 x 1 Z h 0 g sindzdx: (B.29) The momentum inside the control volume is specied by equation B.15. The net momentum eux through the control surfaces is described by equation B.16. The summation of all forces on the control volume is specied by equation B.29. Those 191 three equations are substituted into the Reynolds transport theorem for momentum which is equation B.14 to yield d dt Z x 2 x 1 Z h 0 udzdx + Z h 0 u 2 dz = 1 2 g cosh 2 + Z h 0 t @u @x dz Z x 2 x 1 b dx + Z x 2 x 1 Z h 0 g sindzdx: (B.30) The water density , gravitational acceleration g, and the channel inclination are constant. Similar to the conservation of mass, the length of the control volume x is set to be innitely small. Therefore the ow variables and bottom shear stress is uniform inside the control volume. Approximation B.1 is also adopted for the conservation of momentum. Dividing equation B.30 by and the length of the control volume and using the partial dierentiation operator, we can derive @ @t (hu) + @ @x hu 2 + 1 2 g cosh 2 @ @x t h @u @x =g sinh b : (B.31) We utilize the Ch ezy formula to model the bottom shear stress. The Ch ezy formula species the relationship between the ow velocity and hydraulic radius in steady-uniform open channel ows. It is written as u 0 =C c p R h tan (B.32) where u 0 is the normal ow velocity, R h denotes the hydraulic radius, and C c is the Ch ezy coecient. We can use the Darcy-Weisbach friction factor f to evaluate the 192 Ch ezy coecient: C c = s 8g cos f : (B.33) For water ows in a wide, rectangular open channel, the hydraulic radius closely approximates to the water depth. Hence the Ch ezy formula can be rewritten as u 0 = r 8 f g sinh 0 (B.34) or f 8 u 2 0 =g sinh 0 (B.35) in which h 0 is the normal water depth. On the basis of equation B.35, the bottom shear stress is modeled as b =C d ujuj (B.36) where C d represents the dimensionless drag coecient and C d = f 8 : (B.37) The bottom drag as the ow resistance always acts in the opposite direction of the ow velocity, and this is why one of u has an absolute value sign in equation B.36. After substitution from equation B.36, the momentum conservation equation B.31 becomes @ @t (hu) + @ @x hu 2 + 1 2 g cosh 2 @ @x t h @u @x =g sinhC d ujuj: (B.38) 193 Equation B.8 along with equation B.38 are the conservative form of the Saint Venant equations. We can make further simplications on the momentum equation. From the conservation of mass equation B.8, we know @h @t =h @u @x u @h @x : (B.39) Equation B.38 can be rewritten as h @u @t +u @h @t +2hu @u @x +u 2 @h @x +g cosh @h @x @ @x t h @u @x =g sinhC d ujuj: (B.40) After substitution from equation B.39, equation B.40 is divided by the water depth h to give @u @t +u @u @x +g cos @h @x 1 h @ @x t h @u @x =g sinC d ujuj h : (B.41) Equation B.8 and equation B.41 are the pair of the equations of motion to be numeri- cally solved in the present study. They are equivalent to equation 3.1 and equation 3.6. B.3 Disscusion on Approximation B.1 and B.2 Two approximations are adopted in the derivation of the Saint Venant equations. By adopting Approximation B.1, we are seeking a single value of the depth-averaged ow velocity at each cross-section. The real velocity prole is disregarded. Indepen- dent experiments about the turbulent boundary layer have shown the real velocity prole in channel ows obeys the law of the wall and the velocity defect law (Pantan 194 2005). However the ow velocity almost increases to the surface velocity very near the boundary wall. Approximation B.2 is inferred from the assumption that the vertical component of the particle acceleration has a negligible eect on the water pressure. The vertical component of the Navier-Stokes equations for open channel ows is written as Dw Dt = 1 @p @z +r 2 w +g cos (B.42) wherew is the vertical velocity and is the kinematic viscosity. The local acceleration @w/@t and the rate of strainrw are assumed to be zero. Hence the material deriva- tive, representing the total acceleration, and stresses are omitted from the equation. This is how equation B.19 comes from. Those approximations can be made in open channel ows because the horizontal dimension is much larger than the vertical scale. Therefore they are commonly known as the shallow water theory. In open channel ows, the wavelength and channel width are much larger than the water depth. The Saint Venant equations are capable of describing other shallow water ows, such as tidal bore propagation and wave run-up on a sloping beach.

Abstract (if available)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Wave induced hydrodynamic complexity and transport in the nearshore

PDF

The development of a hydraulic-control wave-maker (HCW) for the study of non-linear surface waves

PDF

Numerical and experimental study on dynamics of unsteady pipe flow involving backflow prevention assemblies

PDF

Modeling and analysis of parallel and spatially-evolving wall-bounded shear flows

PDF

Computational fluid dynamic analysis of highway bridge superstructures exposed to hurricane waves

PDF

Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders

PDF

On the simulation of stratified turbulent flows

PDF

Tsunami-induced turbulent coherent structures

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Understanding the transport potential of nearshore tsunami currents

PDF

Feasibility of using a lined portion of the Los Angeles River for artificial recharge

PDF

Flow and thermal transport at porous interfaces

PDF

Inverse modeling and uncertainty quantification of nonlinear flow in porous media models

PDF

An experimental study of shock wave attenuation

PDF

Numerical analysis of harbor oscillation under effect of fluctuating tidal level and varying harbor layout

PDF

Tracking and evolution of compressible turbulent flow structures

PDF

Modeling and simulation of complex recovery processes

PDF

Reconstruction and estimation of shear flows using physics-based and data-driven models

PDF

Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles

PDF

Wave method for structural system identification and health monitoring of buildings based on layered shear beam model

Asset Metadata

Creator Huang, Ziyi (author)

Core Title Open channel flow instabilities: modeling the spatial evolution of roll waves

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering (Environmental Engineering)

Publication Date 09/10/2013

Defense Date 08/14/2013

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag breaking-slope channel,modeling roll waves,OAI-PMH Harvest,open channel flow,Saint Venant equations,shallow water flow,water flow instability

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Lee, Jiin-Jen (committee chair)

Creator Email huangziyi2013@foxmail.com,ziyihuang2011@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-327707

Unique identifier UC11294071

Identifier etd-HuangZiyi-2030.pdf (filename),usctheses-c3-327707 (legacy record id)

Legacy Identifier etd-HuangZiyi-2030.pdf

Dmrecord 327707

Document Type Dissertation

Format application/pdf (imt)

Rights Huang, Ziyi

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA