• Categories
• Análisis numérico
• Ecuaciones
• Dinámica de fluidos
• Ecuaciones diferenciales
• Dinámica de fluidos computacional

Citation preview

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

Lecture notes

WL |

delft hydraulics

Computational Hydraulics Version 2006-1

Copyright © 2005, Adri Verwey The right of Adri Verwey to be identified as the author of this work has been asserted in accordance with the Copyright, Designs and Patents Act 1988

Lecture notes

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Contents

1

2

WL | Delft Hydraulics

Introduction...............................................................................................................7 1.1

Introductory remarks......................................................................................7

1.2

Unsteady flow hydraulics ..............................................................................8

1.3

Numerical methods......................................................................................11

1.4

Mathematical modelling model ...................................................................11

Introduction to numerical methods.......................................................................17 2.1

Introduction..................................................................................................17

2.2

Introduction to finite difference approximations .........................................18

2.3

The Euler scheme ........................................................................................18

2.4

The improved Euler scheme ........................................................................21

2.5

The implicit scheme.....................................................................................22

2.6

The Newton-Raphson scheme .....................................................................22

2.7

A more formal analysis of accuracy ............................................................23

2.8

Stability........................................................................................................26

2.9

Consistency and convergence ......................................................................27

2.10

The choice of a time step .............................................................................27

2.11

Reservoir routing .........................................................................................29

2.12

Routing of a decayable pollutant through a reservoir..................................33

2.13

Non-uniform steady flow in channels..........................................................34

2.14

An inverse scheme for the backwater curve computation ...........................38

2.15

Non-uniform flow in non-uniform channels................................................40

2.16

What have we learnt?...................................................................................41

2.17

Questions and small assignments ................................................................41

3

Computational Hydraulics Version 2006-1

3

4

5

WL | Delft Hydraulics

© 2005 Adri Verwey

UNESCO - IHE

Basic Unsteady Channel Flow Equations .............................................................43 3.1

Introduction..................................................................................................43

3.2

Continuity equation .....................................................................................43

3.3

Momentum equation....................................................................................44

3.4

Transformation to the characteristic form....................................................45

3.5

The significance of the characteristics.........................................................48

3.6

The method of characteristics ......................................................................51

3.7

More complex boundary conditions ............................................................55

3.8

The formation of positive and negative hydraulic jumps ............................57

3.9

The limited practical importance of the method of characteristics..............59

3.10

Questions and assignments ..........................................................................60

Introducing Numerical Solutions for Partial Differential Equations.................61 4.1

Introduction..................................................................................................61

4.2

The advection equation................................................................................63

4.3

The characteristic solution ...........................................................................64

4.4

Finite difference schemes ............................................................................65

4.5

Characteristic solutions on a fixed grid .......................................................70

4.6

Introducing diffusion ...................................................................................72

4.7

An explicit finite difference scheme ............................................................73

4.8

Explicit schemes for the combination of advection and diffusion...............76

4.9

Implicit schemes for the diffusion equation.................................................77

De Saint Venant equations and their solutions.....................................................81 5.1

Introduction..................................................................................................81

5.2

The continuity equation ...............................................................................82

5.3

The momentum equation .............................................................................84

5.4

Numerical solutions .....................................................................................88

4

Computational Hydraulics Version 2006-1

6

7

8

WL | Delft Hydraulics

© 2005 Adri Verwey

UNESCO - IHE

5.5

Description of hydraulic structures..............................................................92

5.6

Topological model schematisation...............................................................94

5.7

Hydraulic model schematisation..................................................................95

5.8

Boundary- and initial conditions..................................................................97

5.9

Model calibration and validation .................................................................99

Mathematical modelling of floods .......................................................................101 6.1

Introduction................................................................................................101

6.2

Flood model requirements .........................................................................101

6.3

The role of new data collection technologies ............................................103

6.4

The nature of flood wave propagation .......................................................104

6.5

Deformation of flood waves ......................................................................105 6.5.1

The role of varying celerities ........................................................105

6.5.2

The role of the diffusion term .......................................................106

6.5.3

The role of the lateral flow terms..................................................107

6.6

Link to hydrologic flood routing models ...................................................107

6.7

Two-dimensional modelling of floods .......................................................108

6.8

Integrated 1D/2D modelling ......................................................................113

6.9

Exercise......................................................................................................116

Water Hammer .....................................................................................................117 7.1

Introduction................................................................................................117

7.2

Water Hammer Equations ..........................................................................119

7.3

The Method of Characteristics for Water Hammer....................................125

7.4

Exercise......................................................................................................128

References..............................................................................................................129

5

Computational Hydraulics Version 2006-1

WL | Delft Hydraulics

© 2005 Adri Verwey

UNESCO - IHE

6

Computational Hydraulics Version 2006-1

1 1.1

© 2005 Adri Verwey

UNESCO - IHE

Introduction Introductory remarks

The topic of computational hydraulics is dealing with the question how water resources engineers and planners can be assisted in dealing with complex hydraulic problems. When these problems are very local, they can generally be addressed by using empirical relationships and a well trained engineer usually is equipped with methods and tools dealing with it. In these cases, the use of laws and relationships in hydraulics is often limited to steady flow approximations. For larger scale problems, however, the unsteady nature of flows becomes more dominant and methods used will become more complex. Whereas until various decades ago, the focus has been on developing approximations, partly empirical, of the full hydrodynamic behaviour of the water system, the use of computers has made it possible to describe hydraulic systems quite accurately. Over the past decades enormous progress has been made in developing simulators or mathematical models for all kinds of hydraulic systems, such as rivers, drainage networks, irrigation networks, water distribution networks etc. Currently, the level of accuracy of such simulators is primarily limited by the quality of data available to construct and calibrate the models. A variety of good software packages is available to construct such models. However, good schematisations of hydraulic systems in models requires some insight in the laws and techniques behind these modelling systems in order to use them correctly in building one’s own model. In this series of lectures we address this need by providing insight into the nature of one-dimensional unsteady flow. After the introduction of the unsteady flow equations, the link between the equations and the physical system is shown by the characteristic celerities of disturbances propagating along channels. This provides the basis for the numerical schemes developed to solve the equations. Moreover, it shows clearly the effect of boundary variations and, in particular, the effects of control of hydraulic systems. With this understanding as a basis, numerical method is introduced. First, only the so-called ordinary differential equations are treated, enabling us to do simple backwater computations, water quality simulation in well mixed reservoirs etc. However, at the same time numerical concepts and their evaluation are introduced relating to the accuracy, stability, robustness and efficiency of numerical operations. This will serve as a necessary basis for those who want to develop their own models, as well as for those applying existing modelling systems.

WL | Delft Hydraulics

7

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Subsequently, numerical descriptions of more complex problems described by partial differential equations are discussed, primarily to provide enough insight as a basis for good mathematical modelling practice. Also here, the role and choice of numerical parameters in mathematical model simulation is emphasised. Finally, the overall construction of good mathematical models is discussed. Based upon a clear idea about the objective of model use, rules for the most economical and correct schematisation are developed, starting with a topological schematisation of the system. In addition, the derivation of correct physical parameters is discussed, such as those describing channel conveyance, storage and lateral flows. Moreover, the possibilities and limitations of such models is discussed, with special attention given to the use of models under extrapolated conditions, such as often the case in the modelling of floods.

1.2

Unsteady flow hydraulics

In teaching, the topic of hydraulics is generally first introduced from the steady flow concept. In principle, it would be more logical to start with the unsteady flow concept and than introduce steady flow as a special case of unsteady flow. In this case a clear indication has to be made of the underlying assumptions and the simplification of the problem. In each application these assumptions have to be verified. It should be realised that a real steady flow does not exist in nature. There will always be some small variations in the flow distribution, even if there are no observable water level variations. The steady flow is a concept of our engineering mind, which sometimes can simplify engineering without unacceptable differences from reality. However, the danger exists that engineers turn too easily towards the concept of steady flow, even in cases where this is a quite wrong schematization of reality and where this approach may lead to quite wrong conclusions. For this reason, it is important to familiarize oneself with problems which typically show more significant variations both in time and in space. Hydraulic problems in channels are governed by two important concepts: storage and conveyance. For incompressible flow, the first concept deals with water volume conservation. The second concept deals with balance of forces acting upon the water mass and their effect on the momentum balance. Of particular importance in this description is the magnitude of flow momentum losses due to channel friction relative to the gain of momentum due to gravity or other forces. Before we can discuss more in detail the difference between steady and unsteady flow, the concept of boundary conditions has to be introduced. In general, one is only interested in a specific part of a hydraulic system. To illustrate this idea, it is useful to consider the hydrological cycle. Water evaporates from the sea surface, precipitates partly above land and flows via the land surface, or via infiltration through the subsurface, to the rivers and, in most cases, back to the sea. Let us now consider the river part of this cycle, or even a small part of this river system. The link to the upstream part of this river subsystem is specified in the form of a boundary condition and, more specifically, as the upstream boundary condition.

WL | Delft Hydraulics

8

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

The link to the lower part of the river system, or to the sea, is specified as the downstream boundary condition. In the sequel, we will discuss these boundary conditions in more detail, including the question of the real need for such conditions in our computations. The question of the flow type is very much linked to the nature of the changes at the boundaries of the area described (or modelled). When the part of the river system studied adapts its state only slowly to the changes at the boundaries, the system is considered to be unsteady. However, when the river system adapts its state nearly instantaneously to the changes at the boundaries, the system is said to be quasisteady and some steady flow concepts might be used. This is usually the case when local or near-field problems are studied, such as local structures with their typical backwater effects. One of the important parameters influencing this ease to adaptation is the storage in the system. If this storage capacity is large compared to the difference between inflow and outflow the time scale of adaptation is also large. It is very likely, then, that we will observe a strongly unsteady flow phenomenon. If, however, there is little storage capacity available between the boundaries, the adaptation of the state of the system to the new boundary conditions may be fast and the flow may pass through a sequence of nearly-steady states. This adaptation is also dependent on the facility of the flow to accelerate or decelerate. If the adaptation of flow particle velocities to changing boundary conditions is fast, the system will pass through a series of nearly-steady states. Where the adaptation is slow compared to the changes at the boundaries, the result will be clearly an unsteady flow. These concepts are best illustrated by giving some practical examples. Flood waves in rivers

Floods in a river are the result of the surface- and subsurface runoff generated during periods of intense rainfall. This response usually leads to a typical hydrograph shape discharge and water level variation in the river, which is more pronounced when the rain is uniformly distributed over the period of the rain event. However, as rainfall is generally not uniformly distributed in time, the discharge distribution from the catchments into the river is usually less irregular. While the flood wave propagates down the river it undergoes further deformations due to varying storage and conveyance characteristics of the river channel and due to additional lateral inflows from other catchments. These processes together form a typical unsteady flow phenomenon, which can only be studied by simulations on the basis of unsteady flow equations. Flow in an urban drainage or combined sewer system

Drainage from urban catchments is to be seen as a special case of flood routing. Differences with a rural catchment are the faster response to the rainfall, the smaller amount of storage usually available and the more important role of channel conveyance, compared to channel or reservoir storage. In systems with drainage pipes, the open channel flow may temporarily change into a pressurized flow. The unsteadiness of the phenomenon, however, is very pronounced.

WL | Delft Hydraulics

9

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Despite these differences, it will be shown in the sequel that these problems require the same modelling tools as the ones used for the study of flood routing in rivers. Flow in a desert wadi

In wadis, the unsteady nature of the flow is even more pronounced than in many other rivers, due to the infiltration of the water into the river and flood plain bed. This leads to the steepening up of the flood wave front. Also, if often leads to the complete disappearance of the flood wave after some time. Flow in irrigation systems

Flow in irrigation systems is usually controlled for the optimum use of the water resources. This control introduces variations in time. A fast control may even lead to the formation of hydraulic jumps travelling along the channel. Although the design of the irrigation canals is often based on steady flow concepts, the performance under operation usually requires checks on the basis of unsteady flow computations. Flow over or through a hydraulic structure

The description of flow through a hydraulic structure within a river branch is usually based on an assumption of steady flow. The discharge at each moment in time is directly dependent on the water level boundary conditions. Although these water levels may vary rather fast, the discharge generated will respond more or less instantaneously. The immediate adaptation of the flow to changing boundary conditions is the result of the lack of storage between the upstream and downstream section and the presence of a relatively small water mass to be accelerated or decelerated. As will be shown in the sequel, the possibility to link a steady flow channel element to channel elements where the flow has a distinctly unsteady nature, depends on the relative importance of the various terms of the equations describing the flow problem. Flow in pipe networks

Flow in pipe networks for water distribution in a town is often computed as steady flow. For given constant water demands at various places in town and constant water levels in reservoirs, the flow and pressure distribution over the complete network can indeed be computed. These computed pressures can be checked against minimum pressure requirements. However, the assumption of constant demands is an oversimplification of reality. Water demands usually have daily and weekly cycles. Reservoirs may be filled during the night at cheaper electricity rates and emptied during the day, when demands are higher. Fire fighting may suddenly change the water demands over the network. These varying demands and storage lead to unsteady flow phenomena in pipe networks, although these are still often computed as series of quasi-steady states. Complete unsteady flow may occur as a result of sudden changes caused by failures or misoperation of the distribution network. This may cause unacceptable water hammer and cavitation effects. In this case computations are based on the use of the full unsteady flow equations for pipe flow, including storage of water resulting from water compressibility and pipe diameter expansion.

WL | Delft Hydraulics

10

Computational Hydraulics Version 2006-1

1.3

© 2005 Adri Verwey

UNESCO - IHE

Numerical methods

Although numerical methods have been developed already for centuries, the possibilities offered by computers have given a strong impulse to the further development of these methods. With numerical methods one tries to solve sets of differential equations, such as the De Saint Venant equations for unsteady flow in river- and channel systems. In the derivation of the differential equations one starts with finite control volumes for the definition of balance equations and than assumes that these volumes reduce to an infinitesimal small size. Under this assumption, these equations provide a correct and valid description of continuous flows. With numerical methods this process is reversed and balance equations are derived over finite control volumes, starting from the differential equations. For this reason, one cannot expect that this procedure leads to the same results as those which might have been obtained by solving the partial differential equations directly. However, in most cases such solutions do not exist, particularly when the equations are nonlinear. This, unfortunately, is usually the case when solving practical problems in hydraulics. Currently, for nearly all applications in hydraulics, numerical methods have been developed that work well and have the potential of limiting the differences between exact solutions and the approximate solutions. In these lectures we are discussing the differences by dealing with concept such as consistency, stability, robustness and economy of numerical operations. In particular, it will be shown how partial differential equations can be transformed into linear finite difference equations, to which extent these linearization’s require iterations and what sort of algorithms exist to solve the systems of equations in an economical way. It will also be shown, wherever applicable, that numerical behaviour is related to the physical behaviour described. Most obvious is the relation between boundary data requirements and the way changes at boundaries are affected by the hydraulic system. However, also the performance of iteration techniques are influenced by the physical behaviour of the system. The objective of dealing with numerical methods in this lecture series is to provide enough background information to serve as a basis for the correct development of models and for the best choice of modelling systems offered for use in a project.

1.4

Mathematical modelling

Modelling has become a frequently used tool for studies in hydraulic and environmental engineering. Whereas in the past many engineers turned to physically based models or simplified descriptions for the support of engineering studies, the increasing availability of personal computers and the powerful developments in computer graphics, data bases and on-line control, software has brought computer support to the desk of consulting engineers. In line with these developments we also see a strongly increased availability and use of mathematical modelling software tools.

WL | Delft Hydraulics

11

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

For a better understanding of what a model represents, let us look at one of the many possible ways of defining such a tool: a model is a physical or mathematical description of a physical system, including the interaction with its outside world, which can be used to simulate the effect of changes in the system itself or the effect of changes in the conditions imposed upon it. In the development of a mathematical model one may distinguish the following main elements: · · · · · · · · ·

definition of objectives schematisation equations and conditions solution algorithm software choice or development data collection model calibration model verification simulations

Definition of objectives Definition of objectives is a very important and often underestimated element in the decision process which leads to the use of a model. The first question to be posed is whether a model can add important information to what is already understood about a system. It also involves the estimation of the possibility to save project costs by using a model and the economic value that may be attached to the development and use of the model. This process will also lead to the choice of the level of complexity of model description. The use of models is generally associated to their use for simulations. However, another interesting field of application is in the analysis of possible hypothesis about empirical relations. Finally, models define relations between the various variables describing a state of a physical system and consequently models may be used for data consistency checking. As typical examples of model objectives in the fields of environmental, hydraulic and hydrological engineering one may mention: · · · · · · · · · · ·

WL | Delft Hydraulics

effects of hydraulic works; simulation of the impacts of floods; on-line flood forecasting; flood prediction; design of urban drainage systems; simulation of the impacts of dam breaks; study of field irrigation water supply; control of salt intrusion in estuaries; BOD-DO computations along rivers; ecological effects of heat loads from power plants; estimation of sedimentation in reservoirs;

12

Computational Hydraulics Version 2006-1

· · · ·

© 2005 Adri Verwey

UNESCO - IHE

soil remediation studies; effects of sewerage overflows upon the receiving waters; consistency checking of water quality data; consistency checking of hydrological data.

Model schematization The schematisation of the physical system follows from the complexity of the processes and the economical interest in studying these processes in all their details. The use of simplified models, based on a simplified description of the physical processes is only justified if the results of the model can still be used reliably in the design process. In other words, if the results still fall within the reliability range of data used by the designer. All processes in nature are of a three-dimensional and unsteady nature. The choice in the model schematisation is primarily: · ·

choice of the number of spatial dimensions; choice of time variability.

From the spatial dimensions and the time, the modeller selects one or more independent variables x,y,z or t, or other independent variables if certain transformations are applied, e.g. r and in a polar coordinate system. Such a choice is strongly linked to the model objectives. For example, a reservoir can be schematized into a single point, if one wants to study the water level variations and the reservoir outflow as a function of time. The same reservoir, however, will require a three-dimensional schematization in space in a study of wind-induced circulations or velocity patterns following from density stratification. Equations and boundary conditions A mathematical model is based primarily on the choice of equations describing the state of the physical system. Water levels, discharges, velocities, temperatures, salinities etc. are so-called state variables or dependent variables. One equation is needed for each of these variables describing the state of the physical system. Most of these equations are balance equations of mass, or, simplified, of volumes. Other equations are based on balances of other physical quantities, such as momentum, energy or turbulence. Often, also, simplified forms of these equations are used for the description of processes within our computational domain in space and time. In space this domain might represent the axis of a river flowing from town A to town B further downstream and in time the duration of a typical flood wave, including an antecedent period. Additional conditions have to be provided at the model boundaries, to specify the interaction of the outside world with the domain described by our model. These conditions will follow from the nature of the physical processes, translated into mathematical conditions through the equations describing the system.

WL | Delft Hydraulics

13

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Solution algorithms The balance equations in space and time provide us, generally, with partial differential equations. These equations are transformed from a continuous form into a discrete form by writing them out as relations between variables at points in the computational domain. Such formulations are based on finite differences, finite elements or boundary integrals and provide systems of linear equations. Completed with a linearized formulation of boundary conditions and special elements, such as hydraulic structures, the total system of linear equations may be solved with a variety of algorithms, ranging from simple Gaussian elimination as a direct matrix solver, to iterative techniques such as the conjugate gradient method. The solution algorithm usually has to follow a specific sequence of operations, consistent with the physical links in the system. Software For most environmental studies standard software tools are available. For simple problems engineers usually turn to spread sheet packages, whereas for problems related to open channel flow in networks, reservoirs, flow in pipe networks, sewer systems, water hammer, coastal management, short waves computations etc. various software products are available, developed at specialized hydraulic research institutes. The use of these packages assures a flexible user environment and a reliable solution of all sorts of numerical problems. Data collection Over the past years more and more effort has been spent on the collection of all sorts of data and the processing and storage of these data in data bases. However, for model development, data available in data bases is not always sufficient, as the calibration of models often requires data measured over short periods in time, available at various locations simultaneously. For this reason, the models are usually set up with whatever data available in the standard data base, completed with data collected during some specially organised campaigns. Model calibration and verification Some data required for a model can be collected directly in the field. Examples are salinities, channel cross-sections, discharges, water levels and concentrations of dissolved substances. Some data can only be collected with a certain degree of uncertainty, such as the details of the topography in the flood plains. Other types of data cannot be measured at all, such as Manning numbers and diffusion coefficients. Such data can only be estimated on the basis of a sound engineering judgement, based on the interpretation of recorded values of other variables and parameters. The more uncertainty we have in the model parameters, the more we are dependent on a good set of calibration data. The fit between measurements and computations and the knowledge of the processes enables the adjustment of the parameters until an acceptable fit has been obtained. For model calibration one will usually select a number of events, which are complementary to each other in terms of the calibration parameters.

WL | Delft Hydraulics

14

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

A calibration of a river model, for example, will typically start with low flow calibration and finish up with the calibration of some typical flood events. This procedure allows for the calibration of channel roughness parameters, prior to the calibration of flood plain conveyance and storage parameters. After completion of a model calibration the model should be verified on a set of data not yet used for parameter estimation. However, in practise it is not easy to reserve such a set and even if such verification runs are made, the differences between model and prototype performance may lead to lengthy arguments about the quality of the model data. In other words, why should one trust the results of verification more than the results of a model calibration? Simulations Once the model has been accepted it can be used for the typical simulations following from the definition of the model objectives. It should be kept in mind that the use of the model with modified parameters may, in turn, modify other parameters. As an example, the construction of river embankments may also change the bed roughness. The use of such models, therefore, should always be accompanied with sound engineering judgement based on a thorough knowledge of the physical processes.

WL | Delft Hydraulics

15

Computational Hydraulics Version 2006-1

WL | Delft Hydraulics

© 2005 Adri Verwey

UNESCO - IHE

16

Computational Hydraulics Version 2006-1

2

© 2005 Adri Verwey

UNESCO - IHE

Introduction to numerical methods

2.1

Introduction

This chapter deals with the solution of problems described on the basis of one single independent variable, such as x or t. A typical example is the simulation of flow through a reservoir, where water level variations in the horizontal plan are neglected and the point water level fluctuates in time as a result of time-variations in inflow and outflow. The water volume balance leads to an ordinary differential equation of the first order which can conveniently be solved by a finite difference approximation. Another typical example connected to this problem is the description of the variation in time of the concentration of a chemical substance in the reservoir water, assuming that this reservoir is well-mixed. The balance equation for this substance also leads to an ordinary differential equation with time as the independent variable. In a further extension one may consider the interaction of various substances, leading to coupled systems of first-order ordinary differential equations. A wellknown example is given by the Streeter-Phelps equations, describing the decay of organic pollutants and the effects on the oxygen concentration in a well-mixed pond. The same set of equations is found by considering the BOD and DO concentrations in a particle of water moving with the stream velocity in a river and neglecting the influence of exchange by diffusion between various water particles. Although this set of equations may be described as a function of time, it might also be convenient to apply a transformation which takes the channel axis x as the independent variable. Such transformation facilitates the inclusion of typical influences along the river, such as point effluent loads, weirs providing additional reaeration and the effects of variations in the channel topography on stream velocity and reaeration. Although the simplest problems of this category can be studied by exact solutions of the differential equations, the more flexible formulation of model equations and parameters requires solutions formulated by numerical schemes, such as finite difference schemes or finite element schemes. In the sequel, a range of finite difference approximation are introduced and discussed in terms of accuracy, stability and convergence of solutions.

Figure 2.1

WL | Delft Hydraulics

Sketch of a cylindrical groundwater reservoir

17

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

2.2

UNESCO - IHE

Introduction to finite difference approximations

Let us consider one of the simplest ordinary differential equations expressing the volume balance of groundwater in a cylindrical reservoir as shown in Figure (2.1). Assuming a constant reservoir surface area A, a constant porosity n and an outflow defined as a linear function of the hydraulic head h, the water level is described by the continuity equation and Darcy law, respectively, as

nA

dh dt

= -Q ;

Q = - kAh

or, by elimination of Q, Article 2. where

dh dt

= -ah

(2.1)

is a linear reservoir coefficient.

This equation has the exact solution

h = h0 e-a t

(2.2)

where h0 is the initial water level in the reservoir. This exact solution is given here primarily with the purpose of comparing various numerical schemes, solved with a variety of numerical parameters for these schemes. The numerical schemes that will be introduced successively are the · · · ·

Euler scheme; Improved Euler scheme; Implicit scheme; Newton-Raphson scheme.

The numerical parameter in this example is the time step of numerical integration t.

2.3

The Euler scheme

For the definition of a finite difference approximation we will first return to the concept of derivatives and the differential equations that may be constructed from these. The meaning of Equation (2.1) is that at any point along t, the derivative dh/dt to the function h(t) is equal to the value of the right hand side of that equation.

WL | Delft Hydraulics

18

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

Figure 2.2

UNESCO - IHE

Situation sketch for the definition of a derivative (tangent) to a function h(t)

Keeping in mind (Figure 2.2) that this derivative is defined as

dh dt

= L im Dt ® 0

h(t + Dt ) - h(t ) Dt

(2.3)

one may use the inverse of this definition to construct a finite difference scheme of Equation (2-3) on the basis of the approximation Dh dh h n +1 - h n @ = @ - a hn Dt Dt dt or (2.4) h n +1 @ (1 - aDt ) h n where it should be noted that n is introduced as a counter for the time step and written as a superscript to the variable h or to any other quantity defined at a grid point at time t=n t. This notation should not be confused with an exponent. The numerical scheme introduced this way is called an Euler scheme. The right hand side of Equation (2.6) is taken at time t=n t. Assuming that the value of h is known at that point, we call such a scheme a forward difference scheme as we construct a solution forward in time proceeding from a point where the solution is already known.

Figure 2.3

WL | Delft Hydraulics

Sketch of a finite difference approximation as the inverse of the definition of a derivative at a point A at t=n t

19

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Figure (2.3) also shows that for any non-linear solution the finite difference scheme introduces at each time step an error , defined as the difference between the exact solution of the equation and the solution obtained with the scheme. From the figure it may also be concluded that this error usually increases with an increasing time step. It will also be clear that the error depends on the curvature of the function, expressed by the influence of the higher-order derivatives. In other words, for solutions deviating slowly from straight lines one may take larger time steps than for solutions which deviate fast from straight lines, if one wants to obtain the same relative accuracy of the solution. Let us demonstrate the effect of the choice of time step by solving Equations (2.1) and with the Euler scheme of Equation (2.4) for the following data: h0 = 10 mm = 0.1 day-1 Table 2.1 shows results at T=4 days for the exact solution of the equations, compared with numerical solutions obtained with Equation (2.4) for time steps t = 0.5 days, 1.0 day and 2.0 days respectively. Differences between the exact solution and the numerical solution are 1 %, 2 % and 4.5 %, respectively. It is also observed that the errors are larger than those found at T = 2 days, demonstrating the accumulative effect of the errors during the integration of the differential equation along the time axis. Table 2.1

