Strojniški vestnik - Journal of Mechanical Engineering 59(2013)5, 301-310 © 2013 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2012.669	Original Scientific Paper
Received for review: 2012-06-18 Received revised form: 2012-11-23 Accepted for publication: 2013-01-17
Evaluation of Fluid Flow over Stepped Spillways Using the Finite Volume Method as a Novel Approach
Hamid Reza Vosoughifar1* - Azam Dolatshah2 - Seyed Kazem Sadat Shokouhi1 - Seyed Reza Hashemi Nezhad1 1 Islamic Azad University, South Tehran Branch, Department of Civil Engineering, Iran 2 Islamic Azad University, Dezful Branch, Department of Civil Engineering, Iran
This research deals with the development of a Computational Fluid Dynamics (CFD) code called V-Flow using MATLAB for the two-dimensional modelling of unsteady flow over stepped spillways. V-Flow can be coupled with GAMBIT software and different spillway geometries can be modelled using voronoi mesh elements. The governing equations of flow over stepped spillways were discretized using the Finite Volume Method (FVM). The Power-Law scheme, implicit time approximation, Gauss-Seidel method, and SIMPLE algorithm were used for the discretization procedure. The flow was considered to be a laminar fluid flow with no turbulence model. The V-Flow model was validated against velocity vectors, streamlines, static pressures, dynamic pressures, and total pressures over the spillway obtained from the FLUENT model application. The experimental model simulated using both V-Flow and FLUENT for validation according to Gonzalez's experiments. A comparison between the results obtained from the application of the different models showed good agreement according to the mass imbalance, which serves as both a useful indicator of the convergence of the fluid flow solution and as an accuracy criterion to compare the V-Flow and FLUENT results.
Keywords: stepped spillways, finite volume method, discretization, voronoi mesh, numerical modelling, roller compacted concrete
0 INTRODUCTION
A stepped spillway is a hydraulic structure, which is an integrated part of the dam, that allows the safe passage of overtopping flows. Due to increased air entrainment (i.e. aeration flow), stepped spillways dissipate the energy of the flow. This means that they reduce the level of energy in dissipater structures downstream from the spillway and lower the risks of cavitation. In recent years, the development of the Roller Compacted Concrete (RCC) technique has improved stepped spillways due to the low-cost and relatively high-speed construction this method enables.
Most studies investigating the flow over stepped spillways have been undertaken utilizing both laboratory experiments on scaled models and analytical and numerical approaches. For example, Degoutte et al. [1] suggested that the energy dissipation is higher in the jet flow than the skimming flow. They also discovered that there are two types of jet flow: one with fully developed hydraulic jumps and the other one with partially developed hydraulics jumps. Gonzalez and Chanson conducted an experimental study to gain a better understanding of the flow properties in stepped chutes with slopes typical of embankment dams. Their work yielded a new design procedure including some key issues not foreseen in prior studies [2]. Pfister and Hager presented visual observations made with a high-speed camera and air concentration measurements in the vicinity of
the pseudo-bottom air inception point on a stepped model spillway [3]. Felder and Chanson conducted a physical study in a moderate slope-stepped chute and tested five configurations. The results were compared in terms of flow patterns, energy dissipation, and flow resistance [4]. Chanson and Felder also studied the two-phase gas-liquid flow properties of highvelocity open channel flows in a large-size channel experimentally. The results demonstrated high levels of turbulence in the high-speed, highly turbulent free-surface flows [5]. Barani et al. [6] used the feasible direction method and a wooden physical model to determine the optimum slope and step height of a stepped spillway and to calculate the dissipated energy of the flow. Si-ying et al. [7] studied the effect of aerator type on countering cavitation damage in the stepped chute at the Murum Hydraulic Power Station. They optimized the shape of the aerator and observed the aeration effects by hydraulic modelling.
With the advent of high-performance computers and the development of robust Computational Fluid Dynamic (CFD) software, the development of complementary analytical tools for resolving the intricacies of the flow pattern have become feasible, for examples see Chen et al. [8], Tabbara et al. [9], and Naderi Rad and Teimouri[10]. Carvalho and Martins [11] also investigated hydraulic jumps on the steps of stepped spillways analytically, physically, and numerically. They designed a conceptual prototype using classic hydraulic formulae. A large-scale model was adapted and an experimental study was conducted
*Corr. Author's Address: Islamic Azad University, South Tehran Branch, Tehran, Iran, vosoughifar@azad.ac.ir
301
to examine the similarity of hydraulic jumps at each step, minimizing the hydraulic jump length, and maximizing the discharge per unit width. A numerical model based on the Two-Dimensional (2D) Reynolds-averaged Navier-Stokes (RANS) equations was used for evaluating velocities and pressures and for characterizing the hydraulic forces on the baffles and sills. Bombardelli et al. [12] presented and discussed the results of a comprehensive study addressing the non-aerated region of the skimming flow in steep stepped spillways. They used a relatively large physical model of the spillway to acquire data on flow velocities and water levels. Numerical simulations using a commercial code were then performed to reproduce those experimental conditions.
In addition, Hanbay et al. [13] created two intelligent models to predict flow conditions and aeration efficiency in stepped cascades using critical flow depth, step height, and channel slope. Salmasi [14] developed a procedure using an Artificial Neural Network (ANN) for calculating the energy dissipation of flow over a stepped spillway chute. The work was based on a wide variety of experimental facilities in large-size models.
In this paper, the advantages of the CFD software are considered and a program written to investigate the flow pattern problems. In general, the Finite Volume Method (FVM) and voronoi mesh grid are used to solve the governing equations and to estimate the flow characteristics over stepped spillways, especially the velocity vectors, streamlines, static pressures, dynamic pressures, and total pressures over steps.
1 RESEARCH METHODOLOGY 1.1 Governing Equations
The 2D governing equations of a Newtonian, viscous, incompressible, unsteady state and turbulent flow can be described via the Continuity and Navier-Stokes equations, Eqs. (1) and (2), respectively [15], whereas the last term of the Navier-Stokes equation is eliminated when describing the laminar flow state.
V .F = 0,	(1)
p^ + pV-Vy = -VP + pgi +	-Uv'. (2)
dt
where V , p, ft t, P, g., p, u' and v' are the velocity vector, water density, velocity vector components in the x or y direction (^ = u, v), time, pressure, gravity acceleration in the x or y direction (i = x, y),
dynamic viscosity and fluctuating values of velocity vector components in the x ory direction (ft = u', v') respectively.
An additional term called Reynolds stresses (-pV • u'v') in the right side of Eq. (2) arises from the time-averaged equations for turbulent flow compared to laminar flow. This is due to an averaging operation in which it is assumed that there are random fluctuations about the mean value. It means that =	and
P = P + P', while 0 and P are the mean values of the velocity vector components in the x or y direction (0 = u, v ) and the pressure, respectively. In addition, ft and P' are fluctuating values of the velocity vector components and pressure around their mean values, respectively. However, as mentioned before, if the last term is ignored, the Navier-Stokes equation represents the laminar flow state.
Many turbulence models such as the k-e and Renormalization Group (RNG) models employ the concept of a turbulent viscosity to express the turbulent stresses. The result is that the time-averaged equations for turbulent flow have the same appearance at the equations for laminar flow, except that the laminar exchange coefficient such as viscosity is replaced by the effective exchange coefficient.
However, only laminar flow was considered in this research, i.e. no turbulent model was taken into account. Instead, attempts have been made to link the pressure and velocity fields somehow and to correct them during the solution procedure for the discrete governing equations using the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm [16]. The main features of this process are explained in the following section.
1.2 Numerical Modelling Algorithm
In Finite Volume Method (FVM), the integration domain is covered with control volumes. Each control volume surrounds a node, which lies on a grid mesh. Mesh grids are generally divided into structured and unstructured meshes. Due to its generality, the FVM can handle any type of mesh, structured as well as unstructured. Unstructured grid methods employ an arbitrary collection of elements to fill the domain. For instance, the voronoi and Delaunay are types of unstructured meshes. Dirichlet [17] used a special form of the voronoi tessellation in his study of positive quadratic forms. Voronoi later published a generalization of this concept that would apply to higher dimensions and so introduced the concept in its modern form. A Delaunay triangulation is created
by several algorithms such as the Watson-Bowyer and Step-By-Step Algorithms [18].
In this paper, the studied domain was discretized using the FVM for solving the discrete governing equations and the unstructured voronoi mesh grid established using the MATLAB software [19]. In this approach, the domain was divided into several distinct control volumes without any overlapping. By integrating the governing differential equations over every control volume, an algebraic equations system was created for every single control volume, and each equation links the parameter investigated in a control volume node to different numbers of adjacent nodes, consequently the parameter was computed in each node [20]. In order to solve discrete equations the parameter investigated in each node was computed considering its discrete equation and newest adjacent nodes' values. The solution procedure can be expressed as follow:
•	Initial guess for the investigated parameter (i.e. in each node at initial conditions.
•	Calculation of the value in an investigated node considering its discrete equation.
•	Performing previous step for all nodes in the studied domain; one cycle is performed by repetition of this step.
•	Verification of the convergence clause. If this clause is satisfied, the computing will end, otherwise the computation will be repeated from the second step.
The Riemann boundary condition was also utilized for computing values in boundaries, because the exact values of the boundary conditions were not distinct. Therefore, we can assume that a layer which is near the boundary layer, where d / dx = 0 and d / dy = 0, is defined and that the calculated value for the boundary adjacent nodes is transformed into the related boundary nodes. This procedure will continue until the difference in results converges to zero.
1.3 Discretization of Governing Equations
Discretization of the Continuity and Navier-Stokes equations was performed utilizing the FVM. Fig. 1 illustrates the 2D voronoi mesh grid used for describing these equations. The discrete Continuity and the unsteady state Navier-Stocks equations were discretized as shown in Eqs. (3) and (4), respectively.
where;
anb = DfA(\D-\) + [[0, -Ff ]], Df
a =Yj(anb + F )f
f
Df =
HÄ
Si
xnb - x) +( - y0]
Ff =P A)f.
(6)
(7)
(8)
The general discrete governing equation of the unsteady state flow problem was also discretized using implicit time approximation as shown in Eq. (9).
m = ZiWnb) ~Y{pA)f + PSA9 + prr^, (9)
where;
anb = DfA
D,
+ [[0, -Ff ]].	(10)
and are defined in Eqs. (7) and (8), respectively. The two general discrete equations for evaluating the velocity components in two directions (u, v) can be determined from Eq. (9), where the velocity components (u, v) need to be inserted into the equation instead of Furthermore, other changes due to different coordinate directions (x, y) should be also considered.
)f = 0,
(3)
Fig. 1. The 2D schematic voronoi mesh cell used to describe the discretization of the governing equations
) -YiP~A )f +PSi A9> (4)
1.4 Pressure and Velocity Correction Equations and Mass Imbalance
Once the two general discrete velocity equations have been derived, we can attempt to solve the different terms and parameters of these equations. The two obtained equations, which evaluate the velocity components in two directions, should satisfy the continuity equation. Thus, the continuity equation will be obtained as shown in Eq. (11). Fig. 2 shows the velocity components of the 2D schematic voronoi cell used to describe the discretization of the governing equations.
(AX) P2àyf
K )f =
pAf
Ff =pA) =(P0 -p) +
(11)
Fig. 2. Velocity components of the 2D schematic voronoi cell for
describing the discretization of the governing equations
With regard to the SIMPLE algorithm, let u* and v* be the discrete u and v fields resulting from a solution of the discrete u and v momentum equations. In addition, let P* represent the discrete pressure field that is used in the solution of the momentum equations. If the pressure fieldp* is only a guess or a prevailing iterate, the discrete u* and v* obtained by solving the momentum equations will not, in general, satisfy the discrete continuity equation. So therefore we should let p = p* + p'. In addition, to know how the velocity components respond to this change in pressure, we should let u = u* + u' and v = v* + u'. These corrected values satisfy Eq. (11). Finally, Eq. (12) can be obtained for the evaluation of the pressure correction.
a0=Yj{a'nb )f ,
f ,0
b = -! f
f ,0
(13)
(14)
(15)
p0 =Z(£
f ,0
nbpnb
)f + ^
(12)
The mass source, b, serves as a useful indicator of the convergence of the fluid flow solution. If the face flow rates, F*, satisfy the discrete continuity equation, b will be zero. Indeed the b term is a measure of the resulting mass imbalance. The pressure correction equation corrects the pressure and velocity fields to ensure that the resulting field annihilates this mass imbalance [21] and [16]. In this research, the mass imbalance was considered to be a factor indicating the accuracy of the results. The results were finally studied and compared with each other according to the mass imbalance factor at the end of the modelling practice.
1.5 Numerical Model
In this study, a finite volume program called V-Flow was developed using the MATLAB programming for 2D modelling of unsteady flows over stepped spillways. V-Flow can discretize the domain using the unstructured voronoi mesh grid with different geometries, boundary conditions, and arbitrary vertices. V-Flow can be connected to the GAMBIT software [22] to model the arbitrary and complicated geometry. The study domain was modelled in GAMBIT using triangular elements. Then, the coordinates of the created nodes (not the created triangular mesh grid) were exported from GAMBIT and imported into V-Flow. V-Flow created the Delaunay triangulation based on the coordination of the grid. Then, it created the voronoi mesh grid according to the relationship between the Delaunay triangulation and voronoi mesh generation. Generally, the Power-Law [16] and fully implicit schemes were used for discretization of the governing algebraic equations.
Although it is somewhat more complicated than other schemes, the Power-Law scheme expressions are not particularly expensive to compute and they provide an extremely good representation of the exponential behaviour. However, when the general discrete equation is related as Eq. (4) and anb is
a
a
described as Eq. (5), then A(\Ff / Df |) can be written as Eq. (16) in terms of the Power-Law scheme.
Ff	Ff 5
A(| D I) = [[0,(1 - 0.1\D~ I)5]]. (16)
In addition, the iterative solution Gauss-Seidel method was employed to solve the algebraic equations set resulting from the discretization. The SIMPLE algorithm was also used for the pressure correction equation. The flow was considered to be a laminar fluid flow and no turbulent model was applied. Consequently, V-Flow can evaluate the velocity and pressure field in both directions and it can display the flow pattern over stepped spillways.
On the other hand, the FLUENT CFD software [23] was implemented as a reference in order to evaluate the precision of V-Flow. FLUENT also uses the FVM to solve the governing equations. It provides physical models on different types of unstructured 2D/3D meshes. Its modelling and grid generation is conducted using GAMBIT. Here the SIMPLE algorithm, Power-Law scheme, and the segregated solver in the laminar flow state were used to investigate the verification model compared with the V-Flow model. The results obtained from both V-Flow and FLUENT were then compared with each other, especially where the mass imbalance was concerned.
1.6 Case Study
To test the precision of the V-Flow code developed in this study, a physical model of a large stepped spillway with flat horizontal steps was considered for simulation of the flow properties using the V-Flow code and FLUENT software. This broad-crested stepped spillway physical model was created by Gonzalez [24] at the hydraulics laboratory of University of Queensland, Australia. The model is 3.15 m long, 0.9 m high and 1 m wide. It has 9 steps and each step is 0.1 m high and 0.35 m long. The slope angle of the chute is 15.94°, the length of the crest is 0.617 m from the upstream face of the spillway to the tip of the first step, and the total length of the model is 4 m including upstream and downstream in addition to the spillway.
In this paper, the broad-crested profile has an angle steepness of 15.94° which is equal to the chute slope angle, because the V-Flow was prepared to model the non-free surface flow. Therefore, it was essential to change the crest of the model.
1.7	Boundary Conditions
Boundary conditions comprise four boundaries. The first boundary condition is the area beneath the steps made from soil. This is a non-slip wall boundary condition. The second one is the inlet flow boundary just at the left face of the spillway. This boundary is defined as the velocity inlet boundary condition. The total velocity, which is parallel to the slopping crest, is considered to be 2 m/s and the water depth is assumed to be 0.208 m. The third boundary condition is the air pressure boundary condition, which is a line parallel to the chute at top of the flow region over the spillway. The pressure was defined as Pa. The fourth one is the outflow boundary condition just at the right face of the spillway and both the velocity components and pressure are undefined.
1.8	Convergence Criterion
In statistics, the Mean Absolute Error (MAE) is a quantity used to measure how close predictions are to the eventual outcomes. It is a common measure of forecast error in time series analysis. The Mean Absolute Error is given by Eq. (17).
MAE =1 ¿|f -e,|,	(17)
n ,=i
wheref is the prediction and e, is the true value.
In this paper, MAE has been considered equal to 10-6 in order to lead to convergence for both V-Flow and FLUENT.
2 RESULTS
2.1 Prepared Code
A large stepped spillway physical case study was modelled using GAMBIT software. The coordination nodes of the triangular mesh grid were imported into the V-Flow code. A Delaunay triangulation was produced based on these nodes and a voronoi mesh grid was consequently created utilizing the Delaunay triangulation. Fig.3 illustrates the Delaunay triangulation and voronoi mesh generation in V-Flow.
The model was run using the V-Flow code and Figs. 4 to 8 show the results of the analysis procedure, which include velocity vectors, streamlines, static pressure, dynamic pressure and total pressure contours over the spillway.
Fig. 3. The Delaunay triangulation and voronoi mesh generation created by V-Flow
Fig. 4. Velocity vectors over the stepped spillway obtained using V-Flow
Fig. 5. Streamlines and static pressures over the stepped spillway obtained using V-Flow
Fig. 6. Static pressure contours over the stepped spillway obtained using V-Flow
Fig. 7. Dynamic pressure contours over the stepped spillway obtained using V-Flow
Fig. 8. Total pressure contours over the stepped spillway obtained using V-Flow
Fig. 9. Velocity vectors over the stepped spillway obtained using FLUENT
Fig. 10. Streamlines over the stepped spillway obtained using FLUENT
2.2 CFD Software for Validation
2.3 Comparison
The large stepped spillway physical case study was simulated using FLUENT software for validation of the V-Flow code. The GAMBIT software was used for triangular mesh generation which is similar to the Delaunay triangulation for V-Flow. Figs. 9 to 13 show velocity vectors, streamlines, static pressure, dynamic pressure and total pressure contours over the spillway.
Nodes, mesh grids, velocity vectors, streamlines, static pressures, dynamic pressures, and total pressures over the spillway are shown in Figs.3 to 13 for both V-Flow and FLUENT. Total pressure contours indicate a decline in pressure of 0.1043 to 0.1008 MPa from the fluid surface to the tip and vertical face of steps. During the solution procedure, as the obtained
I
1.03e+0s [Pa]
I
-^T g ^			
Y (cm) I— □-1> X(cm)	_ w c	1	
Fig. 11. Static pressure contours over the stepped spillway obtained using FLUENT
Fig. 12. Dynamic pressure contours over the stepped spillway obtained using FLUENT
Fig. 13. Total pressure contours over the stepped spillway obtained using FLUENT
responses become closer to final responses, the mass imbalance at each control volume will be at its lowest. So, the lowest mass imbalance leads to achieving the best result. Table 1 shows the mass imbalance values for both V-Flow and FLUENT. As it can be seen, the relative mass imbalance on average for an individual cell is 0.0000213 or 0.0000329 for V-Flow or FLUENT, respectively. The difference between these two codes is about 0.0000129, and V-Flow provides a mean mass imbalance less than that of FLUENT. A small difference is observed between the maximum and minimum mass imbalance values of V-Flow compared to FLUENT. The sum total of the mass imbalance at all control volumes can be used as a criterion to determine whether the continuity equation is satisfied. This relative sum total is 0.04331 and 0.0626 in the laminar flow state for V-Flow and FLUENT, respectively. Finally, Table 2 shows a comparison between some other properties of V-Flow and FLUENT, such as mesh type, number of elements, and iterations.
Table 1. Comparison between mass imbalance values
Mass Imbalance (b)			
Minimum	Maximum	Mean	Sum Total for all cells
V-Flow -0.00997	0.00571	0.0000213	0.04331
FLUENT -0.02157	0.00961	0.0000324	0.0626
Table 2. Comparison between some basic properties of V-Flow and			
FLUENT			
		V-Flow	FLUENT
Mesh grid type		Voronoi	Triangular
Number of elements [-]		2032	3692
Maximum time [s]		11	14
Time step[s]		0.0001	0.000178
Iterations after several time steps [-]		600	3000
Mean Absolute Error (MAE)		10-6	10-6
3 CONCLUSIONS
The FVM with the Power-Law scheme and implicit time approximation can discretize the governing equations of the 2D unsteady laminar incompressible flow problem over stepped spillways to prepare a CFD numerical code named V-Flow using MATLAB programming. V-Flow can be coupled with the GAMBIT software and used to model different geometries with voronoi mesh elements.
According to the results, velocity vectors and streamlines illustrate re-circulating vortexes at the corner of steps in both V-Flow and FLUENT, as
they were observed in experiments performed by Gonzalez [24]. This is an important property as a result of establishing the skimming flow regime over stepped spillways, which can be simulated using V-Flow precisely. On the other hand, V-Flow was prepared using a no turbulence model, whereas the recirculating vortices clearly appeared at the corner of steps. An important property of V-Flow is to reduce the procedure solution. Furthermore, the results demonstrate that the static pressure significantly increased where the streamlines are in contact with the horizontal face of steps. Static pressure contours show an increase in pressure from the fluid surface to the bottom. Dynamic pressure contours show that the pressure is higher in regions where the velocity is high. Increasing the total pressure at the tip and vertical face of steps raises the risk of cavitation in this area compared to others.
Despite the fact that the number of elements used by the V-Flow code is less than that used by FLUENT, the V-Flow results are very precise. This may be due to the efficiency of the voronoi mesh grid in modelling the geometry and to the discretized equations. The node that is being studied in the voronoi mesh grid is at the centre of the cell and it is considered to be the central node of the computational control volume. This is not the case for the triangular mesh grid used by FLUENT due to the placement of the nodes at the tip of the triangular control volumes. It is essential to transfer the nodes that are on boundaries from the tip of triangles to the middle of the sides using some approximations. Therefore, voronoi elements can model the geometries better than triangular elements.
The mass imbalance criterion can be used to investigate the precision of the results and the convergence during the iterative solution. Both V-Flow and FLUENT give precise results due to a very small mass imbalance and similar output results due to a very small mean mass imbalance difference. A small difference was observed between the maximum and minimum mass imbalance values of V-Flow compared with FLUENT. This means that V-Flow can make a more uniform mass balance than FLUENT. The sum total of the mass imbalance of V-Flow is less than that of FLUENT. This may indicate that the velocity field obtained using V-Flow satisfies the continuity equation approximately three times better than the velocity field obtained using FLUENT.
4 REFERENCES
[1] Degoutte, G., Peyras, L., Royet, P. (1992). Disscussion
of skimming flow in stepped spillways by Rajaratnam.
ASCE Journal of Hydraulic Engineering, vol. 118, no. 1, p. 111-114., D0I:10.1061/(ASCE)0733-9429(1992)118:1(111.2).
[2]	Gonzalez, C.A., Chanson, H. (2007). Hydraulic design of stepped spillways and downstream energy dissipaters for embankment dams. Dam Engineering, vol. 17, no. 4, p. 223-244.
[3]	Pfister, M., Hager, W.H. (2011). Self-entrainment of air on stepped spillways. International Journal of Multiphase Flow, vol. 37, no. 2, p. 99-107, DOL10.1016/j.ijmultiphaseflow2010.10.007.
[4]	Felder, S., Chanson, H. (2011). Energy dissipation down a stepped spillway with non-uniform step heights. Journal of Hydraulic Engineering, vol. 137, no. 11, p. 1543-1548, D0I:10.1061/(ASCE)HY.1943-7900.0000455.
[5]	Chanson, H., Felder, S. (2010). Turbulence measurements in air-water self-aerated flows: basic analysis and results. 7th International Conference on Multiphase Flow, Tampa.
[6]	Barani, G.A., Rahnama, M.B., Bagheri, H. (2005). Optimization of stepped spillway dimensions and investigation of flow energy dissipation over a physical model. Journal of Applied Sciences, vol. 5, no. 5, p. 878-882, DOI:10.3923/jas.2005.878.882.
[7]	Si-ying, W., Dong-mei, H., Cai-huan, W. (2012). Aerator of Stepped Chute in Murum Hydropower Station. Procedia Engineering, vol. 28, p. 803-807, DOI:10.1016/j.proeng.2012.01.813.
[8]	Chen, Q., Dai, G., Liu, H. (2002). Volume of fluid model for turbulence numerical simulation of stepped spillway overflow. ASCE Journal of Hydraulic Engineering, ,vol. 128, no. 7, p. 683-688, D0I:10.1061/ (ASCE)0733-9429(2002)128:7(683).
[9]	Tabbara, M., Chatila, J., Awwad, R. (2005). Computational simulation of flow over stepped spillways. Computer & Structures, vol. 83, no. 27, p. 2215-2224, DOI:10.1016/j.compstruc.2005.04.005.
[10]	Naderi Rad, I., Teimouri, M. (2010). An investigation of flow energy dissipation in simple stepped spillways by numerical model. European Journal of Scientific Research, vol. 45, no. 4, p. 544-553.
[11]	Carvalho, R., Martins, R. (2009). Stepped spillway with hydraulic jumps: application of a numerical model to a scale model of a conceptual prototype. ASCE Journal
of Hydraulic Engineering, vol. 135, no. 7, p. 615-619, D01:10.1061/(ASCE)HY. 1943-7900.0000042.
[12]	Bombardelli, F.A., Meireles, I., Matos, J. (2011). Laboratory measurements and multi-block numerical simulations of the mean flow and turbulence in the non-aerated skimming flow region of steep stepped spillways. Environmental Fluid Mechanics, vol. 11, no. 3, p. 263-288, D0I:10.1007/s10652-010-9188-6.
[13]	Hanbay, D., Baylar, A., Batan, M. (2009). Prediction of aeration efficiency on stepped cascades by using least square support vector machines. Express System whit Application, vol. 36, no. 3, p. 4248-4252, D0I:10.1016/j.eswa.2008.03.003.
[14]	Salmasi, F. (2010). An artificial neural network (ANN) for hydraulics of flow on stepped chutes. European Journal of Scientific Research, vol. 145, no. 3, p. 450457.
[15]	Shaughnessy, E.J., Katz, I.M., Schaffer, J.P. (2005). Introduction to Fluid Mechanics. Oxford University Press, Oxford.
[16]	Patankar, S.V. (1980). Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York.
[17]	Dirichlet, G.L. (1850). Über die Reduktion der positiven quadratischen Formen mit drei unbestimmten ganzen Zahlen. Journal für die Reine und Angewandte Mathematik, vol. 40, p. 209-227, D0I:10.1515/ crll.1850.40.209. (in German)
[18]	Delaunay, B. (1934). On empty sphere; A memory of Georges Voronoi. Bulletin of the Academy of Sciences of the U.S.S.R., vol. 7, no. 6, p. 793-800.
[19]	MATLAB. (2011). MATLAB R2011a (7.12), User's Guide, Mathworks, San Francisco.
[20]	Fletcher, C.A.J. (1991). Computational Techniques for Fluid Dynamics, Volume 2, 2nd ed., Springer Verlag, New York.
[21]	Murthy, J.Y. (2002). Numerical Methods in Heat, Mass, and Momentum Transfer. School of Mechanical Engineering, Purdue University, Purdue.
[22]	Gambit. (2007). Version 2.4.6 User's Guide, Fluent, Lebanon, New Hampshire.
[23]	Fluent. (2009). Fluent User's Guide. Ansys Inc., Lebanon.
[24]	Gonzalez, C.A. (2005). An Experimental Study of Free-Surface Aeration on Embankment Stepped Chutes. Ph.D. thesis, Department of Civil Engineering, The University of Queensland, Brisbane.