- Email: [email protected]
Dam-break flow simulation: some results for one-dimensional models of real cases
Journal of Hydrology 216 (1999) 227–247
Dam-break flow simulation: some results for one-dimensional models of real cases P. Garcia-Navarro*, A. Fras, I. Villanueva Department of Fluid Mechanics C.P.S., University of Zaragoza, 50015 Zaragoza, Spain Received 20 July 1998; accepted 12 January 1999
Abstract In many countries, the determination of the parameters of the wave, likely to be produced after the failure of a dam, is required by law, and systematic studies are mandatory. There is a necessity to develop adequate numerical solvers which are able to reproduce situations originated from the irregularities of a non-prismatic bed and to model the complete equations that progress despite the irregular character of the data. Many hydraulic situations can be described by means of a one-dimensional (1D) model, either because a more detailed resolution is unnecessary or because the flow is markedly 1D. Many techniques have been developed recently for systems of conservation laws in 1D (in the context of gas dynamics). Some years after their adoption for solving problems in gas dynamics, upwind and total variation diminishing (TVD) numerical schemes have been successfully used for the solution of the shallow water equations, with similar advantages. Their use is nevertheless only gradually gaining acceptance in this sector. The performance of some of these techniques for practical applications in river flow is reported in this work. 䉷 1999 Elsevier Science B.V. All rights reserved. Keywords: Flood routing; Explicit TVD methods; Unsteady shallow water flows; Dam break flow
1. Introduction Many problems of river management and civil protection consist of evaluation of the maximum water levels and discharges that may be attained at particular locations during the development of an exceptional meteorological event. There is another category of events of catastrophic nature whose effects also fall into the civil protection area. It is the prevision of the scenario subsequent to the almost instantaneous release of a great volume of liquid. The situation is that of the breaking of a man made dam. In many countries, the determination of the parameters of the wave, likely to be produced after the * Corresponding author. Fax: ⫹34-976761882. E-mail address: [email protected] (P. Garcia-Navarro)
failure of a dam, is required by law (Molinaro,1991; Betamio de Almeida and Bento Franco, 1993), and systematic studies are mandatory. There are works based on scaled physical models of natural valleys, but they represent expensive efforts that are not devoid of difficulties. Hence, there is a necessity to develop adequate numerical models which are able to reproduce the situations that originate from the irregularities of a non-prismatic bed. It is also necessary to trace their applicability considering the difficulty of developing a model capable of producing solutions of the complete equations despite the irregular character of the river bed. For many practical applications, it is accepted that the unsteady flow of water in a one-dimensional (1D) approach is governed by the shallow water or St. Venant equations. These represent the conservation
0022-1694/99/$ - see front matter 䉷 1999 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(99)00007-4
228
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 1. Dam break problem. Lax–Friedrichs versus exact solution.
of mass and momentum along the direction of the main flow. The prediction of water profiles and discharges for unsteady as well as steady situations of hydraulic systems is often linked to the use of numerical methods able to solve the system of equations. As in any other case of integration of a system of non-linear partial differential equations, the first difficulty is the choice of the numerical scheme. There is, at present, a wide range of techniques suitable for subcritical shallow water flows. However, it has traditionally been difficult to have a single method able to reproduce automatically any general situation.
The numerical modelling of unsteady flow in rivers is a complicated task and the difficulties grow as the pretensions to obtain better quality or more general solutions do (Cunge et al., 1980) There exist good methods, developed to deal with Gas Dynamics problems (Euler equations), able to cope with complex systems of discontinuities and shock waves (Hirsch, 1990; MacCormack, 1971; Roe, 1989). The shallow water, or St. Venant equations, being a hyperbolic partial differential system, represent a good candidate for the application of many of these techniques. Many new schemes have been
Fig. 2. Dam break problem. MacCormack versus exact solution.
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
229
Fig. 3. Dam break problem. TVD MacCormack versus exact solution.
reported successful for the flow in channels (Glaister, 1988; Alcrudo and Garcia-Navarro, 1992; GarciaNavarro and Alcrudo, 1992). However, their application to the river flow is not so common in the literature (Jin and Fread, 1997) and their adaptation to river hydraulics is not straightforward because of the irregular topography. The presence of extreme slopes, high roughness and strong changes in the irregular geometry represent a great difficulty that can lead to important numerical errors presumably arising from the source terms of the equations. In this work, the performance of three finite difference techniques is reported and compared so that their suitability for this kind of problems can be discussed. First, the system of equations modelling 1D unsteady river flow will be presented in Section 2. The finite difference numerical techniques applied to solve the equations are the Lax–Friedrichs scheme and the explicit MacCormack scheme with and without TVD corrections. TVD corrections produce a variation of the classical predictor–corrector scheme based on the theory of total variation diminishing methods and introduce upwind effects in the solution while maintaining central difference accuracy; They are described in Section 3. Section 4 starts with a validation of the numerical results provided by the three schemes against an exact solution for a problem of an ideal dam break. Then, the numerical treatment applied to the practical simulation of river flow will be discussed. Finally, in Section 5, the results
obtained for the simulation using data from real rivers are presented and the performance of the numerical schemes is compared. 2. Mathematical model Many hydraulic situations can be described by means of a 1D model, either because a more detailed resolution is unnecessary or because the flow is markedly 1D. The fundamental hypothesis implied in the numerical modelling of river flows are formalised in the equations of unsteady open channel flow. Under the St. Venant hypotheses, 1D unsteady flow can be written in the form: 2U 2F ⫹ R; 2t 2x
1
with U A; QT ; Q2 ⫹ gI1 F Q; A
!T ;
R 0; gI2 ⫹ gA S0 ⫺ Sf T which emphasises the conservative character of the system in the absence of source terms. The effects of the wind as well as those of the Coriolis force have been neglected and no lateral
230
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 4. Cross-section profiles at different locations. Aragon river.
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
231
Fig. 5. Water profiles at different times, Aragon river.
inflow/outflow is considered. In Eq. (1), A is the wetted cross-sectional area, Q, the discharge and g, the acceleration due to gravity. I1 represents a hydrostatic pressure force term as described in Cunge et al. (1980). The pressure force can have a component in the direction of the main stream due to the reaction of the walls in case of variations in shape along this direction. The amount of this force depends on the cross-sectional variation for constant depth. It is important to note that the validity of this approach is linked to the hypothesis of gradual variation. If sudden expansions or contractions take place, the approach is not valid. I2 accounts for the pressure force in a volume of constant depth h due to longitudinal width variations. The bed slope is S0. Sf stands for the energy grade line and is defined, for example, in terms of the Manning’s roughness coefficient n: Sf
Q兩Q兩n2 ; A2 R4=3
with R A=P; P being the wetted perimeter.
The Jacobian matrix of the system (Eq. (1)) is ! 0 1 2F A ; 2U c2 ⫺ u2 2u p where u Q=A and c gA=b. The eigenvalues and eigenvectors of A are: a1;2 u ^ c e1;2 1; u ^ cT : Other formulations are possible and frequently encountered. 3. Numerical model There is a long list of finite difference techniques suitable for the numerical solution of the equations presented (Cunge et al., 1980; Glaister, 1998; Garcia-Navarro and Saviron, 1992; Alcrudo and Garcia-Navarro, 1992). If there is an interest to deal with unsteady problems that may contain rapidly varied flows with
232
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 6. Water depth history at two points, Aragon river.
Fig. 7. Discharge history at two points, Aragon river.
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
233
Fig. 8. Mass balance error, Aragon river.
the possible presence of transcritical flow or formation of surges, attention must be focused on a class of these techniques: those able to progress automatically towards a weak or discontinuous solution of the equations. Among these, there is still a variety of methods. The explicit finite difference schemes used for the numerical computations in this work are presented next.
3.1. Lax–Friedrich’s scheme This is a first order in space and time finite difference technique. For a general homogeneous conservation system
updating: Un⫹1 uUni ⫹ i ⫺
1⫺u n Ui⫹1 ⫹ Uni⫺1 2
Dt Fn ⫺ Fni⫺1 : 2Dx i⫹1
0 ⱕ u ⬍ 1: The value u 1 renders the scheme unstable. More stability and more numerical diffusion is introduced as u approaches zero. A value u 0.1 is usually adopted. It is easy to implement and very robust. The interest of considering its performance was to establish the extent of the numerical error introduced for problems highly dominated by the source terms. 3.2. MacCormack’s scheme
2U 2F ⫹ 0; 2t 2x the procedure to update one time step Dt the interior points 2,…, N-1 of a regular grid is based on a nodal
A classical finite difference scheme suitable for the discretisation of (Eq. (1)) is the explicit two step predictor–corrector MacCormack method (MacCormack, 1971). It is a technique of proved efficiency for unsteady open channel flow modelling
234
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 9. Bed level profile, Neila river.
(Garcia-Navarro and Saviron, 1992; Garcia Navarro and Zorraquino, 1993; etc). For linear problems it reduces to the Lax–Wendroff scheme so that, the Courant–Friedrich–Lewy criterion Dt ⱕ
Dx max{ai }
ensures numerical stability. Being second order accurate in space and time, it offers good resolution and has a great conceptual simplicity. It is composed of a sequence of two sub-steps in which the spatial derivatives are taken in alternate directions as follows: Predictor step ! Upi Uni ⫺
Dt n F ⫺ Fni Dx i⫹1
Corrector step ! Uci Upi ⫺
Dt p F ⫺ Fpi⫺1 Dx i
so that the updated solution is given by Uin⫹1
1 p U ⫹ Uci : 2 i
It is well known, however, that classical second order schemes show oscillatory behaviour near discontinuities, and McCormack’s scheme is no exception. This is why the introduction of artificial viscosity became usual in reported numerical simulation practice (Garcia and Kahawita, 1986). The difficulty with the artificial viscosity approach is that it is hard to determine an appropriate form that introduces just enough dissipation to preserve monotonicity without causing unnecessary smearing. 3.3. Total variation diminishing MacCormack In the 1D theory of TVD methods, it is shown that some classical schemes can be reformulated in a way that improves their performance without great effort. In particular it is possible to add a conservative and non-linear modification to the Lax–Wendroff scheme to render it non-oscillatory at points of extrema without losing its accuracy in other regions of the flow. The application of this idea to the MacCormack version of the Lax–Wendroff scheme for the 1D
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
235
Fig. 10. Initial top width, Neila river.
shallow water equations is reported in (GarciaNavarro and Alcrudo, 1992). The TVD correction can be cast as a modification of the last step of MacCormack scheme with inclusion of new terms to eliminate the oscillations. Un⫹1 i
1 p Dt n U ⫹ Uci ⫹ D ⫺ Dni⫺1=2 : 2 i Dx i⫹1=2
The form of the new terms is based on Roe’s linearisation of the system (Eq. (1)). According to Roe (1989) a matrix A can approximate locally the Jacobian A of the system. The eigenvalues a k and eigenvectors ek of A satisfy X k k dUi⫹1=2 Ui⫹1 ⫺ Ui a e i⫹1=2 k
dFi⫹1=2 Fi⫹1 ⫺ Fi
X
a k ak ek i⫹1=2
1D shallow water equations (k 2) in Alcrudo and Garcia-Navarro (1992). The particular expressions used in the present work are ~ a 1;2 u~ ^ c; ~ T; e1;2 1; u~ ^ c
a1
a2
~ dA ⫺ dQ c~ ⫹ u : 2c~
The average velocity and celerity at the interface i ⫹ 1/2 are p p Qi⫹1 Ai ⫹ Qi Ai⫹1 p p ; u~ p Ai Ai⫹1 Ai⫹1 ⫹ Ai
k
inside every computational cell (i,i ⫹ 1). Expressions for ak ; ak and ek can be found, for the 1D Euler equation (k 3) in Roe (1989) and for the
~ d A ⫹ dQ c~ ⫺ u 2c~
c~
r g A=bi ⫹ A=bi⫹1 2
and dA Ai⫹1 ⫺ Ai , dQ Ai⫹1 ⫺ Ai .
236
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 11. Initial Froude number distribution, Neila river.
The third step is then made using: 1X k Dt k 兩a i⫹1=2 兩 ai⫹1=2 兩a ki⫹1=2 兩 1 ⫺ Di⫹1=2 2 k Dx 1 ⫺ Fki⫹1=2 eki⫹1=2 : The function F is a flux limiter responsible for the smoothing near fronts. Mathematical expressions for this function can be found in Sweby (1984). 4. Applications 4.1. Dam break test case The idealised dam break problem was chosen because it is a classical example of non-linear flow with shocks to test conservation in numerical schemes and, at the same time, has an analytical solution (Stoker (1957). This problem is generated by the homogeneous 1D shallow water equations for the ideal case of a flat and
frictionless channel of unit width and rectangular cross-section, with the initial conditions. 8 L > > if x ⱕ < hL 2 h x; 0 : > L > : hR if x ⬎ 2 Q x; 0 0: If the calculation times used are so as to avoid interaction with the extremities of the channel, the boundary conditions are trivial. Figs. 1–3 show a sample of the solution for a value of the initial height ratio hL:hR 20. They represent the profiles of the water surface 30 s after the dam break as computed on a grid of 200 cells. The analytic solution appears as a continuous line. The dotted lines represent the numerical solutions obtained with Lax– Friedrich’s scheme in Fig. 1, MacCormack’s scheme in Fig. 2 and TVD MacCormack’s scheme in Fig. 3. The excessive numerical diffusion introduced by the first order scheme and the lack of robustness of the
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
237
Fig. 12. Discharge history at Neila village, Neila river.
second order scheme are shown. The TVD step provides the necessary numerical viscosity to cope well both with the shock front and with the critical point at the former dam location. The size of the time step is chosen from the condition CFL 0.9. 4.2. Application to river flow The numerical scheme applied must also be able to reproduce the surface discontinuities that may appear as a result of the irregular geometry itself or because of the hydraulic structures existing along the river. Independently of the particular problem, a preprocessing and storage of the geometrical and hydraulic data is always required. The resolution is compatible with their usage in the form of numerical functions (data tables). They furnish the information required by the equations through a process of interpolation. 4.3. Geometrical data The 1D model requires the selection of a series of
river cross-sections each one representing the geometry of a reach of length Dx. The approximation of every cross-sectional shape as a juxtaposition of trapezoids simplifies the systematic calculation of the values of cross-section area, surface width or wetted perimeter as increasing functions of the water depth for a sequence of discrete values. At the same time, the consistency of the data is ensured as independent calculations of these functions may lead to numerical instabilities (Bento Franco and Betamio, 1991). The simulation starts from a definition of the topographic features in the form of numerical functions (data tables) of the depth and the distance along the river. These are matrices of dimensions (number of cross-sections x number of water depth levels). The data tables do not correspond in general to equidistributed points along the river and never to the computational grid positions. One option is the numerical generation of intermediate sections by the interpolation of the surveyed cross-sections. However, the particular value of any of the functions at a nodal
238
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 13. Water depth history at Neila village, Neila river.
position for an arbitrary water depth can be obtained each time via interpolation from the data tables. More or less sophisticated interpolations can be performed. Considering that any interpolation may introduce numerical errors in the results and the irregular character of the river, a liner interpolation between data tables was used in this work. The bed slopes were determined using the values of the bottom level at each cross-section. Another difficulty arises from the presence of functions of these variables such as the integrals I1 and I2 in the system of equations. Their evaluation in a case of simple geometry is direct and analytical. The irregular character of the river cross-sections requires some kind of approximation. Considering a cross-section well approximated by a polygonal line of j vertices, the integral of the hydrostatic pressure at a section is estimated as I1 h
Zh 0
h ⫺ hb hdh
jX ⫺1
Zhk ⫹ 1
k0
hk
h ⫺ hb hdh;
where inside every interval (k, k ⫹ 1) b h b hk ⫹
h b hk⫹1 ⫺ b hk ; hk⫹1 ⫺ hk
hk ⱕ h ⱕ hk⫹1 : The pressure force from the walls in an interval Dx is estimated from I2
1 I h; x ⫹ Dx ⫺ I1 h; x: Dx 1
The numerical schemes presented are just defined for the homogeneous part of the system of equations. It still remains the question of the best discretisation of the source terms. The simplest pointwise option for a generic source term function R A; Q can be expressed as ÿ Ri ⬇ R Ani ; Qni ; which consists of a simple evaluation of the source term functions at the nodal point using the already known dependent variables. This is the procedure
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
239
Fig. 14. Discharge history at Villarvelayo village, Neila river.
applied to the difference schemes tried. Other possibilities include some form of average between neighbouring nodal points or an upwind treatment (Glaister, 1998). Problems may arise when the source terms contain spatial derivatives and a definition of their discrete representation must be established. This difficulty is met for the discretisation of the bottom slope (S0) and variable width (I2) source terms. When dealing with the Lax–Friedrichs scheme, a standard central difference was applied. The option selected in the case of the MacCormack scheme was 2b b ⫺ bi ⬇ i⫹1 Dx 2x 2b bi ⫺ bi⫺1 ⬇ Corrector step ! Dx 2x Predictor step !
and similarly for 2z=2x: The limitation imposed by the friction term can be relaxed if this term is not discretised using an explicit approach. A pointwise semi-implicit discretisation of
the friction term (Sf) proved efficient in reducing instabilities, Huang and Song (1985). Q兩Q兩 AR2=3 ; K 2 n K n 兩Q兩 uQni ⫹ 1 ⫺ uQn⫹1 ; Sfi ⬇ i K2 i Sf
0 ⱕ u ⱕ 1:
This strategy can be very easily implemented in an explicit procedure. The term containing Q n⫹1 joins the left hand side so that the updating expression for the discharge is modified but the method remains explicit. The weighting parameter u controls the degree of implicitness. It must also be signaled that in this 1D model a single value of the roughness coefficient was used. The mathematical model also considers a single average water velocity at every cross-section and no distinction was made between river bed and flood plain. A uniform distribution of nodal points was considered. All these limitations in the model require
240
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 15. Water depth history at Villarvelayo village, Neila river.
a careful analysis of the results and invite further improvements.
5. Results 5.1. Case 1 Taking the topographical data of an upper reach of about 14 440 m in the river Aragon in the Spanish Pyrenees, an application of the numerical simulation is performed for the estimation of the propagation of a dambreak wave. A 15 m high dam is placed at the location x 8 000 m and the hypothesis of total and sudden failure is made. The interest is in the performance of the numerical techniques. The roughness coefficient is set to n 0.05. The initial discharge is 50 m 3/s and a water profile upstream and downstream of the dam provides the initial conditions for depth. The reach is discretised in 500 intervals (Dx 28 m). Fig. 4 displays some of the measured transverse
cross-sections along the river reach. Fig. 5 represents a series of computed longitudinal river profiles for subsequent instants of time starting with the initial condition. The time step size was the maximum allowed by numerical stability. The bottom level profile is also displayed in the figure. The profiles correspond to the results obtained with the MacCormack scheme. They do not differ very much from those supplied by the other schemes in this case. Figs. 6 and 7 show depth and discharge histories, respectively, at two locations along the reach (x 8500 and 1200 m). The predictions provided by MacCormack are compared to those from Lax–Friedrich’s scheme. Some differences between them are noticeable, in particular the smearing introduced by the latter scheme. The consequence of this smearing is quantified in the mass balance error. As we are dealing with a fluid of constant density, mass balance is equivalent to volume balance. It is calculated as the difference between the volume variation along the reach and the net volume inflow through the boundaries. In order to deal with a dimensionless quantity,
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
241
Fig. 16. Froude number distribution, Neila river.
the difference is calculated relative to a meaningful volume of water, the initial volume in this case. The estimation of the error in mass conservation introduced by the numerical calculations is shown in Fig. 8 using the two classical schemes. No noticeable improvement was introduced by the TVD correction in this case. It is worth noting here that the considerable difference in the error produced by the two numerical schemes as shown in Fig. 8 may seem excessive if Fig. 7, for instance is taken as a reference. The reason that justifies this discrepancy, is the different nature of the two plots. Fig. 7 shows the discharge history at a point whilst Fig. 8 represents the global volume error evolution in time. Spurious discharge oscillations at a single point along the channel make this global error accumulate though the solution is well behaved in the rest of the domain. 5.2. Case 2 Hydraulic phenomena of marked unsteady character can be included under the form of a more or less
steep inflow hydrograph Q Q t. Predicting the evolution of the flooding wave that it represents is the aim of the method. From a mathematical point of view, the transit from a value of the discharge to a higher one can happen in different ways. The parameter T (characteristic time of breach formation), with T 0 in the case of a sudden release of water, can be useful and may help to analyse aspects and consequences of the influence of this variable. Our case 2 deals with such a kind of event at the upstream inflow point of Neila river (Burgos, Spain). A reach of 18000 m of this river measured from its upstream dam was studied to analyse the flood risk at Neila village, located 5350 m downstream from the dam. Several dam break and friction roughness hypotheses were made and the corresponding discharge hydrographs were calculated. From the set of results obtained only a few were selected for inclusion in the present work; those corresponding to the more extreme breaking process with T 0.94 h and a peak discharge of 160 m 3/s.
242
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 17. Froude number distribution, Neila river.
Manning roughness was assumed to be uniform and equal to 0.05. The topographic characteristics of this river are represented in Fig. 9 (bottom level variation) and Fig. 10 (plan view of the breadth variations along the axis of the river). These, together with the small value of its base discharge (Q 1 m 3/s) and the corresponding initial values of the Froude number (Fig. 11) rendered the numerical computation more challenging than the former case and accentuated the differences among the capabilities of the numerical schemes. Fig. 12 displays the discharge hydrographs at Neila village location as calculated from the three schemes. Fig. 13 shows the equivalent plot of the depth histories. In the same order, Figs. 14 and 15 show the discharge and water depth changes in time at Villa velayo, a village located 16 875 m downstream from the dam site. Figs. 16 and 17 show the distribution of the Froude number along the river at four different times and clearly indicate the transcritical character of the flow. It is interesting to point out here that the upstream
boundary condition consisted of the inflow hydrograph Q Q t. This value was complemented with the value of the water depth at the upstream point obtained from the characteristics theory whenever there was subcritical flow at the point. The value of the Froude number at the second nodal point was continuously controlled for this purpose. In case of supercritical conditions, the equation S0 Sf supplied the required water depth, that is, normal depth was imposed. The TVD correction of MacCormack scheme was crucial in stabilizing the calculation in this example and was the only way to return to the correct base flow (initial steady state) after the passage of the flood wave. The simple Lax–Friedrichs scheme proved very inefficient in dealing with the dominating source terms present in this river reach. In Fig. 18 the improvement provided by the TVD scheme is clear. The discussion concerning the important discrepancy between the magnitude of this global error and that of the discharge hydrographs at the particular locations is again as in the previous test case.
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
243
Fig. 18. Mass balance errors, Neila river.
5.3. Case 3 The third example is based on one of the cases presented by Jin and Fread (1997). These authors use a numerical formulation based on flux splitting and pointwise evaluation of the source terms. In an attempt to validate our models against those present in the literature, one of their dam break problems was calculated using the TVD McCormack scheme. This is a 36.6 m high dam located at 11.3 km in a 32.2 km reach of nearly uniform trapezoidal crosssection shape. A constant Manning friction factor n 0.04 is used. The dam fails according to a breach law that depends on a failure time Tf. Tf 0.02 h in this example as this case was the least favourable leading to the formation of supercritical regions near the dam. The breaching dam is included in the calculation as an internal boundary separating two reaches upstream and downstream from it. One mass continuity equation and a dynamic weir type equation are used together with two numerical boundary conditions (characteristics).
Fig. 19 displays several longitudinal water profiles and the river bed profile. They correspond to the times indicated on the figure. The mass error balance, using as reference volume the initial volume of water, gives a value bounded by 5%. Fig. 20 shows the discharge hydrographs at three locations, one near the dam (x 11.3 km) and two downstream from it (x 16.1 km, x 24.1 km). The results from Jin and Fread are plotted in lines and compared with the numerical results from our model. Figs. 21 and 22 show the water depth and the Froude number distributions in space downstream from the dam at 4 different times after the beginning of the breach formation (t 0.05, 0.10, 0.15, and 0.2 h). There are some differences arising from a different water depth value at the right hand side of the dam. They might be due to a difference in the numerical boundary conditions at that section. Whilst the reference authors use a discretisation of the mass conservation equation, we apply a semi-Lagrangian treatment, that is, the trajectory of the backwards characteristic is used to bring information from the old to the new time level. This test case
244
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 19. Longitudinal profiles at different times.
is presented by Jin and Fread as an example in which, because of the short time of breach formation, regions of supercritical flow appear near the dam location. Unfortunately, no experimental data are supplied.
6. Conclusions A first order (Lax–Friedrich’s), a second order (MacCormack) and a high resolution (TVD MacCormack) explicit finite difference schemes are applied to the simulation of 1D unsteady flow in rivers. During the unsteady processes, some of the selected examples contain zones of supercritical flow that the schemes are able to cope with automatically. The presence of extreme slopes and strong changes in the irregular geometry represent a great difficulty that can lead to important numerical errors presumably arising from the source terms of the equations. The TVD version of the MacCormack scheme appears as superior in this kind of problems. The conceptual and programming simplicity of MacCormack scheme
is not greatly altered by the inclusion of the TVD correction. For smoother situations, the three schemes performed well. There is a difference in mass conservation ability between Lax–Friedrichs and MacCormack for a given grid size, with a clear superiority of the latter. TVD MacCormack does not present a remarkable superiority in these cases. In what concerns the particular problem of the dam break, given the initial conditions of the system and a temporal law for the breach formation (usual treatment of the problem in the literature) the simulation furnishes simultaneously the characteristics of the propagation of the wave downstream and the evolution of the water profile upstream of the dam site. In contrast to other numerical techniques incorporated in commercial codes (DAMBRK Fread, 1985), the procedure applied in the present work deals very efficiently with the presence of locally supercritical conditions anywhere along the river. The time step size is rather restricted in general to ensure stability. In this sense, irregular grids and
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
245
Fig. 20. Discharge hydrographs at x 11.3 km (dam site), x 16.1 km and x 24.1 km. Jin and Fread: Lines. Present model: Symbols.
Fig. 21. Depth profiles at t 0.05, 0.10, 0.15 and 0.020 h after the beginning of the breach formation. Jin an Fread: Lines. Present model: Symbols.
246
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 22. Froude profiles at t 0.05, 0.10, 0.15 and 0.20 h after the beginning of the breach formation. Jin and Fread: Lines. Present model: Symbols.
implicit versions of the schemes are envisaged as future work. Acknowledgements We would like to thank SERS S.A. for the partial funding of the second author and Mr. V. Zorraquino from SERS S.A. for many useful and interesting discussions during the course of this work. Thanks also to Dr. Jin from NOAA for supplying the necessary data from his own work. References Alcrudo, F., Garcia-Navarro, P., 1992. Flux difference splitting for 1D open channel flow equations. Int. J. for Numerical Methods in Fluids 14, 00. Bento Franco, A., Betamio, A., 1991. Simulaco uni e bidimensional de cheias provocadas por roturas de barragens em planicies de inundaco, V Simposio Luso-Brasileiro de Hidraulica e Recursos Hidricos SILUSB.
Betamio de Almeida, A., Bento Franco, A., 1993. Modelling of dam break flows. In: Chaudhry, M.H., Mays, L.W. (Eds.). Computer Modelling of Free Surface and Pressurized Flows, NATO ASI Series, 343–375. Cunge, J.A., Holly, F.M., Verwey, A., 1980. Practical Aspects of Computational River Hydraulics. Pitman, London. Fread, D.L., 1985. Channel routing. In: Anderson, M.G., Burt, T.P. (Eds.). Hydrological Forecasting, Wiley, New York. Garcia, R., Kahawita, R.A., 1986. Numerical solution of the St.Venant equations with the MacCormack finite difference scheme. Int. J. for Numerical Methods in Fluids 6, 259– 274. Garcia-Navarro, P., Alcrudo, F., Saviron, J.M., 1992. 1D open channel flow simulation using TVD McCormack scheme. J. of Hydraulic Engineering, ASME 118 (10). Garcia-Navarro, P., Saviron, J.M., 1992. McCormack’s method for the numerical simulation of one-dimensional discontinuous unsteady open channel flow. J. of Hydraulic Research 30, 95– 105. Garcia Navarro, P, Zorraquino, V., 1993. Numerical modelling of flood propagation through a system of reservoirs. J. of the Hydraulic Division, Am. Soc. of Civil Engineering 119 (3), 380–389. Glaister, P., 1998. Approximate Riemann solutions of the shallow water equations. J. of Hydraulic Research 26 (3), 293–306.
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247 Hirsch, C., 1990. Numerical Computation of Internal and External Flows, 2, Wiley, New York. Huang, J., Song, C.C.S., 1985. Stability of dynamic flood routing schemes. J. of Hydraulic Engineering, ASCE 111 (12), 1497– 1505. Jin, M., Fread, D.L., 1997. Dynamic flood routing with explicit and implicit numerical solution schemes. J. of Hydraulic Engineering, ASCE 123 (3). MacCormack, R.W., 1971. Numerical solution of the interaction of a shock wave with a laminar boundary layer. In: Holt, M. (Ed.). Proceedings 2nd International Conference on Numerical Methods in Fluid Dynamics, Springer, Berlin, pp. 151.
247
Molinaro, P., 1991. Dam break analysis: a state of the art. In: Bensar, D., Brebbia, C.A., Ovazar, D. (Eds.). Computational Water Resources, CMP. Roe, P.L., 1989. A survey on upwind differencing techniques, Lecture Series in CFD, Von Karman Institute for Fluid Dynamics, Belgium. Stoker, J.J., 1957. Water Waves. Interscience, New York. Sweby, P.K., 1984. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J. of Numerical Analysis 21, 995–1011.
Dam-break flow simulation: some results for one-dimensional models of real cases P. Garcia-Navarro*, A. Fras, I. Villanueva Department of Fluid Mechanics C.P.S., University of Zaragoza, 50015 Zaragoza, Spain Received 20 July 1998; accepted 12 January 1999
Abstract In many countries, the determination of the parameters of the wave, likely to be produced after the failure of a dam, is required by law, and systematic studies are mandatory. There is a necessity to develop adequate numerical solvers which are able to reproduce situations originated from the irregularities of a non-prismatic bed and to model the complete equations that progress despite the irregular character of the data. Many hydraulic situations can be described by means of a one-dimensional (1D) model, either because a more detailed resolution is unnecessary or because the flow is markedly 1D. Many techniques have been developed recently for systems of conservation laws in 1D (in the context of gas dynamics). Some years after their adoption for solving problems in gas dynamics, upwind and total variation diminishing (TVD) numerical schemes have been successfully used for the solution of the shallow water equations, with similar advantages. Their use is nevertheless only gradually gaining acceptance in this sector. The performance of some of these techniques for practical applications in river flow is reported in this work. 䉷 1999 Elsevier Science B.V. All rights reserved. Keywords: Flood routing; Explicit TVD methods; Unsteady shallow water flows; Dam break flow
1. Introduction Many problems of river management and civil protection consist of evaluation of the maximum water levels and discharges that may be attained at particular locations during the development of an exceptional meteorological event. There is another category of events of catastrophic nature whose effects also fall into the civil protection area. It is the prevision of the scenario subsequent to the almost instantaneous release of a great volume of liquid. The situation is that of the breaking of a man made dam. In many countries, the determination of the parameters of the wave, likely to be produced after the * Corresponding author. Fax: ⫹34-976761882. E-mail address: [email protected] (P. Garcia-Navarro)
failure of a dam, is required by law (Molinaro,1991; Betamio de Almeida and Bento Franco, 1993), and systematic studies are mandatory. There are works based on scaled physical models of natural valleys, but they represent expensive efforts that are not devoid of difficulties. Hence, there is a necessity to develop adequate numerical models which are able to reproduce the situations that originate from the irregularities of a non-prismatic bed. It is also necessary to trace their applicability considering the difficulty of developing a model capable of producing solutions of the complete equations despite the irregular character of the river bed. For many practical applications, it is accepted that the unsteady flow of water in a one-dimensional (1D) approach is governed by the shallow water or St. Venant equations. These represent the conservation
0022-1694/99/$ - see front matter 䉷 1999 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(99)00007-4
228
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 1. Dam break problem. Lax–Friedrichs versus exact solution.
of mass and momentum along the direction of the main flow. The prediction of water profiles and discharges for unsteady as well as steady situations of hydraulic systems is often linked to the use of numerical methods able to solve the system of equations. As in any other case of integration of a system of non-linear partial differential equations, the first difficulty is the choice of the numerical scheme. There is, at present, a wide range of techniques suitable for subcritical shallow water flows. However, it has traditionally been difficult to have a single method able to reproduce automatically any general situation.
The numerical modelling of unsteady flow in rivers is a complicated task and the difficulties grow as the pretensions to obtain better quality or more general solutions do (Cunge et al., 1980) There exist good methods, developed to deal with Gas Dynamics problems (Euler equations), able to cope with complex systems of discontinuities and shock waves (Hirsch, 1990; MacCormack, 1971; Roe, 1989). The shallow water, or St. Venant equations, being a hyperbolic partial differential system, represent a good candidate for the application of many of these techniques. Many new schemes have been
Fig. 2. Dam break problem. MacCormack versus exact solution.
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
229
Fig. 3. Dam break problem. TVD MacCormack versus exact solution.
reported successful for the flow in channels (Glaister, 1988; Alcrudo and Garcia-Navarro, 1992; GarciaNavarro and Alcrudo, 1992). However, their application to the river flow is not so common in the literature (Jin and Fread, 1997) and their adaptation to river hydraulics is not straightforward because of the irregular topography. The presence of extreme slopes, high roughness and strong changes in the irregular geometry represent a great difficulty that can lead to important numerical errors presumably arising from the source terms of the equations. In this work, the performance of three finite difference techniques is reported and compared so that their suitability for this kind of problems can be discussed. First, the system of equations modelling 1D unsteady river flow will be presented in Section 2. The finite difference numerical techniques applied to solve the equations are the Lax–Friedrichs scheme and the explicit MacCormack scheme with and without TVD corrections. TVD corrections produce a variation of the classical predictor–corrector scheme based on the theory of total variation diminishing methods and introduce upwind effects in the solution while maintaining central difference accuracy; They are described in Section 3. Section 4 starts with a validation of the numerical results provided by the three schemes against an exact solution for a problem of an ideal dam break. Then, the numerical treatment applied to the practical simulation of river flow will be discussed. Finally, in Section 5, the results
obtained for the simulation using data from real rivers are presented and the performance of the numerical schemes is compared. 2. Mathematical model Many hydraulic situations can be described by means of a 1D model, either because a more detailed resolution is unnecessary or because the flow is markedly 1D. The fundamental hypothesis implied in the numerical modelling of river flows are formalised in the equations of unsteady open channel flow. Under the St. Venant hypotheses, 1D unsteady flow can be written in the form: 2U 2F ⫹ R; 2t 2x
1
with U A; QT ; Q2 ⫹ gI1 F Q; A
!T ;
R 0; gI2 ⫹ gA S0 ⫺ Sf T which emphasises the conservative character of the system in the absence of source terms. The effects of the wind as well as those of the Coriolis force have been neglected and no lateral
230
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 4. Cross-section profiles at different locations. Aragon river.
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
231
Fig. 5. Water profiles at different times, Aragon river.
inflow/outflow is considered. In Eq. (1), A is the wetted cross-sectional area, Q, the discharge and g, the acceleration due to gravity. I1 represents a hydrostatic pressure force term as described in Cunge et al. (1980). The pressure force can have a component in the direction of the main stream due to the reaction of the walls in case of variations in shape along this direction. The amount of this force depends on the cross-sectional variation for constant depth. It is important to note that the validity of this approach is linked to the hypothesis of gradual variation. If sudden expansions or contractions take place, the approach is not valid. I2 accounts for the pressure force in a volume of constant depth h due to longitudinal width variations. The bed slope is S0. Sf stands for the energy grade line and is defined, for example, in terms of the Manning’s roughness coefficient n: Sf
Q兩Q兩n2 ; A2 R4=3
with R A=P; P being the wetted perimeter.
The Jacobian matrix of the system (Eq. (1)) is ! 0 1 2F A ; 2U c2 ⫺ u2 2u p where u Q=A and c gA=b. The eigenvalues and eigenvectors of A are: a1;2 u ^ c e1;2 1; u ^ cT : Other formulations are possible and frequently encountered. 3. Numerical model There is a long list of finite difference techniques suitable for the numerical solution of the equations presented (Cunge et al., 1980; Glaister, 1998; Garcia-Navarro and Saviron, 1992; Alcrudo and Garcia-Navarro, 1992). If there is an interest to deal with unsteady problems that may contain rapidly varied flows with
232
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 6. Water depth history at two points, Aragon river.
Fig. 7. Discharge history at two points, Aragon river.
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
233
Fig. 8. Mass balance error, Aragon river.
the possible presence of transcritical flow or formation of surges, attention must be focused on a class of these techniques: those able to progress automatically towards a weak or discontinuous solution of the equations. Among these, there is still a variety of methods. The explicit finite difference schemes used for the numerical computations in this work are presented next.
3.1. Lax–Friedrich’s scheme This is a first order in space and time finite difference technique. For a general homogeneous conservation system
updating: Un⫹1 uUni ⫹ i ⫺
1⫺u n Ui⫹1 ⫹ Uni⫺1 2
Dt Fn ⫺ Fni⫺1 : 2Dx i⫹1
0 ⱕ u ⬍ 1: The value u 1 renders the scheme unstable. More stability and more numerical diffusion is introduced as u approaches zero. A value u 0.1 is usually adopted. It is easy to implement and very robust. The interest of considering its performance was to establish the extent of the numerical error introduced for problems highly dominated by the source terms. 3.2. MacCormack’s scheme
2U 2F ⫹ 0; 2t 2x the procedure to update one time step Dt the interior points 2,…, N-1 of a regular grid is based on a nodal
A classical finite difference scheme suitable for the discretisation of (Eq. (1)) is the explicit two step predictor–corrector MacCormack method (MacCormack, 1971). It is a technique of proved efficiency for unsteady open channel flow modelling
234
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 9. Bed level profile, Neila river.
(Garcia-Navarro and Saviron, 1992; Garcia Navarro and Zorraquino, 1993; etc). For linear problems it reduces to the Lax–Wendroff scheme so that, the Courant–Friedrich–Lewy criterion Dt ⱕ
Dx max{ai }
ensures numerical stability. Being second order accurate in space and time, it offers good resolution and has a great conceptual simplicity. It is composed of a sequence of two sub-steps in which the spatial derivatives are taken in alternate directions as follows: Predictor step ! Upi Uni ⫺
Dt n F ⫺ Fni Dx i⫹1
Corrector step ! Uci Upi ⫺
Dt p F ⫺ Fpi⫺1 Dx i
so that the updated solution is given by Uin⫹1
1 p U ⫹ Uci : 2 i
It is well known, however, that classical second order schemes show oscillatory behaviour near discontinuities, and McCormack’s scheme is no exception. This is why the introduction of artificial viscosity became usual in reported numerical simulation practice (Garcia and Kahawita, 1986). The difficulty with the artificial viscosity approach is that it is hard to determine an appropriate form that introduces just enough dissipation to preserve monotonicity without causing unnecessary smearing. 3.3. Total variation diminishing MacCormack In the 1D theory of TVD methods, it is shown that some classical schemes can be reformulated in a way that improves their performance without great effort. In particular it is possible to add a conservative and non-linear modification to the Lax–Wendroff scheme to render it non-oscillatory at points of extrema without losing its accuracy in other regions of the flow. The application of this idea to the MacCormack version of the Lax–Wendroff scheme for the 1D
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
235
Fig. 10. Initial top width, Neila river.
shallow water equations is reported in (GarciaNavarro and Alcrudo, 1992). The TVD correction can be cast as a modification of the last step of MacCormack scheme with inclusion of new terms to eliminate the oscillations. Un⫹1 i
1 p Dt n U ⫹ Uci ⫹ D ⫺ Dni⫺1=2 : 2 i Dx i⫹1=2
The form of the new terms is based on Roe’s linearisation of the system (Eq. (1)). According to Roe (1989) a matrix A can approximate locally the Jacobian A of the system. The eigenvalues a k and eigenvectors ek of A satisfy X k k dUi⫹1=2 Ui⫹1 ⫺ Ui a e i⫹1=2 k
dFi⫹1=2 Fi⫹1 ⫺ Fi
X
a k ak ek i⫹1=2
1D shallow water equations (k 2) in Alcrudo and Garcia-Navarro (1992). The particular expressions used in the present work are ~ a 1;2 u~ ^ c; ~ T; e1;2 1; u~ ^ c
a1
a2
~ dA ⫺ dQ c~ ⫹ u : 2c~
The average velocity and celerity at the interface i ⫹ 1/2 are p p Qi⫹1 Ai ⫹ Qi Ai⫹1 p p ; u~ p Ai Ai⫹1 Ai⫹1 ⫹ Ai
k
inside every computational cell (i,i ⫹ 1). Expressions for ak ; ak and ek can be found, for the 1D Euler equation (k 3) in Roe (1989) and for the
~ d A ⫹ dQ c~ ⫺ u 2c~
c~
r g A=bi ⫹ A=bi⫹1 2
and dA Ai⫹1 ⫺ Ai , dQ Ai⫹1 ⫺ Ai .
236
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 11. Initial Froude number distribution, Neila river.
The third step is then made using: 1X k Dt k 兩a i⫹1=2 兩 ai⫹1=2 兩a ki⫹1=2 兩 1 ⫺ Di⫹1=2 2 k Dx 1 ⫺ Fki⫹1=2 eki⫹1=2 : The function F is a flux limiter responsible for the smoothing near fronts. Mathematical expressions for this function can be found in Sweby (1984). 4. Applications 4.1. Dam break test case The idealised dam break problem was chosen because it is a classical example of non-linear flow with shocks to test conservation in numerical schemes and, at the same time, has an analytical solution (Stoker (1957). This problem is generated by the homogeneous 1D shallow water equations for the ideal case of a flat and
frictionless channel of unit width and rectangular cross-section, with the initial conditions. 8 L > > if x ⱕ < hL 2 h x; 0 : > L > : hR if x ⬎ 2 Q x; 0 0: If the calculation times used are so as to avoid interaction with the extremities of the channel, the boundary conditions are trivial. Figs. 1–3 show a sample of the solution for a value of the initial height ratio hL:hR 20. They represent the profiles of the water surface 30 s after the dam break as computed on a grid of 200 cells. The analytic solution appears as a continuous line. The dotted lines represent the numerical solutions obtained with Lax– Friedrich’s scheme in Fig. 1, MacCormack’s scheme in Fig. 2 and TVD MacCormack’s scheme in Fig. 3. The excessive numerical diffusion introduced by the first order scheme and the lack of robustness of the
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
237
Fig. 12. Discharge history at Neila village, Neila river.
second order scheme are shown. The TVD step provides the necessary numerical viscosity to cope well both with the shock front and with the critical point at the former dam location. The size of the time step is chosen from the condition CFL 0.9. 4.2. Application to river flow The numerical scheme applied must also be able to reproduce the surface discontinuities that may appear as a result of the irregular geometry itself or because of the hydraulic structures existing along the river. Independently of the particular problem, a preprocessing and storage of the geometrical and hydraulic data is always required. The resolution is compatible with their usage in the form of numerical functions (data tables). They furnish the information required by the equations through a process of interpolation. 4.3. Geometrical data The 1D model requires the selection of a series of
river cross-sections each one representing the geometry of a reach of length Dx. The approximation of every cross-sectional shape as a juxtaposition of trapezoids simplifies the systematic calculation of the values of cross-section area, surface width or wetted perimeter as increasing functions of the water depth for a sequence of discrete values. At the same time, the consistency of the data is ensured as independent calculations of these functions may lead to numerical instabilities (Bento Franco and Betamio, 1991). The simulation starts from a definition of the topographic features in the form of numerical functions (data tables) of the depth and the distance along the river. These are matrices of dimensions (number of cross-sections x number of water depth levels). The data tables do not correspond in general to equidistributed points along the river and never to the computational grid positions. One option is the numerical generation of intermediate sections by the interpolation of the surveyed cross-sections. However, the particular value of any of the functions at a nodal
238
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 13. Water depth history at Neila village, Neila river.
position for an arbitrary water depth can be obtained each time via interpolation from the data tables. More or less sophisticated interpolations can be performed. Considering that any interpolation may introduce numerical errors in the results and the irregular character of the river, a liner interpolation between data tables was used in this work. The bed slopes were determined using the values of the bottom level at each cross-section. Another difficulty arises from the presence of functions of these variables such as the integrals I1 and I2 in the system of equations. Their evaluation in a case of simple geometry is direct and analytical. The irregular character of the river cross-sections requires some kind of approximation. Considering a cross-section well approximated by a polygonal line of j vertices, the integral of the hydrostatic pressure at a section is estimated as I1 h
Zh 0
h ⫺ hb hdh
jX ⫺1
Zhk ⫹ 1
k0
hk
h ⫺ hb hdh;
where inside every interval (k, k ⫹ 1) b h b hk ⫹
h b hk⫹1 ⫺ b hk ; hk⫹1 ⫺ hk
hk ⱕ h ⱕ hk⫹1 : The pressure force from the walls in an interval Dx is estimated from I2
1 I h; x ⫹ Dx ⫺ I1 h; x: Dx 1
The numerical schemes presented are just defined for the homogeneous part of the system of equations. It still remains the question of the best discretisation of the source terms. The simplest pointwise option for a generic source term function R A; Q can be expressed as ÿ Ri ⬇ R Ani ; Qni ; which consists of a simple evaluation of the source term functions at the nodal point using the already known dependent variables. This is the procedure
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
239
Fig. 14. Discharge history at Villarvelayo village, Neila river.
applied to the difference schemes tried. Other possibilities include some form of average between neighbouring nodal points or an upwind treatment (Glaister, 1998). Problems may arise when the source terms contain spatial derivatives and a definition of their discrete representation must be established. This difficulty is met for the discretisation of the bottom slope (S0) and variable width (I2) source terms. When dealing with the Lax–Friedrichs scheme, a standard central difference was applied. The option selected in the case of the MacCormack scheme was 2b b ⫺ bi ⬇ i⫹1 Dx 2x 2b bi ⫺ bi⫺1 ⬇ Corrector step ! Dx 2x Predictor step !
and similarly for 2z=2x: The limitation imposed by the friction term can be relaxed if this term is not discretised using an explicit approach. A pointwise semi-implicit discretisation of
the friction term (Sf) proved efficient in reducing instabilities, Huang and Song (1985). Q兩Q兩 AR2=3 ; K 2 n K n 兩Q兩 uQni ⫹ 1 ⫺ uQn⫹1 ; Sfi ⬇ i K2 i Sf
0 ⱕ u ⱕ 1:
This strategy can be very easily implemented in an explicit procedure. The term containing Q n⫹1 joins the left hand side so that the updating expression for the discharge is modified but the method remains explicit. The weighting parameter u controls the degree of implicitness. It must also be signaled that in this 1D model a single value of the roughness coefficient was used. The mathematical model also considers a single average water velocity at every cross-section and no distinction was made between river bed and flood plain. A uniform distribution of nodal points was considered. All these limitations in the model require
240
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 15. Water depth history at Villarvelayo village, Neila river.
a careful analysis of the results and invite further improvements.
5. Results 5.1. Case 1 Taking the topographical data of an upper reach of about 14 440 m in the river Aragon in the Spanish Pyrenees, an application of the numerical simulation is performed for the estimation of the propagation of a dambreak wave. A 15 m high dam is placed at the location x 8 000 m and the hypothesis of total and sudden failure is made. The interest is in the performance of the numerical techniques. The roughness coefficient is set to n 0.05. The initial discharge is 50 m 3/s and a water profile upstream and downstream of the dam provides the initial conditions for depth. The reach is discretised in 500 intervals (Dx 28 m). Fig. 4 displays some of the measured transverse
cross-sections along the river reach. Fig. 5 represents a series of computed longitudinal river profiles for subsequent instants of time starting with the initial condition. The time step size was the maximum allowed by numerical stability. The bottom level profile is also displayed in the figure. The profiles correspond to the results obtained with the MacCormack scheme. They do not differ very much from those supplied by the other schemes in this case. Figs. 6 and 7 show depth and discharge histories, respectively, at two locations along the reach (x 8500 and 1200 m). The predictions provided by MacCormack are compared to those from Lax–Friedrich’s scheme. Some differences between them are noticeable, in particular the smearing introduced by the latter scheme. The consequence of this smearing is quantified in the mass balance error. As we are dealing with a fluid of constant density, mass balance is equivalent to volume balance. It is calculated as the difference between the volume variation along the reach and the net volume inflow through the boundaries. In order to deal with a dimensionless quantity,
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
241
Fig. 16. Froude number distribution, Neila river.
the difference is calculated relative to a meaningful volume of water, the initial volume in this case. The estimation of the error in mass conservation introduced by the numerical calculations is shown in Fig. 8 using the two classical schemes. No noticeable improvement was introduced by the TVD correction in this case. It is worth noting here that the considerable difference in the error produced by the two numerical schemes as shown in Fig. 8 may seem excessive if Fig. 7, for instance is taken as a reference. The reason that justifies this discrepancy, is the different nature of the two plots. Fig. 7 shows the discharge history at a point whilst Fig. 8 represents the global volume error evolution in time. Spurious discharge oscillations at a single point along the channel make this global error accumulate though the solution is well behaved in the rest of the domain. 5.2. Case 2 Hydraulic phenomena of marked unsteady character can be included under the form of a more or less
steep inflow hydrograph Q Q t. Predicting the evolution of the flooding wave that it represents is the aim of the method. From a mathematical point of view, the transit from a value of the discharge to a higher one can happen in different ways. The parameter T (characteristic time of breach formation), with T 0 in the case of a sudden release of water, can be useful and may help to analyse aspects and consequences of the influence of this variable. Our case 2 deals with such a kind of event at the upstream inflow point of Neila river (Burgos, Spain). A reach of 18000 m of this river measured from its upstream dam was studied to analyse the flood risk at Neila village, located 5350 m downstream from the dam. Several dam break and friction roughness hypotheses were made and the corresponding discharge hydrographs were calculated. From the set of results obtained only a few were selected for inclusion in the present work; those corresponding to the more extreme breaking process with T 0.94 h and a peak discharge of 160 m 3/s.
242
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 17. Froude number distribution, Neila river.
Manning roughness was assumed to be uniform and equal to 0.05. The topographic characteristics of this river are represented in Fig. 9 (bottom level variation) and Fig. 10 (plan view of the breadth variations along the axis of the river). These, together with the small value of its base discharge (Q 1 m 3/s) and the corresponding initial values of the Froude number (Fig. 11) rendered the numerical computation more challenging than the former case and accentuated the differences among the capabilities of the numerical schemes. Fig. 12 displays the discharge hydrographs at Neila village location as calculated from the three schemes. Fig. 13 shows the equivalent plot of the depth histories. In the same order, Figs. 14 and 15 show the discharge and water depth changes in time at Villa velayo, a village located 16 875 m downstream from the dam site. Figs. 16 and 17 show the distribution of the Froude number along the river at four different times and clearly indicate the transcritical character of the flow. It is interesting to point out here that the upstream
boundary condition consisted of the inflow hydrograph Q Q t. This value was complemented with the value of the water depth at the upstream point obtained from the characteristics theory whenever there was subcritical flow at the point. The value of the Froude number at the second nodal point was continuously controlled for this purpose. In case of supercritical conditions, the equation S0 Sf supplied the required water depth, that is, normal depth was imposed. The TVD correction of MacCormack scheme was crucial in stabilizing the calculation in this example and was the only way to return to the correct base flow (initial steady state) after the passage of the flood wave. The simple Lax–Friedrichs scheme proved very inefficient in dealing with the dominating source terms present in this river reach. In Fig. 18 the improvement provided by the TVD scheme is clear. The discussion concerning the important discrepancy between the magnitude of this global error and that of the discharge hydrographs at the particular locations is again as in the previous test case.
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
243
Fig. 18. Mass balance errors, Neila river.
5.3. Case 3 The third example is based on one of the cases presented by Jin and Fread (1997). These authors use a numerical formulation based on flux splitting and pointwise evaluation of the source terms. In an attempt to validate our models against those present in the literature, one of their dam break problems was calculated using the TVD McCormack scheme. This is a 36.6 m high dam located at 11.3 km in a 32.2 km reach of nearly uniform trapezoidal crosssection shape. A constant Manning friction factor n 0.04 is used. The dam fails according to a breach law that depends on a failure time Tf. Tf 0.02 h in this example as this case was the least favourable leading to the formation of supercritical regions near the dam. The breaching dam is included in the calculation as an internal boundary separating two reaches upstream and downstream from it. One mass continuity equation and a dynamic weir type equation are used together with two numerical boundary conditions (characteristics).
Fig. 19 displays several longitudinal water profiles and the river bed profile. They correspond to the times indicated on the figure. The mass error balance, using as reference volume the initial volume of water, gives a value bounded by 5%. Fig. 20 shows the discharge hydrographs at three locations, one near the dam (x 11.3 km) and two downstream from it (x 16.1 km, x 24.1 km). The results from Jin and Fread are plotted in lines and compared with the numerical results from our model. Figs. 21 and 22 show the water depth and the Froude number distributions in space downstream from the dam at 4 different times after the beginning of the breach formation (t 0.05, 0.10, 0.15, and 0.2 h). There are some differences arising from a different water depth value at the right hand side of the dam. They might be due to a difference in the numerical boundary conditions at that section. Whilst the reference authors use a discretisation of the mass conservation equation, we apply a semi-Lagrangian treatment, that is, the trajectory of the backwards characteristic is used to bring information from the old to the new time level. This test case
244
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 19. Longitudinal profiles at different times.
is presented by Jin and Fread as an example in which, because of the short time of breach formation, regions of supercritical flow appear near the dam location. Unfortunately, no experimental data are supplied.
6. Conclusions A first order (Lax–Friedrich’s), a second order (MacCormack) and a high resolution (TVD MacCormack) explicit finite difference schemes are applied to the simulation of 1D unsteady flow in rivers. During the unsteady processes, some of the selected examples contain zones of supercritical flow that the schemes are able to cope with automatically. The presence of extreme slopes and strong changes in the irregular geometry represent a great difficulty that can lead to important numerical errors presumably arising from the source terms of the equations. The TVD version of the MacCormack scheme appears as superior in this kind of problems. The conceptual and programming simplicity of MacCormack scheme
is not greatly altered by the inclusion of the TVD correction. For smoother situations, the three schemes performed well. There is a difference in mass conservation ability between Lax–Friedrichs and MacCormack for a given grid size, with a clear superiority of the latter. TVD MacCormack does not present a remarkable superiority in these cases. In what concerns the particular problem of the dam break, given the initial conditions of the system and a temporal law for the breach formation (usual treatment of the problem in the literature) the simulation furnishes simultaneously the characteristics of the propagation of the wave downstream and the evolution of the water profile upstream of the dam site. In contrast to other numerical techniques incorporated in commercial codes (DAMBRK Fread, 1985), the procedure applied in the present work deals very efficiently with the presence of locally supercritical conditions anywhere along the river. The time step size is rather restricted in general to ensure stability. In this sense, irregular grids and
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
245
Fig. 20. Discharge hydrographs at x 11.3 km (dam site), x 16.1 km and x 24.1 km. Jin and Fread: Lines. Present model: Symbols.
Fig. 21. Depth profiles at t 0.05, 0.10, 0.15 and 0.020 h after the beginning of the breach formation. Jin an Fread: Lines. Present model: Symbols.
246
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247
Fig. 22. Froude profiles at t 0.05, 0.10, 0.15 and 0.20 h after the beginning of the breach formation. Jin and Fread: Lines. Present model: Symbols.
implicit versions of the schemes are envisaged as future work. Acknowledgements We would like to thank SERS S.A. for the partial funding of the second author and Mr. V. Zorraquino from SERS S.A. for many useful and interesting discussions during the course of this work. Thanks also to Dr. Jin from NOAA for supplying the necessary data from his own work. References Alcrudo, F., Garcia-Navarro, P., 1992. Flux difference splitting for 1D open channel flow equations. Int. J. for Numerical Methods in Fluids 14, 00. Bento Franco, A., Betamio, A., 1991. Simulaco uni e bidimensional de cheias provocadas por roturas de barragens em planicies de inundaco, V Simposio Luso-Brasileiro de Hidraulica e Recursos Hidricos SILUSB.
Betamio de Almeida, A., Bento Franco, A., 1993. Modelling of dam break flows. In: Chaudhry, M.H., Mays, L.W. (Eds.). Computer Modelling of Free Surface and Pressurized Flows, NATO ASI Series, 343–375. Cunge, J.A., Holly, F.M., Verwey, A., 1980. Practical Aspects of Computational River Hydraulics. Pitman, London. Fread, D.L., 1985. Channel routing. In: Anderson, M.G., Burt, T.P. (Eds.). Hydrological Forecasting, Wiley, New York. Garcia, R., Kahawita, R.A., 1986. Numerical solution of the St.Venant equations with the MacCormack finite difference scheme. Int. J. for Numerical Methods in Fluids 6, 259– 274. Garcia-Navarro, P., Alcrudo, F., Saviron, J.M., 1992. 1D open channel flow simulation using TVD McCormack scheme. J. of Hydraulic Engineering, ASME 118 (10). Garcia-Navarro, P., Saviron, J.M., 1992. McCormack’s method for the numerical simulation of one-dimensional discontinuous unsteady open channel flow. J. of Hydraulic Research 30, 95– 105. Garcia Navarro, P, Zorraquino, V., 1993. Numerical modelling of flood propagation through a system of reservoirs. J. of the Hydraulic Division, Am. Soc. of Civil Engineering 119 (3), 380–389. Glaister, P., 1998. Approximate Riemann solutions of the shallow water equations. J. of Hydraulic Research 26 (3), 293–306.
P. Garcia-Navarro et al. / Journal of Hydrology 216 (1999) 227–247 Hirsch, C., 1990. Numerical Computation of Internal and External Flows, 2, Wiley, New York. Huang, J., Song, C.C.S., 1985. Stability of dynamic flood routing schemes. J. of Hydraulic Engineering, ASCE 111 (12), 1497– 1505. Jin, M., Fread, D.L., 1997. Dynamic flood routing with explicit and implicit numerical solution schemes. J. of Hydraulic Engineering, ASCE 123 (3). MacCormack, R.W., 1971. Numerical solution of the interaction of a shock wave with a laminar boundary layer. In: Holt, M. (Ed.). Proceedings 2nd International Conference on Numerical Methods in Fluid Dynamics, Springer, Berlin, pp. 151.
247
Molinaro, P., 1991. Dam break analysis: a state of the art. In: Bensar, D., Brebbia, C.A., Ovazar, D. (Eds.). Computational Water Resources, CMP. Roe, P.L., 1989. A survey on upwind differencing techniques, Lecture Series in CFD, Von Karman Institute for Fluid Dynamics, Belgium. Stoker, J.J., 1957. Water Waves. Interscience, New York. Sweby, P.K., 1984. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J. of Numerical Analysis 21, 995–1011.