WL | Delft Hydraulics

Influence of the numerical scheme and the time step on the accuracy of results

time (days)

Exact

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

10.00 9.05 8.19 7.41 6.70

t=½ (days) 10.00 9.50 9.03 8.57 8.15 7.74 7.35 6.98 6.63

Euler t=1 (days) 10.00

t=2 (days) 10.00

Improved t=1 (days) 10.00

Euler t=2 (days) 10.00

Implicit t=1 (days) 10.00

t=2 (days) 10.00

9.50

9.00

9.00

9.05

9.05

8.60

8.10

8.00

8.19

8.20

8.19

8.18

7.78

7.29

7.38

7.41

7.41

7.04

6.56

6.40

6.71

6.72

6.70

6.69

20

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

Figure 2.4

2.4

UNESCO - IHE

Various derivatives used in the finite difference approximation: (a) derivative used in the Euler scheme; (b) ideal choice of derivative bringing the solution from point A to point B; (c) approximate derivative used in the centred schemes

The improved Euler scheme

From Figure (2.4) it is observed that one possible way of reducing the errors is by selecting a better location for the derivative along the time axis. A good location would be a place C, approximately halfway the time levels n t and (n+1) t. However, the exact location of the point is not known, due to the influence of the higher-order derivatives on the shape of the solution function. Moreover, even if this exact location would be known, the solution of the function is not known at this point and, hence the value of the derivative. The principle of a better positioning of the derivative leads us to the improved Euler scheme based on an additional iteration of the solution. As a first step in the improved Euler method the value of hn+½, halfway the time step, is approximated by h n +½

@

(1 - ½aDt ) h n

(2.5)

Substitution of this approximate value in the expression of the derivative gives the finite difference approximation over the total time step from grid point n to (n+1) as h n +1 - h n Dt

@ - a h n +½

or h n +1

@ h n - aDt h n +½

(2.6)

Results of this scheme, also given in Table (2.1), show considerable improvements in accuracy, with errors of 0.15 % and 0.3 % for time steps of 1.0 and 2.0 days, respectively. Even considering that the amount of computational work done with the improved Euler scheme is approximately twice the amount of work done with the normal Euler scheme, the improvement in terms of efficiency is still remarkable. For an equivalent computational effort the improved Euler scheme produces only 15 % of the error of the normal Euler scheme. Although similar improvements in accuracy are not always obtained for all problems, the example demonstrates the potential of efficiency improvements by using higher-accurate numerical schemes.

WL | Delft Hydraulics

21

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

2.5

UNESCO - IHE

The implicit scheme

The numerical schemes introduced so far are based on the concept that the derivative of the function is known at the initial point in the computation over a time step and does not change during that time step or, at best, is adapted in its value and location during an additional iteration. A next logical step for improvement, then, is to include the variables composing the derivative in the expression for the yet unknown variable h at time level (n + 1) t. For Equation (2.1) this leads to a centred implicit finite difference scheme of the form h n +1 - h n Dt or

h n +1

=

@ - (½a h n + ½a h n +1 )

1 - ½aDt n h 1 + ½aDt

For the given value of

(2.7)

and a time step of 1 day this leads to the simple relation

h n +1 = 0.9048 h n . Results of this computation are shown in Table (2.1). The errors are 0.017% and 0.13% for time steps of 1.0 and 2.0 days, respectively. In terms of accuracy, this approach gives another considerable improvement over the earlier introduced socalled explicit schemes. Although this conclusion may not be generalized to other numerical schemes without exceptions, the implicit schemes, in general, have the potential of providing more accurate results, at lower computations cost, due to the better centring of the finite difference equations. However, for problems involving more unknown variables, the implicit schemes lead to systems of linear equations, which have to be solved simultaneously through matrix operations. In most cases the overall solution algorithm leads to more numerical operations per time step. However, this is generally compensated by the much larger time steps that can be taken. As a consequence, currently most numerical algorithms are based upon implicit schemes. Apart from the higher accuracy, another important advantage of implicit schemes is their improved stability or robustness behaviour, as discussed in § 2.8.

2.6

The Newton-Raphson scheme

One special form of implicit schemes is the Newton-Raphson scheme. Although the Newton- Raphson approach is generally presented as a method for solving nonlinear equations, we introduce it here as an approach to the formulation of linear finite difference schemes and notably the linearization of individual terms and coefficients in the scheme.

WL | Delft Hydraulics

22

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

Assume, for example, that the coefficient function of h by the relation

a

UNESCO - IHE

in Equation (2.1) is given as a linear

= a0 + a1h

where a0 and a1 are constants. Substitution into Equation (2.1) gives dh dt

= - a0 h - a1h 2

= F

(2.8)

where F is a general function of h. An implicit finite difference scheme, centred halfway the time step, could be written as

Dh Dt

= F n + ½DF

As F is a function of h, the Newton-Raphson approach gives dF DF = Dh = - ( a0 + 2a1h ) Dh dh leading to the finite difference scheme Dh =

F n Dt 1 + ½Dt ( a0 + 2a1h n )

(2.9)

The advantage of this approach is that the change in the value of the coefficient is already partly taken into account during the integration over the time step. From Equation (2.13) one may conclude that the process is not yet ideal, as the value of h in the right hand side of the expression is taken at grid point n, whereas the substitution of a value at grid point (n+½) would be more precise. Referring, again, to the implicit scheme of § 2.5, the value of would be taken at grid point n, whereas an improved centring would require an additional iteration, similar to the approach followed in the improved Euler method. The Newton- Raphson implicit scheme is generally used without such additional iteration as in most cases this extra step is hardly cost effective.

2.7

A more formal analysis of accuracy

Although the concept of accuracy was discussed on the basis of a common sense approach and such a reasoning should always accompany the use of numerical and physical concepts, it is also useful to introduce more fundamental analysis techniques. One of the most frequently used tools for the analysis of the accuracy of numerical schemes is the Taylor's series expansion.

WL | Delft Hydraulics

23

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Referring to Equation (2.2) and assuming a continuous behaviour of the function h and all its first and higher-order derivatives in time, the value of h at grid point n+1 can be expanded from its value at grid point n through the infinite Taylor's series

h

n +1

Dt k æ d k h ö = h +å ç k ÷ k =1 k ! è dt ø n

¥

n

(2.10)

with all derivatives taken at grid point n. In Equation (2.10) k! should be read as “k factorial”, defined as the product 1*2*3*… .*k. It will be useful to discuss the meaning of the term under the summation sign in a pragmatic way. The term tk refers to a step in time raised to the power k and cancels, at least in magnitude, against the contribution dtk. The notation dkh has the meaning of expressing differences in function values at points in the vicinity of the point where the Taylor's series is expanded upon and has a value comparable in order of magnitude to the average function value at these points. As the value of k! increases rapidly with increasing k, it may be expected then that the contribution of the higher-order derivative terms in this expanded series decreases rapidly with increasing k. Referring to Equation (2.10), one may divide all terms by t to give h n +1 - h n Dt

n

n

n

2 2 3 æ dh ö Dt æ d h ö Dt æ d h ö = ç ÷ + ç 2÷ + ç ÷ + h.o.t . 2 è dt ø 6 è dt 3 ø è dt ø

(2.11)

where h.o.t. refers to all higher-order terms in this series expansion. Substitution of the Euler scheme of Equation (2.4) into Equation (2.11) and dropping the superscript n provides

dh dt

Dt d 2 h Dt 2 d 3 h = -ah + h.o.t . 2 3 2 dt 6 dt 14444244443 -TE

(2.12)

Comparison of Equation (2.12)with Equation (2.1) leads to a difference TE between the differential equation and the Euler scheme defined with the intention to solve this equation. This difference TE is called the truncation error of the finite difference scheme. For the improved Euler scheme, Figure (2.3) visualizes this truncation error by the magnitude . For the special case of Equation (2.1), the magnitude of all higherorder derivatives can be expressed in terms of the value of h at those points, as shown for constant by successive differentiation of the equation with respect to t.

WL | Delft Hydraulics

24

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

At all points in time this gives the truncation error in the form TE =

a2 h a3 h 2 Dt Dt + h.o.t. 2 6

(2.13)

From Equation (2.13) it is readily seen that by decreasing t in the integration of Equation (2.1), the total error in h decreases linearly with t with respect to the first term in the right hand side of Equation (2.13) and decreases even faster than linearly with respect to the remaining terms in the truncation error. A truncation error of this type is said to be of first order in t or simply O( t), while the numerical scheme is said to be first order accurate. As the numerical integration proceeds in time, the errors of each individual time step are accumulated. It should be remarked here that the accumulative error for any t is subject to a decay by virtue of the meaning of Equation (2.1) and for any realistic time step t the accumulated error will tend to zero for t à . Another interesting observation regarding the truncation error in the numerical integration of Equation (2.1) is the nearly-linear decrease of the error when reducing the time step t from 2.0 days to 0.5 days, as shown in Table 2. 1. This behaviour points at a rapidly decreasing influence of the higher-order derivatives in the truncation error, for this application. The much smaller truncation error in the implicit scheme of Equation (2.7) can also conveniently be demonstrated by a Taylor's series expansion. This derivation follows a more common introduction of Taylor's series in numerical schemes. After the selection of the appropriate centre point of the scheme all values of the dependent variables introduced at neighbouring grid points are expanded from that centre point. For the implicit scheme, centred at point n+½, the Taylor's series expansion gives 2 3 æ ö dh ( ½Dt ) d 2 h (½Dt ) d 3 h + + h o t = . . . ÷ (1 + ½aDt ) çç h + ½Dt + 2 3 ÷ dt dt dt 2 6 è ø 2 3 æ ö dh ( -½Dt ) d 2 h ( -½Dt ) d 3 h h o t + + . . . ÷ (1 - ½aDt ) çç h + ( -½Dt ) + 2 3 ÷ dt dt dt 2 6 è ø

where all values of h and the derivatives with respect to time are taken at grid point n+½. Expanding these expressions further leads to the equation dh dt

d 2h 1 d 3h a = - a h - Dt 2 2 - Dt 2 3 + h.o.t. 8 dt 424 dt 144444 2444444 3

(2.14)

-TE

which is the equation in reality solved by the implicit scheme.

WL | Delft Hydraulics

25

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Considering, again, that for this application all third- and higher-order derivatives are small, the remaining part of the truncation error is indeed much smaller than that of Equation (2.12) for any realistic choice of t. A realistic time step in this context is a choice which relates t to the value of as discussed in § 2.10.

Figure 2.5

Oscillations and instabilities generated for large time steps: (a) exact solution; (b) stable oscillatory solution for t=18 days; (c) unstable solution for t=22 days

2.8

Stability

Despite the truncation error in the computation presented in Table (2.1) most results are quite acceptable for practical purposes. However, it is interesting to observe what kind of results would have been produced if the time step had been taken much larger. As an example, let us consider a time step of 25 days in an application of the Euler scheme for the same equation and data as used in Table (2.1). For successive time steps the sequence of results would be 10, -15, 22.5, -33.75, 50.63, -75.94 etc., leading to infinity or an exponent overflow message on a digital computer after a sufficiently large number of time steps. In any hand computation the sequence of operations would have been interrupted after one or a few time steps as the results would appear to be unrealistic for any practical interpretation. A computation of this kind is called an instability. In an unstable computation results will always exceed a limit which has been set by the engineer as a realistic maximum or minimum value, for which exceedance is not to be expected (Figure 2.5). A frequently used analysis for the definition of stability criteria is based on the notion of amplification factors between results at successive time steps. If the absolute value of the amplification exceeds unity at each and every step in time, one has sufficient proof of the unstable nature of the computation. The application of this analysis to the Euler scheme of Equation (2.4) gives | A| =

h n +1 hn

=

1 - aDt

£ 1

as a stability condition for the scheme. Since condition for the Euler scheme is derived as t

WL | Delft Hydraulics

(2.15) t is definite positive the stability 2.

26

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Applying the same reasoning to the implicit scheme of Equation (2.7) leads to

1 - ½aDt 1 + ½aDt

£ 1

(2.16)

and to the conclusion that the implicit scheme is unconditionally stable. However, the stability limit of t=20 days obtained for this application of the Euler scheme is far beyond any time step that would be set as a maximum from the accuracy point of view. Even for the implicit scheme the results obtained with this time step are very inaccurate, as h drops to zero over the first time step and remains zero over all successive time steps, whereas in the exact solution the value of h decreases exponentially from the given initial value and only approaches the value of zero for tà . It may be included that the unconditional stability of the implicit scheme does not have special advantages in this application to the solution of ordinary differential equations. When moving to applications on partial differential equations, however, the increased stability of the implicit schemes will turn out to be of such great importance that currently, nearly all finite difference schemes are based upon implicit formulations.

2.9

Consistency and convergence

The Taylor's series expansion of the numerical scheme visualizes the difference between the differential equation as a continuous description of the physical system and the finite difference equation as a discrete description on a set of grid points. As t decreases, the difference between both equations reduces to the extent that they become equivalent as tà0. In such case the difference equation is said to be consistent With the differential equation. In the case of ordinary differential equations this generally implies that the results also converge towards the exact solution as tà0. An additional condition for such convergence is that the results are computed under stable conditions.

2.10

The choice of a time step

The choice of a time step in the numerical integration of the differential equation is a balance between the maximum tolerable error and the economy of the numerical operations. In any model application, errors are introduced from the following sources : · · · · ·

WL | Delft Hydraulics

accuracy of basic data; choice of additional parameter values; model schematization; choice of simulation data; numerical errors.

27

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Of all these sources of errors the numerical error is easiest controlled and a general requirement in model simulations is that the numerical error does not add to the uncertainties in the results introduced by the other sources. To satisfy this condition, the numerical error should be of an order of magnitude smaller than the overall expected error. So, if the overall error is expected to be some centimetres in level, the admissible numerical error should not exceed a few millimetres. Whereas it is already difficult to estimate the magnitude of the overall error, it is even more difficult to estimate the numerical error generated during one single time step and even more so the accumulated effect over various steps. As the truncation error contains higher-order derivatives, a first estimate of the time step is based on an idea about the curvature of the solution function. Strongly curved solution functions require smaller computational steps than solutions with more gentle variations.

Figure 2.6

Process of diminishing errors by reducing the time step t

A better estimate of the time step is based on a sensitivity analysis. By taking successively smaller grid steps, solutions are compared and if the differences appear to be sufficiently small, similar computations can be made with that acceptable grid step. It should be noted that the difference obtained by successively halving the time step is less than half the total numerical error , as shown in Figure (2.6), where this total numerical error at a given point in time is plotted as a function of the number of computational time steps N over a given period of integration. This leads to the conclusion that it is more efficient to refine coarse grids than refining further already fine grids. In a pragmatic approach one might also limit the allowable change in the function derivative over a single time step. For the simple Equation (2.4), this condition is equivalent to setting a maximum to the change in the value of h from one time step to the other. Setting, for example, as a criterion that over a single time step a change of 5% is allowed, this criterion leads to t 0.05, or t 0.5 days. Even in this simple case, however, it remains difficult to estimate the accumulated effect of this error over various time steps. In an attempt to do so and keeping in mind the decaying nature of the error, the accumulated error E at time step n=N, is approximated as

E =

N

å ½ Dt a

2

h n (1 - aDt )

N -n

(2.17)

n =1

WL | Delft Hydraulics

28

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

where in this notation N-n represents an exponent. In principle, it is possible to follow the magnitude of E during the computational process, assuming that the effect of third- and higher-order derivative contributions to the truncation error are negligible. Such an approach, however, is not practical for more complex problems and it is better to turn to sensitivity analysis for determining an acceptable time step.

2.11

Reservoir routing

Consider a reservoir with a level-dependent storage area A, an inflow Qi given as a hydrograph and a level-dependent spillway outflow Q0, as shown in Figure (2.7).

Figure 2.7

Sketch of a reservoir with spillway flow

Reservoir volume balance and spillway flow are given by the following set of equations

A

dh dt

= Qi ( t ) - Qo

(2.18)

Q0 = 1.71 m L ( h - hcr )

1.5

(2.19)

where A h

- level dependent surface area of the reservoir; - reservoir water level above a general reference (e.g. mean sea level

MSL); hcr - level of the spillway crest; L - length of the crest; m - discharge coefficient. In principle, Equation (2.19) may be substituted into Equation (2.18) to give one single equation with h as the only dependent variable. However, in general it is preferred to do such substitutions at the level of the numerical formulation after linearization of the equation and/or the finite difference formulation. In a simple Euler approach the equation reads as follows: Q0 n

WL | Delft Hydraulics

= 1.71 m L ( h n - hcr )

1.5

(2.20)

29

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

h n +1

= hn +

Dt Qi ( t n ) - Qo n An

{

UNESCO - IHE

}

(2.21)

where n, again, has the meaning of a superscript indicating the grid point along the time axis. To demonstrate the computational algorithm, a small example is worked out with the following data:

hinitial= 24.70 m hcr= 24.00 m L = 30 m t = 1 hour m = 1.1

t (hours) 0 1 2 3

Qi(t) (m3/s) 50 150 360 340

h (m) 24 25 26 27

A(h) m2 0.4*106 0.8*106 1.0*106 1.1*106

The results of the computation are given in Table (2.2). Table 2.2 time (hrs) 0.0 1.0 2.0 3.0

Results of the reservoir routing simulation with the Euler method grid point 0 1 2 3

hn (m) 24.700 24.790 25.345 26.472

Qon (m3/s) 33.05 39.62 88.02 219.32

Qin (m3/s) 50 150 360 340

An (m2) 6.80*105 7.16*105 8.69*105

h (m) 0.090 0.555 1.127

With reference to Figure (2.8) it is readily seen that the time step of 1 hour is too large, as the error introduced by using a reservoir surface area at time n t is significant. Moreover, the use of the discharge at time n t contributes to the error in the time derivative, although it does not affect the volume balance in a direct way.

Figure 2.8

Volume error at successive time steps in reservoir routine using the Euler scheme

A much improved formulation is based on the Newton-Raphson approach. Apart from the better centring of the derivative of Equation (2.23) the variation in the surface area is included implicitly. Although this variation is only introduced as a linear function, it is a great improvement over including A as a constant over a computational step, defined at time level n t.

WL | Delft Hydraulics

30

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

For the Newton-Raphson formulation the equations are rewritten in the form

dh dt

1 ( Qi ( t ) - Qo ) = F ( h, Qi , Qo ) A

=

(2.22)

and discretized, jointly with Equation (2.20), as

Dh Dt

= F n + ½DF

(2.23) n

Qo + DQo n

where DF

æ dQ ö = Qo ( h ) + ç o ÷ Dh è dh ø n

(2.24)

=

¶F ¶F ¶F DQi + DQo + Dh ¶Qi ¶Qo ¶h

=

