Advances in River Sediment Research – Fukuoka et al. (eds)© 2013 Taylor & Francis Group, London, ISBN 978-1-138-00062-9

Numerical study of flushing half-cone formation due to pressurized
sediment flushing
W.J. Pringle & M.E. Meshkati Shahmirzadi
Urban Management & Environmental Engineering Department, Kyoto University, Kyoto, Japan

N. Yoneyama & T. Sumi
Disaster Prevention Research Institute, Kyoto University, Kyoto, Japan

S. Emamgholizadeh
Department of Water and Soil, Agriculture College, Shahrood University of Technology, Iran

ABSTRACT: Pressure flushing cleans up the hydropower plant entrance from the sediment deposition.
This technique results in localized scouring, that is called flushing cone. In this study, the flushing cone
formation is simulated by implementing a RANS based 3D numerical model that consists of two successive stages: hydraulic flow and sediment erosion simulation. The velocity distribution profiles andthe final geometry of the flushing cone were used to verify the numerical model outcomes. The results
obtained by the hydraulic flow simulation revealed that the calculated velocity distribution profiles werecomparable with the experimental data and the universal wall equations with roughness, ks  10 mm as asediment boundary condition was adequate for the application. However, when considering the simulation of the sediment erosion process using the moving boundary method the model should be improvedfor better accuracy in predicting flushing cone geometry.

1

INTRODUCTION

One of most pressing problems of river and dam
management in current times is the excess deposition of sediment. Severe sediment deposition
which occurs in dam reservoirs decreases the water
storage capacity and thus reduces the efficiencyof the dam operation related to flood control andflow regulation (Graf, 1984). On top of the issues
directly related to sedimentation in dam reservoirs,
including intake clogging and machinery abrasion,a number of environmental problems also arise
downstream of the dam. These include riverbed
and bank erosion, armouring of the river bed, anda reduction in the quality of habitat for aquatic
species, thereby threatening reductions in river system biodiversity (Kondolof, 1997).
To negate the issues related to sedimentation,
a number of techniques have been proposed suchas density current venting, flushing, by passing,dredging and sluicing (Brandt, 2000). However,flushing has proven to be one of the most effective and efficient in removing reservoir sedimentthrough the bottom outlet (Shen and Lai, 1996).
This technique can be categorized into two types;
free flushing and pressure flushing. In free flushing
the water level is drawn down to a relatively low

level which causes large scale scouring in the formof a channel along the reservoir. Whereas, in pressure flushing the water level is maintained at a nearconstant level high above the bottom outlet suchthat the loss of water, an important dam resource,
is conserved. Here, the scour formation shape issimilar to that of a half-cone, often simply referred
to as a flushing cone. This paper is concerned withpressure flushing and the formation of the cone.
Previously a few studies have experimentallyinvestigated the flushing cone formation and proposed empirical equations to estimate the final
volume and length of flushing cone (Scheuerlein
et al., 2004, Emamgholizadeh, 2005, Meshkati
et al., 2010, and Fathi-Moghadam et al., 2010).
While empirical regression-based relations canallow for a quick and simple estimate of final cone
geometry, these equations are often flawed becauseof the high nonlinearity between the cone geometry and the significant parameters.
Instead of using experimental and analyticalmethods, a less costly and potentially more accurate alternative is to use a numerical model (Olsen,
1999). A number of numerical studies have been
conducted to simulate the flushing operation in
dam reservoirs. Most of these studies are devoted
to either 1D or 2D numerical models. Holly and1249

