Iterative explicit simulation of 1D surges and dam-break ows
Chih-Tsung Hsu
∗;†and Keh-Chia Yeh
Department of Civil Engineering; National Chiao Tung University; Hsinchu; 30050 Taiwan; Republic of China
SUMMARY
The one-dimensional Saint Venant equations for shallow-water ows are used to simulate the ood wave resulting from the sudden opening (or closure) of a gate or collapse of a dam. An iterative explicit characteristics-based <nite-di=erence method, based on the explicit <nite analytic method, is proposed to discretize the dynamic equation, and the conservative control volume method is used for the discretization of the continuity equation. Surge and dam-break ows in a horizontal, rectangular and frictionless channel were <rst considered, under such conditions the analytic solutions exist. For the surge simulation, numerical results of the proposed scheme are nearly identical to those obtained from the Preissmann scheme. For the dam-break simulations addressing three ratios of tailwater depth to water depth in the reservoir, the proposed scheme, as compared with the analytic solutions, yields better results than those obtained by the MacCormack scheme, the Gabutti scheme, and Jha et al.’s ux splitting scheme (J. Hydraul. Res. 1996; 34(5):605–621). As the depth ratio approaches zero, the accuracy of the proposed scheme is still satisfactory, even with the dry-bed condition. Investigations then were made for more realistic dam-break ow waves propogating in a sloped and frictional channel. Lacking analytic solutions, the simulating results from the proposed scheme as well as those from Chen’s scheme (J. Hydraul. Div. 1980; 106(HY4):535–556) were compared with the laboratory data collected in 1960– 1961 at the United States Army Engineer Waterways Experiments Station (WES). An assumed initial ow was required for the computer-simulated condition in Chen’s model. However, this is not the case in the proposed model, i.e., a real dry-bed condition was set as the initial condition in the downstream channel of the dam. The consistency between the two simulated results is obvious compared with the experimental data. Copyright ? 2002 John Wiley & Sons, Ltd.
KEY WORDS: iterative explicit scheme; surge; dam-break wave; explicit <nite analytic method; staggered-grid
1. INTRODUCTION
The task of estimating the movement of a surge (or shock) or a dam-break wave, resulting from the sudden upstream opening (or the sudden downstream closure) of a sluice gate for emergencies or dam failures, has occupied the attention of researchers as well as practicing engineers for several decades. The determination of the surge height at di=erent locations along the channel provides important information for the design of the bank height. The
∗Correspondence to: Chih-Tsung Hsu, Department of Civil Engineering, National Chiao Tung University, 1001 Ta
Hsueh Road, Hsinchu, Taiwan 30050.
†E-mail: [email protected]
dreadful disaster due to the dam-break ood wave reminds the decision- makers to pay more attention to the dam-safety problem.
The rapid varied ows in open channels, with steep fronts resulting from surges or dam-break waves, are often simulated as one-dimensional (1D) ows. Such ows can be described by the Saint Venant equations (shallow-water equations), which are a set of non-linear hy-perbolic partial di=erential equations. However, the analytic solutions of these equations are not available, except for certain special simpli<ed conditions. Therefore, consistent e=orts have been made to develop numerical schemes to solve the equations. There are three basic approaches to compute the equations of ows with steep fronts: the shock <tting, pseudo-viscosity, and ‘through’ (shock capturing) methods [1]. The shock <tting method requires internal boundary conditions to determine the position and the associated ow characteristics of the shock, yielding a complicated solution algorithm. The pseudo-viscosity method, which is a common practice adopted in the non-dissipative <nite-di=erence schemes, operates by adding an arti<cial di=usion term to suppress the oscillation near the front. The ‘through’ method is most often used in practice because it solves the Saint Venant equations directly without any particular arrangements of the algorithm. The proposed scheme in this paper belongs to the ‘through’ method.
Many shock capturing schemes are available as described in the literature below. To check their applicability to the real conditions of wave propagations, the simulated results were usually compared to the analytical solutions of the shocks propagating in a 1D horizon-tal, frictionless channel <rst, then to the experimental data in a sloped, frictional channel. The numerical schemes using the <nite-di=erence method include the Preissmann four-point scheme [2], Holly–Preissmann two-point together with the reach-back characteristics scheme [3], MacCormack scheme [4–7], Lambda scheme [7], Gabutti scheme [7], Beam-Warming scheme and its modi<cations [8–10], ux di=erence splitting scheme [11–14], Godunov-type upwind scheme with a Riemann solver [13; 15–19], modi<ed Godunov scheme [20], total vari-ation diminishing (TVD) scheme [21–23], and semi-Lagrangian scheme [24]. Additionally, the numerical schemes using the <nite-element method include the Eulerian–Lagrangian linked Galerkin scheme [25], the dissipative Galerkin scheme [26–28], Petrov–Galerkin scheme [18; 29], Taylor–Galerkin scheme [30], and Addcollocation scheme [31]. Additionally, the explicit scheme by Cheng-lung Chen mentioned above used the characteristics method.
The numerical schemes used for solving the dam-break problem face severe challenges when they are applied to cases characterized by large upstream water depths (hu) in the
reservoirs and small downstream water depths (hd) in the channels. A mixed- ow regime
(i.e., supercritical and subcritical ows co-exist) occurs for ows in the horizontal, frictionless channels when the value of hd=hu is smaller than 0.138 [32]. This ow regime is diMcult to
be handled by the numerical schemes. Furthermore, the height, shape, and celerity of the simulated front under such a condition may signi<cantly deviate from the exact solution [11]. Under the extreme condition of dry bed (hd= 0); few of the above mentioned numerical
schemes [15; 33] can work because it is not easy to deal with the moving boundary condition of zero water depth. To cope with this diMculty, a commonly used technique in numerical computation is to assume that a minimum water depth or discharge exists on the dry bed of the channel. As pointed out by Zhang et al. [34], the accumulated error due to this assumption as the computation proceeds might eventually distort the solutions and hence render the numerical solution near the wave front doubtful. Therefore, Zhang et al. [34] proposed an implicit staggered-grid scheme, analogous to the SIMPLER method [35], to analyze the dam-break
problems on dry beds. Bellos and Sakkas [36] employed the explicit MacCormack scheme with an expedient treatment to determine the computational grid at the boundary, wet or dry, to solve the same problem. Moreover, Di Monaco and Molinaro [37] developed a Lagrangian model for simulating the dam-break waves on dry beds using the <nite-element method.
The purpose of this paper is to propose an iterative explicit characteristics-based <nite-di=erence scheme suitable for the computation of surges and dam-break waves propagating both in rectangular, horizontal, and frictionless channels and in sloped, and frictional channels. The proposed scheme is a modi<cation of the explicit <nite analytic (FA) method, which was <rst employed by Dai [38] for the computation of 2D and 3D cavity ows. In his scheme, based on the ow conditions of the previous time step, the convective transport equation is solved with a local analytic solution and the viscous di=usion and source terms are approximated by the <nite di=erences. In this paper, the time dependent variables, including the linearized characteristic curve, are updated by the newly calculated ow conditions in the process of iterative computation during a time step. The superiorities of accuracy and simpli<cation of the proposed scheme to other numerical schemes for the dam-break problem can be seen under various cases as well as the most critical dry-bed case.
2. GOVERNING EQUATIONS
The 1D unsteady free-surface ow in a wide rectangular channel without lateral ow can be described by the well-known Saint Venant equations, those based on the conservation of mass and momentum of the conservative form:
@h @t + @q@x= 0 (1) @q @t + @@x uq + 1 2gh2 = gh(S0−Sf) (2)
where h is ow depth; u is ow velocity; q is volumetric ow rate per unit width; t is time; x is distance along the channel; g is gravitational acceleration; S0 is bed slope; and Sf is
frictional slope. The key assumptions in the governing equations are the hydrostatic pressure distribution and small bed slope. The applicability of the Saint Venant equations to surges and dam-break ows with steep fronts, i.e., those with large surface curvatures, has been examined by Basco [39]. He concluded that the Saint Venant equations are valid if the wave period of the input hydrograph is greater than about 100 seconds. Fortunately, most ood waves resulting from dam breaks of practical interest have periods much larger than this value.
For open channel ows with steep fronts, the convective terms in Equation (2) dominate. Thus, Equation (2) can be written in the total derivative form based on the concept of the characteristics method Dq Dt =−gh @h@x−q @u@x+ gh(So−Sf) (3 ) along dx dt = u (4)
Figure 1. Sketch of characteristic curve and staggered-grid system.
where Dq=Dt is the total derivative of q along the characteristic curve described by Equation (4). The explicit FA method [38] assumes that the variation of u in Equation (4) and the source terms on the right-hand side of Equation (3) are relatively small in a given time increment Stn (Figure 1), and can be simpli<ed as known constants. Consequently,
Equation (3) reduces to a linear hyperbolic partial di=erential equation and has the analytic solution (q 2) at time step n if the initial condition (q 1) at the previous time step (n−1)
is speci<ed. In such a case, the characteristic is a straight line. In this paper, however, the time dependence of u and the source terms are further considered. Hence, the integration of Equations (3) and (4) can be solved numerically and can be expressed as
q 2−q 1= t 2 t 1 − gh @h@x+ q @u@x + gh(So−Sf) dt (5) x 2−x 1= t 2 t 1 u dt (6)
where the subscripts 1 and 2 denote the foot and the head of the characteristic, respectively (see Figure 1).
3. ITERATIVE EXPLICIT SCHEME 3.1. Discretization
To avoid the so-called ‘checkerboard’ phenomenon of the computational results, a staggered-grid system (Figure 1) is used in the proposed iterative explicit scheme, which is also adopted
in Dai’s approach [38]. Discretization of Equations (1) and (5) results in hn i+1=2= hi+1=2n−1 − Stn Sx(qn i+1−qn i ) + (1−)StSxn(qi+1n−1−qn−1i ) (7) qn i = q 1− {[ghin(S∇hn ) + qn i (S∇un )−ghn i (So−Sf)ni] + (1−)[gh 1(S∇hn−1) + q 1(S∇un−1)−gh 1(So−Sf) 1]}Stn (8) where Sf= N 2u2 R4=3
and N is Manning coeMcient; R is hydraulic radius; n
( = q; h or u) denotes the newly up-dated value of the variables during the iteration at the present time step n; is the weighting factor in time; Stn is the time interval, which is determined by the Courant condition for
ensuring the stability of the explicit scheme; Sx is the subreach length; and 1 is
approxi-mated by the linear interpolation of n−1
i−1 and n−1i
1= n−1i x 1Sx−xi−1+ n−1i−1xiSx−x 1 (9)
where x 1 is obtained from Equation (6) with the following approximation:
x 1= xi−uni−1=2 Stn (10)
Note that the source terms @h=@x and @u=@x in Equation (8) are discretized by the di=erence scheme, which is represented by the symbol S∇. In the subcritical- ow regime, a central di=erence scheme is adopted to consider the physical phenomenon that the disturbance waves travel both downstream and upstream. On the other hand, due to the unique direction of the traveling disturbance waves in the supercritical- ow regime, the above two terms are discretized by the upwind scheme, i.e., forward or backward di=erence scheme according to the direction of ow velocity.
3.2. Boundary conditions
Unique discrete solutions at each cross-section of a channel at each time step, based on Equations (7) and (8), exist when the initial and boundary conditions are properly speci<ed. The initial conditions include the values of q and h at t = 0. As for the boundary conditions, it should be recalled that the only general technique available for dealing with the boundary conditions is the method of characteristics [3; 5; 40] based on the equations:
@q @t + c−@q@x + a−@h @t + c−@h@x = gh(So−Sf) (11) @q @t + c+@q@x + a+@h @t + c+ @h@x = gh(So−Sf) (12)
where
c−= u−gh; a−= −u−gh;
c+= u +gh; a+= −u +gh:
When the ow is supercritical, two boundary conditions related to q(t) and h(t), or q(h) should be provided at the upstream boundary of the channel and Equations (11) and (12) are used to determine the ow conditions at the downstream boundary. When the ow is subcritical, one boundary condition should be provided at the upstream boundary and Equa-tion (11) is used to determine the other ow condiEqua-tion; same to the downstream boundary but Equation (12) is used. The method has been applied to solve the seaward boundary in predicting the ow characteristics on rough slopes for speci<ed, normally incident wave trains by Kobayashi et al. [40], and the external and internal boundaries for rapidly varying open channel ows by Garcia-Navarro and Saviron [5].
However, under subcritical ow condition, it is usually diMcult to specify the downstream boundary condition, especially for unsteady ow, except that experimental or <eld measured data can be obtained. The assumption of uniform ow condition or the method of extrapolation from the neighbouring points in the computational domain is taken to provide the boundary condition in some numerical methods. In this paper, as will be described below, an approach based on the characteristic concept is introduced to approximate the ow variables, both q and h, at downstream boundary.
In order to determine the two variables q and h at the downstream boundary, one more equation is needed to be solved simultaneously with Equation (12). The continuity equation is rewritten in the following form:
@h
@t + c@h@x= −h@u@x (c = u) (13) Figure 2 shows the grid points and the two corresponding characteristic lines of Equation (12) (labeled c+) and Equation (13) (labeled c) near the downstream end are shown. The values
of variables at points L and R have to be determined to make the solution progress to point P at the present time step n. The algorithm is an iterative approximation technique and is described as following:
(i) The position and corresponding variables of point R are calculated in the same way as Equations (9) and (10) with the characteristic speed c.
(ii) Equation (13) is used to approximate the water depth at point P. hP= hR−hR un−1 M −uM−1n−1 Sx St (14)
(iii) The position and corresponding variables of point L are calculated in the same way as Equations (9) and (10) with the characteristic speed c+.
(iv) With the result of Step (ii), Equation (12) is used to approximate the discharge at point P.
Figure 2. Sketch of characteristic curves near the downstream boundary.
where
f = gh(So−Sf)St
The values of c; c+; a+ and f are calculated by variables at point M for the <rst
approximation, and should be adjusted in the proceeding iterations to obtain the correct ow conditions. Hence the following steps:
(v) The values of c+; a+ and f are recalculated by averaging those at points L and P with
the last approximated dependent variables. c+= 1 2[(u + gh)L+ (u + gh)P] a+= 1 2[(−u + gh)L+ (−u + gh)P] f = 12g{[h(So−Sf)]L+ [h(So−Sf)]P}
(vi) Repeat Steps (iii) to (v) until discharge at point P converges.
(vii) The value of c is recalculated by averaging the last approximated dependent variables at points R and P.
(viii) Repeat Steps (i) to (vii) until the water depth at point P converges. 3.3. Algorithm
Given the proper initial conditions, the solution procedure of the proposed iterative explicit scheme involves the following steps:
(ii) Estimate q 1 using Equations (9) and (10). (iii) Compute qn 2 ; qn 3 ; : : : ; from Equation (8). (iv) Compute hn 3=2; hn
5=2; : : : ; from Equation (7), using the updated qn
i (i = 2; 3; : : :) from Step
(iii).
(v) Repeat Steps (ii) to (iv) until the following criterion is satis<ed:
|(qn
i )new−(qin)old|
(qn
i )old 6; for i = 2; 3; : : : (16) |(hn
i+1=2)new−(hi+1=2n )old|
(hn i+1=2)old
6; for i = 1; 2; : : : (17) where the subscripts ‘new’ and ‘old’ denote the present and the previous iterations, respectively; and is the allowable error, which was set to 10−6 in this study.
(vi) Let qn i+1= (qn
i+1)new; hi+1=2n = (hn
i+1=2)new, and ui+1n = (un
i+1)new; for i = 1; 2; : : :
(vii) Repeat Steps (i) to (vi) for the next time Step (n + 1). 3.4. Stability analysis
From Equations (9) and (10), it can be seen that the foot of the characteristics line, x 1,
must lie between xi and xi−1 to ensure the stability of the convection-dominated momentum
equation as shown in Equation (8). Therefore |u|St=Sx61 is a necessary condition for the stability of Equation (8), which can be proved by the Fourier analysis with the Von Neumann condition (see Appendix A).
Similar to other explicit numerical method applied to the open-channel ow, the Courant condition, de<ned as Cr = (|u|+gh)St=Sx61 by the larger absolute value of the charac-teristic speeds, u±gh, is the limit on determining the time step. When the Courant condition is satis<ed, the stability condition, |u|St=Sx61, for the proposed scheme will be satis<ed, too. Hence the proposed scheme is proved to be valuable under the Courant condition.
4. ANALYTICAL SOLUTIONS
To examine the accuracy of the proposed scheme, the analytical solutions of surges and dam-break ows in horizontal, rectangular, and frictionless channels are used. The exact solution of surge celerity (us) and surge height (h1) can be obtained by simultaneously solving the
following two equations based on the conservation of mass and momentum: us= qh2−q1 2−h1 (18) us= q2 2 h2 + 12gh 2 2 − q2 1 h1 + 12gh 2 1 q2−q1 (19)
where the subscripts 1 and 2 represent the upstream and downstream sides of the surge, respectively. For the dam-break problem, one needs an additional equation describing the movement of the negative simple wave:
u−2c = u0−2c0= const: (20)
where u0 and c0 are the initial velocity and the surface-disturbance wave celerity in the
reservoir, respectively; and c =gh. Thus, the water surface pro<le of the negative simple wave can be determined by
dx
dt= u + c = 3c + u0−2c0 (21) 5. EXAMPLES AND SIMULATION RESULTS
The advantage of the proposed iterative explicit scheme is its simplicity in formulation and fast convergence in numerical iterations. Nevertheless, two parameters, i.e., the Courant number Cr and weighting factor in Equations (7) and (8), in the proposed scheme need to be examined for further applications. Suggestion of the two values is possible by comparing the relative accuracy of the simulated results for cases to analytic solutions. To simplify the analysis and to compare the numerical solutions with the exact solutions of surges and dam-break ows, bed slope and resistance are neglected <rst. The proposed scheme performs well in the simulations of (1) downstream surge propagation due to the sudden opening of a gate, (2) upstream surge propagation due to the sudden closure of a gate, and (3) dam-break ow in a wet or even dry downstream channel bed, when compared with the analytic solutions and numerical solutions obtained by other schemes. Furthermore, to account for more realistic ow situation and to compare with the experimental data of dam-break ow waves, bed slope and resistance e=ects are taken into account.
5.1. Flows in a horizontal and frictionless channel
5.1.1. Sudden opening of a gate. The sudden opening of a gate to convey water into a chan-nel from still waters at a constant depth results in the formation of a surge that propagates downstream. Further opening of the gate will convey even more water and, consequently, create a higher surge that will pass over the previous surge. Figure 3sketches a gate or a dam located somewhere in a channel. The example assumes that the gate is initially closed with a still tailwater of 1 m in depth (h2) and a surge of q1= 20 m2s−1. According to
Equa-tions (18) and (19), one can <nd that us= 8:48 m s−1 and h1= 3:36 m. The traveling distance
of the surge front can be determined by multiplying us by the elapsed time after the sudden
opening of the gate. Note that the surge propagation in this example is in the supercritical- ow regime and is more diMcult to solve from the viewpoint of numerical stability and accuracy.
Numerical testing of various Cr’s and ’s based on the above example shows that Cr = 0:5 and = 0:9 yield the best accuracy in comparison with the analytic solution. Therefore, these two calibrated values are recommended for use in future analyzes on surge problems. Figure 4 demonstrates the simulated results for = 0:5; 0:7, and 0.9 given Cr = 0:5 after 10 seconds of gate opening; it can be seen that the phase error of the surge front slightly decreases
Figure 3. Sketch of surge or dam-break problem.
Figure 4. Simulated results of surge propagation for Cr = 0:5.
as increases. On the other hand, Figure 5, as an illustration, shows the simulated re-sults for Cr = 0:4; 0:5, and 0.6 given = 0:9 at the same time as in Figure 4; the simulated surge celerity is a little slower for the case of Cr = 0:4, and is a little faster for the case of Cr = 0:6.
One of the other tested cases was used to demonstrate the validation of using Cr = 0:5 and = 0:9. Reducing q1 from 20 m2s−1 to 2 m2s−1 and keeping other conditions unchanged as
in the above example, exact values of us and h1 was determined to be 4:23m s−1 and 1:47 m,
respectively. Figure 6 shows the simulated results of the proposed scheme and the Preissmann scheme after 10 seconds of gate opening. The Preissmann scheme used a time-weighting factor
Figure 5. Simulated results of surge propagation for = 0:9.
Figure 6. Comparison of simulated results for downstream propagating surge.
of 0.7 to avoid the numerical oscillations near the surge front. Time increment St = 0:2 s was used in the Preissmann scheme. Even though the Preissmann scheme is slightly more accurate than the proposed scheme in the surge celerity, the simulated result by the proposed scheme is still satisfactory compared with the analytic solution.
Figure 7. Comparison of simulated results for upstream propagating surge.
5.1.2. Sudden closure of a gate. A uniform ow is assumed to have q1= 2 m2s−1 and
h1= 2 m in a channel. The sudden closure of a gate may occur at the downstream end of the
channel due to some unexpected reasons, e.g., an abrupt shutdown of a hydro-power plant connected to the channel. Such an emergent operation of the gate will result in a upstream propagating surge. The exact value us and h2 is 4:23m s−1 and 2:47 m, respectively. The
cal-ibrated values of Cr = 0:5 and = 0:9 for surge problems are used for the proposed scheme. On the other hand, the time increment St and the time-weighting factor for the Preissmann scheme is 0:15 s and 0.7, respectively. The space increment Sx is 1 m for both schemes. Figure 7 shows the simulated results by the proposed scheme and the Preissmann scheme after 10 seconds of gate closure. It can be seen that the proposed scheme performs slightly better than the Preissmann scheme with less phase error of the surge front. In this example, the satisfactory simulated result validates the above calibrated parameters.
5.1.3. Dam-break /ow on wet bed. Computation of dam-break ood waves resulting from the sudden collapse of a dam is necessary for the dam-safety evaluation. The ratio of hd=hu
(Figure 3) is an important index to judge the applicability and accuracy of the numerical schemes to the dam-break problem. As mentioned earlier, both subcritical- and supercritical- ow regimes exist simultaneously in a rectangular, horizontal, and frictionless channel when hd=hu is smaller than 0.138 [32]. In such a mixed ow regime, some of the existing schemes
will encounter the numerical stability problem. In particular, as can be seen later, the simu-lation accuracy will usually decrease signi<cantly as hd=hu approaches zero for most of the
existing numerical schemes.
To compare the simulated results of the proposed scheme with that of other schemes, the example in Jha et al. [11] will be used. In addition to Jha et al.’s ux splitting scheme for comparison, the results obtained by the MacCormack and the Gabutti schemes mentioned in Jha et al. [11] will also be included. In the example, a 2000-m-long channel is assumed to
Figure 8. Simulated results of dam-break ow for Cr = 0:7.
be separated by a dam located in the middle of the channel, as shown in Figure 3. The initial depth (hu) in the reservoir is assumed to be 10 m, and the tailwater depth (hd) is varied
with the ratios of hd=hu= 0:5; 0:05, and 0.005. As was pointed out by Jha et al. [11], the
MacCormack and the Gabutti schemes failed to execute when hd=hu was less than 0.05. If a
proper arti<cial di=usion term was added, both schemes could work, even for smaller ratios of hd=hu. However, the Gabutti scheme would still fail when hd=hu was equal to or less than
0.005. Jha et al.’s ux splitting scheme, on the other hand, has no limitations on hd=hu except
for the critical case of hd= 0.
Before comparing the performances of various numerical schemes, the proper values of Cr and in the proposed scheme for the dam-break ows need to be examined <rst. The case of hd=hu= 0:05 was used for calibrating the values of Cr and . The distance increment (Sx)
used in the simulation was 5 m, and this value was also adopted by the proposed scheme in later analysis. By varying Cr and for test cases, one can conclude that the simulated results under Cr = 0:7 and = 0:7 accurately <t the analytic solution. Hence, these two calibrated values are recommended for use in further analyzes on dam-break problems. Figures 8 and 9 demonstrate that the shape and celerity of the dam-break wave are insensitive to both Cr and values. Figure 8 shows the variation of the solutions due to using di=erent ’s with Cr = 0:7. One can see that the simulated surge height and celerity are more accurate when = 0:7. On the other hand, Figure 9 shows the variations of the solutions due to using di=er-ent Cr’s with = 0:7. The best resolution of the simulated surge height and celerity occurs when Cr = 0:7.
Figures 10–12 illustrate the comparisons among the simulated results obtained by the pro-posed scheme and the other schemes, as well as the analytic solutions, for hd=hu= 0:5; 0:05,
and 0.005, respectively. The superiority of the proposed scheme can be realized from these comparisons. As shown in Figure 10(a), when hd=hu= 0:5, the Gabutti scheme yields the
Figure 9. Simulated results of dam-break ow for = 0:7.
numerical oscillations at the surge front and the slower celerity of the surge front com-pared to the exact solution. Figure 10(b) shows the corresponding velocity distribution along the channel when hd=hu= 0:5. The oscillations of the velocity near the surge front by the
Gabutti scheme are clearly observed. When hd=hu= 0:05, the errors in the location and the
height of the surge front obtained by the MacCormack and the Gabutti schemes, as can be seen in Figure 11, are signi<cant. The results obtained by Jha et al.’s ux splitting scheme are similar to those obtained by the MacCormack scheme when hd=hu= 0:5 and
0.05, and thus are not shown in Figures 10 and 11. Figure 12 shows that Jha et al.’s scheme performs less satisfactory than the MacCormack scheme when hd=hu= 0:005.
Fur-thermore, Figure 13demonstrates the accuracy of the proposed scheme when hd=hu= 0:001.
Instead of using the ux splitting scheme, Jha et al. [9] proposed another scheme based on the modi<cations of the Beam-Warming scheme for the same dam-break example, except that hu was changed to 15 m. The simulated results were improved, as can be seen in their
paper.
5.1.4. Dam-break /ow on dry bed. As mentioned earlier, most of the existing numerical schemes fail when the downstream channel of the dam is initially an dry bed, i.e., hd= 0.
However, the proposed scheme does not have such a limitation. We used the same dam-break example above, but with the initial condition hd= 0, to demonstrate the capability of the
proposed scheme. Figure 14(a) and 14(b) show the simulated water-surface pro<le and the corresponding velocity pro<le along the channel, respectively, after 50 s of dam failure. The results obtained by the proposed scheme are still accurate compared with the analytic solution. In order to compare with the results obtained by Zhang et al. [34] using an ana-log of the SIMPLER method, we need to non-dimensionalize the ordinate and abscissa of Figure 14(b) by dividing them by ghu and t
Figure 10. (a) Comparison of simulated water-surface pro<les for dam-break ow with hd=hu= 0:5; (b) comparison of simulated velocity pro<les for dam-break ow with hd=hu= 0:5.
after the sudden dam break. Figure 15 shows the comparison of the simulated results of the two schemes. Again, the excellent performance of the proposed scheme under the dry-bed condition is observed.
5.2. Dam-break /ow on dry bed in a sloped and frictional channel
The examples demonstrated in the previous sections veri<ed the proposed scheme under the condition of neglecting the e=ects of the bed and friction slope in the Saint Venant equation.
Figure 11. (a) Comparison of simulated water-surface pro<les for dam-break ow with hd=hu= 0:05; (b) comparison of simulated velocity pro<les for dam-break ow with hd=hu= 0:05.
By comparison with the analytic solutions, two important parameters Cr and were examined. With the suggested values of them, it has been found that the proposed scheme simulated most the illustrated cases with excellent accuracy, even the critical ow conditions. However, more important, the capability of simulating ow conditions in a sloped and frictional channel makes the extension of the proposed model to practical applications possible and will be investigated in the section.
Figure 12. (a) Comparison of simulated water-surface pro<les for dam-break ow with hd=hu= 0:005; (b) comparison of simulated velocity pro<les for dam-break ow with hd=hu= 0:005.
The laboratory data under WES test condition 1.1 in 1960–1961 [41] are used to verify the simulated results of the proposed scheme, as well as Cheng-lung Chen’s scheme [43] devel-oped on the basis of an explicit scheme of the characteristics method. WES conducted dam-break ood-wave experiments in a 400-ft (122-m) long, 4-ft (1.22-m) wide, sloped (0.005), and rectangular ume with a 1-ft (0.305-m) high model dam located midway of its length. The model dam was lifted almost instantaneously by a pulley-weight system to simulate the
Figure 13. (a) Simulated water-surface pro<le for dam-break ow with hd=hu= 0:001; (b) simulated velocity pro<le for dam-break ow with hd=hu= 0:001.
sudden dam failure. Under WES test condition 1.1, width of breach is 4-ft (1.22-m), depth of breach is 1-ft (0.305-m), and the channel roughness, Manning coeMcient, is 0.009. The values of parameters Cr and for the proposed scheme maintained 0.7 that was recommended in the previous section.
Figure 16(a) and 16(b) show the water surface pro<les and their corresponding velocity distributions at 0, 10, 20, 30, 50, 90, and 160 s after dam break, respectively. The initial
Figure 14. (a) Simulated water-surface pro<le for dam-break ow with hd= 0; (b) simulated velocity pro<le for dam-break ow with hd= 0.
condition in the channel downstream of the dam for the proposed method can be realized from the <gures that real dry bed condition, i.e. u = 0 and d = 0, was given. This is not the case for Chen’s scheme, in which a computer-simulated condition Qo for no initial ow
was assumed to be 0.001 cfs. At the trailing edge in the receding part of the ood near the upstream end, the same negligibly small in ow discharge Qo was also used to
Figure 15. Comparison of dimensionless velocity pro<les for dam-break ow with hd= 0.
the front of the wave. The assumptions of small initial ows are not needed in the pro-posed scheme and numerical noise does not happen either, as can be seen in Figure 16(a) and 16(b).
According to the classi<cation mentioned previously, Chen’s cheme is a ‘shock <tting’ method. The shock equations (often referred to as the Rankine–Hugoniot equations) are adopted to determine the propogation velocity of the shock front, and the conjugate water depths in back and front of the shock. By imposing shock equations, it has been assumed that the ow downstream of the shock front is always uniform under a given ow rate. Therefore, as imposed on a dry bed, shock equations become valid only under some assump-tions of ow immediately behind the leading edge. For example, Whitham [42] circumvented this problem by assuming uniform ow velocity in the direction of ow behind the leading edge.
However, the proposed scheme belongs to the ‘through’ method which solves the Saint Venant equations directly without any particular treatment of the shock front. Figure 17 shows the water depth hydrographs at six stations; three upstream the dam recording the negative wave at x = 100; 150 and 200 ft; three downstream the dam recording the positive wave at x = 225; 280 and 350 ft. The water depths obtained from the proposed and Chen’s schemes are quite accurate at the upstream stations, and are also close at the downstream stations with slightly small peak depth by the proposed scheme. The arrival time of the positive wave by the proposed scheme at the downstream stations are a little closer to those of the measured data. Figures 18 and 19 show the velocity and discharge hydrographs at the downstream stations, respectively. It can be seen that the peak discharges and corresponding velocities by the proposed scheme are slightly larger than those by Chen’s scheme. Of course, the arrival time of the positive wave are the same as for depth hydrographs and still a little closer to those of the measured data. In general, the two numerical schemes have similar performance and have considerably reasonable results compared with the measured data.
Figure 16. (a) Water-surface pro<les at di=erent time for ood under WES test condition 1.1; (b) velocity distributions at di=erent time for ood under WES test condition 1.1.
6. CONCLUSIONS
The objective of this paper is to propose an iterative explicit scheme suitable for solving 1D unsteady open-channel ow problems based on the Saint Venant equations, especially for the surges and dam-break ows. For each time step, the proposed scheme <rst solves the momentum equation for discharge per unit width and then the continuity equation for water depth, and the procedure is iterated until convergence is reached. A staggered-grid system was adopted in the scheme to avoid the ‘checkerboard’ phenomenon.
Figure 17. (a) Stage hydrographs at station x = 100 ft for ood under WES test condition 1.1; (b) stage hydrographs at station x = 150 ft for ood under WES test condition 1.1; (c) stage hydrographs at station x = 200 ft for ood under WES test condition 1.1; (d) stage hydrographs at station x = 225 ft for ood under WES test condition 1.1; (e) stage hydrographs at station x = 280 ft for ood under WES test condition 1.1; (f) stage hydrographs at station x = 350 ft for ood under WES test condition 1.1.
Figure 17. (Continued).
To determine the parameters Cr and , comparisons were made between the proposed scheme and other schemes against the analytic solution for surge and dam-break ows in the horizontal, rectangular, and frictionless channel. For the surge problem, the suggested values
Figure 18. Discharge hydrographs at stations (a) x = 225 ft, (b) x = 280 ft, (c) x = 350 ft downstream of breach for ood under WES test condition 1.1.
for Cr and in the proposed scheme were 0.5 and 0.9, respectively. For the dam-break problem, the suggested Cr and were both 0.7. Three ratios of hd=hu, i.e., 0.5, 0.05, and
Figure 19. Velocity hydrographs at stations (a) x = 225 ft, (b) x = 280 ft, (c) x = 350 ft downstream of breach for ood under WES test condition 1.1.
height and celerity increased, except by the proposed scheme. In particular, when compared with the results obtained by Zhang et al.’s [34] modi<ed SIMPLER scheme under the dry bed condition, the proposed scheme still has satisfactory accuracy.
Further investigations of the capability of the proposed scheme to more realistic channels with bed slope and frictional e=ects were carried out to demonstrate the extension of this scheme to other potential practical applications are possible. The laboratory data, WES test condition 1.1 in 1960–1961 recording the ow conditions of the dam-break ood waves on a dry bed, were used to verify the numerical schemes. Performance of the proposed scheme is satisfactory compared with the measured data. Because of its algorithmic simplicity and its accuracy for modeling surge and dam-break problems, the proposed scheme may serve as an attractive alternative to other practical, unsteady open-channel problems.
APPENDIX A
The iterative, explicit, characteristics-based <nite-di=erence scheme, Equation (8), for the Saint Venant dynamic equation is stable if [(|u|St)=(Sx)]61:
[Proof] Equation (8) is rewritten in the expression with the discharge q being the unknown variable only. qn j = q 1+ [qn j + (1−)q 1]C1·St + [qn2 j + (1−)q2 1]C2·St + S·St (A1)
where C1 and C2 are regarded as constants, and S is the hydrostatic pressure term and is
neglected hereinafter in the proof. The Fourier transform taking r as the so-called ampli<cation factor is de<ned as
qn
j = rn( )eijS x (A2)
where is the wave number, and r may be complex. Consider u being positive and ! = uSt=Sx61, the following inequality is obtained.
061−4!(1−!) sin2 Sx
2 61 (A3)
Taking the Fourier transforms of Equation (A1) without the term SSt yields r = (1−!) + !e−iS x +{r+ (1−)[(1−!) + !e−iS x ]}C 1·St +{r2eiS x + (1−)[(1−!)eiS x + !]2·e−iS x } C 2·St = (1−!) + !e−iS x +{r+ (1−)[(1−!) + !e−iS x ]} C 1·St
+{r2eiSx + (1−)(1−!)·[(1−!)eiSx + !]
+ (1−)!·[(1−!) + !e−iSx ]} C 2·St
where r is the ampli<cation factor between the discharge newly updated during iterations
relations
|! + (1−!)eiSx |=|(1−!) + !e−iSx |=1−4!(1−!) sin2 Sx
2
1=2
the magnitude of the ampli<cation factor is
|r| 6 1−4!(1−!) sin2 Sx 2 1=2 + |r|+ (1−)1−4!(1−!) sin2 Sx 2 1=2 C1·St + |r|2+ (1−)(1−!)·1−4!(1−!) sin2 Sx 2 1=2 + (1−)!· 1−4!(1−!) sin2 Sx 2 1=2 C2·St 61 + C1·St + C2·St = 1 + O(St)
By the Von Neumann condition, the scheme is stable. The other case of negative u can be proved to be stable likewise.
ACKNOWLEDGEMENTS
This research is supported by the National Science Council of the Republic of China, under grant No. NSC84-2211-E-009-034. The writers are grateful to W. F. Tsai and C. Y. Shen for their valuable advice during the course of this study.
REFERENCES
1. Cunge JA, Holly FM, Verway A. Practical Aspects of Computational River Hydraulics. Pitman: London, 1980. 2. Cunge JA. Rapidly Varied Flow in Power and Pumping Canals in Unsteady Flow in Open Channel. Water
Resources Publications: Fort Collins, CO, 1975: Ch. 14.
3. Yang JC, Chen KN, Lee HY. An accurate computation for rapidly varied ow in an open channel. International Journal for Numerical Methods in Fluids 1992; 14:361–374.
4. Aguirre-Pe J, Quisca S, Plachco FP. Tests and numerical one-dimensional modelling of a high-viscosity uid dam-break wave. Journal of Hydraulics Research 1995; 33(1):17–26.
5. Garcia-Navarro P, Saviron JM. McCormack’s method for the numerical simulation of one-dimensional discontinuous unsteady open channel ow. Journal of Hydraulics Research 1992; 30(1):95–105.
6. Garcia R, Kahawita RA. Numerical solution of the St. Venant equations with the MacCormack <nite-di=erence scheme. International Journal for Numerical Methods in Fluids 1986; 6:259–274.
7. Fennema RJ, Chaudhry MH. Explicit numerical schemes for unsteady free-surface ows with shocks. Water Resources Research 1986; 22(13):1923–1930.
8. Jha AK, Akiyama J, Ura M. A fully conservative Beam and Warming scheme for transient open channel ows. Journal of Hydraulics Research 1996; 34(5):605–621.
9. Jha AK, Akiyama J, Ura M. Modeling unsteady open-channel ows—modi<cation to Beam and Warming scheme. Journal of Hydraulics Engineering, ASCE 1994; 120(4):461–476.
10. Fennema RJ, Chaudhry MH. Simulation of one-dimensional dam-break ows. Journal of Hydraulics Research 1987; 25(1):41–51.
11. Jha AK, Akiyama J, Ura M. An implicit model based on conservative ux splitting technique for one dimensional unsteady ow. Journal of Hydroscience Hydraulics Engineering, JSCE 1994; 11(2):69–82.
12. Alcrudo F, Garcia-Vavarro P, Saviron JM. Flux di=erence splitting 1D open channel ow equations. International Journal for Numerical Methods in Fluids 1992; 14:1009–1018.
13. Glaister P. Approximate Riemann solutions of the shallow water equations. Journal of Hydraulics Research 1988; 26(3):293–306.
14. Glaister P. Flux di=erence splitting for open-channel ows. International Journal for Numerical Methods in Fluids 1993; 16:629–654.
15. Mingham CG, Causon DM. A fully conservative Beam and Warming scheme for transient open channel ows. Journal of Hydraulics Engineering, ASCE 1998; 124(6):605–614.
16. Mingham CG, Causon DM. Calculation of unsteady bore di=raction using a high resolution <nite volume method. Journal of Hydraulics Research 2000; 38(1):49–56.
17. Toro EF. Riemann problems and the WAF method for solving the two-dimensional shallow water equations. Philosophical Transactions of the Royal Society of London, A 1992; 338:43–68.
18. Yang JY, Hsu CA, Chang SH. Computations of free surface ows, part 1: one-dimensional dam-break ow. Journal of Hydraulics Research 1993; 31(1):19–34.
19. Hu K, Mingham CG, Causon, DM. A bore-capturing <nite volume method for open-channel ows. International Journal for Numerical Methods in Fluids 1998; 28:1241–1261.
20. Savic L, Holly FM. Dambreak ood waves computed by modi<ed Godunov method. Journal of Hydraulics Research 1993; 31(2):187–204.
21. Baines MJ, MaMo A, Di Filippo A. Unsteady 1D ows with steep waves in plant channels: the use of Roe’s upwind TVD di=erence scheme. Advances in Water Resources 1992; 15:89–94.
22. Louaked M, Hanich L. TVD schemes for the shallow water equations. Journal of Hydraulics Research 1998; 36(3):363–378.
23. Garcia-Navarro P, Priestley A, Alcrudo F. An implicit method for water ow modelling in channels and pipes. Journal of Hydraulics Research 1994; 32(5):721–742.
24. Garcia-Navarro P, Priestley A. A conservative and shape-preserving semi-Lagrangian method for the solution of the shallow water equations. International Journal for Numerical Methods in Fluids 1994; 18:273–294. 25. Moin SMA, Lam DCL, Smith AA. Eularian–Lagrangian linked algorithm for simulating discontinuous open
channel ows. Proceedings of the VII International Conference on Computer Methods in Water Resources, MIT, USA, 1988; 2:363–368.
26. Hicks FE, SteVer PM. Characteristic dissipative Galerkin scheme for open-channel ow. Journal of Hydraulics Engineering, ASCE 1992; 118(2):337–352.
27. Hicks FE, SteVer PM, Yasmin N. One-dimensional dam-break solutions for variable width channels. Journal of Hydraulics Engineering, ASCE 1997; 123(5):464–468.
28. Katopodes ND. A dissipative Galerkin scheme for open-channel ow. Journal of Hydraulics Engineering, ASCE 1984; 110(4):450–466.
29. Berger RC, Stockstill RL. Finite-element model for high-velocity channels. Journal of Hydraulics Engineering, ASCE 1995; 121(10):710–716.
30. Katopodes ND, Wu CT. Explicit computation of discontinuous channel ow. Journal of Hydraulics Engineering, ASCE 1986; 112(6):456–475.
31. Alam MM, Bhuiyan MA. Collocation <nite-element simulation of dam-break ows. Journal of Hydraulics Engineering, ASCE 1995; 121(2):118–128.
32. Henderson FM. Open Channel Flow. Prentice-Hall: New Jersey, 1966.
33. Fraccarollo L, Toro EF. Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type. Journal of Hydraulics Research 1995; 33(6):843–864.
34. Zhang H, Youssef H, Long ND, Kahawita R. A 1D numerical model applied to dam-break ows on dry beds. Journal of Hydraulics Research 1992; 30(2):211–224.
35. Patankar SV. Numerical Heat Transfer and Fluid Flow. McGraw-Hill Inc.: New York, NY, 1980.
36. Bellos CV, Sakkas FG. 1D dam-break ood wave propagation on dry bed. Journal of Hydraulics Engineering, ASCE 1987; 113(12):1510–1523.
37. Di Monaco A, Molinaro P. Finite element solution of the Lagrangian equations of unsteady free-surface ows on dry river beds. Finite Elements in Water Resources, K.P. Holz et al. (eds). Springer: Berlin, Germany, 1982: 4.25–4.35.
38. Dai W. Numerical Solutions of Unsteady Navier–Stokes Equations Using Explicit Finite Analytic Scheme, PhD Thesis, The University of Iowa, Iowa, 1994.
39. Basco DR. Limitations of de Saint Venant equations in dam-break analysis. Journal of Hydraulics Engineering, ASCE 1989; 115(7):950–965.
40. Kobayashi N, Otta AK, Roy I. Wave re ection and run-up on rough slopes. Journal of Waterway, Port, Coastal and Ocean Engineering, ASCE 1987; 113(3):282–298.
41. Floods resulting from suddenly breached dams. Miscellaneous paper No. 2-374. United State Army Corps of Engineers, Waterways Experiment Station, Vicksburg, Miss., Report 1, Conditions of minimum resistance, 1960; Report 2, Conditions of high resistance, 1961.
42. Whitham GB. The e=ects of hydraulic resistance in dam-break problems. Proceedings of the Royal Society of London, Series A 1995; 227:399–407.
43. Chen CL. Laboratory veri<cation of a dam-break ood model. Journal of Hydraulics Division 1980; 106(HY4):535–556.