1 1 n æ dA ö DQi - DQo ) Qi ( t ) - Qo n ç ÷ Dh 2 n ( A è dh ø ( An )

(

)

n

(2.25)

and dQo dh

= 1.5*1.71* m L ( h - hcr )

0.5

(2.26)

Special attention is drawn to the use of Q0(hn) in Equation (2.24) instead of substituting the already known value Q0n. As shown in Figure (2-9) the value Q0n was obtained from a linearized Q-h relation at time level (n-1) t (line a at point A). The solution should proceed from point B along the line b during the next time step. If in the right hand side of Equation (2.24) Qon had been used instead of Qo(hn), the solution would proceed from point C along the dotted line c and the accumulation of errors over various time steps would bring us further and further away from the correct Q-h relation. Note that for consistency reasons this correction is not included in the volume balance equation.

Figure 2.9

WL | Delft Hydraulics

Linearization of the Q0-h relation

31

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

The coefficients in the two linear equations can be collected to the form a11Dh + a12 DQo

= b1

a21Dh + a22 DQo

= b2

(2.27)

or -

-

Ax = b

(2.28) -

-

with matrix A and vectors x and b given as éa A = ê 11 ë a21

a12 ù ; a22 úû

é Dh ù x = ê ú; ë DQ0 û

éb ù b = ê 1ú ëb2 û

(2.29)

where a11 a12 a21 a22

(

)

2 An 1 æ dA ö = + n Qi ( t n ) - Qo n ç ÷ Dt A è dh ø = 1

æ dQ ö = ç 0÷ è dh ø = -1

n

n

(

b1

= 2 Qi ( t n ) + ½DQi - Qo n

b2

= Qo n - Qo ( h n )

(2.30)

)

Elimination of Q leads to

Dh =

b1 + b2 a11 + a21

DQ = b1 - a11Dh

(2.31) (2.32)

The computation over the same time steps as taken for the demonstration of the Euler method gives results as shown in Table (2.3). The results, indeed, are more realistic than those of Table 2.2. However, the differences between Q0n and Q0(hn) indicate that it is advisable to use a smaller time step of integration.

WL | Delft Hydraulics

32

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

Table 2.3

UNESCO - IHE

Results of the reservoir routing simulation with the implicit Newton-Raphson method

time (hrs) 0.0 1.0 2.0 3.0 time (hrs) 0.0 1.0 2.0 3.0

grid point 0 1 2 3

hn (m)

Q0n (m3/s)

Q0(hn) (m3/s)

24.700 24.992 25.688 26.379

33.05 53.73 114.45 199.71

33.05 55.75 123.78 207.02

Qi(t) (m3/s) 50 150 360 340

grid point 0 1 2 3

a11

a21

b1

b2

388 491 573

70.8 84.3 110.0

133.9 402.5 481.1

0.00 -2.03 -9.33

2.12

An*105 (m2) 6.80 7.97 9.38

h (m) 0.292 0.696 0.690

dA/dh *105 (m) 4.00 4.00 2.00

dQ0/dh (m2/s) 70.82 84.31 109.98

Q0 (m3/s) 20.68 60.72 85.27

Routing of a decayable pollutant through a reservoir

Let us consider next a similar reservoir which is polluted by a decayable substance with concentration c. Assuming that the reservoir volume can be schematized as well mixed and introducing a first order decay, with reaction coefficient k, gives d ( cV ) = ci Qi - c Qo - k cV dt

(2.33)

or

V

dc dh ( c ) + cA dt dt

= ci Qi - c Qo - k cV

(2.34)

Substitution of Equation (2.18) into Equation (2.34) and division of all terms by the reservoir volume V leads to

dc dt

=

Qi ( ci - c ) - kc V

(2.35)

where the volume V, as a function of the reservoir level, follows from the integration of the surface area A along h. Such integration is best based upon a simple trapezium rule as the surface area, obtained from planimetering a topography from a map, is never accurate enough to justify higher-accuracy integrations of the A-h relation, such as the Simpson's rule or even integrations based on cubic spline functions. Moreover, linear A-h relations are used in the Newton-Raphson method and the use of the trapezium rule for the integration of the V-h relation is consistent with this approach. Although in a quick analysis the dilution coefficient Qi D = (2.36) V is often set as a constant, giving for Equation (2.35) the exact solution c = c0 e - ( k + D ) t -

WL | Delft Hydraulics

D ci ( e - ( k + D )t - 1) k+D

(2.37)

33

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

We will focus here on the general case where D is a function of time and the exact solution of Equation (2.35) does not exist. The simplest numerical solution of Equation (2.35) is, again, based on the Euler method, giving n æQ ö = c n + Dt ç i ÷ ( ci - c ) - kc n Dt èV ø

c n +1

(2.38)

Numerical solutions have the advantage over analytical solutions that all parameters can easily be made a function of the dependent variables c and h and of the independent variable t. Examples of such further generalization of relations both in the water quantity and quality part are: · · · · · · · ·

outflow partly given as a user demand, possibly affected by reservoir operation levels; inclusion of the effects of precipitation to and evapotranspiration from the reservoir surface in a longer term water balance simulation; stage dependent width of the spillway crest; controlled spillway crest level as a function of h or t; the use of m, calibrated as a function of the reservoir level; various sources of pollutant inflows into the reservoir; decay described as a function of the time-varying water temperature; inclusion of additional substances, such as dissolved oxygen with dependence on wind-generated reaeration, photosynthesis, bottom sediment processes etc.

Again, in this example more accurate results are obtained for a given time step with the improved Euler method and implicit methods, including the Newton-Raphson formulation. As the water balance is not affected by the pollutant concentration, any implicit technique leads to a system of linear equations with the unknown concentration at the new time level decoupled from the unknown level and out flowing discharge. It is also rather common to simulate the water quantity part first, write results to a data base and subsequently retrieve the necessary information in a separate water quantity computation. We will return to this approach when discussing water quality studies in rivers, coastal areas and in reservoirs where the assumption of a well mixed estuary is not correct.

2.13

Non-uniform steady flow in channels

As an example of computations which describe steady processes with variations in the spatial x-direction let us consider the backwater computation in a uniform channel. A typical channel cross-section is shown in Figure (2.10), including the conveyance K plotted as a function of the water level, where K

=

jj

å j =1

WL | Delft Hydraulics

1 Aj R j 2 / 3 nj

(2.39)

34

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

with nj Aj Rj jj

= Manning roughness coefficient for the sub section j; = sub area of the cross-section j; = local hydraulic radius of the sub area j; = number of vertical slices used in the integration across the section.

Figure 2.10 Cross-section and cross-sectional parameters used in the non-uniform flow computation

Assuming the absence of lateral flow, the steady state De Saint Venant equations reduce to the form Q = Qcons tan t

(2.40)

Q 1 d æ Q ö dh + I0 + 2 ç ÷+ gA dx è A ø dx K 2

2

= 0

(2.41)

Further differentiation of the first term of Equation (2.41) gives, for constant Q, æ Q 2 dA ö dh Q2 ç1 - 3 ÷ + I0 + 2 K è gA dh ø dx

= 0

(2.42)

or dh dx

=

æ Q2 ö I + ç 0 ÷ Fr 2 - 1 è K2 ø 1

(2.43)

with the dimensionless Froude number Fr defined as a function of the flow velocity u, the cross-sectional area A and the storage width b, as Fr 2

=

u2 A g bs

=

Q 2 dA gA3 dh

(2.44)

For a given discharge, the right hand side of Equation (2.43) is a function of the water level above the channel bottom, giving

dh dx

WL | Delft Hydraulics

= F ( h)

(2.45)

35

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Table (2.4) gives an example of such function F(h), where it is readily seen at which depth the equilibrium slope is reached. This depth, for which F(h) = 0, is called the normal depth. The solution can easily be found in a Newton-Raphson approach, where, following the usual notation, a subscript i has been introduced as a grid point counter along the x-axis. The implicit finite difference scheme reads as Dhi Dx

æ dF ö = Fi + ½ ç ÷ Dhi è dh øi

(2.46)

or

Fi Dx (2.47) æ dF ö 1- ½ ç ÷ Dx è dh øi For the cross-sectional data given in Figure (2.10), a discharge Q=3.40 m3/s, a bed slope I0=-10-4, the uniform flow depth is found for dh/dx=0, or Q2=-K2I0, giving K=340 m3/s and h=1.00 m from the conveyance function shown in Table (2.4). At the downstream end of the channel the flow is controlled by a weir with a free overflow given by Equation (2.19) with parameter values: Dhi

=

m = 1.1 L = 11.0 m hcr = 2.00 m

(weir coefficient); (weir crest length); (weir crest level),

giving a water level above the bottom just upstream of the weir of h0=2.30 m. The value serves as the downstream boundary condition for the non-uniform flow computation. In the situation sketch of Figure (2.11) attention is drawn to the fact that the numerical integration of Equation (2.43) proceeds in the negative xdirection, requiring a negative value for the distance step x. A suitable choice for this distance step follows from the consideration that in the first integration step the change in h should be limited to a reasonably small fraction of the difference between the downstream boundary water depth and the uniform flow depth. Considering that the Newton-Raphson approach gives fairly accurate results, a first step h = 0.10 m, covering approximately 10 % of the total difference, should be acceptable. Linear interpolation in Table 2.3 gives a value K = 1498 m3/s at the downstream boundary. For a friction slope Q2/K2 = 0.052 * 10-4, giving a difference of 0.95 * 10-4 with the bed slope, a choice x = -1000 m satisfies the criterion selected. Completing the table with the water level gradient function F of Equation (2.43) for the given discharge (Table 2.4), the upstream water depths can be computed by applying Equation (2.47). Results are shown in Table (2.5) over a distance of 12 km upstream of the weir.

WL | Delft Hydraulics

36

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Figure 2.11 Positive backwater upstream of a hydraulic structure Table 2.4

Cross-sectional data given as a tabulated function of the depth, including discharge dependent functions for Q = 3.4 m3/s (see also Figure 2.10)

water depth h (m) 1.0

bs (m)

A (m2)

K (m3/s)

Q2/K2 *10-4

Fr

F *10-4

10.0

9.00

340

1.000

0.126

0.000

1.5

11.0

14.25

688

0.244

0.066

0.756

2.0

12.0

20.00

1143

0.088

0.042

0.912

2.5

20.0

28.13

1734

0.038

0.032

0.962

(dF/dh) *10-4 (m-1) 1.512 0.312 0.100

Table 2.5

Non-uniform flow computation based on the Newton-Raphson implicit approach (a) and compared with results of the Euler method (b) (a)

Newton-

Raphson

(b)

Euler

h (m)

hi (m)

Fi *10-4

h (m)

0.369

2.300 2.206 2.113 2.021 1.930 1.841 1.755 1.671 1.590 1.512 1.436 1.370 1.314

0.942 0.933 0.923 0.914 0.890 0.862 0.836 0.809 0.784 0.760 0.659 0.559

-0.094 -0.093 -0.092 -0.091 -0.089 -0.086 -0.084 -0.081 -0.078 -0.076 -0.066 -0.056

x (m)

grid point i

hi (m)

Fi *10-4

0 -1000 -2000 -3000 -4000 -5000 -6000 -7000 -8000 -9000 -10000 -11000 -12000

0 1 2 3 4 5 6 7 8 9 10 11 12

2.300

0.942

(dF/dh)i *10-4 (m-1) 0.100

1.930

0.890

0.312

0.335

1.595

0.786

0.312

0.296

1.327

The results are also computed by using the Euler scheme which differs from the Newton- Raphson scheme by the second term in the denominator of Equation (2.47). Assuming that the Newton-Raphson results are rather accurate, for this distance step, the Euler scheme produces results which might be considered acceptable for practical use. Although the difference of 1.3 cm at x=-12000 m will still increase further upstream, the cross-section schematization, the constant bed slope and the choice of roughness coefficients will give rise to significantly larger uncertainties about the correct water levels.

WL | Delft Hydraulics

37

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

For another observation we refer to Table (2.4). The number of levels for which the cross- sectional parameters are given is fairly small and produce unacceptable errors in the linear interpolation of these values. The errors introduced are most likely larger than those caused by using a lower accuracy numerical integration scheme. Although this might be acceptable in a quick analysis by hand computation, any computer code would be programmed to set up the table with much smaller steps. This table step size should not exceed substantially the step h in the numerical integration. For our problem it would have been advisable to present the functions bs, A, K, Fr and F at intervals of 10 cm. The advantage of the Newton-Raphson scheme over the Euler scheme is the higher accuracy obtained at the cost of the simple additional computation of the denominator of Equation (2.47). For a comparison of its efficiency with that of the normal Euler scheme let us recompute the solution with a four times larger distance step, as shown in Table (2.6). Table 2.6

x (m) 0 -4,000 -8,000 -12,000

Recomputation of the backwater curve with x=-4,000 m

grid point i 0 1 2 3

(a)

Newton-

Raphson

hi (m) 2.300 1.931 1.596 1.300

Fi *10-4 0.942 0.890 0.786

(dF/dh)i *10-4 0.100 0.312 0.312

h (m) 0.369 0.335 0.296

(b)

Euler

hi (m) 2.300 1.923 1.568 1.257

Fi *10-4 0.942 0.888 0.777

h (m) -0.377 -0.355 -0.311

The Newton-Raphson approach still maintains its high accuracy during the first two steps. During the third step the effect of the large table step in felt in the correction through the term dF/dh. The difference of 1.7 cm with the computation taking x =1,000 m, would certainly have been much smaller if Table (2.4) had been set up with smaller steps h. The result obtained with the Euler method for this case differs 7.0 cm from the most accurate result and is quite unacceptable. The advantage of the Newton-Raphson scheme becomes clearer when negative backwaters with much stronger water profile curvatures are computed. Such computations start from downstream boundary conditions with depths below the normal depth h=1.00m. There is no doubt that the Newton-Raphson scheme is far more efficient in its application than the Euler scheme.

2.14

An inverse scheme for the backwater curve computation

It has been demonstrated that once the table of the function F of Equation (2.45) has been constructed for a uniform cross-section, the numerical operations for the nonuniform flow computation are fairly simple. For uniform channels with a constant bed slope we will now introduce a scheme which combines a high accuracy with an extreme simplicity. The computation of depths with the Newton-Raphson method introduced in the previous section is an asymptotic process, requiring continuous computations with ever decreasing changes in the depth.

WL | Delft Hydraulics

38

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Of course, one could have decided to increase the distance step as the depth comes closer and closer to the normal depth. However, returning to Equation (2.45), one could also introduce a numerical scheme taking equal steps h in the dependent variable h, a so-called inverse scheme. The advantage of this approach is that the right hand side of the equation, being only a function of h, can be centred without requiring the iteration of the improved Euler scheme. The simple inverse scheme looks as follows:

Dxi

=

Dh Fi +1/ 2

(2.48)

and results are shown in Table (2.7) for steps h = 0.10 m. These results are based on the same interval for presenting topographic data as given in Table (2.4). Interpolation at x = - 8,000 m in Table (2.7) gives a depth of 1.596 m, which only differs 0.1 cm from the result obtained with the Newton-Raphson scheme. At x = 12,000 m, however, the difference is 0.5 cm, which difference must mainly be attributed to the Newton-Raphson integration based on too large table intervals (grid point 9-10). The normal depth h=1.00 m is found at x = -33 km, whereas theoretically this should be at x = - . However, this deviation is of no practical importance at all. Although the inverse scheme is very handy for uniform channel computations, it looses its advantage or even does not work at all in applications where the channel parameters vary along the channel axis. Table 2.7

WL | Delft Hydraulics

Results of the inverse scheme for h = 0.10 m

i 0

xi (m) 0

1

-1,067

2

-2,146

3

-3,237

4

-4,353

5

-5,509

6

-6,708

7

-7,954

8

-9,250

9

-10,720

10

-12,610

11

-15,256

hi (m) 2.300 2.25 2.200 2.15 2.100 2.05 2.000 1.95 1.900 1.85 1.800 1.75 1.700 1.65 1.600 1.55 1.500 1.45 1.400 1.35 1.300 1.25 1.200 1.15

Fi+½ * 10-4

xi (m)

0.9370

-1,067

0.9270

-1,079

0.9170

-1,091

0.8964

-1,116

0.8652

-1,156

0.8340

-1,199

0.8028

-1,246

0.7716

-1,296

0.6804

-1,470

0.5292

-1,890

0.3780

-2,646

0.2268

-4,409

39

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

12

-19,665

13

-32,893

2.15

1.100 1.05 1.000

0.0756

UNESCO - IHE

-13,228

Non-uniform flow in non-uniform channels

Let us now consider the somewhat more complex case of steady flow in nonuniform channels. For such channels the cross-sectional area is a function of both h and x, giving the variation in the x-t space

dA =

¶A ¶A dx + dh ¶x ¶h

(2.49)

dA = dx

¶A ¶A dh + ¶x ¶h dx

(2.50)

or

Equation (2.43) should now read as dh dx

=

æ Q 2 Q 2 ¶A ö I + ç ÷ 0 Fr 2 - 1 è K 2 gA3 ¶x ø 1

(2.51)

For the solution of this equation it is most practical to process the function F (x, h) as a tabulated function at equal depth steps h, starting from a reference point near the channel bottom. It is advisable to choose a h that has the same value at all cross-sections given along the river. For each river section, between two successive cross-section points, a set of parallel lines is found along which F (x, h) reflects the change in F along x, while keeping h constant. This approach will require at each cross-section the processing of one table reflecting the values dA/dx applicable to the left hand channel section and another table reflecting their values applicable to the right hand section. The Newton-Raphson integration now proceeds with the scheme Dhi Dx

æ ¶F ö æ ¶F ö = Fi + ½ ç ÷ Dx + ½ ç ÷ Dhi è ¶x øi è ¶h øi

(2.52)

or Dhi

=

æ ¶F ö 2 Fi Dx + ½ ç ÷ Dx x ¶ è øi æ ¶F ö 1- ½ ç ÷ Dx è ¶h øi

(2.53)

in an integration performed over an integer number of distance steps along each successive channel section marked by its end cross-sections and its local reference slope. Again here, the procedure is most efficient if only one Newton-Raphson iteration is made per time step and care is taken that the distance steps are not too large.

WL | Delft Hydraulics

40

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Some final remarks: 1. Both methods applied to uniform and to non-uniform channels can easily be 2. 3. 4. 5.

extended to include the effects of lateral flows. Both methods apply to fully irregular cross-section geometry and allow for a variation of roughness coefficients across the channel. The methods are applicable for any formulation of the friction laws entered into the integration of the channel conveyance K. The efficient inverse scheme of Equation (2.48) can, in general, not be applied to non-uniform channel flow computations. Both methods can be applied to sub critical flow (Fr1), if the appropriate direction of integration (algorithmic structure) is followed.

2.16

What have we learnt?

In this Chapter we introduced various numerical methods for the solution of (sets of) ordinary differential equations. Many of these problems can simply be programmed in Excel. This is very straightforward in cases where equation parameters are constants. Cases with non-linear parameters can still be programmed easily by using macros for the function descriptions. When programming solutions for incidental application it is recommended to use the simple Euler method and to make sure that a sufficiently small grid step is used. Higher accuracy methods only justify the additional effort of programming and testing of the code when repetitive use of this code is made and a single computation is not finalized instantaneously. With the current computer speed this is no longer the case for the relatively simple problems programmed in Excel. The higher accuracy integration methods have been introduced primarily, as the principles demonstrated are commonly applied in the numerical methods for solving partial differential equations. Applications of iterative methods as applied in the improved Euler method, implicit schemes and Newton Raphson principles, will be discussed in Chapters 4, 5 and 6.

2.17

Questions and small assignments

1. Extend the computation of all cases shown in Table (2.1) up to T = 6.0 days. 2. Explain in your own words why the improved Euler method is more attractive in

use than the standard Euler method. 3. Determine the type of numerical scheme used in the well-known reservoir

routing equation

DS Dt

WL | Delft Hydraulics

=

1 ( I1 + I 2 - O1 - O2 ) 2

41

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

where S = reservoir storage ; I = reservoir inflow ; O = reservoir outflow, with subscripts 1 and 2 referring to time levels n t and (n+l) t respectively. 4. Explain why the resulting functions in Table (2.1), both obtained with the Euler

method and the implicit scheme, is below the exact solution. Also explain why the resulting function obtained with the improved Euler method lies above the exact solution. 5. Find with a computation carrying more decimal positions than shown in Table

(2.1) the ratio between the errors obtained with the Euler scheme and the implicit scheme for a time step of 1 day. Discuss this ratio in the light of the truncation errors for both schemes. 6. Verify Equation (2.14). 7. Show with a more precise computation of the results of the implicit scheme

computation in Table (2.1) that at T = 4 days the error obtained with t = 2 days is four time larger than the error obtained with t = 1 day. Explain this difference from the truncation error of the scheme as shown in Equation (2.14). 8. Give the fourth order derivative term of the truncation error of Equation (2.14). 9. Define an Euler scheme for Equation (2.1) based on a numerical integration in

the negative time direction and show that the scheme is unconditionally unstable. 10. Extend the computation shown in Table (2.2) with one additional time step. 11. Extend the computation shown in Table (2.3) with one additional distance step. 12. Extend the computation shown in Table (2.5) with one additional distance step. 13. Explain why the depths calculated with the Euler scheme in Table (2.5) are

smaller than those obtained with the Newton-Raphson scheme. 14. Explain why in Table (2.6) the depth obtained in a computation with a distance

step of 4 km is smaller than the depth obtained with a 1 km step. 15. Recompute the solution of Table (2.7) with h = 0.05 m over the range 2.1

h

2.3 m. Comment on the results. 16. Explain what sort of problems may result from the application of Equation (2.48)

for the computation of non-uniform channel flow.

WL | Delft Hydraulics

42

Computational Hydraulics Version 2006-1

3

© 2005 Adri Verwey

UNESCO - IHE

Basic Unsteady Channel Flow Equations

3.1

Introduction

Unsteady channel flow will be first derived for a simple uniform channel of unit width, a horizontal bottom with friction along the walls and the bottom neglected. Referring to Figure (3.1), the water depth is symbolized by h and the average flow velocity by u. These variables u and h are usually called the dependent variables, defined as functions of the independent variables x (space coordinate) and t (time coordinate). The definition of two dependent variables requires the formulation of two equations for their solution. In channel flow these equations are usually based upon volume and momentum conservation. The volume conservation approach assumes incompressible flow and constant water density. This is a quite realistic assumption for free surface fresh water channel systems. The equation following from this volume conservation principle is called the continuity equation.

Figure 3.1

Control volume along a channel axis

3.2

Continuity equation

The derivation of the equation is based on the concept that over a given time the difference between inflow to and outflow from a control volume along the channel balances the change in the storage in this control volume over this same time. Referring, again, to Figure (3.1), it is readily seen that the inflow per unit time into the control volume, or the flux of water volume through the upstream cross-section, is given by

inflow = uh Assuming that both u and h are continuous functions of x, the outflow at the downstream section of the control volume can then be defined as ¶ outflow = uh + (uh) dx ¶x

WL | Delft Hydraulics

43

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

The storage in the control section per unit time is expressed in terms of the change in water level as ¶ storage = (h dx) ¶t based on the concept that for continuous water level changes with initial depth h0, the depth h1 after a step in time dt equals

¶h dt ¶t Balancing terms gives the continuity equation in the form h1 = h0 +

¶h ¶ + (uh) = 0 ¶t ¶x

(3.1)

Differentiating out the second term leads to another, often used, form

¶h ¶h ¶u +u +h =0 ¶t ¶x ¶x

3.3

(3.2)

Momentum equation

It should be recalled that momentum is defined as the product of velocity and mass. The total momentum M of the fluid contained in our control volume of unit width, therefore, is defined as M = u r h dx where M and u are both vectors acting along the channel axis. Over a certain step dt in time, the change of M is balanced by the net amount of momentum carried by the fluid velocity through the boundaries of the control volume and the impulse generated by the forces acting upon it. Both quantities are to be taken with their component acting in the x-direction. Following the same principles as used in the derivation of the continuity equation, the net amount of momentum carried with the flow velocity is defined by

M a = u r uh - [u r uh +

¶ (u r uh) dx ] ¶x

where the first term refers to the inflow of momentum and the term within braces to the outflow of momentum per unit time. We will furthermore recall Newton's second law

F dt = d(mu)

WL | Delft Hydraulics

44

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

where F is a force acting on the fluid with mass m and velocity u. The product of the force and the time period over which it is acting is called an impulse. Newton's second law states that this impulse generates a change in momentum.

Figure 3.2

Forces acting on a horizontal uniform channel section

Under the simplifications introduced above, the only forces acting on the control volume are the hydrostatic forces at the upstream and downstream ends of the control section. Integrating the hydrostatic pressures over the verticals (Figure 3.2) gives the net force ¶ dF = ½r g h 2 - [ ½r g h 2 + ( ½r g h 2 ) dx ] ¶x Balancing all terms over a unit time step leads, after division of all terms by the constant dx, to the equation

¶ ¶ ¶ ( r uh ) + ( ru 2 h ) + (½r g h 2 ) = 0 ¶t ¶x ¶x Dividing, furthermore, all terms by the constant , leads to

¶ ¶ ¶h (uh) + ( u 2 h ) + gh = 0 ¶t ¶x ¶x

(3.3)

Differentiating out the product terms in the derivatives, substituting the continuity equation and dividing all the remaining terms by h, leads to

¶u ¶u ¶h +u + g =0 ¶t ¶x ¶x

3.4

(3.4)

Transformation to the characteristic form

Although the form in which the set of continuity and momentum Equations (3.1) and (3.3) is shown is very convenient for the formulation of solutions on the basis of finite differences, it is useful to introduce a transformation which enables us to understand better the meaning of these equations. For this purpose we will introduce the auxiliary variable

WL | Delft Hydraulics

45

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

c = gh

UNESCO - IHE

(3.5)

Multiplication of Equation (3.1), substitution of Equation (3.3) it and division of all terms by c gives

¶c ¶c ¶u + 2u + c = 0 ¶t ¶x ¶x Similarly, the substitution of Equation (3.4)into Equation (3.3)s 2

(3.6)

¶u ¶u ¶c +u + 2c = 0 ¶t ¶x ¶x

(3.7)

The set of Equations (3.5), (3.6) represents the characteristic form of the flow equations. The special significance of this form will become clear when respectively Equation (3.5) is added to and subtracted from Equation (3.6), to give the set of equations in the form

and

¶ ¶ (u + 2c) + (u + c) (u + 2c) = 0 ¶t ¶x

(3.8)

¶ ¶ (u - 2c) + (u - c) (u - 2c) = 0 ¶t ¶x

(3.9)

Before drawing any conclusions from these relations, we will introduce the general mathematical formulation expressing the change of the value of a function f when moving from a given point P to a neighbouring point Q.

Figure 3.3

Relation between function values at neighbouring points

As shown in Figure (3.3), this change results from the partial derivatives of the function in t-direction and x-direction respectively as

df =

¶f ¶f dt + dx ¶t ¶x

(3.10)

or, after division of both sides by dt

WL | Delft Hydraulics

46

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

df dt

=

¶f dx ¶f + ¶t dt ¶x

UNESCO - IHE

(3.11)

Taking for f the function u+2c, readily shows by comparison of Equations (3.7) and (3.10) that along the line with direction

dx = u+c dt the expression d (u + 2c) = 0 dt

(3.12)

(3.13)

holds, while the comparison of Equations (3.8) and (3.10) shows that along the line with direction

dx = u-c dt

(3.14)

the expression

d (u - 2c) = 0 dt

(3.15)

holds (Figure 3.4).

Figure 3.4

Characteristic directions and their Riemannn invariants

The line with the direction expressed by Equation (3.11) is called the c+ characteristic, whereas the line with the direction expressed by Equation (3.13) is called the c- characteristic. The integration of Equation (3.12) gives

u + 2c = J + = constant

(3.16)

whereas the integration of Equation (3.14) yields the expression

WL | Delft Hydraulics

47

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

u - 2c = J - = constant

UNESCO - IHE

(3.17)

In this formulation J+ and J- are the so-called Riemann invariants. For the simplified channel flow equations these invariants have constant values along the so-called c+ and c- characteristics respectively. When introducing additional terms in the equation, such as lateral flow, bottom slope and bottom friction, these Riemann invariants will contain correction terms and will no longer be "invariant".

3.5

The significance of the characteristics

The conditions derived along the characteristics enable us to compute solutions from known conditions at an earlier point in time. Referring to Figure (3.5), it is seen that the state of the fluid at point P can be computed from given states at points A and B. Given, for example, that at point A, u=1.0 m/s and h=5 m, while at point B, u=1.2 m/s and h=4.8 m, the Riemann invariants are computed as J+ = u+2Ögh = 15 m/s at the points A and P, and as J- = u-2Ögh = -12.52 m/s at the points B and P. Solving these two equations representing the Riemann invariants at point P gives u=1.24 m/s and h=4.83 m.

Figure 3.5

Characterisitics leaving from points A and B, defining the solution of the flow state at a point P

The significance of the characteristics, however, goes far beyond the possibility to make such computations. The fact that along the characteristics a special condition is valid, implies that they pass on some information on the state of the fluid. In other words, the characteristics are to be seen as lines along which information on the state of the fluid propagates. The uniqueness of the equivalence between Equation (3.10) on the one hand and Equations (3.7) and (3.8) on the other hand, implies that these lines along which information propagates are unique lines and that information only travels along just these characteristics. The auxiliary variable c introduced in Equation (3.4) appears to get a special meaning as part of the direction of the characteristics. For the special case where u=0, both characteristic directions have the same magnitude, though their signs are opposite (Figure 3.6a). As c is a function of the water depth, this means that the effect of any disturbance or control at a point along the river, propagates in both channel directions with the same speed, which is a function of the square root of the water depth. The variable c gets here the meaning of a celerity with which disturbances propagate in stagnant water.

WL | Delft Hydraulics

48

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

Figure 3.6

UNESCO - IHE

Characteristic directions for: (a) and (b) subcritical flow; (c) critical flow and (d) supercritical flow

For velocities u>0, the characteristic directions no longer have the same magnitude and, at a particular point, information on the fluid state travels faster in the downstream then in the upstream direction (Figure 3.6b). For increasing water velocities this leads to the particular case where the celerity c equals the water velocity u (Figure 3.6c) and at this point information is no longer able to travel in the upstream direction. Introducing the Froude number, defined as |u| (3.18) gh it is readily seen that for the situation of Figure (3.6c), where u=c, the Froude number equals unity. This special case is called critical flow, whereas the cases of Figures (3.6a, 3.6b), where the Froude numbers are less than unity, represent sub critical flow. For water velocities exceeding the celerity c, the Froude numbers exceed unity and the flow is called super critical (Figure 3.6d). Fr =

From Figure (3.6) it is readily seen that for the case of sub critical flow, the state of the fluid at any point of space and at any point in time is controlled by both the upstream and the downstream conditions. In the case of super critical flow, however, the state of the fluid is only controlled by the upstream conditions. So, if at any point A at a channel, we have supercritical flow, downstream control of the conditions at A is impossible, as long as the super critical state at A is maintained. It should be noted, however, that this supercritical state may change into a subcritical state when downstream control, by a weir for example, generates a region of subcritical flow upstream of the control structure to the extent that this subcritical flow region reaches point A. From Figure (3.6) it may also be concluded that it takes a certain time before the effect of the opening of the gate at the upstream end of an irrigation channel reaches the point where the water is required. Taking, for example, a channel with stagnant water at 1 meter depth, giving a celerity of 3.13 m/s, it lasts nearly 1 hour before additional water reaches a demand point after the further opening of a gate 10 km upstream.

WL | Delft Hydraulics

49

Computational Hydraulics Version 2006-1

Figure 3.7

© 2005 Adri Verwey

UNESCO - IHE

Regions of determinacy of the state at a point P for various flow conditions

To make these principles even more clear, Figure (3.7) shows various cases of subcritical and supercritical flow, where the shaded areas show the region of determinacy of the state at point P. Any disturbance or control within the shaded region will have some effect on the water velocity and depth at point P, while any disturbance or control outside the shaded region will not be able to affect conditions at point P. The propagation of information along the characteristics will also determine the initial- and boundary condition requirements for the solution inside a computational domain, bounded in the x-t plane by the lines t=0 and t=T and the lines x=0 and x=L. Referring first to Figure (3.8a), point A at time t=0 receives two characteristics coming from outside the computational domain. If point A would have been located inside the computational domain, the solution would have followed from the Riemann invariants along the two characteristics passing through it, much similar to the situation sketched in Figure (3.5). With A located at the boundary of the computational domain, it is can be concluded that the information on the Riemann invariants has to be replaced by initial conditions, as we have no information on the state of the fluid outside our computational domain and the values of the Riemann invariants associated with it.

Figure 3.8

Relation between characterisitic directions and the boundary- and initial data requirements

Similarly, at point B in Figure (3.8a), one characteristic arrives from inside the computational domain and is supposed to carry information from another point within the computational domain. The other characteristic arrives from outside the domain, where the computation is not made and the information, which would normally be passed on via the Riemann invariant, has to be replaced by a boundary condition. Such boundary condition is usually provided in the form of a given velocity, water level, discharge or stage-discharge relationship.

WL | Delft Hydraulics

50

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

The same applies to point C, where one characteristic carries information from inside the computational domain, while the other characteristic comes from outside and requires a boundary condition to replace the information which would otherwise be provided by its Riemann invariant. Figure (3.8b) shows the case of supercritical flow with positive u. The problem is completely upstream controlled and two upstream boundary conditions are needed for a solution. We define this need as a two-point boundary condition against a onepoint boundary condition for the case of subcritical flow, given in Figure (3.8a). Initial data requirements are the same as for the subcritical flow case and two-point initial data are required for all types of flow. Usually, values for u and h are to be given. Figure (3.8c), finally, shows the case of a supercritical flow with velocity u in the negative x-direction. For the same reasons as pointed out earlier, this system requires two-point initial data and two-point boundary data at the upstream end, which in this case is located at the point x=L. From the discussion above the general conclusion can be derived: one initial- or boundary condition is required for every characteristic entering the computational domain, where entering is defined by following the characteristic forward in time.

3.6

The method of characteristics

In § 3.3 we have introduced the Riemann invariants and the computation of the state variables u and h with these Riemann invariants. Where necessary, these Riemann invariants are complemented with the initial- and boundary conditions. We may now proceed to the description of a complete solution algorithm, defined as the method of characteristics. This method is based on the construction of a network of points where the state variables are computed. As shown in Figure (3.5), from a given set of two points on the network a new point may be constructed by drawing a characteristic of one family from one point, which intersects with a characteristic of the other family drawn from the other point. In general, these directions are not known initially, as they depend on the solution of the state variables at the intersecting point. However, the solution of these state variables can be found through the Riemann invariants, even if the exact location of the point is not known yet. The nearly exact direction of these characteristics can then be drawn based on the average characteristic direction at the connecting points. The procedure is best explained by considering a practical example. Figure (3.9) shows the application of the method on the problem of closing a gate in an irrigation canal. Initially, the water flows at steady state with a velocity of 1.0 m/s at a depth of 1.60 meter. The gate is closed over a period of 100 seconds, controlled to give a linearly decreasing water velocity immediately downstream of the gate.

WL | Delft Hydraulics

51

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

This assumption serves our exercise, as most gate control operations are based upon water level sensoring. For a further simplification in this example it is also assumed that the downstream channel is very long and that the wave, generated by the gate closure, is not reflected at a downstream side.

Figure 3.9

Hodograph and physical planes for a graphical solution of the computations along characteristics

For ease of computation the acceleration due to gravity is taken as g=10 m2/s. For the initial steady state, the following data can be derived: c+ = u+ Ögh = 5.00 m/s c- = u- Ögh = -3.00 m/s J+ = u+2Ögh = 9.00 m/s J- = u-2Ögh = -7.00 m/s The disturbance generated by the gate operation will travel in the downstream direction with a celerity c+=5.00 m/s. This implies, for example, that at a point L meter downstream of the gate, the flow remains undisturbed over a period of L/c+ seconds.

WL | Delft Hydraulics

52

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

As a result, one may expect a steady state region downstream of the gate, shown as region A in Figure (3.9). This is the region enclosed by the lines (t=0; x 0) and (x=5.00 t; x 0). The solution will first be defined at points along the t-axis, immediately downstream of the gate. As a first step a set of points will be selected and numbered from 1 onward at intervals of 25 seconds. Subscripts u and d will be used for the locations of these points upstream and downstream of the gate respectively. Let us first consider the solution of the problem downstream of the gate. The solution at point 1d of Figure (3.9) is found by applying the boundary condition u=0.75 m/s at t=25 s (Figure 4.9b), and the condition that the Riemann invariant J-=7.00 m/s at point 1d, as it keeps its constant value with which it arrives along the negative characteristic from the steady state region A. The solution of this set of two equations gives h=1.50 m. Similarly, the solution at points 2d,3d and 4d is found by applying the velocity boundary condition combined with the Riemann invariant carried along the negative characteristic arriving from the steady state region A as: point 2d: u=0.50 m/s & J-=-7.00 m/s --> h=1.41 m; point 3d: u=0.25 m/s & J-=-7.00 m/s --> h=1.31 m; point 4d: u=0.00 m/s & J-=-7.00 m/s --> h=1.22 m. A summary of these data is shown in Table 3.1. Table 3.1 Data at various characteristic net points of the physical plane Point

c+ (m/s)

c(m/s)

2 gh

4.00

5.00

-3.00

1.125 0.705 0.328 0.000

3.88 3.75 3.62 3.50

4.63 4.25 3.88 3.50

1.125 0.705 0.328 0.000

4.20 4.31 4.39 4.50

4.84 4.69 4.56 4.50

u (m/s)

h (m)

q (m2/s)

0

1.0

1.60

1.600

1d 2d 3d 4d

0.75 0.50 0.25 0.00

1.50 1.41 1.31 1.22

1u 2u 3u 4u

0.64 0.38 0.17 0.00

1.76 1.86 1.93 2.03

gh (m/s)

J+ (m/s)

J(m/s)

8.00

9.00

-7.00

-3.13 -3.25 -3.37 -3.50

7.75 7.50 7.25 7.00

8.50 8.00 7.50 7.00

-7.00 -7.00 -7.00 -7.00

-3.56 -3.93 -4.22 -4.50

8.40 8.62 8.79 9.00

9.00 9.00 9.00 9.00

-7.76 -8.24 -8.62 -9.00

(m/s)

Let us now turn to points x 0 away from the boundary and draw the positive characteristic leaving from point 1d. Conditions at all points along this line are given by the same Riemann invariants: J+ = u+2Ögh = 9.00 m/s and J- = u-2Ögh = -7.00 m/s. Consequently, the same solution is found at all points along this line, giving u=0.25 m/s and h=1.50 m. This characteristic, then, is a straight line defined by a celerity c+ = u1+Ögh1 = 4.63 m/s.

WL | Delft Hydraulics

53

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Similarly, the positive characteristics leaving from points 2,3 and 4 must be straight lines with celerities of 4.25, 3.88 and 3.50 m/s respectively. The region between the positive characteristics leaving from points 1d and 4d is defined as a simple wave region, characterised by straight (positive) characteristics of one family and curved characteristics of the other family. However, it is not necessary here to draw the curved c- characteristics, as the visualisation of these lines does not add any essential information on the state of the flow in such a simple wave region. For the same reason, only the characteristics bordering a steady state region are drawn and none of the characteristics inside the steady state region. Having defined a simple wave region one might define a complex wave region as one where both families of characteristics are curved lines, indicating that velocities and depths vary from point to point in both directions. In such a complex wave region a network of both families of characteristics is drawn. This situation would occur if the simple wave generated from points 1d to 4d would be reflected against a downstream boundary as could easily be verified by introducing, for example, a constant water level condition at a point, say, 400 m downstream of the gate. As will be explained later, however, we do not want to carry the presentation of the method of characteristics too far and refer the interested reader to standard text books on the topic (e.g. Abbott, 1979). However, as an example of a curved characteristic and for demonstration purpose only, consider the negative characteristic between point 1d and the steady state region A shown in Figure (3.9). The direction of this characteristic is varying from point to point, as follows from the varying J+-values, while moving along this characteristic from point 1d to region A. As a reasonable approximation the average direction is taken as c - 1d + c - A = -3.13 - 3.00 = - 3.07 m/s 2 2 and the approximate characteristic is drawn backward in time, leaving from point 1 until its arrival at the steady state region. c - 1d ® A =

From this procedure it may be concluded that, although the method of characteristics is generally considered a quite accurate solution of the simplified partial differential equations describing channel flow, it does not represent an exact solution of these differential equations. The approximation is in the location of the points of the network. The solution of the simplified equations at the network points, however, is exact. The more generally applied finite difference methods, to the contrary, have exact locations of the network- or grid points, while the solution of the equations at these points, as a rule, are approximations of the true solutions.

WL | Delft Hydraulics

54

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

3.7

UNESCO - IHE

More complex boundary conditions

Let us now consider the upstream part of the channel, where the water velocity is reduced and the flux of momentum from upstream is likely to lead to a piling up of the water against the control structure. The effect of any control will be travelling in upstream direction with the celerity c- = -3.00 m/s. As also in this case some time will elapse before the disturbance reaches an upstream point along the canal, we will also find here a steady state region shown as D in Figure (3.9). The solution upstream of the gate requires gate boundary conditions. These are not the same as given for the downstream part, as the velocities given are related to cross-sections and water levels just downstream of the gate. The solution at this downstream section, however, enables the combination of velocity and depth to give the specific discharges, which are the same just upstream and downstream of the structure. The solution of the state at point 1u then follows from the computed specific discharge q=udhd=1.125 m2/s and the positive Riemann invariant arriving from steady state region D, J+ = u+2Ögh = 9.00 m/s. Substitution of one equation into the other leads to a third degree polynomial in h, which requires a Newton-Raphson iteration algorithm for their solution. For non-computerized solutions it is therefore easier to turn to a graphical method. Graphically, the solution of any set of two equations is found at the intersection of the lines representing these equations in an appropriate coordinate system. For the given problem, a suitable choice of the coordinates is a plot of values 2Ögh linearly along the vertical axis and values u along the horizontal axis. The plane formed by these two variables is called the hodograph plane. The choice of the vertical scale leads to the representation of the Riemann invariants as lines under 45o (Figure 3.9c). By plotting the 2Ögh-axis vertically down, for subcritical flow the sign of the directions coincides with that of the characteristic directions in the physical plane. This convention is quite convenient in its practical use. Table 3.2

WL | Delft Hydraulics

Data for drawing q-relations in the hodograph plane

q m2/s

u m/s

h m

1.125

.60 .70

0.705 0.328

gh

2 gh

m/s

m/s

1.87 1.61

4.32 4.01

8.65 8.02

.35 .40

2.01 1.76

4.48 4.20

8.97 8.39

.18 .16

1.82 2.05

4.27 4.53

8.53 9.06

55

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Figure (3.10) also shows various boundary conditions, which may all easily be presented in graphical form, such as u given (Figure 4.10c), h given (Figure 3.10d), a given specific discharge q (Figure 3.10e) and a critical flow condition (Figure 3.10f). This last relation is a special form of a rating curve, or Q-h relation, which, for the simplified channel flow equations, would normally be presented as a u-h relation. The plot of a given specific discharge, as applied at point 2u of Figure (3.9c), is easily made by selecting a number of u-values and calculating the associated h- and 2Ögh-values. It should be realised that q should always be plotted with the appropriate sign. Also, the line is curved and requires a certain number of computed points in the expected neighbourhood of intersection for a sufficient accuracy.

Figure 3.10 Various boundary conditions represented graphically in the hodograph plane

We return again to Figure (3.9c), which shows the hodograph plane belonging to the physical plane of Figure (3.9a). The boundary conditions q=1.125, 0.705 and 0.328 m2/s at points 1, 2 and 3 respectively, are presented by three different lines in the hodograph plane, whereas the discharge q=0 is represented by the line u=0. The solutions at points 1u,2u,3u and 4u are shown in the hodograph plane as the intersections of the lines representing the specific discharge at those points and the condition J+=9.00 m/s. The characteristics leaving the t-axis over the period of the closure, again, define a simple wave region, with characteristics now converging towards each other. Further upstream this will lead to the intersection of these ccharacteristics. Such intersection of characteristics of the same family leads to twovalued solutions for both u and h at such intersection, or, in other words, to the formation of a hydraulic jump. It should be realised that the intersection of characteristics is found here on the basis of the equations which do not contain bed slope and bed friction. Where these terms are playing a significant role the computations become more complex and are the subject of another chapter.

WL | Delft Hydraulics

56

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Although not frequent in nature or in man-made systems, such travelling hydraulic jumps are known phenomena, for example, as tidal bores in estuaries or as undular hydraulic jumps in tail race channels of peak power generating hydropower systems.

3.8

The formation of positive and negative hydraulic jumps

The discussion on the formation of hydraulic jumps is simplest presented on the basis of an extreme operation mode of the gate. Figure (3.11) shows the case of an instantaneous closure of otherwise the same problem as discussed in Figure (3.9). Considering first the upstream side of the gate, the positive Riemann invariant carried from the steady state region D has to be combined with the condition u=0 at x=0, for t>0. As shown in the hodograph plane the solution of these two conditions is u=0 m/s; h=2.03 m. The negative characteristic celerity is c-=-4.50 m/s, which is greater than c-=-3.00 m/s found over the steady state region D. This difference in celerities leads to the immediate intersection of the negative characteristic directions just upstream of the gate at t=0 and to the formation of a positive hydraulic jump.

WL | Delft Hydraulics

57

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Figure 3.11 Method of characteristics for the case of a sudden gate closure

We have assumed that despite this intersection of characteristics, the positive Riemann invariant remains constant across the jump. We will, furthermore, assume that the celerity of the jump is approximated by the average of the two celerities of the intersecting characteristics forming this jump. Both hypotheses have been investigated and appear to be approximately true for small jumps (Abbott 1979). Based on these assumptions the celerity of the jump is computed at c-=-3.75 m/s. The line representing this propagation in the physical plane separates the steady state region D from another steady state region C bordered by the line x=0, for t>0 and the path of the hydraulic jump.

WL | Delft Hydraulics

58

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

A similar construction at the downstream side of the gate gives a steady state region B, where h=1.22 m. Steady state regions A and B give characteristic celerities c+=5.00 m/s and c+=3.50 m/s respectively, which show diverging directions originating from point (x=0, t=0). Only at the very origin, a discontinuity is present, which rapidly spreads along x while travelling downstream. We call such special form of a discontinuity a negative hydraulic jump and the subsequent simple wave region a centred simple wave. In this simple wave region the characteristics drawn can be selected freely and Figure (3.11) shows the choice of some characteristics representing steps of 10 cm in h. The results of the computation along x, at T=125 s are shown in Figure (3.11c). The discontinuity formed at the positive jump and the spreading negative hydraulic jump in the downstream part are clearly recognised.

3.9

The limited practical importance of the method of characteristics

The method of characteristics has the principal advantage of visualising the way in which flow disturbances or the effects of flow control travel through the hydraulic system. As will be shown in the Chapters dealing with mathematical modelling, the structure of the net of characteristics is also very important in understanding the numerical procedures required for the practical solution of hydraulic phenomena. In other words, a good understanding of the physics of a hydraulic phenomena guides us in the choice of appropriate numerical solution techniques and their numerical parameters. As already discussed in § 3.4, the characteristics show us the initial- and boundary data requirements for any hydraulic engineering problem simulation. However, where the concept of characteristics is very important to us, the method of characteristics provides us only for simple practical cases with an algorithm for the computation of unsteady flow systems. For this reason we do not expand here on the inclusion of energy generating or dissipating terms, such as gravity terms and bottom friction, especially not for non-uniform channels with irregular topographies. In this text we also do not show how to include lateral flows. Inclusion of such terms is much easier and flexibly handled in programmable solution algorithms based on implicit finite difference schemes, as discussed in Chapter 5. We will return, however, to the practical use of characteristics in the case of water hammer and flood propagation simulations.

WL | Delft Hydraulics

59

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

3.10

Questions and assignments

3.10.1

Derive Equations (3.1) and (3.3) again for the case where a channel width b is introduced.

3.10.2

The simplified model is extended with a lateral flow term, by introducing the lateral flow per unit width and length of the channel q (m/s). a. Derive the new form of Equations (3.1) and (3.3). b. Give the transformation to the characteristic form. c. Does this term give any change in the characteristic directions? d. Derive the new expressions for the Riemann invariants.

3.10.3

Calculate u and h at point P of Figure (3.5) for the following data: Point A: u=0.5 m/s; h=3.00 m; Point B: u=0.6 m/s; h=2.80 m.

3.10.4

A rectangular channel of 10 meter width carries a discharge of 6 m3/s. For a water velocity u=2m/s, calculate the state of the flow. What initial and boundary conditions do we need for modelling a certain reach of the channel?

3.10.5

At the downstream end of a channel the flow is critical. A positive characteristic arrives at this boundary carrying a Riemann invariant J+=7 m/s. Calculate the outflow velocity and depth.

3.10.6.

A gate in an irrigation canal of unit width is closed over a period of 4 minutes. The initial velocity of the water in the canal is 1.0 m/s at a depth of 0.9 m. Assume a linear decrease of the discharge during the closure. a. Sketch the structure of the characteristics upstream and downstream of the gate. Take a total time of 8 minutes and a channel length of 2000 meter upstream and downstream of the gate. Take one characteristic every 2 minutes (only for the purpose of this exercise, as for practical solutions a somewhat denser net will be needed). b. Draw the hodograph plane and solve the depths and velocities at the intersection points of the characteristics. c. Draw a more precise network of characteristics. d. Draw the water level and velocity functions along the canal at 8 minutes after the initiation of the closure. e. What is going to happen at a later moment in time?

3.10.7.

WL | Delft Hydraulics

Repeat the same exercise for an instantaneous closure of the gate.

60

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Introducing Numerical Solutions for Partial Differential Equations

4

4.1

Introduction

Partial differential equations are introduced when the dependent variables of differential equations are a function of more than one independent variable. In many practical applications these independent variables are x for space and t for time. In this case, derivatives of the dependent variables have fractions both in space and in time. For this reason we speak about partial derivatives, when dealing with a derivative along one of these two or more independent variables. Consequently, we speak about partial differential equations when rates of changes of these dependent variables are related along the various dependent variables. Depending of the nature of the physical processes, these partial differential equations can be classified as belonging to one of the following types: · · ·

hyperbolic equations; parabolic equations; and elliptic equations.

The classification of the Equation (s) for a certain problem leads to distinct differences in which these equations are solved. Typical application areas for these various classes are: · · ·

hyperbolic equations: unsteady surface water flow; unsteady flow in the unsaturated soil zone; flood wave propagation; transport of pollutants, channel morphology; parabolic type: dispersion of pollutants; flood wave peak dampening; dispersion of heat or energy; elliptic equations: steady surface water flow; steady ground-water flow.

The difference in type lies in the way in which disturbances propagate from one point of a given (physical) system to other points. In the case of hyperbolic problems it takes a definite time before a change of the state of the system at one point influences the state at another point. The effect of these changes propagates along so-called characteristics, or lines along which disturbances propagate. An obvious example is the propagation of flood waves along rivers. Generally, it takes a few hours or a number of days before a rain event in the upper catchments leads to floods at locations downstream. In parabolic problems, the effect of a change somewhere in the system, immediately affects the state everywhere else. In this case, the characteristics travel at an infinite speed.

WL | Delft Hydraulics

61

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

In a step further, elliptic problems even have imaginary characteristics only and these lead to an even more pronounced effect of the influence of the change of the state at one single point on the state of the total system. Examples are state descriptions of steady surface water or ground-water flow. It should be realised that these are rather artificial descriptions, as true steady states do not exist in nature. For this reason, most parabolic or elliptic problems are solved as if one is dealing with a hyperbolic problem. The simplest of the equations of hyperbolic type is the so-called advection equation (see Equation (4.3)). Advection describes how some quantity, such as mass, momentum, energy, heat, a pollutant etc. is transported by a carrier, such as water or air. When dealing with the water sector, there is a variety of problems which can be described by equations where advection processes play an important role. Examples are: ·

·

·

·

WL | Delft Hydraulics

convective momentum term in the unsteady flow equations. Although often not of much importance in most gradually varied flow situations, the advection term (or convective momentum term) plays a dominant role in the formation of hydraulic jumps, including the special appearances of moving hydraulic jumps, such as tidal bores; flood wave propagation, as a special case of the application of unsteady flow equations. Unsteady river- or channel flow is described by the De Saint Venant equations. These form a set of second order partial differential equations of hyperbolic type with two characteristic directions. In the most common case of sub critical flow, these characteristics have opposite directions. In the special case of flood wave propagation along rivers, gravity and friction play a dominant role. This results in the emergence of an underlying simplified advectiondiffusion equation, of which the advection part is only first order hyperbolic. The characteristic direction of this part is shown to us by the propagation celerity of the flood wave peak; in hydrology also the movement of ground-water in the unsaturated zone is based upon a combination of advection and diffusion processes, generally described by the Richards equation. The advection part of this equation is driven by gravity and by soil suction forces (matric forces). Due to the complex relation between soil permeability and soil moisture, the process is extremely non-linear, leading to the propagation of soil moisture shock fronts, much similar to the propagation of hydraulic jumps in free surface flow; transport and dispersion of substances in water. A good description of these processes is essential for the assessment of water quality processes, such as the spreading of pollutants and heat, growth of algae, sedimentation processes etc. An important element of this description, again, is the advection part of the advection-dispersion equation.

62

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

4.2

UNESCO - IHE

The advection equation

As an example, consider the advection of a conservative matter (e.g. a chemical pollutant) in a channel. The derivation of the equation describing this process proceeds via the introduction of a control element along the channel axis as shown in Figure (4.1).

Figure 4.1

Control element for the derivation of the advection equation.

As a first step, the equation of mass conservation for water is derived. referring to Figure (4.1), the mass balance for a control element with length dx along the x-axis is set up by defining the change in storage to be equal to the difference between inflow and outflow, or over a time period dt, ¶ ( rV ) dt ¶t

ì ¶ ü é = í rQ - ê rQ + ( rQ ) dx ùú ý dt ¶x ë ûþ î

Replacing the volume V by the product of cross-sectional area A and distance step dx it is found, for constant density , that

¶A ¶Q + =0 ¶t ¶x

(4.1)

In a similar way, the mass balance for a conservative matter in the water with concentration c can be derived as ¶ ¶ (4.2) ( Ac ) + ( Qc ) = 0 ¶t ¶x or, after differentiating out the terms of this equation and substituting Equation (4.1),

¶c ¶c +u ¶t ¶x

= 0

(4.3)

where the velocity u is defined as Q/A. In this derivation it is assumed that the concentration c is uniformly distributed over the cross-section. For the time being, it is also assumed that the cross-sectional area

WL | Delft Hydraulics

63

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

for flow is equal to the cross-sectional area representing storage. Equation (4.3) is valid at every point of the x-t plane and shows the relation at the solution surface between the tangent in x-direction (to the function which shows the distribution of c along the river at a fixed time) and the tangent in t-direction (to the function showing how c varies in time at a fixed point along the channel).

4.3

The characteristic solution

To get more insight into the meaning of Equation (4.3) let us consider the total derivative

dc =

¶c ¶c dt + dx ¶t ¶x

(4.4)

This equation simply states how the variation of c in any direction in the x-t plane is composed of variations of c along the t- and x-axis, respectively. Division by dt shows that, by equivalence of the Equations (4.3) and (4.4) along the line

dx dt

= u

(4.5)

the condition holds that

dc dt

= 0

(4.6)

or, after integration of Equation (4.6), that c remains constant along this line. The line defined by Equation (4.5) is called the characteristic of Equation (4.3). In Figure (4.2) it is shown how a solution of Equation (4.3) is obtained, for the simple case where u=constant and where, consequently, the characteristics are parallel lines in the x-t plane.

Figure 4.2

Solution of the advection equation with the method of characteristics.

It is also seen that a solution is only obtained when in the computational domain initial data are defined along the line t=0 and, where an upstream boundary is

WL | Delft Hydraulics

64

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

present, also boundary data along the line x=0. As will be shown later, these characteristics play an important role in defining where initial- and boundary data are required, as well as the number and the type of these conditions. For practical applications the method of characteristics is not often applied in this form, mainly because it has the disadvantage that it is not based upon a fixed grid in space. This makes it more complicated to include in the solution such processes as sources and sinks, decay and diffusion, which will be discussed later.

4.4

Finite difference schemes

On the grid, shown in Figure (4.3), the derivative c/ x can be written in several ways in terms of finite differences and, at grid point (i,n) more in particular as ¶c ¶x ¶c ¶x ¶c ¶x

Figure 4.3

@ @ @

cin - cin-1 Dx n ci +1 - cin-1 2Dx n ci +1 - cin Dx

backward difference approximation; centred difference approximation; forward difference approximation.

Grid for defining finite difference schemes of the advection equation.

In a similar way finite difference approximations for c/ t can be defined. For the advection scheme these approximations can be combined to nine possible finite difference schemes, already. This number even increases when combinations of 1), 2) and 3) are made by introducing, for example, other time levels than n t. Not all these schemes will give satisfactory results and it is needed, therefore, to define criteria for judging which scheme can be selected as the best out of the large number of possible schemes. Let us consider, in the first instance, a scheme approximated by a backward difference in the x-direction and a forward difference in the t-direction. Equation (4.3) than reads as