Rahuel (1990); Lai and Shen (1996) and Basson
and Rooseboom (1997), individually developed a
1-D model to simulate the bed topography change
during the drawdown phase in free flushing. Thesemodels may only be used to study flushing processes in long and narrow reservoirs. Ruland andRouve (1992) utilized a 2D finite element model
to simulate the erosion process in the reservoir.
Olsen (1999) solved the depth-averaged NavierStokes equation on a 2D mesh and extrapolated
the corresponding flow field onto a 3D mesh tosolve sediment concentration in flushing processes.
Although, in addition to its difficult analyticalcharacterization, the formation of the flushingcone is a complex three-dimensional phenomenon(Scheuerlein et al., 2004), where turbulent eddies,secondary currents and vertical velocities are significant, particularly in the vicinity of the bottom
outlet. However, as can be found in the literature
only a few studies have considered a 3D numerical model and the majority of these are devoted
to only free flushing phenomena, e.g.: Haun and
Olsen (2012a, b) where a RANS model that solves
the transient-convection equation for suspended
sediment transport and Van Rijn’s (1984) empiricalformula for bed load transport was utilized.
In this study, the formation of the flushing cone
is simulated utilizing a RANS based 3D numericalmodel verified by Yoneyama (2010) and Yoneyama
et al. (2012a, b) and adapted for this analysis. The
simulation process consists of two stages. In thefirst stage the velocity distribution and equilibrium bed topography (final shape) of the flushing half-cone, from the experimental investigation
conducted by Emamgholizadeh (2005), were usedfor verification and calibration of the hydraulicmodel. In the second stage the numerical model
commences from the initial sediment condition
incorporating the moving boundary method tosimulate the sediment erosion process. The moving boundary method, while primitive, was chosen because of its simplicity and since only the
final bed shape is of interest in this analysis. The
results of the simulation in the second phase arecompared to the geometry of the flushing cone asfound experimentally.

2

EXPERIMENTATION

2.1

Experimental setup

A hexahedral shallow basin of size 3.0 m length,
1.5 m width and 1.5 m height was used as thereservoir model. A schematic side and plan viewof the model as well as the scoured flushing conein the vicinity of bottom outlet are shown in
Figure 1. A gate valve was installed in the center

Figure 1. Schematic side and plan view of the reservoir model, the flushing cone and its relevant hydraulicparameters.

of the front wall of the model (dam wall) as arepresentation of the bottom outlet. A 50.8 mm
diameter gate valve was employed. Inside thereservoir, sediment was placed and paved in
the reach covering the entire width and 2.0 mof the length from the dam wall. The sediment
thickness was 0.32 m, measured from the lower
invert of the outlet to the sediment surface. Sand
particles with a uniform size distribution wereused as sediment. The particles had a median
diameter, d50 of 1.2 mm and specific gravity,
Gs of 2.65. Inflow discharges of 1.0, 3.0, 4.5, 6.0
and 8.0 L/s were used at the upstream end. The
discharge at the outlet was kept the same as thatof the inflow so that the reservoir water level wasmaintained constant at the required elevation.
Three constant water levels of 0.425, 0.825 and1.125 m were examined.
A digital point gauge device was utilized to measure the scour cone configuration after each experiment. The vertical velocity profile was measured
along six different axis locations by Laser-Doppler(LDV). The axis locations were selected along thex-axis (centerline of the bottom outlet), separated
by equal 0.1 m intervals and numbered I–VI in
Figure 1.
2.2

Experimental procedure and design

The following general procedure was enacted
for all tests conducted in the experimental study.
Before each test the sediment was wetted, flattened

1250

and leveled to a specific level above the lower invertof the bottom outlet, Hs. To begin the test run, themodel was slowly filled with water by adjusting
the incoming flow rate into the reservoir. When thewater surface elevation in the reservoir above thelower invert of the outlet, Hw reached the desired
level, the bottom outlet was carefully opened until
the outflow discharge became equal to that of the
inflow discharge. Subsequently, the sediment flushing began, flowing out from the reservoir into thesettling basin. Experiments were continued until
equilibrium conditions prevailed in which the sediment concentration released from the reservoir was
negligible, signaling the end of the experiment. The
development of the flushing cone was fairly quickwith the bulk of the formation finishing within a
few minutes or less.
At this stage, the vertical velocity profiles at thesix axes were measured using LDV. Afterwards, atthe end of each test case as soon as the pump was
turned off and the incoming discharge was set tozero, the outlet was closed immediately. Thereafter,
water was carefully and slowly drained from thereservoir via a pump. The final scour contour wasmeasured using digital point gauges. The volumeof the flushing cone was calculated by Surfer (8.0)
Software.

3

NUMERICAL METHOD

3.1

Basic equations

