717
Technical Committee 103 /
Comité technique 103
l
i
ij
i
ij
vh
vh
vh
) (
) ( ) (
r
(12)
ij
l
i
ij
i
ij
ij
ij
z
z
h
) (
r
(13)
where
z
ij
is the height of the flow bed at the midpoint of the
j
th
side of the
i
th cell, which is obtained as the average of the bed
height at the two vertices making the side, and
r
ij
is the position
vector from the cell centre to the midpoint of the side. As seen
in equations (11) to (13), momentums (
uh
)
ij
and (
vh
)
ij
on the cell
face are directly evaluated with their unlimited gradients, while
flow depth
h
ij
is calculated by subtracting
z
ij
from
ij
after the
height of water surface
ij
is evaluated on the cell face with its
unlimited gradient.
3.3 Treatment of source term
The source terms in equation (1) are divided into the slope and
the friction terms seen as
S
0
and
S
f
in equation (2). The
treatment of these terms has a great influence on the accuracy of
the numerical scheme. When the surface gradient method is
employed, adequate attention needs to be paid to balance the
first numerical flux term and the second source term on the
right-hand side of equation (7) under a steady state. For this
purpose, the second and third components of slope term
S
0
are
rewritten in the following form:
2
2
)
(
z g z g z z
g z gh
(14)
Then, these components of the
i
th cell are numerically evaluated
as
i
i
i
i
z g z
g z gh
2
2
(15)
in which
3
1
1
j
ij ij
ij
i
i
z
A
z
t
,
3
1
2
2
1
j
ij ij
ij
i
i
z
A
z
t
(16)
where
t
ij
denote the normal unit outward vector at the
j
th side of
the
i
th cell. The slope term of
S
0
, calculated from equations (15)
and (16), can be balanced with the numerical flux term, and the
right-hand side of equation (7) vanishes under the steady
condition of no-flow velocity and a constant water level.
In treating the friction source terms, a simple explicit
method may induce numerical instabilities when the water
depth is very small. To overcome this problem, the friction
terms are treated in a fully implicit way with the operator-
splitting technique proposed by Yoon & Kang (2004).
i f
i
dt
d
,
S U
(17)
i
j
ij
ij
i
i
A dt
d
,0
3
1
*
1
S
E
U
(18)
The right-hand side of equation (17) includes only friction
source terms. Equations (17) and (18) are solved in implicit and
explicit ways, respectively.
3.4 Alteration of flow bed elevation
The values of
u
,
v
and
h
are stored at the centroids of the cells,
while the values of
z
are placed at the vertices. The height of
flow bed
z
changes in accordance with the erosion rate of the
bed material, which is a function of the bed shear stress related
to the flow velocity and the flow depth. Therefore, the values of
u
,
v
and
h
at the vertices must be known and the linearity-
preserving interpolation method proposed by Holmes & Connel
(1989) is used to calculate their vertex values. From these
values at the vertices, the temporal changes in the flow bed are
computed by the following equation:
k
k
k
E
dt
dz
1
(19)
which is the ordinary differential equation with respect to
t
at
the
k
th vertex based on equation (4). The TVD Runge-Kutta
scheme was applied for the time integration of the spatially
discretized equations (18) and (19), which simultaneously
solved the shallow water equations and the temporal changes in
the bed height.
4 NUMERICAL SIMULATION
While the governing equations of equations (1) and (4) are
solved over the two-dimensional computational domain, i.e.,
x-y
plane, the elevation changes of water surfaces and erosion beds
can be succesively computed, which implies that the numerical
methods can produce quasi-three-dimensional results over two-
dimensional computational domain. This is a great advantage
for reducing computational load of three-dimensional numerical
analyses which usually need enormous computational effort and
time. An example of three-dimensional numerical simulation of
embankment breaching induced by the concentration of
overtopping water flow is presented herein.
Figure 2 shows the initial profiles of embankment and water
surface and the imposed boundary conditions. The embankment
has the dimensions of 30cm in height, 130cm in bottom length
and 60cm in thickness. The centre of the crest is 2cm lower than
the other part to induce the concentration of the overflow. 2,050
finite volume cells with 1096 nodes were used for the spatial
discretization. As the boundary conditions for the water flow,
the inflow per unit width at a rate of 0.029 m
3
/s/m was given
from the extreme upstream (the right extreme in Figure 4), and
the free outfall condition and the free slip condition were
imposed onto the downstream end and the sides of the
calculation domain, respectively. The steady flow velocity and
flow depth under these boundary conditions were adopted as the
initial conditions for the water flow on the embankment. The
flow velocity vector field is shown in the
x-y
plane of Figure 4.
For the embankment material properties, a porosity of 0.395, a
critical shear stress of
c
=0.1Pa and erodibility constants of
=8.42×10
-5
m/s/Pa
3/2
and
=1.5 were given. The value of the
Manning’s roughness coefficient was assumed to be 0.0158.
Figure 3 shows the computed embankment profiles and
water surfaces 100 and 600 seconds after the initiation of the
embank-ment erosion. As seen in the figure, the central part was
dominantly eroded, the overflowing water concentrated there
and a flow channel passing through the embankment appeared.
The grey cells on
x-y
plane indicate the dry surface which
appeared because of the flow concentration. The numerical
results shown in Figure 3 reflect the stability and the feasibility
of the proposed method for three-dimensional analysis.
5 CONCLUSIONS
This paper has presented the numerical simulation of
embankment erosion caused by overflow. As the governing
equations, the two-dimensional shallow water equations were
adopted as the governing equations for describing the water
flow onto embankments, and the temporal changes in the flow
bed were formulated by the erosion rates. The finite volume
method was employed for the spatial discretization of the flow
domain, and the HLL Riemann solver was used to evaluate the
flux through the cell interfaces. The surface gradient method
(SGM) was incorporated into the finite volume approach; this
enabled the stable computation of the flow field on the erosion
surface having complex undulation. The three-dimensional
analysis has shown the natural and stable results of embankment
breaching which includes the concentration of water flow and
embankment erosion, and the creation of the breach channel.