WL | Delft Hydraulics

65

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

cin +1 - cin c n - cin-1 +u i Dt Dx

UNESCO - IHE

= 0

or, expressing cin+1 in terms of the other variables and defining uDt r = Dx

cin +1

= r cin-1 + (1 - r ) cin

(4.7) (4.8)

This relation can be represented graphically by the operator shown in Figure (4.4).

Figure 4.4

Operator for the scheme of Equation Error! Reference source not found..

A solution using this scheme is shown in Figure (4.5) for the following data: u

= 0.0185 m/s; x = 100 m; t = 0.75 hours.

Upstream boundary data are set as c=0 and the initial data are given by the Gaussian distribution of the form c = e

-½ (

x-m 2 ) s

(4.9)

with mean value =500 m and standard deviation =100 m.

WL | Delft Hydraulics

66

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

Figure 4.5

UNESCO - IHE

Numerical solution of the advection equation using the scheme given by Equation (4.8) with r=½ and various values for x.

Figure (4.5) shows the effect of varying numerical parameters on the solution. The exact solution is defined by a shift of the initial Gaussian distribution over 400 m along the x-axis. Numerically, the same solution is found for r=l. For r=½ the wave is strongly dampened. This effect is called numerical diffusion of the wave. The figure also shows that a smaller value for x, while maintaining the same value r=½, reduces the numerical diffusion. For a value r=2 the solution will not look realistic. Negative values for c are found as well as values greater than 1. These values cannot possibly be correct as the method of characteristics shows that in the absence of source- and sink terms extreme values for c are found along the boundaries where characteristics enter the computational domain. The problems faced here are related to the stability of the scheme. It is clear that not only the choice of a specific scheme is important for obtaining good results in a computation, but also the choice of the numerical parameters, such as x and t, for that scheme. The accuracy aspects of a particular scheme can be judged by writing each term in a Taylor's series expansion from the selected centre point of the scheme. At the centre point (i,n) we obtain for the scheme given by Equation (4.8): n +1 i

c and

Dt = c + 1! n i

n

n

n

Dt 2 æ ¶ 2 c ö Dt 3 æ ¶ 3 c ö æ ¶c ö + + ç ÷ ç ÷ + h.o.t . ç ÷ 2! è ¶t 2 øi 3! è ¶t 3 ø i è ¶t øi