The numerical method consists of the RANS equations combined with the FAVOR method (Hirt &Sicilian, 1985) for modelling the interface between
fluid and obstacles, and the VOF model (Hirt &Nichols, 1985) for modelling the free surface. Then,the continuity, momentum and advection of fluid
equations 1–3 in this form are written as;Continuity:
tc v
tt

tc aj u jtx j

c aj u j tui
c

v

(1)

tx j

 gi

1 tp
l txi

t ¤ tui¥v
x j ¦ tx j
ttx³
u` i u` j ´µ(2)

Advection of fluid:
tc v Fttu` i u` j 

tc aj Fu
F j
tx j

0

(3)

C * k 2 ¤ tui¥a ¦ tx j

tu j ³´txi µ

2k` i j
3

ta
tt

(4)

a

1 tc j ku j
c v tx j2t §¤ C * k ³ tk ¶ ` ` tui¨¥ v· ui u j´
tx j
tx j ¨ ¦m k a µ tx j ·¸
©

tk
tt

a

(5)

2
a
t § ¤ C * k ³ ta ¶1 tc j a u j
¨
·

v
¥
´
tx j ¨ ¦tx
c v tx j
m a a µ tx j ·¸©

Ca 1a ` ` tuiui u j
tx j
k

Ca 2

a2
k

(6)

where k x u` i u` j / 2 is the turbulent kinetic energy,
a x vu` i j u` i j is the turbulent energy dissipation,and `i,j is the Kronecker delta. The constants inequations 3–5 are mk  1.0, ma  1.3, C*  0.09,
Ca1  1.45, and Ca2  1.92. The discretized form ofequations 1 and 2 are solved using the SIMPLE
method (Patankar et al., 1972) on a staggered grid
system. Finally, while in this study the free surface is not such a significant parameter, the VOF
method is used for its calculation.
3.2

0

Momentum (i  1, 2, 3):tui
tt

