Zhang Xiuzhong and Wang Guangqian
Key laboratory for Water and Sediment Sciences of Ministry of Education, Tsinghua University, Beijing, 100084, China
Abstract:
Based on the 2D shallow water equations and the nonequilibrium suspended
sediment transport equation, a mathematical model is presented to predict flow
motion and scour hole evolution of an instantaneous dyke- or dam-breach. To
simulate the complex flow phenomenon
accurately and conveniently, a high-resolution finite
element scheme is applied to capture the dam-break bore. Flow analysis and scour
computation are conducted for a variety of breach widths, discharges and soil
composition. The results obtained are used to provide hydrodynamic conditions
for structural and stability analysis of the steel, timber and earth-rock
combined dam.
Keywords:
dike breach, scour hole, high-resolution finite element, numerical simulation
Levee or dam failures during floods often cause tremendous loss of lives and properties. For example, a 164-m long breach of the Putuo River in Hebei province of China resulted in flooding of 23 villages and 18 thousand hectares of farmland by the released 230 million cubic meter water in August 1996 [1]. Another failure of the anti-flood wall of the Jiujiang City by the Yangtze River in 1998 not only aggravated flood disaster but drove the flood control work into a severer situation. Whether the breach can be closed up quickly or not, and how high cost will be paid all depend on the hydraulic characteristics and the scour development of the breach. Therefore, further researches on problems of the dyke-breach play a vital role in preventing, reducing, and harnessing the flood disaster.
There are two major causes for a dike failure [2]. One is piping which causes collapse of levee and gives rise to overflow. The other is that flood stage exceeds the top of the levee to lead to overtopping. Both of the two kinds of flows are open channel flow, quite similar to the dam-break flow. Thus the dike-break flow can be described by the shallow water equations (also named the Saint-Venant equations for the 1D case) mathematically. After the dike fails with breaches, the bed at or near the breach will be continuously eroded to form scour hole. As the scour hole becomes bigger and deeper, it will further endanger the safety of the levee, and might cause a wider range of dyke breaching. So exactly predicting the flood routing and evolution of the scour hole downstream the breach will undoubtedly direct the implementation of flood prevention.
With the development of computer powers and the advance in numerical techniques, numerical simulation of dike-break problems is not only feasible but efficient. However, due to the convection-dominated and degenerative hyperbolic features of the shallow water equations, especially for the complex structure of the dyke-breach flow, with subcritical, supercritical flows, and even discontinuities either presenting simultaneously in different parts or occurring in sequence at different time. Thus how to overcome these numerical difficulties is the key to simulating the dyke-break flow.
So far, calculating results of the eroding process near the breach are rarely available in the literature. This paper is devoted to the study of flow and scour of the dyke breach. For that, a high-resolution FEM is adopted to track the dam-break bore [3]. Assuming that the bed material is noncohesive, and its motion satisfies the sediment continuity equations, so that the nonequilibrium sediment transport model can be applied to predict the motion of the sediment and the evolution of the scour hole.
It is commonly accepted to consider the
dam-break flow as a shallow water phenomenon, and its solution with jump or
shock discontinuities is said to constitute a weak solution of the conservation
law. Under
the hydrostatic approximation, the simplified N-S equations to describe
incompressible free surface flow are as follows
(1)
(2)
The nonequilibrium sediment transportation and bed deformation are described by
(3)
(4)
where
;
and
= depth averaged velocity and unit
discharge in
direction, respectively; C = the
Chezy resistance coefficient; Eddy viscosity is estimated by
, in which
= Von Karman’s constant,
= friction velocity;
;
;
;
,
and
= depth integrated concentration,
sediment-carrying capability and setting velocity of the kth
grain size fraction;
= recovery coefficient;
= dry bulk density of the bed material;
= sediment diffusive factor assumed
to be the same as
in this paper; water surface
elevation
is defined by water depth
and bed surface level
.
In view of sediment either from upstream flow-carrying, or turbulence diffusion entering or the both, it is reasonable to calculate the gradation of sediment carrying capacity by considering both the flow condition and size distribution of the bed material simultaneously. So the treatment approach proposed by Li [4] is adopted.
Given the erosion or deposition depth
of the kth fraction sediment
, and the total erosion or deposition depth
, the computation of the bed material gradation can be divided into two cases:
(1)
, which means deposition occurring, the gradation of the mixed layer can be
estimated by
where
,
= gradations of the mixed layer at the time level
and
, respectively;
and
are the corresponding thickness of
the mixed layer.
(2)
, which means that erosion happened, the gradation of the mixed layer is then
computed by
where
= the average gradation in some sub-layers of the memory layer.
To simulate the armoring phenomenon of
the bed material in degradation, two layers are incorporated in the model. The
upper layer of the bed material is called mixed layer and the lower is named
memory layer. The thickness of the mixed layer is
, and its gradation is
. The memory layer is divided again into many sub-layers, say n sub-layers, and
the thickness and the gradation of each sub-layer are
,
, respectively. The number of the sub-layers of the memory layer will increase
to n+1, and the gradation of the newly added sub-layer equals to that of the
mixed layer in time level t, if the riverbed aggrades. Otherwise, if the
riverbed degrades, the number of the sub-layers will subtract some layers, and
the gradation of the memory layer should be adjusted correspondingly.
The FEM is adopted due to its advantage to treat
practical problem with complex geometries.
Discretization for spatial variables and integration by part leads to
(5)
(6)
,
where M = lumped
mass matrix; C = convective matrix; D = diffusive matrix; F = source term. They
are expressed as follows.
;
;
;
Where N, W are interpolation and weight function, respectively.
When W=N,
it constitutes the Galerkin FEM which equals to central difference, so numerical oscillation will be created when solving convection-dominated problems.
In this study, a high-resolution finite element scheme is incorporated to
rearrange the convective term. Its main characteristic is that the scheme is
obtained through the introduction of almost equal amount of diffusion and
anti-diffusion, with limiter to preserve monotunity condition, so that both the
requirement of high-order of accuracy and non-oscillatory condition are
satisfied.
For the sake of convergence or fasting convergence, the source term F is transposed and linearized. After discretization for spatial variables, Eqs.(5), (6) are a system of ordinary differential equations, which can be solved by various explicit and implicit schemes. The Runge-Kutta time integration is used herein for unsteady flow. A sequential procedure is employed to solve the discretized equations, where Eq.(5) is a 2D convective equation for computing h; Eq.(6) is a 2D convection-diffusion transport equation to calculate u, v and s; Erosion depth can be obtained by Eq. (4).
To verify accuracy, and allow comparisons with analytical solutions, the scheme is used to simulate a dam-break flow. The problem considered involves sudden and complete failure in a horizontal rectangular channel without friction. Initially the water is rest with an upstream depth of 10 m and a downstream of 1 m. Comparisons with stoker’s solution [5] at 120 second are displayed in Fig 1 and Fig 2.
Calculations
concerning both dam-break and dike-breach are carried out. For the dam-break
case, attention is primarily paid to considering flow and scour development near
the breach under conditions that the inlet water surface elevation remain
unchanged and all flow pass the breach. As for the case of dike-burst, focus is
mainly put to observe flow and scour evolution under conditions that the channel
diverts the flood and the upstream water stage is changeable.
The physical region
is shown in Fig 3. The horizontal channel is 2500m in length, 1000m in width.
The dam is assumed to be 30m thick in the direction of flow. The Manning
roughness is set to be 0.015, and the median size of the bed material is 0.3 mm.
Initially the water is rest with an upstream depth of 7.5m and a downstream of
0.5m. As to the boundary conditions, water stage is specified to be 7.5m at the
upstream, non-slip conditions are satisfied at and upstream the breach; at the
downstream Neumann condition is applied to the variable u.
Three cases with the breach width of 50m, 100m, and 200m are simulated, respectively. Confined to the length of the paper, only some representative result figures are incorporated. More figures are available in the report [6]. With the breach width of 50m, Fig. 5 shows the varying process of the scour depth at the central line of the breach. At the same wide breach, Fig. 6 and 7 give the scour depth contours at different time. Fig. 8 and 9 present the velocity vector when the breach is 100m wide.
To demonstrate effects of the breach width on scour depth, and to investigate the relation between the sour depth and the breach width, the maximum scour depth and its location are listed in Table 1.
Table 1 Erosion results of different widths of the breach
|
|
maximum erosion depth (m) |
location of maximum erosion (m) |
||||
|
width (m) |
50 |
100 |
200 |
50 |
100 |
200 |
|
100s |
0.091 |
0.16 |
0.14 |
96.0 |
120.0 |
165.0 |
|
600s |
0.633 |
1.15 |
1.30 |
96.8 |
119.0 |
165.0 |
|
1200s |
1.181 |
2.01 |
2.53 |
78.0 |
98.0 |
142.0 |
|
1800s |
1.640 |
2.67 |
3.50 |
60.0 |
80.0 |
120.0 |
|
3600s |
2.75 |
4.1 |
5.60 |
45.0 |
80.0 |
99.0 |
|
7200s |
4.375 |
6.1 |
|
45.0 |
62.0 |
|
|
Notes: D50=0.3 mm, numbers in the first column is time. |
||||||
It can be seen from Tab.1 that as the breach widens, the scour depth increases, and the distance between location of maximum erosion and the beach increases as well. The location of the maximum erosion depth gradually approaches the breach with time going on, till at a stable position at last.
From the scour depth contour we can see that the plane view of the scour hole is elliptical, and this trend is more obvious as time becomes longer. From Tab.1 and figure of scour hole elution we can also draw the following conclusions. The scour rate increased firstly, then dropped, and got zero at last. Compared to its downstream slope, the scour hole has a steeper upstream slope.
As Fig.4 shows, the computational domain consists of the upper stream reach, the flood diversion area below, and the breach zone. The stream is 800m wide with a bed slope and a levee width being 0.0002 and 50m, respectively. Both the breach and the flood diversion area are horizontal. The Manning factor of the three portions is set to be 0.002. As boundary conditions of the stream, discharge process is given at the inlet, and the relation between water stage and discharge, which can be derived by the Chezy formula, is specified at the outlet. The initial conditions are determined by the computation of the channel flow.
At the same discharge and bed material, four cases with the breach width of 100m, 200m, 300m, and 500m are computed. With the same breach width, calculations of different discharges and bed material composition are also included. Similarly, only typical result figures are attached to conserve space. More figures can refer to [6]. With the breach width of 300m, Figs.10, 11 give the velocity vectors, Fig.12, 13 describe the evolution of the scour depth at the central line of the breach. The velocity process at the central line of a 300-meter wide breach is demonstrated in Fig. 14.
To compare the different cases and to explore some useful laws, the maximum scour depths and its locations at time of 7200s are listed in Tab. 2 and 3.
Table 2 Maximum erosion depth of different widths of the breach
|
Time(s) |
7200 |
|||
|
Breach width (m) |
100 |
200 |
300 |
500 |
|
Maximum erosion depth (m) |
2.01 |
2.14 |
2.0 |
0.95 |
|
Location of maximum erosion depth (m) |
50. |
85. |
86. |
86. |
|
Notes: D50=0.25 mm, Q=10000.0 m3/s. |
||||
Table 3 Effects of different discharge and sediment on scour depth
|
D50(mm) |
0.25 |
0.1 |
|
|
discharge(m3/s) |
10000. |
15000. |
15000. |
|
Maximum erosion depth (m) |
2.0 |
2.85 |
3.1 |
|
Location of maximum erosion depth (m) |
86. |
90. |
120. |
|
Notes: the breach width is 300 m, and the time is 7200 s. |
|||
From Tabs. 2, 3 as well as the scour depth evolution figure, It is not
difficult to see that the maximum scour depth increases firstly, then decreases
when the breach widens. At about 200m of the breach, the scour depth is deepest.
At the same breach width, the larger the inlet discharge of the channel is, the
greater will be the scour depth, and vice versa. Moreover for the same inlet
discharge and breach width, the scour hole depth will be reduced if coarse bed
material is used. In addition, the scour length is proportional to the scour
depth. Fig.14 indicates that the velocity at the central line of the breach
increases initially and then dropped.
The numerical scheme is validated by comparison
with Stoker’s analytic solutions. As the present model is based on
unstructured grid, it is very flexible to deal with geometrically complex problems that
will occur in practical engineering.
The plane view of the scour hole is elliptic. Contrast to its downstream slope, the scour hole has a steeper upstream slope. There is a linear relation between the scour distance and the scour depth, which manifests that the geometrical similarity exists in the evolution process of the scour hole. The scour rate increases firstly, then drops, and gets zero at last.
For the dam-break scour with upstream water stage remaining unchanged, the wider the breach is, the deeper the scour hole gets. For dyke-break scour with the stream inflow discharge being constant, the scour depth will have a maximum value when the breach has a certain width. In addition, the evolution characteristics of the scour hole is closely related to the stream discharge and the grain size of the bed material. The velocity is of magnitude near the breach, and it almost approached 10. m/s in the computation. It is a fact that the velocity increases in the beginning and then decreases.
Owing to the complexity of problems of the dike burst, some simplifications or assumptions are made in the model, so improvements and further researches are essentially required.
References
[1] The 27th Army Group headquarters of Chinese People's Liberation Army. 1998, Report on application research of the breach-stopping technique for steel, timber and earth-rock combined dam. Rep. (in Chinese).
[2] Gui, S. X., Zhang, R. D., and Xue, X. Z. 1998, Overtopping reliability models for river levee. J. Hydr. Engrg., ASCE, 124(12), 1227-1234.
[3] Zhang, X. Z., Wang, G. Q., and Jin, S. 2001, Finite-element analysis for shallow water flows and high-resolution FEM scheme. J. Yangtze River Scientific Res. Inst., No.1. (in Chinese).
[4] Li, Y. T. 1987, Preliminary study on the gradation of bed load in an equilibrium state. J. Sediment Res., No. 1. (in Chinese).
[5] Stoker, J. J. 1957, Water waves. Interscience Publications, Wiley, New York.
[6]
Tsinghua University, the 27th Army Group. 2000, Hydraulic characteristics and
scour hole computation of the dike-break. Rep., Beijing (in Chinese).

Fig. 1 Comparison of water surface elevation Fig. 2 Comparison of velocity
Fig. 3 Definition sketch for partial dam break Fig. 4 Definition sketch for Dyke breach
Fig.
5 Depth evolution of scour
hole
Fig. 6 Contour of scour hole at t=600
s
Fig. 7 Contour of scour hole at t=7200 s Fig. 8 Velocity vectors at t=600 s
Fig.
9 Velocity vectors at t=3600 s
Fig. 10 Velocity vectors at t=100 s

Fig.
11 Velocity vectors at t=7200 s
Fig. 12 Depth evolution
of scour hole

Fig.
13 Depth evolution of scour
hole
Fig. 14 Evolution of the velocity at the
central line of the breach