n -Dx 3 ) æ ¶ 3c ö ( -Dx æ ¶c ö Dx 2 æ ¶ 2 c ö = c + ç ÷ + ç ÷ + h.o.t . ç ÷ + 1! è ¶x øi 2! è ¶x 2 øi 3! è ¶x3 øi n

n i -1

c

WL | Delft Hydraulics

n

n i

67

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

where h.o.t. stands for higher order terms. Collecting terms in Equation (4.8) and dividing by t gives at point (i,n): ¶c ¶c +u ¶t ¶x

=

-Dt ¶ 2 c r Dx 2 ¶ 2c Dt 2 ¶ 3c r Dx 3 ¶ 3c + + h.o.t. (4.10) 2 ¶t 2 2Dt ¶x 2 6 ¶t 3 6Dt ¶x 3

The right hand side of this equation is called the truncation error of the finite difference scheme. It can be concluded that by using the scheme of Equation (4.8), another equation, containing higher order derivatives, is solved than the differential equation (4.3). For looking more closely at the lowest order terms of the truncation error, Equation (4.3) is differentiated with respect to t to give, for constant u, ¶ 2c ¶t 2

= -u

¶ 2c ¶x¶t

= -u

¶ 2c ¶t ¶x

= u2

¶ 2c ¶x 2

Substitution of this relation into Equation (4.10) gives ¶c ¶c +u ¶t ¶x

= ½u Dx (1 - r )

¶ 2c + h.o.t. ¶x 2

(4.11)

Referring already to Equation (4.17), it will be shown in Chapter 4 that this equation adds diffusion to the advection equation. For the example of Figure (4.5), with r = ½ and x = 100 m, the numerical diffusion coefficient is 0.463 m2/s. Comparison of Figure (4.5) with Figure (4.9) of § 4.7 confirms the magnitude of the diffusive effect of the truncation error. Equation (4.11) shows, furthermore, that the lowest order term of the truncation error goes to zero for r=1. The computed results for r = l indicate that a further substitution will cancel out all the higher order terms in the truncation error for this specific value of r, to give the exact solution of the differential equation. For values r = ½ the numerical diffusion can be reduced by taking smaller and smaller values for x. In the limit, for x --> 0, the truncation error will go to zero (see Equation (4.11)) so that the numerical solution converges towards the true solution of Equation (4.3). This convergence will not be found for the solutions with r=2. It can be shown that for this r-value subsequent reduction of x will give a solution deviating more and more from the true solution. The results become more and more unstable. There are several ways in which this stability problem can be investigated: 1. 2. 3. 4.

by interpreting the dominant effects of the truncation error of the scheme; by showing that the values in the computational domain will always remain between the extreme values specified at the boundaries; by an analysis based on the role of characteristics on the solution (CourantFriedrich-Lewy or CFL-condition); by a Fourier series analysis of the solution function and its transformation from one time level to another.

Ad 1) For positive values of the numerical diffusion coefficient in Equation (4.10), the solution will be dampened within the computational domain.

WL | Delft Hydraulics

68

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Negative diffusion coefficients are physically unrealistic and would lead to amplification of the solution. This would affect most strongly the shorter wave lengths as those have the largest values for 2c/ x2. In this way, small errors will rapidly be blown up to large values. This leads to the conclusion that only positive diffusion coefficients will give stable results and thus, from Equation (4.11), that r should not exceed unity. Ad 2) The exact solution of Equation (4.3) shows that extremes in the computational domain should always remain between minimum and maximum values given at the boundaries of this domain. It can be shown that the scheme, indeed, gives values which are limited by the extremes at the boundaries for r 1 (Abbott, 1979). Ad 3) Referring to Figure (4.7) and following the characteristics through point (i,n+l), the scheme will be stable when the characteristic intersects the line t=n t between the grid points (i-1,n) and (i,n), as in this case cin+1 is found from an interpolation of earlier computed results. When the intersection of the characteristic and the line t=n t is outside the line element between the points (i-l,n) and (i,n), the value cin+1 follows from an extrapolation of data and the scheme will become unstable. The time step can be compared with the time tc required to travel along the characteristic over a distance x. The ratio t/ tc is defined as the Courant number Cr. For Cr 1 the scheme will be stable and it will be unstable for Cr>1. This case is shown in Figure (4.7) when, for the characteristic leaving the point (i-1,n+1) backward in time, the numerical scheme of Equation (4.8) is applied at the grid point (i-1,n). Ad 4) A solution at time level n t is decomposed into the Fourier series

c

n j

åx

=

n k

e

k

x ik p l

=

åx

n k

e

ik

j p jj

(4.12)

k

where im = indicator for imaginary numbers; k = wave component number; n k = wave amplitude at time level n t; L = length of the computational domain along x; ii = highest grid point number over domain L. Stability of the solution is guaranteed if, in the transformation of the solution from time level n t to time level (n+l) t by the scheme, none of the wave amplitudes kn are amplified. Substitution of Equation (4.12) into Equation (4.8) gives

x kn +1 e

i im k p ii

i im k p ü ì im k i -1p = x kn íre ii + (1 - r ) e ii ý î þ

or Ak =

WL | Delft Hydraulics

x kn +1 x kn

= 1- r + e

k im p ii

(4.13)

69

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

For a specific wave component the amplification |Ak| is defined as the length of the vector composed of the real and imaginary parts of Equation (4.13). From Figure (4.6) it can be concluded that |Ak| 1 for r 1.

Figure 4.6

Argand wave amplification diagram for the scheme of Equation (4.8)

Table 4.1

Example of unstable solution of Equation (4.8) for r=2

n 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

4.5

i => x (m) => t (s) 0 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000

0 0

1 100

2 200

3 300

4 400

5 500

6 600

7 700

8 800

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1

0 2 -4 6 -8 10 -12 14 -16 18 -20 22 -24 26 -28 30

0 0 4 -12 24 -40 60 -84 112 -144 180 -220 264 -312 364 -420

0 0 0 8 -32 80 -160 280 -448 672 -960 1320 -1760 2288 -2912 3640

0 0 0 0 16 -80 240 -560 1120 -2016 3360 -5280 7920 #### #### ####

0 0 0 0 0 32 -192 672 -1792 4032 -8064 #### #### #### #### ####

0 0 0 0 0 0 64 -448 1792 -5376 #### #### #### #### #### ####

0 0 0 0 0 0 0 128 -1024 4608 #### #### #### #### #### ####

9 10 900 1000 0 0 0 0 0 0 0 0 256 -2304 #### #### #### #### #### ####

0 0 0 0 0 0 0 0 0 512 -5120 #### #### #### #### ####

Characteristic solutions on a fixed grid

The accuracy of the scheme given by Equation (4.8) is rather low for many practical applications. It may be possible to choose r-values equal or close to unity on some parts of the grid, but varying velocities may make it necessary to work with lower values elsewhere. Reduction of x increases rapidly the number of operations and may demand much computer time.

WL | Delft Hydraulics

70

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Higher accuracy schemes introduce the need to define relations between more grid points. Programming wise this may lead to a rather complicated organisation. For a relatively clear and accurate procedure we will return now to the method of characteristics. For the fixed grid, as shown in Figure (4.7), one may follow the characteristic line from grid point (i,n+l) back to the time level n t.

Figure 4.7.

Characteristic solution on a fixed grid.

As the value for c remains constant along this characteristic, the accuracy of the solution technique depends on the accurate positioning of the intersection point at time level n t and on an accurate interpolation. For a constant velocity u and a linear interpolation between c-values at points (i-l,n) and (i,n) it is found that

cin +1

=

uDt n Dx - u Dt n ci -1 + ci Dx Dx

This relation is equivalent to the finite difference scheme given by Equation (4.8). The advantage of this method is that it can be extended to Courant numbers greater than unity (r>1; see Figure (4.7). An improvement of the accuracy is obtained by using a cubic spline function interpolation (Preissmann and Holly, 1977). In this interpolation not only the values for c are used at the grid points (i-l,n) and (i,n), but also their first derivatives in x. The interpolation is than given by n

n +1 i

c

= ac

n 1 i -1

n

æ ¶c ö æ ¶c ö + a c + a3 ç ÷ + a 4 ç ÷ è ¶x øi -1 è ¶x øi n 2 i

(4.14)

with

a1 = r 2 (3 - 2r ) a 2 = 1 - a1 a 3 = r 2 (1 - r )Dx

(4.15)

a 4 = - r (1 - r ) 2 Dx Preissmann and Holly describe the computation of c/ x values via an advection scheme. As under simplified conditions the advection does not alter the shape of the translated function, it must be possible to find an equation for the advection of both the first and the higher order derivatives. For the fractioned step approaches described later, this procedure, however, cannot be followed. A simpler, but also less accurate procedure, is found by approximating the derivatives from earlier computed c-values as

WL | Delft Hydraulics

71

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

¶c ¶x

@

cin+1 - cin-1 2Dx

UNESCO - IHE

(4.16)

The results obtained with this method of characteristics, using cubic spline functions for interpolation, show a rather low numerical diffusion. However, some parasitic waves still appear in the tail of the function.

4.6

Introducing diffusion

The derivation of the equation for advection-diffusion processes proceeds very similar to the one for advection alone. Referring again to Figure (4.1), the advection through the element boundary has to be increased with the contribution of molecular and turbulent exchange processes that can be described by Fick's law as

Tdispersion

= - DA

¶c ¶x

where D is a diffusion coefficient. The balance equation now reads

¶ ¶ ¶ ¶c ( Ac ) + ( Qc ) = - æç - DA ö÷ ¶t ¶x ¶x è ¶x ø or, assuming that the effects of spatial variation in cross-section and diffusion coefficients can be neglected, ¶c ¶c +u ¶t ¶x

= D

¶ 2c ¶x 2

(4.17)

The contribution of molecular diffusion to D is negligible in practical applications. Furthermore, in a schematization where u and t are averaged over the cross-section, the value of D is much larger than following from turbulent diffusion alone. Much of the spreading of the matter along the channel axis is caused by differential advection. Differences in advection are caused by velocity variations in the vertical profile, in the direction perpendicular to the flow lines, as a result of spiral flow in bends, due to storage elements along the channel banks etc. The combined effect of turbulent diffusion and differential advection is called dispersion, so that D in Equation (4.17) generally stands for a dispersion coefficient. Equation (4.17) has an exact solution defined as c ( x, t ) =

M0

2 A (p Dt )

½

e

-

( x -ut )2 4 Dt

(4.18)

where M0 is the total mass of a tracer which is injected instantaneously and uniformly distributed over the complete cross-section. The initial distribution is

WL | Delft Hydraulics

72

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

defined as a Dirac delta distribution, which tends to be dispersed along the channel axis as a Gaussian distribution with mean value =x+ut and with a standard deviation = 4Dt . For an initial Gaussian distribution the value of 0 can be transformed into a value, defined as the time which it would have taken to come from a delta distribution to the given Gaussian form. This leads to s 02 t0 = (4.19) 4D and the exact solution defined as c ( x, t ) =

M0

2 A {p D ( t0 + t )}

½

e

-

{x -( x0 +ut )}

2

4 D ( t0 + t )

(4.20)

where M0, in this case, is the total mass of the tracer integrated along the channel section at time t=0.

4.7

An explicit finite difference scheme

Exact solutions are defined only for very particular cases. The following limitations of the applicability can be mentioned: · ·

the solution is defined only for single point instantaneous sources ; the channel characteristics are often not constant along x.

More general, again, practical solutions are found with finite difference approximations, which will be demonstrate here, in the first instance, on a scheme for the diffusion part of the equation only, given by ¶c ¶t

= D

¶ 2c ¶x 2

(4.21)

Consider a grid, as shown in Figure (3.8).

Figure 4.8

WL | Delft Hydraulics

Operator for an explicit scheme for the diffusion equation.

73

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

The time derivative can be given as a forward finite difference approximation in the form ¶c ¶t

@

cin +1 - cin Dt

The second derivative in x is approximated as n

¶ 2c ¶x 2

=

¶ æ ¶c ö ç ÷ @ ¶x è ¶x ø

æ cin+1 - cin ö æ cin - cin-1 ö ç ÷-ç ÷ è Dx ø è Dx ø Dx

n

æ ¶c ö æ ¶c ö ç ÷ -ç ÷ è ¶x øi +½ è ¶x øi -½ Dx

=

@

cin-1 - 2cin + cin+1 Dx 2

Substitution of these finite difference approximations into Equation (4.21) gives the scheme, with cjn+1 expressed in terms of variables at time level n t,

cin +1

= rcin-1 + (1 - 2r )cin + rcin+1

where

r = D

Dt Dx 2

(4.22)

(4.23)

As the value at grid point (i,n+l) is computed explicitly from known values at time level n t, it is customary to call this an explicit scheme. At a later stage we will introduce a somewhat refined definition for the types of finite difference schemes. A Taylor's series development shows that, in reality, we are solving with the finite difference scheme (4.22) the equation ¶c ¶ 2c -D 2 ¶t ¶x

=

4 -Dt ¶ 2 c 1 2 ¶ c + D D x + h.o.t. ¶x 4 2 ¶t 2 12

An analysis similar to the one for the advection equation leads to the cancellation of the lowest order terms in the truncation error for r=1/6. However, for this same value of r, cancellation of the higher order terms does not occur. A complete analysis will lead to the conclusion that, for this scheme, it is impossible to find any r-value which will produce the exact solution. As an example, consider a channel section of length L, where initially a conservative pollutant with concentration c is present, given as a Gaussian distribution with standard deviation x. Selecting a distance step x=100 m and a dispersion coefficient D=0.463 m2/s, a time step t=l hour will give an r-value 1/6.

WL | Delft Hydraulics

74

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

Figure 4.9

UNESCO - IHE

Results with the numerical scheme of Equation (4.22) for various values of r.

Consistent with the theory of partial differential equations, the solution requires onepoint initial data at t=0 and one point boundary data at x=0 and x=L. In the case of this example, constant values c=0 have been set at both boundaries. For the second order space derivative of this parabolic partial differential equation, both boundary data could also have been given as first order derivatives of c along x. In Figure (4.9), a comparison is shown between the exact solution for this problem and numerical solutions obtained for different grid steps. With respect to the choice of the grid steps it is interesting to note that, unlike for the advection scheme, the highest accuracy is not found at the limit of stability. For practical applications it is best to choose a grid which gives r-values around 1/6. However, Figure (4.9) also shows that the differences are minor. This is certainly true when compared with the sensitivity of numerical parameters in advection schemes. In all cases, it is important to make sure that the stability limit r=½ is never exceeded. The existance of such limit can be investigated again with a Fourier series analysis, to give an amplification factor

A = 1 - 4r sin 2

kmDx 2

(4.24)

so that, since the sin2-term is definite positive, the stability condition |A| 1 implies that

r £ ½

(4.25)

or Dt £

WL | Delft Hydraulics

Dx 2 2D

(4.26)

75

Computational Hydraulics Version 2006-1

4.8

© 2005 Adri Verwey

UNESCO - IHE

Explicit schemes for the combination of advection and diffusion

For a first discussion, let us return to the advection scheme given by Equation (4.8). Comparing the Taylor's series expansion of this scheme shown by Equation (4.11), to the advection-diffusion Equation (4.17) it is seen that physical and numerical diffusion are equal when D = ½u Dx (1 - r )

(4.27)

In principle, advection-diffusion processes could be computed with a scheme for the advection equation alone, by selecting numerical parameters which just give the correct amount of diffusion by satisfying Equation (4.27). In practice this cannot easily be realised as u varies in general with x and t and the grid steps cannot be adjusted from one to the other grid point. Furthermore, the numerical diffusion coefficient has been derived from the idealized case where u and D are constant. It is, therefore, more practical to combine a diffusion scheme with a scheme for advection, which has a numerical diffusion which is small as compared to the physical diffusion. A convenient approach is the fractioned step method where the contributions to the variations in c by advection and dispersion are solved sequentially. As a first step, the equation ¶c ¶c +u ¶t advection ¶x

= 0

(4.28)

is solved, so as to give the c-values at (n+l) t resulting from transport alone, while as a second step the equation

æ ¶c ö ç ÷ è ¶t ø dispersion

=

¶c æ ¶c ö -ç ÷ ¶t è ¶t ø advection

= D

¶ 2c ¶x 2

(4.29)

is solved, giving the complete solution for advection and dispersion together. Equation (4.28) can, again, be solved with the fixed grid characteristic method using the spline function interpolation (see Equation (4.14)), giving intermediate values cj* defined by advection alone. Subsequently, the dispersion step can be given in the finite difference form as cin +1 - ci* Dt

@

ci*-1 - 2ci* + ci*+1 D Dx 2

or

cin +1

WL | Delft Hydraulics

= rci*-1 + (1 - 2r )ci* + rci*+1

(4.30)

76

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

The advection-diffusion equation can easily be extended to the two-dimensional case, giving ¶c ¶c ¶c +u +v ¶t ¶x ¶y

= Dx

¶ 2c ¶ 2c + D y ¶x 2 ¶y 2

(4.31)

where y is the direction in the horizontal plane orthogonal to x, v is the water velocity in y-direction and Dx and Dy are dispersion coefficients in x- and ydirection, respectively. This equation can be solved, following the same principles, in the sequence of fractioned steps: 1) 2) 3) 4)

advection along x; advection along y; dispersion along x; dispersion along y.

4.9

Implicit schemes for the diffusion equation

The stability criterion for the diffusion scheme of Equation (4.22) limits the time step of the solution method and this is not always desirable. For this reason implicit schemes can be introduced, which have less stringent stability conditions. As shown in Figure (4.10), in implicit schemes the unknown variables at the new time level (n+1) t are connected by the partial definition of the space derivative of Equation (4.21) at the new time level as cin +1 - cin Dt

{

}

D (1 - q ) ( cin-1 - 2cin + cin+1 ) + q ( cin-+11 - 2cin+1 + cin++11 ) (4.32) Dx 2

@

where defines a weighting of the space derivative between the time levels n t and (n+1) t. This equation can be rewritten to the form

a i cin-+11 + bi cin +1 + g i cin-+11 = d i where

(4.33)

ai = - q r bi = 1 + 2q r g i = -q r di

(4.34)

(1 - q ) r cin-1 + (1 - 2r (1 - q ) ) cin + (1 - q ) r cin+1

=

As an alternative, the implicit finite difference scheme may be written in terms of changes c from time level n t to (n+1) t as

Dci Dt

@

D Dx 2

{( c

n i -1

}

- 2cin + cin+1 ) + q ( Dci -1 - 2Dci + Dci +1 )

(4.35)

which can be written in a form similar to that of Equation (4.33).

WL | Delft Hydraulics

77

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

Figure 4.10

UNESCO - IHE

Grid for the implicit scheme of the diffusion equation

The system of Equations (4.33) can be solved with the double sweep algorithm, a special form of the Gauss elimination, based upon tri-diagonal materices. In this algorithm, in a first sweep the coefficients in the lower diagonal are eliminated, starting with the left hand boundary condition. In a second sweep values for the unknows are computed, by subsequently eliminating the upper diagonal coefficients.

Figure 4.11

Schematic presentation of the tri-diagonal double sweep matrix solution algorithm

Schematically, the algorithm is presented in Figure (4.10). Assuming the left hand and right hand boundary conditions, respectively, to be of the form

b 0 c0 + g 0 c1 a ii cii -1 + bii cii

= d0 = d ii

(4.36) (4.37)

Gaussian elimination of the lower diagonal leads to the set of equations

cin +1 + g i*cin-+11 = d i*

(4.38)

with the new coefficients in the matrix (now marked with a superscript *), computed from grid point 1 until ii-1 by the relations

g i* =

WL | Delft Hydraulics

g ; b - ag i*-1

d i* =

d i - ad i*-1 b - ag i*-1

(4.39)

78

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

For the more general case of non-linear equations, with coefficients derived as a function of x, also the coeficients , and in Equation (4.33) will get a subscript i. These subscripts will than also appear in the relations given by Equation (4.39). The elimination is started by assigning the starting values for the new coefficients at grid point 0 by substitution of the boundary relation given by Equation (4.36) into Equation (4.33) applied at gid point 1, giving

g 0* =

g0 ; b0

d 0* =

d0 b0

(4.40)

The same result would also have been obtained simply by setting the value of 0 to zero in Equation (4.39). In this case the assignment of new values for i* and i* could have started at grid point 0 instead of 1. After completion of the first sweep, the right hand boundary value is obtained by substitution of the boundary relation given by Equation (4.37) into Equation (4.39) at point ii-1. Successive back substitution (back sweep) of the values for c leads to the overall solution of the system of equations.

WL | Delft Hydraulics

79

Computational Hydraulics Version 2006-1

WL | Delft Hydraulics

© 2005 Adri Verwey

UNESCO - IHE

80

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

De Saint Venant equations and their solutions

5

5.1

Introduction

In Chapter 3 the unsteady channel flow equations were introduced for canals with simple rectangular cross-sections. These equations will now be extended to the more general case of irregular cross-sections. In this approach, the flow in rivers and other type of channels is assumed to be correctly defined by the state variables or dependent variables discharge, Q and water level, as a function of the independent variables t for time and x for space. Basic assumptions and/or limiting conditions are: · ·

·

the discharge is sufficiently well defined as the integral of the velocities through a cross-section, perpendicular to the x-axis and perpendicular to the flow velocity vectors in the flood plain; the water level is constant along the cross-section. This implies that at any time the water level at all points along a given cross-section should rise or fall at the same rate. This assumption is generally justified when the widths of river and flood plain are of the same scale and free of obstacles such as natural levees or embankments; the water level slope or gradient in the x-direction is constant along the crosssection. However, it should be noted that in case this condition is not satisfied, correction factors may be applied to reduce significantly errors in the model parameters, as discussed in this contribution.

Quantitative analysis of the two state variables requires two independent equations. Usually, the following equations are used: · ·

the continuity equation, based upon volume conservation in a control volume defined between two successive points along the channel axis; the momentum equation, based upon the conservation of momentum, including the effect of impulses generated by forces acting upon the water contained in the control volume.

Making, furthermore, the following assumptions: · · ·

WL | Delft Hydraulics

the pressure distribution in the vertical is hydrostatic; the resistance relationship for steady flow is also applicable for unsteady flow; and the bed slope is moderately steep, so that the cosine of the slope can be replaced by unity,

81

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

De Saint Venant (1871) derived the following equations (presented in a slightly adapted form here):