where c v is the void ratio of a cell, c ai is the aperture ratio of a cell boundary, ui is the flow velocity
component, gi is the external force per unit volume,
p is the pressure, l is the density of the fluid, i isthe dynamic viscosity, F is the fluid volume fraction of a cell, ¯ is the Reynolds averaging quantityand ` is the fluctuation in the Reynolds average.
The k-a turbulence closure model is utilized toform equations 4–6 to model the Reynolds stressterm, u` i u` j :Moving boundary method

As shown in equations 1 and 2, the FAVOR
method is used to model the interface between thefluid and obstacles; in this case sediment. Here thevoid ratio, cv and aperture ratios, cai are calculated
to indicate the proportion and orientation of fluidin a cell. Thus, we can view the interface between
the sediment and fluid (in one plane) as consistingof two nodes and a straight boundary line as illustrated in Figure 2, with the corresponding void andaperture ratios indicated.
The magnitude of the tangential velocity justabove the boundary node, Utang is compared withthe critical velocity for sediment transport, Ucr.
When Utang is larger than Ucr, the boundary node is
moved down at a specified velocity, Vd, as illustrated

1251

Figure 2. Calculation cell containing both fluid andsediment, with its corresponding void ratio, cv and aperture ratios, cax.

The length is 3.0 m, width 1.5 m and maximum domain height of 1.0 m. The cell size in thex-direction ranges from 0.05 m at the inlet end, to0.02 m in the area close to the outlet. The size of thecell in the y-direction is 0.025 m at the sides and as
small as 0.01 m at the outlet. Similarly the cell size inthe vertical direction is 0.01 m near the outlet and up
to 0.02 m elsewhere. A total of around 390,000 cells
are used to make up the domain but approximately
220,000 cells are used in calculation. Figure 5 provides a snapshot of the calculation domain with anillustration of the parameters Q, Hw and Hs.
4.2

Bottom outlet representation

The bottom outlet was represented as a polygon
made up of the rectangular cells of the numericalmesh as illustrated in Figure 6. The maximum crosssectional width of the numerical outlet was 50 mm(compared with a diameter of 50.8 mm for theexperiments) and had a total cross-sectional areaof 2.0 r 10 3 m2 (compared with 2.03 r 10 3 m2 forthe experiments). The bottom outlet was modeledby manipulating the streamwise aperture ratio,c ax into the last cell in the calculation domain as
indicated in Figure 6.
Figure 3. Movement of the boundary node when thecritical velocity for sediment transport is exceeded near
the node.

Figure 4. Movement of the boundary node towards theangle of repose in submerged conditions, ]s.

Figure 5. Side-view snapshot of the velocity profile inexperiment Case D (1st stage) with parameters Q, Hw andHs shown.

in Figure 3. Furthermore, the angle of the boundary line cannot exceed the angle of repose in submerged conditions, ]s. In this case the higher of thetwo boundary nodes are lowered at the velocity, Vd,
until the angle of repose is reached (Fig. 4).
4
4.1

NUMERICAL CONDITIONS
DomainThe domain size and conditions are the same as theexperiments conducted by Emamgholizadeh (2005).

Figure 6. Diagram of the bottom outlet representationin the YZ plane, and its aperture ratios, cax.

1252

4.3

5

Boundary conditions

Both the upstream inlet boundary and the downstream bottom outlet boundary were subjected to
Dirichlet boundary conditions.
At the inlet, both the constant water elevation,
Hw, and the resulting ux velocity component calculated from the input flow-rate, Q, over the crosssectional area were specified. Since a single cell
contains two cell boundaries in the x-direction, thevelocity, ux was specified for both cell boundaries.
The velocity distribution was uniform.
The ux velocity component at the bottom outlet
was specified by the setting the outflow equal to theinlet flow rate, Q. This flow rate was divided by thecross-sectional area of the bottom outlet, found insection 4.2, to calculate and set the same ux velocity
component at each cell in the bottom outlet. In addition, since the inflow takes some time to reach andaffect the area downstream, a “warming-up period”
of 2.5s was specified to avoid imparting an initial
rapid outflow velocity with no inflow into the control volume. Within this period, the outflow velocity
was linearly increased from zero at 0s, to the calculated velocity based on the inflow, Q at 2.5s.
4.4

Numerical cases

5.1

Rigid bed analysis

In this stage of the study various sediment boundary conditions were tried to identify the best representation of the experimental work. The simplest
boundary conditions to apply are the no-slip, halfslip or free-slip conditions. The other sedimentboundary conditions applied were the universal
wall equations. These were applied in empiricalform and as suggested by Le Roux (2004) which
gives equations for all three hydrodynamic stages(rough, smooth and transitional). These equations
were adjusted by changing the roughness of therigid bed, ks. The sediment boundary conditions
are shown in Table 2.
The selection of ks  1.0 mm is because it isvery close to the mean sediment diameter whileks  10 mm represents a much rougher condition as
a comparison. Figures 7–10 show velocity profiles
at axis 1 and axis 2 for Cases B and D.
Typically, the velocities at low heights above thesediment surface were either; overestimated near
the bottom outlet (axis 1) except for the free-slip
case or underestimated further from the bottom
outlet (axis 2) in the calculation. For larger heights
above the sediment surface better agreement was

Table 1 shows the experimental conditions thatwere used for simulation.
Each case was made up of two stages. The firststage consisted of solely a hydraulic study withthe sediment, represented as a rigid bed, at the
final stage as measured from the experiments. The
resulting velocity profile was compared with that of
the experiments. The simulation was calibrated byadjusting the boundary condition at the bed. The
second stage took the calibrated hydraulic modeland incorporated the moving boundary method.
The initial condition was that of the sedimentplaced flat at a uniform depth throughout theentire domain. The simulation was run until themodel predicted no further change in scour shape,and compared with the measured experimental
scour shape (the flushing cone).
Table 1.
study.

NUMERICAL RESULTS

Table 2.
applied.

Sediment boundary conditions

ConditionTypeks (mm)

I
II
III
IV
V

Free-slip
No-slip
Half-slip
Wall equation
Wall equation

–
–
–
1.0
10.0

Experimental conditions applied to numerical

Case

Q
(L/s)

Hw
(m)

Hs
(m)

D
(cm)

d50
(mm)

A
B
C
D
E
F

6.0
8.0
6.0
8.0
6.0
8.0

0.425
0.425
0.825
0.825
1.125
1.125

0.32
0.32
0.32
0.32
0.32
0.32

5.08
5.08
5.08
5.08
5.08
5.08

1.2
1.2
1.2
1.2
1.2
1.2

Figure 7.
Case B.

1253

Velocity profile comparison at axis 1 forFigure 8.
Case D.

Velocity profile comparison at axis 1 for

Figure 9.
Case B.

Velocity profile comparison at axis 2 forcondition that sets the pressure of the outflow
equal to atmospheric the outflow will depend on
the hydrostatic pressure in the reservoir rather than
the numerically stipulated outflow velocity calculated using Q. Here, as the outflow increases sothat it exceeds Q, the opening of the bottom outlet,
as defined by the aperture ratios, can be decreasedso that the outflow is equivalent to Q.
It was found that by comparing RMSE and R2
curve fit values between the experimental data andthe numerical data for all axes, considering onlyboundary conditions I–III, the no slip sedimentboundary condition (II) produced more favourable
results overall. When comparing with the wall equations (IV and V), only very tiny differences were evident between each other and the no-slip condition.
This is because the flow is generally hydrodynamically smooth everywhere. Here, the roughness, ks
does not influence on the calculation of the velocity
at the boundary. In fact only at the peak velocitiesof flow near the bottom outlet where the flow could
be in a hydrodynamically transitional or roughstage are any notable differences detected between
the wall equations and the no slip condition. But as
noted above, the peak velocities at axis 1 are overestimated and can be reduced slightly by imposing
the wall equations and a higher value of ks. Moreinvestigation is required to determine the optimum
value of ks for this study but ks  10 mm produced
slightly improved results from ks  1 mm.
5.2 Movable bed analysis
In this stage the calculation is conducted with an
initially uniform sediment bed level (0.32 m for allcases). The wall equations are applied as the sediment boundary condition determined from stage 1.
The calculation is run until no notable change wasevident between a sufficient time gap. Similar to
the experiment this generally occurs after about
60–120s depending on the calculation case.
Figures 11 and 12 show the final bed topography
for Cases B and D respectively. As can be seen a cone

Figure 10.
Case D.

Velocity profile comparison at axis 2 forfound between the experimental and calculation
data, however slightly underestimated in general.
To obtain better results the boundary condition at
the bottom outlet may need to be modelled differently. For example, by using a pressure boundary

Figure 11.

1254

Final flushing cone geometry for Case B.

Figure 12.

Final flushing cone geometry for Case D.

shape is predicted however the simulation could not
yet lead to robust results in terms of flushing volume
and length. It is possible that an analysis considering bed load transport would allow for better results;
however it is assumed that in pressure flushing mostof the sediment is quickly removed through the bottom outlet without significant redeposition. The
simulations are ongoing and by improving the bottom outlet and the sediment boundary conditions,better results using the moving boundary method
are foreseeable in the near future.
6

CONCLUSION

A numerical study to simulate the pressure flushingin dam reservoirs was conducted. An experimental
study by Emamgholizadeh (2005) was used to verify
the numerical analysis. Velocity profiles were shownto have moderate agreement with the experiments.
By defining the bottom outlet velocity as a functionof hydrostatic pressure, the simulation results may
improve. Furthermore, it was found that the wall
equations where ks  10 mm as the sediment boundary condition lead to better agreement in general.
The simulation by the moving boundary method
to calculate the sediment removal is still undergoing and the calculation of its geometrical properties
could not yet be obtained accurately. The authors
are optimistic that work in the near future willattain amenable results by making improvementsin the bottom outlet boundary condition, cell sizes
and sediment boundary condition.
REFERENCES
Basson, G.R. & Rooseboom, A. 1997. Dealing withreservoir sedimentation, guideline and case studies.
ICOLD report. Bulletin 115. Paris.
Brandt, S.A. 2000. A review of reservoir de-siltation.
International Journal of Sediment Research. 15(2):321–342.

Emamgholizadeh, S. 2005. Pressure flushing of sediment through storage reservoir: Laboratory testing.
Journal of the Institution of Eng. (India), CE Div.
89(5): 23–27.
Fathi-Moghadam, M. Emamgholizadeh, S. Bina, M. &Ghomeshi, M. 2010. Physical modeling of pressure
flushing for de-silting of non-cohesive sediment.
Journal of Hydraulic Research. 48(4): 509–514.
Graf, W.H. 1984. Storage losses in reservoirs. International Water Power & Dams Construction. 36(4):37–40.
Haun, S. & Olsen, N.R.B. 2012a. Three-dimensional
numerical modelling of reservoir flushing in a prototype scale. International Journal of River Basin Management. 10(4): 341–349.
Haun, S. & Olsen, N.R.B. 2012b. Three-dimensional
numerical modelling of the flushing process of the
Kali Gandaki hydropower reservoir. Lakes and Reservoirs: Research and Management. 17(1): 25–33.
Hirt, C.W. & Sicilian, J.M. 1985. A Porosity Technique
for the Definition of Obstacles in Rectangular Cell
Meshes. Proc. 4th Int. Conf. Ship Hydro. National
Academy of Science. Washington, D.C. 1–19.
Holly, F.M. & Rahuel, J.L. 1990. New numerical/physical framework for mobile-bed modeling. Journal of
Hydraulic Research. 27(4): 401–416.
Kondolf, G.M. 1997. Hungry water: effects of dams
and gravel mining on river channels. Environmental
Management. 21: 533–551.
Lai, J.S. & Shen, H.W. 1996. Flushing sediment throughreservoirs. Journal of Hydraulic Research. 34(2):237–255.
Le Roux, J.P. 2004. An integrated law of the wall forhydrodynamically transitional flow over plane beds.
Sedimentary Geology. 163(3–4): 311–321.
Meshkati Shahmirzadi, M.E. Dehghani, A.A. Sumi, T.
Naser, Gh. & Ahadpour, A. 2010. Experimental investigation of local half-cone scouring against dam. Proc.
Int. Conf. on Fluvial Hydraulics River. Braunschweig.
Germany. 1267–1273.
Olsen, N.R.B. 1999. Two-dimensional numerical modeling of flushing processes in water reservoirs. Journalof Hydraulic Research. 37(1): 3–13.
Patanker, S.V. & Spalding, D.B. 1972. A Calculation Procedure for Heat, Mass and Momentum Transfer inThree-Dimensional Parabolic Flow. Journal of Heat
Mass Transfer. 15: 1787–1806.
Ruland, P. & Rouve, G. 1992. Risk assessment of sediment erosion in river basins. Proc. of the Fifth International Symposium on River Sedimentation. Karlsruhe.
Germany. 357–364.
Scheuerlein, H. Tritthart, M. & Nunez Gonzalez, F. 2004.
Numerical and physical modeling concerning theremoval of sediment deposits from reservoirs. Proc.
of Hydraulic of Dams and River Structures. Tehran.
Iran. 245–254.
Shen, H.W. & Lai, J.S. 1996. Sustain reservoir useful life
by flushing sediment. International Journal of sediment Research, IRTCES (International Research andTraining Center on Erosion and Sedimentation). 11(3):
10–17.
Van Rijn L.C. 1984. Sediment transport. Part I: bed loadtransport. Journal of Hydraulic Engineering. 110(10):1431–1456.

1255

Yoneyama, N. 2010. Optimal design of a sand discharge
pipe determined by fluid analysis using the moving
boundary method. Proc. of the 11th International
Symposium on River Sedimentation. Stellenbosch,South Africa. CD-ROM.
Yoneyama, N. Saito, A. & Iguchi, K. 2012. Developmentof numerical analysis code for predicting behavior ofremoval sediment from reservoir bed by using discharge pipe based on siphon principle. Proc. of International Symposium on Dams for A Changing World.
24th ICOLD congress. Kyoto, Japan. CD-ROM.

Yoneyama, N. Nagashima, H. & Toda, K. 2012. Threedimensional numerical analysis to predict behaviorof driftage carried by tsunami. Earth Planets Space.
64: 965–972.

1256