¶As ¶Q + = ql ¶t ¶x ¶z Q Q 1 ¶Q ¶ Q 2 + 2 ) }+ { + ( ¶x gA ¶t ¶x A K

(5.1)

= 0

(5.2)

where As is the cross-sectional area representative for storage over a control volume (m2), t the time (s), Q the discharge (m3/s), x the position along the channel axis (m), ql the lateral discharge per unit length of channel (m2/s), A the flow conveying crosssectional area (m2), the stage or water level above a selected horizontal reference plane (m) and K the channel conveyance (m3/s). The meaning of the parameters in these equations becomes clear when we follow the derivation of these equations.

5.2

The continuity equation

The continuity equation for channels with arbitrary cross-section can be derived on the basis of the same principles as applied to the case of the simplified rectangular profile, as described in Chapter 3. It should be realised that in channels and rivers with arbitrary cross-sections even the topography between two successive crosssections is irregular. For this reason, we will base the derivation of the continuity equation upon a storage volume V( ), as a function of stage .

Figure 5.1

Control volume for the derivation of the continuity equation

From Figure (5.1) it is readily seen that the following volume conservation equation holds:

¶V ¶Q + dx = ql dx ¶t ¶x As for each control volume between two successive points along x, the volume is a function of the water level, the following substitution can be made

WL | Delft Hydraulics

82

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

¶V ¶t

=

dV ¶V d V ¶t

=

Asurf

¶V ¶t

= b s dx

UNESCO - IHE

¶V ¶t

(5.3)

where Asurf is the storage area in the horizontal plane defined between two successive control sections. It is common practice, however, to define the storage as a cross-section property. For this reason, we usually transform the parameter Asurf to a parameter linked to a cross-section, by introducing a storage width bs of the channel. This level dependent channel cross-section storage parameter bs is defined as Asurf (for a given level) divided by the distance between two successive control sections (Equation (5.3)). This definition leads to the continuity equation of the commonly used form

bs

¶V ¶Q + ¶t ¶x

= ql

(5.4)

Integration of bs over the height of the cross-section leads to the definition of the storage cross-section As (V ) =

V

ò b dz s

(5.5)

zb

as introduced earlier in Equation (5.1). It should be realised that the storage crosssection parameters bs and As are often different from the parameters B and A used in the momentum equation (see § 5.3). These differences will be discussed in more detail in § 5.7. The second important parameter in the continuity equation is the lateral flow. Examples of lateral flow contributions, either positive or negative, are: · · · · · ·

WL | Delft Hydraulics

catchment runoff hydrographs, either given as point sources or as distributed sources. These consists both of surface water and ground-water components; drainage releases, for example, via pumps; irrigation water extraction; flow over side weirs; water exchange at the river bed, possible linked to a ground-water storage description; rainfall and evapotranspiration directly linked to the water surface.

83

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

5.3

UNESCO - IHE

The momentum equation

For the derivation of the momentum equation we will limit ourselves here to a pragmatic derivation. Recalling § 3.3, with the derivation of the momentum equation for a uniform rectangular cross-section, one of the intermediate forms has been

¶ ¶ ¶h (uh) + ( u 2h) + gh = 0 ¶t ¶x ¶x

(5.6)

Referring to Figure (5.2), the momentum equation of the overall cross-section is obtained by integrating each term of Equation over the width.

Figure 5.2

Cross-section for the integration of the momentum equation

Integration of the first component leads to the acceleration term B

¶ òo ¶t ( u j h j )dy =

B

¶ u j h j dy = ¶t òo

¶Q ¶t

Integration of the second component leads to the convective momentum term B

¶ 2 òo ¶x ( u i hi ) dy =

B

¶ ( u 2j h j ) dy = ò ¶x o

¶ ( b u Q) = ¶x

2 ¶ Q (b ) ¶x A

where is the Boussinesq coefficient of velocity distribution. It corrects for the fact that we have used the average flow velocity u in the integration of the momentum carried by the flow through the cross-section, instead of the local velocities. Consequently, the value of is greater than unity. It should be noted that the importance of this convective momentum term is generally limited and in most practical applications the value of ß is simply set to unity. In some practical applications the value of ß is even set to zero, which means that the complete term is neglected (see diffusive wave approximation of Chapter 6).

WL | Delft Hydraulics

84

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

Figure 5.3

UNESCO - IHE

Definition of the gravity component

Integration of the second component leads to the combined hydrostatic pressure and gravity term B

¶V ¶V òO g h j ¶x dy = g ¶x

B

ò h dy j

O

= gA

¶V ¶x

In the integration it has been assumed that the term (zb+h)/ x, or x, is constant across the section. In order to separate the effects of hydrostatic pressure gradients and gravity, the term can be split up again on the basis of a kind of representative depth hr above a representative bottom level zbr as follows

gA

¶V ¶x

¶h ö ¶h ö æ ¶z æ = gA ç br + r ÷ = gA ç I 0 r + r ÷ ¶x ø ¶x ø è è ¶x

(5.7)

It can be concluded from Equation (5.7) that the way in which we define this bottom level does not affect the results of computations. This definition is a product of our mind and not that of the river itself. The river feels for its flow only the effect of the water level slope (combined with cross-sectional parameters) and not that of a bed slope, as such. Finally, a friction term has to be added. This friction description is rather empirical and it is common practice to follow the same definition as applying to steady flow. In steady flow, the bottom friction for sub section j of the cross-section is supposed to balance the gravity term, or

( t bottom) j = r g h j I o Furthermore, t =t (u2) for turbulent flow, leading to the Chezy relation u j = Cj h j

1/ 2

Io

or to the Manning relation

WL | Delft Hydraulics

85

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

u j=

1 nj

UNESCO - IHE

hj2/3 I0

Proceeding with the Manning relationship and taking into account the side slope of the bottom profile, the local depth is generally replaced by a local hydraulic radius Rj, giving u j=

1 nj

R j 2/3 I0

For unsteady flow, the slope expressed by this shear force relation deviates from the bottom slope and for each sub section a friction slope may be introduced defined as If

=

2

u 2j A2j æ1 ö 2/3 çç Aj R j ÷÷ è nj ø

2

=

Qj

2

Kj

where Kj is the channel conveyance for subsection j. As per definition of a one-dimensional schematisation If is constant over the crosssection, integration of the friction term over this cross-section leads to a contribution to the momentum equation is given as

If

=

Q

2

K

2

(5.8)

with the total channel conveyance K defined as

K

=

jj

1

ån j =1

Aj R j 2 / 3

(5.9)

j

where jj is the number of sub sections in the channel cross-section. The concept of the total conveyance of a cross-section is based on the integration or summation of the contributions Kj for individual vertical slices of the cross-section. For each individual slice (Figure 5.2) the concept of the hydraulic radius is used where Aj (5.10) Rj= Pj with the wet perimeter Pj representing the contact length of water and channel bed.

WL | Delft Hydraulics

86

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

It should be remarked that the concept of a wet perimeter is based on the assumption that the shear force is equal at all points along this perimeter. This is approximately true for a small sub-section or slice of the channel cross-section. However, it is not true for the complete cross-section, as the local shear force is a linear function of the local depth. It is, therefore, strongly advised not to use the concept of the hydraulic radius for rather irregular cross-sections. For such sections one should always work with the concept of the conveyance, or discharge capacity K, obtained through the integration of the local conveyance across the section. In the special cases where the hydraulic radius concept can be used for the oveall cross-section, the Manning and Chezy coefficients are related through

Q = K I

=

1 AR 2 / 3 I n

(5.11)

Q = K I

= CAR 1/ 2 I

(5.12)

and

giving the relation between the Manning and Chezy coefficients as

C =

1 1/ 6 R ; or n = n

1 1/ 6 R C

(5.13)

It is quite well accepted that the Manning coefficient is more constant with depth than the Chezy coefficient. Especially at relatively low flow depth, such as occurring on flood plains, it is recommended to use the Manning coefficient for the friction description. The terms of the momentum equation can now be collected to give the form 2 Q|Q| ¶Q ¶ ¶h Q + (b ) + gA r + gAI 0 r + gA = 0 2 ¶t ¶x A ¶x K

(5.14)

Division of all terms by gA gives the equation in dimensionless form as

QQ 1 ¶Q ¶ Q 2 ¶h { + ( ) } + r + I0r + 2 = 0 ¶x 1424 gA ¶t ¶x A K3 kinematic wave 144 42444 3 diffusive wave 14444444 4244444444 3

(5.15)

dynamic wave

which is equivalent to Equation (5.2). Equation (5.12) also shows how the momentum term can be simplified when some parts of this equation are of lesser importance. These kinematic wave and diffusive wave approximations will be discussed in § 6.4.

WL | Delft Hydraulics

87

Computational Hydraulics Version 2006-1

5.4

© 2005 Adri Verwey

UNESCO - IHE

Numerical solutions

Hydraulic modelling techniques for 1D are currently based upon the numerical solution of the de Saint Venant equations, including the full convective momentum term. These numerical solutions are nearly exclusively based upon implicit finite difference methods. These offer the advantage of unconditional numerical stability, while the various robustness problems of the past, relating to non-linear effects and flooding and drying of channels and flood plains, have been solved satisfactorily. The application of finite elements and finite volume techniques does not provide specific advantages as, in 1D modelling, most of these techniques lead to equivalent forms derived through finite difference formulations. Numerical methods for the de Saint Venant equations may be based upon so-called staggered and non-staggered schemes. The first category represents formulations where the dependent variables Q and are defined alternatingly at successive grid points along the x-axis. For non-staggered schemes, however, the variables Q and are defined at the same grid points. At first sight this last definition offers advantages through the availability of the state variables discharge and water level at the same points along the channel axis. It has been shown, however, that the staggered grid approach offers distinct advantages over non-staggered grids by guaranteeing the convergence of numerical solutions and the better ability to handle flooding and drying of grid sections, as shown by Stelling et al. (1998). For the numerical solution of the de Saint Venant Equations (5.1) and (5.2), we will consider their Eulerian form per unit width of channel by first neglecting the lateral flow and further simplifying the equations to

¶V ¶ (uh) + =0 ¶t ¶x

(5.16)

uu ¶u ¶u ¶V +u +g + cf =0 h ¶t ¶x ¶x

(5.17)

and

where is the water level defined as = h+zb with h defined as the local water depth (m) and zb as the local bottom level (m), u the flow velocity (m/s) and cf the dimensionless bottom friction coefficient.

Figure 5.4

WL | Delft Hydraulics

Staggered grid for unsteady channel flow

88

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Referring to Figure 1 and to Stelling and Duinmeijer (2003) for further details, the staggered grid approach requires that alternatingly at - and u-points, Equations (5.16) and (5.17) are transformed into finite difference form (see also Abbott (1979), Cunge et al. (1986), Hirsch (1990) and Toro (1999)). Taking a -point as the only feasible choice for the transformation of the continuity equation, the finite difference form relates three successive unknowns uin-+1/1 2 , V in +1 and uin++1/1 2 defined at the time level (n+1) t to known values at time level n t as

where

V in +1 - V in * hin+1/ 2 uin++1/q 2 - *hin-1/ 2 uin-+1/q 2 + =0 Dt Dx

(5.18)

u n +q = (1 - q ) u n + q u n +1

(5.19)

and t x n i

= the time step along the t-axis (s); = the space step along the x-axis (m); = a superscript denoting the time step number along t; = a subscript denoting the space step number along x; = the time step weighting coefficient.

The symbol * in Equation (5.18) indicates that the value of h at this grid point has to be approximated by its nearest available upstream value along the grid. For a positive flow direction, for example, this means that hi+½ is approximated by the value of hi. Equation (5.18) can be reformulated as

a1i uin-+1/1 2 + b 1i V in +1 + g 1i uin++1/1 2 = d 1i

(5.20)

where 1i, 1i, 1i and 1i are the coefficients of the linearized implicit finite difference scheme. Equation (5.18) can also be rewritten to provide a choice of time step that guarantees the computation of positive water depths, as follows

Dt ö n æ n +q Dt hin +1 = ç1 - uin++1/q 2 hin-1 ÷ hi + ui -1/ 2 Dx ø Dx è

(5.21)

uin-+1/q 2 ³ 0 and

(5.22)

For

uin++1/q 2 ³ 0

Equation (5.21) shows that for these positive water velocities and for positive water depths, the velocity Courant number condition

WL | Delft Hydraulics

89

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

uin++1/q 2

Dt £ 1; Dx

or Dt £

UNESCO - IHE

Dx uin++1/q 2

(5.23)

can be applied to provide a sufficient condition for obtaining positive water depths at the new time level. This implies that for the time step limitation of Equation (5.23) newly computed water levels can never fall below the bottom of the channel. As a consequence, no artificial bottom slots or other arrangements are required to avoid numerical robustness problems for small water depths. A similar condition can be derived for negative flow directions. Following the same procedure as applied to the continuity equation, the momentum Equation (5.17) can now be defined at a velocity point, by the finite difference form uin+1/ 2 uin++1/1 2 uin++1/1 2 - uin+1/ 2 V in++1q - V in +q n n + a (u , u ) + g + cf =0 * n hi +1/ 2 Dt Dx

(5.24)

where the symbol * has the same meaning as in Equation (5.18) and a(un,un) is a generalization of the discretization of the convective momentum term. Appropriate formulations for a(un,un) follow from the physical conditions, as detailed by Stelling and Duinmeijer (2003). In the paper it is shown that the correct formulation of the convective momentum term depends on the way in which the convective speed of momentum is interpolated on the grid. As a rule, the discretization is based upon a transformation of the convective momentum term to

u

¶u 1 ì ¶ (uq) ¶q ü = í -u ý ¶x h î ¶x ¶x þ

(5.25)

and its discretization * * q -q ¶u 1 æ ui +1 q i +1 - ui q i u ; - ui +1/ 2 i +1 i ç ¶x hi +1/ 2 çè Dx Dx

ö ÷ ÷ ø

(5.26)

where

hi +1/ 2 =

hi + hi +1 q +q ; q i = i -1/ 2 i +1/ 2 2 2

and qi +1/ 2 = *hi +1/ 2ui +1/ 2

Also here the value of *h has the same meaning as in Equation (5.18) The values of * ui are missing at -points on the grid and are approximated by first order unwinding as

ui = ui -1/ 2 , if

*

qi -1/ 2 + qi +1/ 2 ³ 0 and 2

ui = ui +1/ 2 , if

*

qi -1/ 2 + qi +1/ 2 < 0 (5.27) 2

For positive flow direction this yields the simple expression

WL | Delft Hydraulics

90

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

_

q æu -u ¶u ö u ; _ i ç i +1/ 2 i -1/ 2 ÷ ¶x h Dx è ø i +1/ 2

(5.28)

This momentum conservative approximation is applied in cases of gradually varied flow or in flow expansions. Referring to Stelling and Duinmeijer (2003) again, this formulation should not be used in strong flow contractions, as this will breed energy. Instead, the energy head conservation formulation should be used, given by u

ö ¶u ¶V ¶ æ u2 +g =g ç +V ÷ ¶x ¶x ¶x è 2 g ø

(5.29)

and, for positive values of u, the energy head conservative upwind discretization -u ¶u ui2+1/ 2 - ui2-1/ 2 1 æu ö ; = ( ui -1/ 2 + ui +1/ 2 ) ç i +1/ 2 i -1/ 2 ÷ u ¶x Dx 2Dx 2 è ø

(5.30)

For negative flow velocity, this complete expression is shifted one grid point in the positive x-direction. Comparison of Equations (5.28) and (5.30) leads to the conclusion that the difference in obtaining momentum or energy conservation lies in the way in which the convective velocity of momentum is interpolated from the flow field. Terms in Equation (5.24) including (5.26) or (5.30) can be collected to give the generalized relation

a 2i V in +1 + b 2i uin++1/1 2 + g 2i V in++11 = d 2i

(5.31)

By successive elimination of all equations at u-points, the remaining set is given by

a i V in-+11 + bi V in +1 + g i V in++11 = d i

(5.32)

and is solved by applying elimination or conjugate gradient techniques, as discussed further down in this contribution. The numerical discretization is second order accurate for all terms of the de Saint Venant equations, except for the convective momentum term, which is first order accurate. For many practical applications the convective momentum term is of lesser importance and this lower discretization accuracy is quite acceptable then. However, nearly second order accuracy can be achieved by up winded second order extrapolation of u-values, combined with slope limiters, as described by Stelling and Duinmeijer (2003), again. In rapidly varied flow the convective momentum term becomes locally dominant. Modelling these types of flow requires an appropriate implementation of the convective momentum term. By doing so, Delft Hydraulics’SOBEK package has

WL | Delft Hydraulics

91

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

been enabled the highly accurate and robust modelling of phenomena such as supercritical flow in steep channels and moving hydraulic jumps. The numerical scheme given by Equations (5.20) and (5.31) only provides linear equations at internal grid points. At channel boundaries, additional equations are required. As the de Saint Venant equations are of second order hyperbolic type, one boundary condition is required for each of its characteristic lines entering the computational domain. Boundary condition requirements for 1D river models were discussed extensively by Cunge, Holly and Verwey (1986). Here we will limit ourselves by stating that, in most practical applications, inflowing discharges are specified at the upstream ends of channels entering the flood model domain and water levels or rating curves at channels leaving the model domain. At internal boundaries, such as channel junctions, usually a modified continuity equation is applied, jointly with water level compatibility at all channel boundaries at that junction.

5.5

Description of hydraulic structures

Nearly all flood simulation models are dealing with hydraulic structures, such as dams, weirs, bridges etc. The length of the hydraulic structure along the x-axis is usually supposed to be negligible on the scale of the river or channel and therefore the structure can be seen as one single point along x, where both water- and energy levels are discontinuous and the de Saint Venant equations do not apply. At this point a relation is established between the upstream water level, the discharge through the structure and the downstream water level. In this relation it is assumed that the upstream water level is taken at the nearest point just upstream of the structure, where vertical accelerations can still be neglected. Similarly, the downstream water level is taken at the nearest location just downstream of the structure, where the flow can be seen as nearly horizontal and the structure flow no longer contributes to the energy dissipation. Usually, it is also assumed that storage of water around the structure is negligible, so that the relationship is applicable at any moment in time. On this basis, there is a variety of ways in which the flow through hydraulic structures can be formulated: ·

in most cases hydraulic structure descriptions are based upon empirical laws relating the discharge through the structure to the upstream and downstream water level. The relation has the general form

Qcrest = Q (V up , V down )

(5.33)

where up is the water level upstream of the structure (m), Qcrest the discharge through the structure (m3/s) and down the water level downstream of the structure (m). This formulation includes the possibility of using energy levels instead of water levels. Equation (5.33) may be linearized to

a V up + b Qcrest + g V down = d

WL | Delft Hydraulics

(5.34)

92

Computational Hydraulics Version 2006-1

·

·

© 2005 Adri Verwey

UNESCO - IHE

where , , and are the local values of non-linear coefficients at a given state of the flow. A wide variety of structure descriptions is available from literature, e.g. the classical books of Chow (1959), Henderson (1966) and, more recently, Chanson (1999); as an alternative, the state of the flow in the section just upstream of the structure, where the flow is contracted and vertical accelerations occur, can be described by an energy conservation principle. Similarly, the state of the flow just downstream of the structure, where the flow expands with energy losses associated to it, can be described by a momentum and impulse balance relationship. By internal elimination of unknowns in these relationships, the water level at the structure crest or at its most contracted section can be eliminated and a relation similar to Equation (5.34) is obtained; in cases where it is difficult to define the state of the flow in terms of equations, laboratory experiments may be set up to define a matrix relating upstream and downstream water levels to the structure discharge. As an alternative, the structure may be modelled in detail by a 3D numerical code, leading to a similar set of matrix coefficients. Conditions are that a fine grid is used, and that the code is based upon an appropriate numerical description of convective momentum terms and the effect of turbulence.

In all cases shown above, the structure equations could be based upon energy levels, instead of water levels. Through suitable transformations, these relations can be reworked to the form of Equation (5.34) Similar descriptions can be made for local energy loss descriptions along a channel. Appropriate linearization of Equation (5.33) or of the other structure flow descriptions, leads to a numerical scheme of the form of Equation (5.31). In the total set of equations for channel flow this structure relationship replaces the momentum equation that is generally applied between successive -points. In the case of a closed structure or a discharge specified by a pump, the coefficients and of Equation (5.31) are both zero and the structure presents itself as an internal boundary condition with a discharge given. In the case of free flow in the positive x-direction, only the coefficient equals zero, showing also in the numerical relationship that the structure discharge is only dependent on the upstream water level. In literature an abundant number of structure equations can be found (Chow, 1959, Henderson, 1966 and, more recently, Chanson, 1999). Implementing these in a numerical code requires attention for discharge compatibility between the various flow states, especially at the transition between free and submerged flow. Incompatibility in the definition of the discharge may lead to oscillations in the computed flow, in particular at the moment of flow reversal. Problems may be suppressed by defining some inertia to the structure flow. This is implemented by adding a small term to the coefficients and . The numerical structure behaviour is generally better when applying an energy relationship in the accelerating flow section upstream of the structure combined with a momentumimpulse relationship for the section downstream of the structure.

WL | Delft Hydraulics

93

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Where appropriate, energy losses may be added to the description of flow in the contracted section. Internal elimination of local variables leads, again, to a relationship as given by Equation (5.31) . In all cases, the different flow states only affect the computation of the coefficients , , and and have no bearing on the solution algorithm of the resulting system of equations.

5.6

Topological model schematisation

Mathematical models for the simulation of flow in river- and channel networks are based upon the De Saint Venant equations, hydraulic structure equations and initialand boundary conditions. The modeller has to simplify the complete system to a choice of schematization objects which will give a satisfactory representation of the physical behaviour of the hydraulic system. Two processes are of primary importance: storage and conveyance. The actual model construction will start with the selection of the network branches: the topological schematisation. In this process it will be decided up to which level of detail network branches are to be included in the schematization. If too many minor channels are included, the amount of work in data preparation may be unnecessary high. In addition, model simulations may become too slow. The practice of good model construction requires a good insight into the relative importance of these contributions. Starting with conveyance, let us discuss the following issues: ·

·

·

what is the relative contribution of a channel to the overall conveyance of the system? For a good judgement, the Manning equation shows that if two channels give flow into the same direction, their relative importance is found by comparing their conveyance. Applying the Manning concept, Equation (5.11) shows, for example, that for similar roughness, a channel with half the width and depth of a reference channel, has only 16 % of the discharge capacity of that reference channel; can a minor channel form a short cut? Not only the conveyance of a channel is important. Certain channels may provide a short cut to the flow from one part of the network to the other. Even for small conveyance values the discharges may become substantial due to a steep water level slope. Special attention has to be given to a possible over bank flow at higher stages; is the water level always the same all over the channel cross-sections? If this is not the case it may be decided to split the channel up in two or more parallel branches and describe cross flow between the parallel channels via crosschannels or hydraulic structure descriptions. In general, such schematisations lead to significantly more complex model networks.

This last point is also important for the correct representation of storage in the model. Natural levees along rivers may delay the storage in the flood plain and hence, affect the celerity of flood wave propagation along the river and the dampening of flood wave peaks.

WL | Delft Hydraulics

94

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Other questions addressed are: · ·

is a storage area well connected to the network? This means, for example, that it can be fed with the correct conveyances by flow from all relevant directions; are there any barriers in the flood plain? If part of the flood plain is schematized as a storage element, it may contain various embankments or levees. Examples are levees raised around paddy fields. During floods, part of the described storage area may not get flooded at all.

Finally, lateral flows may have their impact on the topological schematisation as follows: ·

·

should we extend the schematisation to include part of a tributary? For example, how far does the backwater influence the tributary and does this effect justify the inclusion of part of the tributary in the topological schematisation? If a measuring station is available to provide the lateral flow time series, it may be quite far from the main river. In this case the schematisation may be extended to include this location; should we couple other models to our hydrodynamic model? For example, the exchange of surface water with the sub soil may be substantial or it may at least be important in relation to the objective of the model development. Rainfallrunoff models may become an integrated part of the total description. In flood wave propagation the role of infiltration of surface water into the sub soil may have to be included in the model.

5.7

Hydraulic model schematisation

Once decisions have been made about the channel network schematisation and its various model elements, choices have to be made about equations and their associated data. In most cases, channel flow will be based upon the description of the full De Saint Venant equations, although also mixed models may be constructed. It could make sense, for example, to describe the flow in upstream, steep river branches with hydrologic models or artificial neural networks (ANNs), whereas in the principal rivers the full De Saint Venant equations are applied. In all cases, however, the models must include representative values for channel storage and conveyance, whether these are specified directly or hidden in the parameters of a hydrologic model. In principle, the model must compensate for conveyance or storage which has been neglected in the topological schematisation. One important decision is the choice of distance step x, which should be based upon a good representation of the hydraulic processes. The choice is based upon: · · ·

WL | Delft Hydraulics

The wave length modelled. A rule of thumb is to model waves with at least 50 to 100 grid point along important wave components; The representation of the local variations in hydraulic parameters, such as sudden contractions and expansions; Placement of grid points closely around hydraulic structures or other discontinuites in the system.

95

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Channel conveyance is described at cross-sections. Under all circumstances these cross-sections must be lines perpendicular to the direction of mean velocities. In the case of meandering rivers these cross-section alignments have to be defined on the basis of sound engineering judgement. A complicating factor is the change of velocity vector directions with changing stage. For man-made channels the specification of cross-section parameters is rather straight forward. As the dimensions are usually quite accurately known, the most sensitive input parameter is the roughness coefficient. For simple channel geometries, the use of the hydraulic radius as an approximation of a representative depth is quite acceptable. The hydraulic radius is defined as the surface of the conveying cross-section divided by the wet perimeter (Equation (5.11) or (5.12)) and assumes a uniform distribution of the shear force along this wet perimeter. For cross-sections with varying depths this assumption is incorrect. The local shear force in the cross-section is a linear function of the local depth. For this reason, the conveyance of compound channels has to be based upon a summation of the conveyances of individual sub-sections where for each of them this shear force is more or less uniform (Equation (5.9)). An advantage of such integration is that each sub section can be given its own roughness value. Another advantage is that by keeping track of the individual contribution of the sub sections to the overall conveyance of the channel, computed discharges can be redistributed across the section to provide water velocities for each sub section. This ability may be important in simulating water quality or morphological processes in flood plains where the use of local velocities is required. In the integration, contributions of sub sections with small conveyance are without any loss of accuracy added to those of the other sub sections. There is no need to neglect the conveyance of shallow and highly resistant parts of the cross-sections as these may still give substantial contributions during floods or may be important in the overall assessment of the morphologic or water quality behaviour. One of the basic assumptions in a one-dimensional schematisation is the use of a constant water level slope all along the cross-section. In meandering channels, however, the slope may be quite different for different sub sections. One way to include this effect in the derivation of the conveyance is by giving each sub section a weight in the integration, as a function of the distance to its equivalent part in the next cross-section. Equation (5.11) will then be modified as follows

( DV ) K ½ ( Dx )

½

Q =

=

jj

1

ån j =1

j

( DV )

½

Aj R j

2/3

( Dx )

½

j

or jj

K

=

å j =1

WL | Delft Hydraulics

Aj R j 2 / 3

n j ( Dx j / Dx )

½

(5.35)

96

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

where x is the distance for which the cross-section parameters are representative along the x-axis and xj a similar value for sub section j. Storage parameters usually are given as a function of stage and linked to crosssections. For uniform channels, as often designed in irrigation and drainage systems, the storage width is equivalent to the flow width. Storage width data are simply extracted from cross-section information. In natural rivers, however, meandering and irregular flood plain topography requires more complex procedures. § 5.2 gives an extensive discussion of the correct definition of the storage width parameter bs. It follows that storage parameters have to be extracted from information provided by topographic maps, currently mostly available in the form of digital elevation models (DEMs) in GIS. For successive compartments along the river axis, approximate water level slopes have to be assumed to extract from these maps storage area as a function of stage. Division of these areas by the length of the compartment along the x-axis, provides the parameter values for bs. For models of lesser economic value, procedures may be simplified, if needed. Storage connected to channels which have been neglected in the topological schematisation must be added to the model as additional storage. In the cross-sections the compensation can be made in the form of additional storage width. It can also be introduced in the schematisation in the form of additional storage areas at nodes. It is recommended to keep track of the changes that this correction phase has given in parameter values derived during the primary phase of schematisation. It is always advised to keep good records of all the steps taken in the development of a model and the way the model parameters have been defined. Without such records it is practically impossible to introduce future improvements or modifications correctly. This is even more so when various persons are involved in the model development and/or the development takes place over a considerable length of time. Hydraulic structures may be described by their empirical relationships, by the application of an energy conservation principle upstream and a momentum conservation principle downstream or by specifying the discharge water level relationships in matrix form. In case of empirical relationships, care must be taken that the parameters describing free flow and submerged flow, successively, give a consistent computed discharge at the moment of transition from one flow state to the other. If applicable, additional entrance energy losses must be specified. Composite cross-sections of the structure require an adaptation of the topological schematisation by defining parallel structures. Similar to the definition of conveyance, the total discharge through the composite structure is defined by adding up the contributions of the individual structure components.

5.8

Boundary- and initial conditions

In relation to boundary conditions two issues are important: the location of a model boundary and the condition applied at this point. The location depends on the model objective. For models meant to study hydraulic problems, a strict rule for a boundary location is that the boundary data given at this point should not affect the outcome computations of various scenarios.

WL | Delft Hydraulics

97

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

In other words, the model changes introduced with these scenarios and the resulting changes of the hydraulic state at the boundaries of the models should not reflect back on the study area. Exceptions to this rule are rare, but may be justified in some cases of water quality modelling, for example. There are various ways in which this can be achieved. Water levels at a sea boundary of a river model are generally not affected by the river discharge. Usually, it is safe to use a sea boundary as a boundary location in a river model and provide water levels as boundary condition. In other cases it must be guaranteed that the downstream river boundary is located at a sufficient distance from the area studied. How far, depends on how an error, or a change in the downstream boundary condition, affects the area under study. In the first place, this depends on the magnitude of the error or change. In the second place this depends on the hydraulic characteristics of the river. A good measure for the required distance is the length of a backwater effect. In case of doubt it is recommended to do a sensitivity analysis on the effect of an error in the downstream boundary condition. At upstream boundaries the same reasoning applies. In this case the scenario applied should not reflect back on the data or condition applied at the upstream boundary. Also here, the required distance follows from an analysis of backwater effects. The type of boundary condition follows from the theory of characteristics, stating that one boundary condition has to be given for each characteristic entering the computational domain. In practice, however, this always means that one boundary condition is given at each boundary of the model. For sub critical flow this is usually a discharge hydrograph at an upstream end and water levels or a rating curve at a downstream end. A rating curve at the upstream end will usually lead to model instabilities as such condition will be nearly equivalent to the information contained in the channel topography. In this way, this boundary condition does not provide additional information that is not already available and the state of the flow remains undetermined. A rating curve conflicting with the topographic information leads mathematically to indeterminacy of the solution and hence, to instabilities. At downstream boundaries, usually water levels or rating curves are given. In river applications one could define this as a weak boundary condition, as it will be shown in Chapter 6 that the wave propagating in a river is usually an advection type phenomenon, only requiring an upstream boundary condition. This implies that for the use of the full hydrodynamic equations the downstream boundary condition becomes redundant. For similar reasons, it is sufficient in practice to give only one upstream boundary condition for super critical flow. The other condition is already implied in the relation between topographic data and the water level slopes following from this.

WL | Delft Hydraulics

98

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

At the downstream boundary, any condition given, will at worst lead to incorrect results over a small stretch of the downstream region. As a rule, the downstream boundary condition given will turn the outflow locally into a state of sub critical flow. Initial data are simplest supplied by starting up flow simulations on a dry bed or a low initial water depth. In the case of tidal models this specification is replaced by a constant mean sea level initially and zero discharges throughout. In tidal flow, usually two tidal cycles are sufficient to arrive at a consistent set of initial data, at least in terms of water levels. Most modelling systems provide the option of generating a set of consistent initial data and writing this state to a file to be used as a so-called hot-start file for further simulations.

5.9

Model calibration and validation

Calibration of a model is the process of removing, or at least reducing, the uncertainties in the choice of model parameter values. In principle, nearly all model data can be collected from the field, except for the roughness values. These parameter values can only be obtained indirectly by comparison of measured and computed variables, such as water levels and discharges. In principle, a model calibration could focus only on the optimization of roughness coefficient values. However, in practice there are uncertainties also in relation to cross-section data, flood plain topography, lateral flow data and inflowing hydrographs. This means that in a calibration also the correctness of storage parameters and measured hydrographs will have to be checked and possibly corrected. Part of the uncertainty is also caused by the schematization of the real system into the model. Neglected storage and conveyance and an inappropriate compensation for these missing effects also influence the calibration results. It is possible, then, that corrections are applied to the wrong parameters and still lead to reasonable calibration results. As will be shown in Chapter 6, flood wave propagation celerity and dampening both depend on storage and conveyance. An error in the storage could, therefore, be corrected through the channel resistance and still improve calibration results. It is, however, not certain that the same improvement is maintained under extrapolated conditions. For river channels, calibration of roughness data should proceed from lower discharges to higher discharges. Only after the roughness values in the main channel have been calibrated, those of the flood plain should be determined. Model validation is meant to provide confidence in the quality of a calibration. In principle, model validation should be applied on a case representing an extrapolation of the calibration events. Many models are developed for use beyond the range of data for which they could be calibrated. Flood models, for example, often are used for studies on flood frequencies of several hundreds of years and it cannot be expected that data of such events is available for model calibration. Moreover, one of the principal problems with model calibration is that results that have been measured during low frequency events, are not very reliable.

WL | Delft Hydraulics

99

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Rating curves for these low frequencies are rather unreliable as measurement conditions are difficult and the staff of catchment authorities generally have other worries than getting the most accurate discharge measurements during an extreme storm event. Only water levels have a reasonable chance of being reliable, also during extreme conditions. A possibility, therefore, is the extrapolation of rating curves on the basis of the use of local 2D models, based upon an accurately measured topography and calibration of roughness values during the higher frequency events. However, this approach is not common practice yet. On the other hand, Public Works in The Netherlands uses detailed 2D models of the flood plains of the principal rivers for their calibration of 1D models. These 1D models can than be used with more confidence in the range of extreme events. It should be realised that this procedure is feasible only if these 2D models are based upon very small size grid cells.

WL | Delft Hydraulics

100

Computational Hydraulics Version 2006-1

6 6.1

© 2005 Adri Verwey

UNESCO - IHE

Mathematical modelling of floods Introduction

All over history, floods have challenged scientists and engineers to master this phenomenon, both by developing mathematical descritions and by taking flood protection measures. Over the past decades, these efforts even have increased. In the first place, because a rapidly increasing world population has shown a strong drive to settle in flood prone coastal zones and river flood plains. In addition, there is an increasing awareness about climate changes which appear to bring more and more precipitation to river catchments. There also are numerous intellectual challenges, as scientific and technological developments in a wide range of fields open up many new ways of supporting flood related studies. In particular, important progress has been made in computer speed and in new measuring technologies. Although numerical methods for solving unsteady flow equations were developed already half a century ago, e.g. Arakawa (1966), it was only recently that numerical techniques behind flood simulation models have reached an acceptable level of perfection. Hydraulic engineers and hydrologists are now able to model flow in channels and over flood plains, irrespective of their bathymetric or topographic complexity; irrespective of the number and location of embankments for flood protection, roads and railways and irrespective of the number and complexity of hydraulic structures and the way we control these.

6.2

Flood model requirements

Flood simulation models may have different requirements, depending on their objective. Criteria for the selection of the appropriate tool are often based on: engineering staff time needed for model development, overall consultancy time for product delivery, speed of computation, completion time for a simulation, accuracy level of results, data requirements, numerical robustness, user-friendliness of the software and possibly others, depending on the objective of the model. These objectives may be related to flood risk analysis, flood forecasting, flood control and be based upon a variety of causes, such as storms, dam or dike breaks, hurricanes, typhoons or similar low atmospheric pressure phenomena. Recently (2004) attention has also been drawn once again to the devastating effects of tsunamis. All these application areas of numerical models have their own requirements, as will be discussed briefly in the sequel. In many countries, insurance companies are using flood risk maps, sometimes based upon relatively simple and quick estimates obtained via simple rules in GIS. In these cases it is assumed that the value of insured property does not justify the more complex laws defining the detailed flow of water.

WL | Delft Hydraulics

101

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

If the economic interest is greater, a first improvement is found by applying 1D (one-dimensional) steady flow models with GIS post processing to develop topography based flood frequency contour lines (e.g. FEMA procedures: www.floodmaps.fema.gov/fhm/). However, there is a tendency now to base flood risk analysis on more detailed combined 1D and 2D unsteady flow models for flood prone areas with valuable assets and a complex infrastructure. Federal and local governments have also become more aware of the potential of using such models for evacuation planning. In The Netherlands, for example, more than 60 % of the country is subject to flood risk and for most of these areas integrated 1D and 2D hydrodynamic models have been developed to study the effects of potential dike breaks and to provide guide lines to authorities in setting up evacuation plans. Besides producing flood depths, these models have to be capable of providing accurate estimates of flood wave propagation celerities over dry beds. Flood forecasting sets quite different requirements. Speed of producing a forecast is one of the most important criteria, especially in areas where flash floods occur. For this reason, numerical models behind a river catchment flood forecasting system are usually 1D hydrodynamic models, gradually replacing the simpler hydrological routing techniques. There is a tendency to include partly 2D hydrodynamic models, which is already common practice in flood forecasting systems for coastal areas and seas. Numerical models for flood forecasting are usually embedded in a flood forecasting platform, such as the Delft FEWS system (Werner et al., 2004), which has recently been installed in the UK to provide flood forecasts for nearly all river basins in the country. Important criteria for numerical models supporting flood control are accuracy, flexible schematization options, numerical robustness and consultancy time for model development and use. Currently, state-of-the-art for flood control is the use of combined 1D and 2D models (e.g. Hesselink, 2003). The former use of flood cells has been replaced by complete 2D flow descriptions, whereas sub-grid channel flow is still better described in 1D. Flood control models should be based upon reliable physical descriptions and schematizations, as part of their use is in extrapolation of calibrated models to extreme situations which have never occurred. One of the reasons to build models for flood control is the study of downstream impacts, especially cross border effects. Downstream impacts of flood control are changed flood wave celerity and changed flood peak attenuation. Higher flood wave celerities result from deepening of the river and the construction of embankments. This, in turn, leads to increased peak floods downstream. The construction of flood retention areas has opposite impacts and may be used to compensate the negative impacts. Model selection criteria then follow from the detail in which potential economic, environmental and social impacts have to be studied. The analysis of floods caused by dam- and dike breaks requires extremely robust numerical methods, especially for the description of flooding of dry areas and the correct propagation of the wave front. Moreover, model accuracy, partly based upon the ability to describe the full hydrodynamic equations, is important, as will be discussed in the section on software and model validation. As dam- and dike break simulations are nearly always made for the prediction of their potential effects, data for model calibration is rarely available.

WL | Delft Hydraulics

102

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

The quality of the model fully depends on its descriptive capabilities of the physical system in terms of topographic and roughness data, the representativeness of the equations and the numerical methods applied. However, it has to be kept in mind that the overall model accuracy also follows from the quality of the description of the dam failure mechanism and the assumptions made here. Floods generated by the passage of low atmospheric pressure zones such as hurricanes, typhoons and the geologically induced tsunamis require the modelling of 2D flow in coastal zones, seas and oceans and may set requirements such as the description of Coriolis forces, the use of spherical coordinates and curvilinear grids, the specification of moving atmospheric pressure fields, special ways of handling initial data etc. A possible integrated use of 1D, 2D and 3D models may provide advantages here.

6.3

The role of new data collection technologies

Models need good quality data if one wants to draw reliable conclusion from their use. In addition to developments in numerical methods, new technologies for the collection of data have recently led to complete changes in the selection of the type of numerical models. In particular, the development of GPS and DGPS technology has led to far cheaper methods of collecting bathymetric and topographic data (Moglen and Maidment, hsa025). In turn, this has led to the gradual replacement of the hydrograph based hydrologic models by the bathymetric and topographic based hydrodynamic or hydraulic models. Similarly, the collection of detailed digital terrain data in river and coastal flood plains has led to the replacement of 1D models by 2D models. Let us first consider the impact of LIDAR. The use of this laser technology, scanning the earth surface with laser beams from airplanes or helicopters, has provided the means to generate highly accurate digital elevation models at relative low cost. As an example, the whole area of The Netherlands has been remapped over the past years, with an accuracy of approximately 10 cm in the vertical at a density of 1 point per 16 m2. Total cost of this project was approximately 10 million euro, or approximately 250 euro per km2 (Verwey, 2001). This new topographic information has been essential for the flood risk analyses and the evacuation planning studies mentioned earlier. Also river bathymetries are obtained at relative low cost now by boats equipped with multibeam echo sounders. Also here, the position of the boat is recorded via DPGS, while the spatial sound signals record bottom depths relative to the boat at an accuracy of approximately 10 cm in the vertical. At a typical boat speed of 4 m/s and a swat width of the order of magnitude of the river depth, the bathymetry of quite large river beds can be obtained in a few days. Experiments are also being made with laser beams in the green range, as these are able to pass through clean water and enable the application of LIDAR technology also for river bathymetries. In similar ways other technologies are advancing rapidly enabling large amounts of data to be collected at relatively low cost.

WL | Delft Hydraulics

103

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Worth mentioning are the ADCP (Acoustic Doppler Current Profiler) technology for tidal discharge measurements and the increased precision of spatially distributed precipitation measurements with radar.

6.4

The nature of flood wave propagation

The set of Equations (5.1) and (5.11) forms the so-called kinematic wave approximation for flood propagation, if the slope I is taken as the bed slope. After substitution of Equation (5.11), or the kinematic wave part of Equation (5.15), into Equation (5.1) and neglecting the lateral flow term, a further simplication is achieved of the form

¶Q ¶Q +c =0 ¶t ¶x

(6.1)

with a flood wave celerity c (m/s) expressed as

c=

1 dQ bs dh

(6.2)

It is important to realise that, as an essential assumption in the derivation of this kinematic wave form, the channel bed slope has been used in the conveyance relationship. Equation (6.1) has the form of an advection equation which expresses that flood waves propagate with a celerity c which is inversely proportional to the available channel storage width bs (m) and a linear function of the derivative of the local flow rating curve. The characteristic celerity of this kinematic wave is lower than the celerity of the dynamic wave characteristic in the same direction. As discussed by Abbott (1979) it is this mechanism that leads to roll waves at flood wave fronts, limiting their propagation speed (see also Stoker, 1957). The first order partial differential equation (6.1) also expresses that along its characteristic celerity c the discharge remains constant, and so does the peak of the flood wave. In other words: there is no dampening effect of the flood peak. Though this is approximately true for rivers with steep slopes, the stretches with milder slopes require a lesser simplification of the De Saint Venant equations. Following, for example, Chaudhry (1993) and defining I of Equation (5.11) as the water level slope z, substitution of the last three terms of Equation (5.15) into Equation (5.1) gives the so-called diffusive wave approximation ¶Q ¶Q ¶ 2Q +c =D 2 ¶t ¶x ¶x

(6.3)

with a flood wave diffusion coefficient D (m2/s) derived as D=

WL | Delft Hydraulics

K 2bs I

(6.4)

104

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Including the lateral flow term the diffusive wave approximation reads ¶Q ¶Q ¶ 2Q +c = D 2 + cq ¶t ¶x ¶x

(6.5)

Returning to Equation (5.15) again, the variable I represents the bed slope in the case of the kinematic wave approximation and the water level slope in the case of the diffusive wave approximation. The description based upon the full set of Equations (5.1) and (5.15) is defined as the full dynamic wave description.

6.5 6.5.1

Deformation of flood waves The role of varying celerities

Returning to Equation (6.2), it is readily seen that for constant celerity at all water levels, there is no deformation of the flood waves. The waves are just translated along the river axis. However, in practice the wave celerity c generally increases with increasing river discharge, especially when the river banks are steep (Figure 6.1a). Moreover, when the flood waves arrive at bank level, the storage width often increases more rapidly than the increase in dQ/dh, resulting in a temporary dip in this celerity function (Figure 6.1b). When the water level rises further, the celerities start increasing again.

Figure 6.1

Flood wave celerity for various cross-section shapes

The higher celerities for increasing discharges cause a steepening of the wave at the front and a stretching of the falling limb of the flood hydrograph (Figure 6.2). This is the most common form of flood wave deformation in natural river systems.

Figure 6.2

WL | Delft Hydraulics

Flood wave deformations resulting from celerity varying with stage

105

Computational Hydraulics Version 2006-1

6.5.2

© 2005 Adri Verwey

UNESCO - IHE

The role of the diffusion term

Equation (6.3) shows that along the characteristic line c=dx/dt the flood peaks are dampened by the integral of the diffusion term. As the diffusive wave approximation of the De Saint Venant equations may include other effects as well (e.g. the form of Equation (6.5)), extending the number of terms of the right hand side, we will discuss these various contributions in a fractioned step approach. Returning now to the right hand side terms of Equation (6.3), the influence of the diffusion fraction can be written as (

2 dQ ¶Q )1 = D 2 ¶x dt

(6.6)

which gives the equation for flood peak attenuation as æ ¶ 2Q ö dQ1 = D ç 2 ÷ dt è ¶ x øpeak

(6.7)

The diffusion mechanism is commonly explained by the fact that the water level slope during rising water level is steeper than the bed slope. Consequently, the discharge at the wave front is higher than might be expected from the rating curve and the water balance in that section forces a lowering of the wave peak. This is further enforced by the lower discharges at the rear of the wave as compared to the discharge which would follow from the stage-discharge relationship. Here, the water level slopes are smaller than the bed slope (Figure 6.3).

Figure 6.3

Illustration of the effect of varying water level slopes on flood wave deformation

Returning to Equation (6.7), it is readily seen that this fraction gives a dampening of the wave peak, as the term ¶2Q/¶x2 is negative, while it leads to an increase of the discharge at the front and at the rear of the wave, where ¶2Q/¶x2 is positive. Observation of the diffusion coefficient given by Equation (6.4) leads to an interesting paradox. It is a well known phenomenon that an increase in storage also increases the dampening of flood waves. This, however, does not follow immediately from Equation (6.4), which shows a decrease in the wave diffusion coefficient with increasing storage width of the river. It should be realised that an increase in storage width also decreases the flood wave celerity, leading to a proportional decrease in the flood wave length. Consequently, the contribution of the

WL | Delft Hydraulics

106

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

term ¶2Q/¶x2 shows a quadratic increase in the wave dampening. Furthermore, the decrease in the wave celerity increases the travel time, or the time over which Equation (6.7) is integrated, over a given stretch of the river, compensating for the decrease in the diffusion coefficient. Taking all these effects into account, an increase in the storage width by a factor two, gives a total increase in wave attenuation of a factor four.

6.5.3

The role of the lateral flow terms

Including also the effect of lateral flow, Equation (6.5) shows its influence. Applying, again, the fractioned step approach, the contribution of the lateral flow to the change in peak discharge is given by

(

dQ ) =cq dt 2

(6.8)

or

dQ 2

= cq peak dt = q peak dx

(6.9)

Equation (6.9) simply tells us that the peak value of the flood wave is increased by the total lateral flow added to the flood as its peak travels along the river. Interesting conclusions can now be drawn on the case where flood water infiltrates into the soil as the flood passes the flood plain. As the infiltration is a loss of water, the lateral flow is negative and leads to a reduction of discharges along the paths given by the characteristics. A typical infiltration function along the flood wave would show the highest water losses during the passage of the front of the wave, which shows up as a pronounced effect when the bed material is highly permeable. A high porosity of the bed material, gravel for example, causes a rapid increase in the infiltration when water is made available by the wave. During the further rise of the water level, the stream width increases and increases the storage volume available for infiltration. At the passage of the peak, the infiltration is more or less completed and the function drops down to zero. A typical response of this infiltration is a pronounced steepening up of the flood wave front. This is why such floods may have a very short lead time for warning the population and may be quite destructive. A typical example is given by floods in wadis. However, also in more moderate climates this effect is known. Note that it is necessary in such cases to include in a model the effect of exchange with the ground-water, as the typical wave deformation found in this case cannot be calibrated by modifying friction parameter values.

6.6

Link to hydrologic flood routing models

Equations (6.1) and (6.3) are rarely used in discretized form anymore. However, they are useful for the flood modeler as to provide insight into the physical nature of flood wave propagation. Moreover, they link hydrologic and hydraulic flood routing

WL | Delft Hydraulics

107

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

techniques. By applying a Taylor’s series expansion to the Muskingum equation, Cunge (1969) has demonstrated that the well-known hydrological Muskingum flood routing technique emulates the solution of the advection-diffusion Equation (6.3). This insight provides guidelines for a suitable choice of the Muskingum parameters on the basis of the expressions for c and D (Equations (6.2) and (6.4), respectively). The numerical Muskingum method belongs to the class of hydrological routing techniques. These methods are based upon the notion that flood wave propagation characteristics can be derived from measured flood hydrographs along the river, rather than from detailed topographic information, as discussed in the WMO report on flood forecasting models by Serban et al. (2005). The basis of these methods is that flood wave propagation and diffusion behavior can be represented by a limited number of parameters, which can be calibrated from observed hydrographs. Usually there are only two parameters, as in the case of the Muskingum routing method. The limitation of this assumption is clearly shown by analyzing the stage dependent expressions for celerity and diffusion given by Equations (6.2) and (6.4). Although the current methods for flood routing are far more precise, there are still many situations where the use of hydrological flood routing techniques is still justified. In many river catchments, and especially in tributaries, the collection of detailed information on river bathymetry and flood plain topography is economically not always justified, despite the emergence of relatively cheap new technologies. An alternative to hydrological forecasting methods is provided by newly developed artificial neural network (ANN) concepts (Minns & Hall, hsa018). Also this technology relies on measured hydrographs to derive relationships between inflowing and out flowing hydrographs, though in practice ANN’s are mostly applied to the development of relationships between rainfall and river catchment runoff. Both hydrological routing and ANN techniques provide limited reliability for the range of events outside those used for calibration. However, this is exactly the range one is interested in when dealing with extreme flood events, which rarely occur and for which measurements are even more rarely available. So even though ANNs and other data-driven modelling techniques may prove quite valuable in e.g. determining rainfall-runoff relations (Solomatine, hsa021), physically based descriptions, such as provided by hydraulic routing techniques and based upon the full use of Equations (5.1) and (5.2), offer better extrapolation possibilities than hydrological routing methods or ANN based models. For this reason, we will focus on hydraulic modelling techniques in the sequel. However, it should be realized that simplified equations, such as Equation (6.5), remain useful, as these offer us a good insight into the physical nature of flood propagation, in particular the concepts of flood peak arrival time and the attenuation of peak discharges.

6.7

Two-dimensional modelling of floods

In the modelling of floods, flows often take short cuts through flood plains where the 1D description may become quite inaccurate. This is even more the case for dam- or embankment failures, where the flow may leave the flood plain completely and inundate natural terrains. For this reason the two-dimensional (2D) shallow

WL | Delft Hydraulics

108

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

water equations will be introduced. Following the same principles as for 1D flow, these equations read ¶h ¶ (uh) ¶ (vh) + + =0 ¶t ¶x ¶y

(6.10)

u u 2 + v2 ¶ (h + zb ) ¶u ¶u ¶u +u +v +g + cf =0 h ¶t ¶x ¶y ¶x

(6.11)

v u 2 + v2 ¶ ( h + zb ) ¶v ¶v ¶v +u +v +g + cf =0 ¶t ¶x ¶y ¶y h

(6.12)

where we now also introduce the y-axis, orthogonal to the x-axis, with its flow velocity v (m/s) associated to it. The friction term, with the dimensionless friction coefficient cf, has in both momentum equations a shear force component derived from the quadratic head loss description along a stream line of the 2D flow. Basic assumptions are similar to those given for the 1D equations, as far as applicable in this form of schematization.

Figure 6.4

Staggered grid for 2D flow simulations

Referring to the 2D grid shown in Figure 2, a volume conservative finite difference form of the continuity equation is given by

V in, +j 1 - V in, j * hin+1/ 2, j uin++1/q 2, j - *hin-1/ 2, j uin-+1/q 2, j + + Dt Dx * n hi , j +1/ 2 vin, +j +q1/ 2 - *hin, j -1/ 2 vin, +j -q1/ 2 =0 Dy

WL | Delft Hydraulics

(6.13)

109

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

uin++1/1 2, j - uin+1/ 2, j Dt

+ a11 ( u n , u n ) + a12 ( v n , u n ) + g 2

u cf

Dt

(

æ n ö n ç ui +1/ 2, j ÷ + vi +1/ 2, j è ø * n hi +1/ 2, j

n +1 i +1/ 2, j

vin, +j 1+1/ 2 - vin, j +1/ 2

UNESCO - IHE

+ a21 ( u , v n

(

n +1 i , j +1/ 2

v cf

n

n i , j +1/ 2

u

)

V in++1,qj - V in, +j q + Dx (6.14)

2

=0

)

V in, +j +q1 - V in, +j q + a22 ( v , v ) + g + Dy

)

æ ö + ç vin, j +1/ 2 ÷ è ø

2

n

* n i , j +1/ 2

h

n

(6.15)

2

=0

where the symbol * has the same meaning as in Equation (5.18) again and a11(un,un), a12(vn,un), a21(un,vn) and a22(vn,vn) are generalizations of the discretization of the convective momentum term. The long double bar over the velocity in the friction term means that this velocity is obtained by averaging over values at four surrounding grid points. The friction term requires special treatment in case of flooding of dry terrain. At the wave front the water velocity rapidly accelerates from zero. Overshoot of velocities can be prevented by a predictor corrector approach. The convective momentum terms are subject to the same principles as discussed for the 1D approximations. For example, for positive flow velocities the momentum conservative discretization of the term a12(un,vn) is given by x

a12 (v , u ) ; n

n

v

q1+1/ 2, j -1/ 2 æ ui +1/ 2, j - ui +1/ 2, j -1 ö ç ÷ x Dy ø hi +1/ 2, j è

(6.16)

whereas it is given by x æ ui +1/ 2, j - ui +1/ 2, j -1 ö a12 (v n , u n ) ; v i +1/ 2, j -1/ 2 ç ÷ Dy è ø

(6.17) x

for the energy conservative discretization. In the first expression v q means the specific discharge in y-direction, averaged over two surrounding points along the x

local x-axis. In the last expression v has the same meaning in relation to the velocity v. The treatment of the convective momentum terms shown above is numerically very robust and allows for the correct description of the effects of sudden expansions and contractions and similar changes in the topography, such as steps in the bed level. Moreover, it allows for the 2D simulation of supercritical flows and the propagation of hydraulic jumps.

WL | Delft Hydraulics

110

Computational Hydraulics Version 2006-1

Figure 6.5

Figure 6.6

© 2005 Adri Verwey

UNESCO - IHE

Layout of the instantaneous dam break experiment

Top view of results of the numerical scheme of Equations (6.13), (6.14) and (6.15) (lower half), compared with video monitored measurements in a physical model (upper half)

As described by Stelling and Duinmeijer (2003), the correct modelling of these phenomena has been demonstrated in a software validation study, where results of the numerical scheme of Equations (6.13), (6.14) and (6.15) were compared with video monitored measurements in a physical model (Figure 6.6). The set-up consists of two reservoirs with different water levels, separated by a wall. The wall contains a gate which can be lifted. The width of both reservoirs is 8.30 m, the length of the

WL | Delft Hydraulics

111

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

upper reservoir 2 m and the length of the lower reservoir 29 m. The gate has a width of 0.4 m and has been placed at the middle of the wall. The numerical experiment was made on a grid of 0.1*0.1 m, giving a total of approximately 25000 grid cells. The time step was set to 0.005 s. The gate was lifted at a speed of 0.16 m/s, to produce a flow spreading out into a 2-dimensional plain. Initial data were set at a depth of 0.60 m for the upstream reservoir and a depth of 0.05 m for the downstream reservoir.

Figure 6.7

Comparison of measured and computed wave front position at various times (see case Figure 4)

Figure 4 shows the results of this simulation. The upper half of this figure presents a video recorded view from above. The lower part of the figure presents the computed results. It is clearly seen that the front propagation, the propagation of the hydraulic jump and the side spreading of the wave are represented reasonably well. Figure 5 shows a comparison of the measured and computed position of the wave front at various times. Simulations were made with various Manning roughness coefficients and both for an initially wet and a dry downstream reservoir. For the propagation of the flood on the dry bed the Manning roughness turned out to be a sensitive parameter. If for dam break models a reliable topography is available, the roughness parameter remains the only parameter to be estimated. Currently, researchers are focusing on better descriptions of roughness parameters by deriving depth dependent relationships on

WL | Delft Hydraulics

112

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

the basis of vegetation characteristics, e.g. Uittenbogaard (2003), Rodriguez (2004) and Baptist (2005).

6.8

Integrated 1D/2D modelling

In flood modelling, there are numerous practical examples where flows are best described by combinations of 1D and 2D schematizations. An obvious example is the flooding of deltaic areas, often characterized by a flat topography with complex networks of natural levees, polder dikes, drainage channels, elevated roads and railways and a large variety of hydraulic structures. Flow over the terrain is best described by the 2D equations, whereas channel flow and the role of hydraulic structures are satisfactorily described in 1D. Flow over higher elevated line elements, such as roads and embankments can be reasonably modeled in 2D by raising the bottom of computational cells to embankment level. However, for a higher accuracy of the numerical description adapted formulations have to be applied, such as energy conservation upstream of overtopped embankments. Another example is the flood propagation in a meandering river, with shortcuts via the flood plain when over bank flow occurs. In large scale models, the flow between the river banks is satisfactorily described by the de Saint Venant equations solved with 1D grid steps several times the width of the channel. An equivalent accuracy of description in 2D would require a large number of grid cells, with step sizes being a fraction of the channel width. However, flow in the flood plain may be better described in 2D and may allow for 2D grid steps often exceeding the width of the river. For this reason, hybrid 1D and 2D schematizations are often used. Basically there are two approaches: one with interfaces defined between 1D and 2D along vertical planes and the other approach with schematization interfaces in almost horizontal planes. Coupling along vertical planes, gives a full separation in the horizontal space of the 1D and 2D modeled domains. In the 1D domain the flow is modeled with the de Saint Venant equations applied over the full water depth. The direction of flow in the 1D domain is assumed to follow the channel x-axis and in the model it carries its momentum in this direction, also above bank level. Without special provisions, there is no momentum transfer accounting applied between the 1D and 2D domains. Momentum and volume entering or leaving the 2D domain at these interfaces, are generated by the compatibility condition applied. As a result, the coupling cannot be expected to be momentum conservative. Depending on the numerical solution applied, the linkage may either be on water level or on discharge compatibility. Particular care has to be taken in applying this form of schematization if water quality processes are to be included in the model. In a model coupled along an almost horizontal plane, 2D grid cells are placed above the 1D domain, as shown in Figure 6. In this schematization, the de Saint Venant equations are applied only up to bank level. Above this level, the flow description in the 2D cell takes over. For relatively small channel widths compared to the 2D cell size, errors in neglecting the effect of momentum transfer at the interface are minor. For wider channels it is recommended to modify each 2D cell depth used in the

WL | Delft Hydraulics

113

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

momentum equation by adding a layer defined by the local hydraulic radius for that part of the 1D cross-section which underlies a 2D cell. Further refinements are possible, including terms describing the momentum transfer between the 1D and 2D domains. Numerical solutions are obtained by discretizing separately the 1D and 2D domains. Assuming that for both domains implicit numerical schemes are applied, the interface compatibility conditions can be modeled either as an explicit or an implicit link. Applying explicit links, first the solutions for the 1D and 2D domains are generated sequentially. Subsequently, exchange flows are computed and added as lateral flows at the next time step. Implicit links are based upon water level compatibility. These equations are then added to the complete sets of equations generated separately for the 1D and 2D domains. There are many approaches to solving the complete set of equations. With the current state of the art, it is no longer necessary to apply for the 1D domain different solvers for so-called simply or multiply connected channel networks. Similarly, in 2D there is no real need anymore for alternating direction algorithms, as the efficiency of the conjugate gradient solvers has increased significantly over the past years.

Figure 6.8

Coupling of 1D and 2D domains in SOBEK

As an example, Delft Hydraulics has developed its combined 1D2D package SOBEK for the modelling of integrated fresh water systems (www.sobek.nl; www.wldelft.nl). The 1D part of hybrid models is based upon the numerical scheme of Equation (5.20) and (5.31). The 2D part is described by the single step 2D scheme given by Equations (6.13), (6.14) and (6.15). For efficiency reasons, the continuity equations for the 1D and 2D domains are combined into one single equation at points where 1D grid sections underlie a 2D cell. As a first step in reducing the total number of equations, SOBEK eliminates all equations at velocity grid points. The second step in the solution algorithm is the elimination of a large number of unknowns by applying a minimum connection search between unknown water levels. As a rule, this leads to an efficient elimination of nearly all unknowns of the 1D domain and a substantial number of unknowns in the 2D domain.

WL | Delft Hydraulics

114

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

This direct solver carries its elimination on, until nearly every second equation in the 2D domain has been eliminated. Beyond this point, it is more economical to apply the conjugate gradient solver to solve the remaining set of equations. Apart from its efficiency, an additional advantage of eliminating nearly every second 2D equation is the improved conditioning of the resulting matrix. This follows from the fact that elimination of an unknown water level at a 2D grid point has the effect of increasing the spatial distance between the remaining adjacent points, where water levels are still unknown. This, in turn, reduces Courant numbers and as a consequence, it leads to changed coefficients at the main diagonal of the matrix which is now more dominant in relation to the other diagonals, e.g. Verwey (1994).

Figure 6.9

Flood modelling of the Vallei and Eem area, The Netherlands

An example of a combined 1D2D model is shown in Figure 7. It represents the schematization of a model of the Eem Valley area in The Netherlands applied in a study of the potential effects of a River Rhine dike breach. This model has been used to provide information on warning lead times and flood depths for evacuation planning. The Rhine branch upstream of the breach has been modeled in 1D. At its upstream end a design hydrograph was specified, whereas the downstream boundary condition of this relatively short branch is given by a rating curve. Such a short downstream reach is permissible as the rating curve automatically corrects for most of the effects of flow deviated through the breach at this boundary. The breach itself has been described in 1D as a structure with a velocity dependent breach growth.

WL | Delft Hydraulics

115

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

North of the dike the 1D link discharges into the 2D domain, given by a 100 * 100 m grid with bottom levels derived from a digital elevation model and resistance coefficients derived from land use maps. Elevated roads and railways are presented as flow barriers by raising the underlying cell bottom levels up to the levels of these embankments. The resulting flood depths presented in Figure 7 clearly show the effect of the 1D channel in the schematization. Due to their greater depth, flood waves propagate faster in these channels than over land. Further downstream, this leads to first signs of the progressing flood wave already one or two days before the main flood arrives.

6.9

Exercise

At a discharge measuring station a cross-section is measured as shown below. a) b) c) d)

e)

Compute the conveyance of the channel for levels h=2.00 m and h=3.00 m, respectively. The Manning numbers vary along the cross-section and their values have been indicated in Figure 6.8.. Compute discharges for both levels given an approximate bed slope I0=10-3. Compute approximate flood wave celerities for the given water levels. Explain the difference in the celerities. During the passage of a flood wave, the water level rises from 2.00 m to 3.00 m above the reference level over a period of 10 minutes. Compute the approximate deviation between the water level slope and the bed slope. (Note: transform the time to a length via the wave celerity). Compute the approximate difference in percentage between the expected measured discharges and the discharges based on a rating curve.

Figure 6.10

WL | Delft Hydraulics

Cross-section for the computation of flood wave propagation celerities

116

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Water Hammer

7 7.1

Introduction

Water hammer is the phenomenon of high-amplitude pressure waves travelling in pipes, caused by the rapid changes in the velocity of the transported fluid. These variations usually result from the operation of flow- and pressure control devices, such as valves, pumps and turbines. d slot

D

Figure 7.1

Sketch of pipe with artificial slot at the top to represent storage effects of pipe expansion and water compressibility

In principle, the equations used for the computation of these pressure waves are the same as those derived for open channel flow. As an analogy one may consider a pipe which has an artificial slot at the top, in which a free surface water level may rise after complete filling of the pipe (Figure 7.1). The wetted area of the slot represents the storage of water resulting from changes in the water pressure. The storage capacity consists of three principal contributions: 1. the elasticity of the pipe wall which leads to a pipe cross-section expansion at

increasing water pressure. It is evident that this type of storage depends on the cross-section shape and on the wall material properties and its composition. 2. compressibility of the water, which usually is neglected in free surface flow computations. In the case of water hammer, however, the compressibility plays a dominant role. As the free surface flow equations are based upon a volume balance, the compression of the water represents a virtual water storage; 3. compressibility of air bubbles or vacuum bubbles contained in the water; Water hammer is important in engineering for the following reasons: 1. high pressures may build up in a pipe line to the extent that it may lead to

bursting of the pipe;

WL | Delft Hydraulics

117

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

2. vacuum bubbles may be generated when the pressure drops to values below

approximately 30% of the atmospheric pressure. The resulting cavitation is highly aggressive to the pipe material and may lead to considerable damage; 3. in other cases, air may be entrained in the fluid at moments and locations where the pressures are low. For these reasons, any design related to transport of fluid in pipes should be checked against the risk of water hammer damage. There are various reasons why water hammer may occur. For details, reference is made to text books, e.g. Chaudhry (1987). Typical reasons are: ·

· · ·

instantaneous or very fast valve closure. Many of us know this phenomenon when closing a tab at home. A wrong design of the valve may cause a sudden deceleration of the water flow and a high pressure builds up due to the transformation of momentum of the water into an impulse; sudden energy demand changes in a hydropower station or turbine failure; pump failure, for example, by a power cut. At the upstream end of the pump the pressure will increase, while at the downstream end a negative wave is generated which may cause cavitation or implosion of the pipe; burst of a pipe line, which in turn generates pressure surges along the pipe.

In many cases it is the way of operating the system that causes the water hammer problems. In other cases there may be accidental causes. In the design of the hydraulic system such operation problems may be foreseen and water hammer may be prevented by designing anti-water hammer arrangements. Examples are: ·

· ·

· · ·

strict control of valve manipulations. As will be shown by the computational procedure, slow closure of valves reduce the over pressures. Slow has to be seen in relation to the length of the pipe and the celerity with which pressure surges are travelling along the pipe; design and construction of a surge chamber as closely as possible upstream of the turbine or pump. This surge chamber serves as an escape for the pressure wave and reduces effectively the dangerous operation time of closure; closed and vented air vessels, which have a function similar to surge chambers. However, in this case the surface is formed by a pressurized air chamber. These devices are used, for example, at the downstream end of sewer pumps in order to prevent excessive under pressures; fly wheels, which prevent the sudden change of the pump rotational speed; by-passes with a check-valve placed in parallel to a pump. It opens during pump rundown and supplies water to the pipe downstream of the pump; pressure release valves, which open when the water pressure exceeds a maximum value admitted.

For a more extensive list of options reference is made again to the literature in this field. As water hammer computations may be complex, use is often made of standard software packages, such as the WANDA system of Delft Hydraulics (www.wldelft.nl). However, simple problems may also be investigated by using Excel. The next paragraphs provide the basis for setting up such computations.

WL | Delft Hydraulics

118

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

7.2

UNESCO - IHE

Water Hammer Equations

For a discussion on the water hammer equations reference is made to the open channel flow equations (Chapter 5), given here as

bs

¶H ¶Q + ¶t ¶x

= 0

¶Q ¶ æ Q 2 ö ¶H f + ç ÷ + gA + Q|Q| = 0 ¶t ¶x è A ø ¶x 2DA

(7.1)

(7.2)

where Q H bs A f D

= pipe discharge (m3/s); = piezometric head above any horizontal reference level (m); = storage width due to fluid compressibility and pipe expansion (m); = cross-sectional area of the pipe (m2); = Darcy-Weisbach friction factor; = pipe diameter (m).

In the momentum equation the friction term has been replaced by a term based on the Darcy-Weisbach concept. In this form, the coefficient expresses the energy head loss DH as the fraction of the velocity head which is lost over a pipe length equivalent to its diameter D. For a length L of the pipe the head loss is then defined as

DH = f

L u2 D 2g

(7.3)

In water hammer, the role of the convective momentum is small and usually neglected. Introducing, furthermore, a characteristic celerity a as

a =

g

A

(7.4)

bs

the water hammer equations are best known in the form gA ¶H ¶Q + 2 a ¶t ¶x

= 0

¶Q ¶H f + gA + Q|Q| = 0 ¶t ¶x 2DA

(7.5)

(7.6)

The physical behaviour of water hammer is easily understood by transforming these equations into the form

WL | Delft Hydraulics

119

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

¶ æ gA ö ¶Q H÷ + a ç ¶t è a ¶x ø

UNESCO - IHE

= 0

(7.7)

¶Q ¶ æ gA ö +a ç H÷ = E ¶t ¶x è a ø

(7.8)

where E is to be seen as a sink term representing the effect of friction losses.

E = -

f 2DA

Q|Q|

(7.9)

Successive addition and subtraction of Equations (7.8) and (7.7) gives

¶ æ ¶ æ gA ö gA ö H÷ + a H÷ = E çQ + çQ + ¶t è ¶x è a a ø ø

(7.10)

¶ æ ¶ æ gA ö gA ö H÷ -a H÷ = E çQ çQ ¶t è ¶x è a a ø ø

(7.11)

Referring to the concept of total derivatives, as introduced in Chapter 3, Equations (7.10, 7.11) express the integration of the equation

d æ gA ö H÷ = E çQ + dt è a ø

(7.12)

along the characteristic line, or simply characteristic, defined by

dx dt

= a

(7.13)

and the integration of

d æ gA ö H÷ = E çQ dt è a ø

(7.14)

along the characteristic line

dx dt

= -a

(7.15)

Neglecting the influences of friction, Equation (7.12) represents the condition

Q+

WL | Delft Hydraulics

gA H = constant a

(7.16)

120

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

along the positive characteristic given by Equation (7.13), while Equation (7.14) leads to

Q-

gA H = constant a

(7.17)

along the negative characteristic, given by Equation (7.15). Before discussing the role of these equations in the algorithm leading to the solution of water hammer problems, let us consider again the meaning of the characteristic celerity a in relation to water storage. In analogy with Hooke's law expressing the elastic deformation of a string under the influence of the axial stress

s where s e E

= eE

(7.18)

= axial stress (N/m2); = strain per unit length of the string; = Young's modulus of elasticity.

The compression of water under the influence of changing pressures is seen as an elastic process described by

dp = - K

dV V

(7.19)

where dp K V dV

= change in fluid pressure (N/m2); = bulk modulus of water elasticity (N/m2); = water volume; = change of this water volume under the influence of pressure change dp.

For a control volume of length dx the compression of fluid in the pipe can be made equivalent to a virtual storage of an incompressible fluid in a virtual pipe slot by the relation

dV = -

V dp = - b s dH dx K

(7.20)

where the change in pressure level dH is related to the change in water pressure dp, by

dp = r g dH

(7.21)

assuming that the contribution of the change in fluid density in this relation can be neglected.

WL | Delft Hydraulics

121

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Substitution of (7.21) into (7.20) and replacing dV by Adx, gives

rgA K

bs =

(7.22)

For the derivation of the effect of pipe deformation it will be assumed that the pipe wall is elastic and that the pipe has expansion joints all along its axis. These joints will prevent the development of axial stresses under the influence of changing water pressure in the pipe. The general stress-strain relation for elastic pipe material under the influence of stresses in axial and tangential direction reads as

s 1 = ms 2 = e E where s1 s2 m

(7.23)

= hoop stress (tangential direction) (N/m2); = axial stress (N/m2); = Poisson's ratio.

For the case where the axial stress cannot develop the equation reads, for a change in tangential stress

ds 1 = e E

(7.24)

The expression of e in terms of the change in the circular pipe circumference with radius r

e

d(2p r) 2p r

=

=

ds 1 E

(7.25)

now leads to the relation

dr r

=

ds 1 E

(7.26)

or

dA A

=

2ds 1 E

(7.27)

Referring to Figure (7.2) the change in pipe material stress is related to the change in pressure dp as

2dF = Ddp

(7.28)

or

WL | Delft Hydraulics

122

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

2eds 1 = Ddp

UNESCO - IHE

(7.29)

D

P

e

t

t

Figure 7.2

Forces on semi-cylinder under the influence of water hammer induced pressures

Substitution of Equations (7.29), (7.21) into Equation (7.27) gives

dA = r gA

D dH eE

(7.30)

Relating this compression volume to the water storage in the fictive pipe slot gives

dV = dA dx =

rgA

D dH dx = b s dH dx e

or

bs = r g A

D eE

(7.31)

The combined effects of water compressibility and pipe expansion leads to the virtual pipe slot width

bs =

r g A æ DK ö ç1 + ÷ K è eE ø

(7.32)

Substitution of this relation into Equation (7.41) gives the expression for the pressure wave celerity

a =

K r DK 1+ eE

(7.33)

The effect of axial stresses and composite pipe walls requires a further generalisation of this equation to the form

WL | Delft Hydraulics

123

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

a =

K r K 1 +y E

UNESCO - IHE

(7.34)

with

y

=

D e

(7.35)

for the case of homogeneous elastic pipe material and the presence of expansion joints all along the pipe. Referring to Chaudhry (1987), the following relations can be derived for Y:

y

=

D (1 - m 2) e

(7.36)

for the case of thin-walled elastic conduits anchored against axial movement throughout their length;

y

=

D (1 - 0.5m ) e

(7.37)

for the case of thin-walled elastic conduits anchored against axial movement at the upper end;

y =1 ;

E=G

(7.38)

for the case of an unlined tunnel with a modulus of rock elasticity G, and

y

=

DE GD + Ee

(7.39)

for a steel-lined tunnel through rock. For values of the various material properties reference is made to literature and to Table 7.1 and Table 7.2 for some very common ones.

WL | Delft Hydraulics

124

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Table 7.1 Pipe material properties

Material

Modulus of elasticity E

Poisson’s ratio

(Gpa) Mild steel

200 - 212

0.27

0.8

0.46

Concrete

14 –30

0.10 –0.15

Cast iron

80 –170

0.25

Polyethylene

Table 7.2 Fluid properties

Fluid

Temperature

Density

Bulk modulus of elasticity K

(oC)

(kg/m3)

(GPa)

Fresh water

20

999

2.19

Sea water

15

1025

2.27

Oil

15

900

1.5

7.3

The Method of Characteristics for Water Hammer

The method of characteristics is by far the most commonly used equation solver in water hammer analysis. In the discussion of Chapter 3 it was concluded that the method of characteristics is a rather unpractical approach to solving free surface flow problems. Most of the reasons given there, however, do not hold for water hammer analysis. As shown by Equations (7.13, 7.15 and 7.33), the characteristic directions for the water hammer equations usually are straight lines, which allow for the construction of an equidistant network, as shown by the following example. Consider the hypothetical case of a pipe with frictionless flow supplied at the upstream end from a large reservoir and controlled at the downstream end by a valve, as shown in Figure (7.3). Data for this problem are given as follows: Hres Q A L a g

= 50 m = 1 m3/s = 1 m2 = 1000 m = 1000 m/s = 10 m/s2

The valve is closed over a period of 5 seconds, controlled to provide a linearly decreasing discharge over this time.

WL | Delft Hydraulics

125

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

The characteristic grid for this problem is shown in Figure 7.4. Initially all discharges and pressure levels are equal to 1 m3/s and 50 m respectively. Referring, again, to Equations (7.16 and 7.17), the solution at point A1 is found by solving the boundary condition

H A1 = 50 m and by Equation (7.17) applied along the negative characteristic from point B0 to point A1 Q A1 - gA/a H A1 = Q B0 - gA/a H B0 = 1 - 0.01*50 = 0.50 m3 /s

giving H A1 = 50 m; Q A1 = 1.0 m 3 /s

Note that this solution is still equivalent to the initial steady state condition along the pipeline between A and B. This is explained by the fact that the effect of the valve closure starting at point B0 cannot influence conditions at the upstream boundary before the negative characteristic through point B0 has arrived at point A1. The triangle formed by the points A0, B0 and A1, therefore, is a steady state region. The solution at point B1 is found in a similar way by combining the downstream boundary condition at t=1 s, given as

Q B1 = 0.8 m3/s with the condition along the positive characteristic passing from A0 to B1, given as Q B1 + gA/a H B1 = Q A0 + gA/a H A0 = 1 + 0.01*50 = 1.50 m3 /s

resulting in H B1 = 70 m; Q B1 = 0.8 m3 /s

Similarly, solutions at subsequent points can be computed, with results shown in Table 7.3. The maximum pressure found at point B is 90 m water column. The gradual valve closure gives a considerable operation improvement over the case of a sudden valve closure. Closing the valve instantaneously at time t=0 gives a pressure H=150 m at point B1, as can be verified easily by applying Equation (7.16) with QB1=0 m3/s. Further improvements are found by slowing the closure operation further down. This example, therefore, demonstrates clearly the origin of water hammer problems and the essential approach to water hammer prevention.

WL | Delft Hydraulics

126

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

reservoir

H supply pipe

L Situation sketch of valve closure problem

t (s)

t (s)

Figure 7.3

10

10

9

9

8

8

7

7

6

6

5

5

4

4

3

3

2

2

1

1

0

0

0

100

0.2 0.4 0.6 0.8 1.0

0

1.0 Q (m 3 /s)

x (m) A Figure 7.4

Characteristic lines for the water hammer problem

Table 7.3

Pressures and discharges computed for the problem of valve closure

Point A B

WL | Delft Hydraulics

B

0

1

2

3

4

5

6

7

8

9

10

variable

1.0 50 1.0 50

1.0 50 0.8 70

0.6 50 0.6 90

0.2 50 0.4 70

0.2 50 0.2 50

0.2 50 0.0 70

-0.2 50 0.0 70

-0.2 50 0.0 30

0.2 50 0.0 30

0.2 50 0.0 70

-0.2 50 0.0 70

Q (m3/s) H (m) Q (m3/s) H (m)

127

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

7.4

UNESCO - IHE

Exercise

Water flows from a reservoir through a horizontal pipe with a length L=10 km. The reservoir water level HR is 50 m +MSL. The pipe has a diameter D of 1.00 meter. At the end of the pipe the outflow is controlled by a valve. The steady state water velocity u=2 m/s. The Darcy-Weisbach friction coefficient f is 0.01. Use g=10 m/s2 in your computations. Questions 1. 2. 3. 4. 5. 6.

7.

8.

Compute the pressure distribution along the pipe line for the steady state flow, assuming that the energy loss at the valve can be neglected. For the computation of Questions 3-7 we will assume that the complete energy loss along the pipe is now concentrated at the valve. Recompute the pressure distribution. Compute the characteristic pressure wave celerity for the data given below. Compute, with the rigid water column theory, the pressure increase at the pipe end if the valve is closed over a period of 1 second. Is this pressure increase realistic? Compute the pressure increase with the water hammer theory for fast closure. Compute the pressure as a function of time at the location of the valve by using the method of characteristics. Use a valve closure function which leads to a linear decrease in out flowing discharge over a period of 1 minute. Neglect friction. While keeping the other parameters constant, what will be the effect on the maximum pressure increase if: a) the closure time is increased; b) the steel pipe is replaced by a PVC pipe; c) the pipe wall thickness is increased; d) the pipe diameter is increased; e) the pipe length is decreased; f) air is entrained into the flow. Repeat the computations for the case where pipe wall friction is included and the energy loss at the valve is neglected. Produce a graph of the results.

Additional data: e = 0.01 m; K = 2 GPa; E = 200 GPa. The pipe has frequent expansion joints.

WL | Delft Hydraulics

128

Computational Hydraulics Version 2006-1

8

© 2005 Adri Verwey

UNESCO - IHE

References

Abbott, M.B. (1979). Computational Hydraulics – Elements of the Theory of Free Surface Flows, Pitman, London/San Francisco/Melbourne. Arakawa, A. (1966). Computational design for long-term numerical integration of the equations of fluid motion: two-dimensional incompressible flow, Part I. Journal of Computational Physics. 1(1), 119-143.

Baptist. M.J. (2005). Modelling floodplain biogeomorphology. Ph.D. thesis, ISBN 90-407-2582-9, 193 pp., Delft University of Technology, Faculty of Civil Engineering and Geosciences, Section Hydraulic Engineering. Chaudhry, M.H. (1987), Applied Hydraulic Transients, Van Nostrand Reinhold Company Inc., New York. Chaudhry, M.H. (1993), Open-Channel Flow, Prentice Hall, Anglewood Cliffs, USA. Chow, V.T. (1959). Open Channel Hydraulics, McGraw-Hill, New York. Chanson, H. (1999). The Hydraulics of Open Channel Flow, Arnold Publishers/Wiley, Paris/New York. Cunge, J.A., (1969). On the subject of a flood propagation computation method (Muskingum method). Journal Hydrology Research 7, No.2. Cunge, J.A., Holly, F.M. and Verwey, A. (1986). Practical Aspects of Computational River Hydraulics (420 pages), Pitman Publishing Ltd., London, Great Britain, 1980. Reprinted at Iowa Institute of Hydraulic Research, USA, 1986. (Also translated into Russian). De Saint Venant, B. (1871). Théorie du Mouvement Non-Permanent des Eaux avec Application aux Crues des Rivières et à l’Introduction des Marées dans leur Lit, Comptes Rendus de l’Académie des Sciences, 73, pp. 148-154, 237-240. Henderson, F.M. (1966). Open Channel Flow, McMillan Company, New York. Hesselink, A.W., Stelling, G.S., Kwadijk J.C.J. and Middelkoop, H. (2003). Inundation of a Dutch river polder, sensitivity analysis of a physically based inundation model using historic data. Water Resour. Res., 39(9), 1234. Hirsch, C. (1990). Numerical Computation of Internal and External Flows, Wiley, New York. Newton, I. (1687). Mathematical Principles of Natural Philosophy (trans. A. Motte, ed. F. Cajuri), Berkeley University Press, California (1947). Rodriguez Uthurburu, R, (2004). Evaluation of physically based and evolutionary data mining approaches for modelling resistance due to vegetation in SOBEK 1D-2D, M.Sc. thesis HH 485, UNESCO-IHE, Delft, The Netherlands. Serban, P., Crookshank, N.L. and Willis, D.H. (2005). Intercomparison of Forecast Models for Streamflow Routing in Large Rivers, WMO, Geneva.

WL | Delft Hydraulics

129

Computational Hydraulics Version 2006-1

© 2005 Adri Verwey

UNESCO - IHE

Stelling, G.S., Kernkamp, H.W.J. and Laguzzi, M.M. (1998). Delft Flooding System: a powerful tool for inundation assessment based upon a positive flow simulation, Hydroinformatics ‘98, Babovi & Larsen (Eds.), Balkema, Rotterdam, 449-456. Stelling, G.S. and Duinmeijer, S.P.A. (2003). A staggered conservative scheme for every Froude number in rapidly varied shallow water flows. Int. J. for Numer. Meth. Fluids, 43, 1329-1354. Stoker, J.J. (1957). Water Waves, Pure and Applied Mathematics, Vol. IV, Interscience Publishers, New York. Toro, E.F. (1999). Riemann Solvers and Numerical methods for Fluid Dynamics, Springer, Berlin.

Uittenbogaard, R. (2003). Modelling turbulence in vegetated aquatic flows, presented at the International workshop on RIParian FORest vegetated channels: hydraulic, morphological and ecological aspects, 20-22 February 2003, Trento, Italy. Verwey, A. (1994). Linkage of Physical and Numerical Aspects of Models Applied in Environmental Studies, keynote lecture in: Proceedings of the Conf. on Hydraulics in Civil Engineering, Brisbane, Australia. Verwey, A. (2001). Latest Developments in Floodplain Modelling - 1D/2D Integration, keynote lecture in : Proceedings of the 6th Conf. on Hydraulics in Civil Engineering, Hobart, Australia. Werner, M.G.F., van Dijk, M. and Schellekens, J., (2004). DELFT-FEWS: An open shell flood forecasting system, in Proceedings of the 6th International Conference on Hydroinformatics, Liong, Phoon and Babovi (Eds.), World Scientific Publishing Company, Singapore, 1205-1212.

WL | Delft Hydraulics

130