Flow Simulation and Energy Loss Estimation in the Nappe Flow Regime of Stepped Spillways with Inclined Steps and End Sill : A Numerical Approach

Recently, the usage of stepped spillways, as energy dissipaters, has increased and led to a reduction in the size of the stilling basin. Extensive experimental considerations, plus the high cost and extended time required for laboratory methods, are among the major issues that require precise attention to determine optimal step design. This research deals with comparing the 2-D numerical simulation and experimental description in stepped spillways equipped with inclined steps and end sill together and presents a brisk, reliable, low-cost, and non-experimental approach to designing the steps. In this new type and complicated geometry, simulation is more complicated than horizontal steps, because it needs more accuracy around the end sills. The VOF Method and the k-ε standard turbulence model are proposed to simulate the flow pattern and evaluate the energy loss over stepped spillway. Energy dissipations obtained through the numerical approach have been compared with laboratory measurements and demonstrate reasonable agreement. Also, the flow pattern, velocity vectors and flow direction resulted from numerical simulation is in a good agreement with the experimental results.


Introduction
Selection of an efficient hydraulic structure to dissipate flow energy affects the length and size of the stilling basin.In recent years, numerous studies have addressed different aspects of stepped spillways; these studies have also introduced them as an efficient energy dissipater [1,2].The discharge rate in stepped spillways leads to three different flow regimes.The Nappe Flow is observed in the low-rate flows; the Transition Flow is represented when the discharge rate increases; and the Skimming Flow is caused by high discharge rates [1,2].Chen et al. analyzed the flow on stepped spillways by using the finite volume method and utilized the k-ε turbulent model to determine flow turbulence in spillways [3].Tabbara et al. analyzed the stepped spillway by the finite element method using ADINA software with the k-ε standard turbulence model [4].Musavi-Jahromi et al. simulated the flow over stepped spillway using ANSYS software [5].Their outcomes indicate that numerical results have a six percent differential with experimental ones.The k-ε model was initially formulated by Launder et al. [6].As the appropriate formulas, two equations for the k-ε turbulence model have been recognized to simulate spillway flow.Other work has been conducted on flow condition, aeration, and energy loss in stepped spillways with horizontal steps [7].Zare and Doering considered flow characteristics on stepped channels and spillways [8,9].Moreover, other works on flow condition, and energy dissipation in stepped spillways equipped with inclined steps or an end sill, have been conducted and the results indicate that inclined steps are reduced more energy than horizontal steps [10][11][12][13].Among all the investigations conducted on spillways, the research performed by Chaturabul [14] can be singled out.In his research, the height of the applied end sills was considered at 5, 10, and 15 millimeters on various step slopes.The result of his investigation demonstrated that relative energy loss increased 8% due to the existence of an end sill.Chinnarasri and Wongwises examined inclined steps and the end sill to determine the energy dissipation rate independently [15,16].They noted that using an end sill is the best approach to increase the energy dissipation rate.As mentioned previously, some researcher simulated flow over horizontal stepped spillways, but no one had tried to simulate flow over a stepped spillway equipped with inclined steps and an end sill together.In this case, simulation is more difficult than horizontal steps, because it needs more accuracy around the end sills.In this research, the finite volume method and the k-ε turbulence model in FLUENT 6.3 have been utilized in two-dimensional mode to simulate the flow pattern, find velocity vectors and flow direction also estimate the energy loss over a stepped spillway equipped with inclined steps with an end sill together.Flow depth and velocity are calculated numerically in critical sections, and energy loss is computed at each section.Dimensionless parameters are used in this research also experimental set up and results are presented in the first part of the research then reasons to select grid type, boundary condition type and turbulence model are discussed in the numerical part of the research.In the end, the efficiency and accuracy of the proposed numerical model to simulate the flow pattern, velocity vectors, flow direction and evaluate energy loss rate in a stepped spillway equipped with inclined steps with an end sill together are discussed.

Energy Loss Equations
The experimental formulas for the Nappe Flow regime are discussed in this section.Among the proposed equations, the equations developed by Chanson also Chamani and Rajaratnam [1,17] are considered the best ones for the Nappe Flow regime.The following equation was proposed by Chanson [17] to determine the energy dissipation rate in the Nappe Flow along with the hydraulic jump in the stepped spillways.
Where is total energy ( ); ∆ is energy dissipated in the length of the chute; is the critical depth of the flow (m); Hdam is the height of the dam and h is the height of the step (m).Chamani and Rajaratnam also presented the subsequent equation to obtain the energy dissipation rate in all Nappe Flows in stepped spillways [1]: Where α = coefficient of energy loss for each step and N = number of steps and l s = horizontal step length (m).

Dimensional Analysis
Buckingham's Pi Method, as one of the most appropriate methods in dimensional analysis, has been used.Chinnarasri and Wongwises [16].Considered effective parameters for energy dissipation in sloped stepped spillways, including: discharge per unit width, flow head, step height, step length, slope height, and gravity acceleration.In this research, other factors such as spillway slope and water density are also considered as follows.Moreover, the number of steps is regarded as another effective variable.Effective factors for energy dissipation in a stepped spillway equipped with inclined steps and the end sill together include:  L =  L (,  0 , ℎ, , , , , ,   , ) Dimensionless parameters obtained from Buckingham's Pi theories:  = f( L , ,  0 , ℎ, , , , , ,   , ) The number of variables is n=11.
The performed dimensional analysis is for a stepped spillway equipped with inclined steps and the end sill together.In this research, only the "m/h" dimensionless parameter is used.

Experimental Setup
Current research has been conducted at the Water Research Institute in Iran.Steps and walls were made of Plexiglass and mounted on the steel frame.The spillway in the chute is a broad-crested weir and the number of steps is 60.In the present investigation, only four steps are inclined and have a sill; all are placed after the middle of the chute (steps 39-42).The horizontal length of the steps is 14 centimeters; the step height is 4.66 centimeters and the chute width is 1.33 meters.The height of the broad-crested weir to the first step is 5 centimeters.Measured parameters during the test include both depth and velocity.The experiments have been conducted for one discharge in the Nappe Flow regime (30 liters per second, q=0.0225 m 2 /s).
Flow depth and velocity have been measured using a liminimeter (with a precision of 1 millimeter) and a pitot-tube, respectively.Flow discharge has been measured by the sharp crested weir at the end of the downstream chute.Characteristics of the used end sills are presented in Figure 1 and Table 1.Three different slopes of 7°, 10°, and 12°, respectively, about the horizon were used.The depth along the spillway width was measured by a liminimeter across each line of the piezometers.This measurement included three depths for Jet 1 (only water), Jet 2 (a mixture of water and air), and Jet 3, which is flow spray.Due to insufficient liminimeter accuracy and the possibility of an outlook error, pictures of the steps were taken during the experiments.Then, the pictures were imported into Arc GIS software to accurately determine the depths.The depths obtained from this method were then compared with liminimeter values.Finally, an average of the three depths across each line of the piezometer was used as the step depth.
To measure velocity, the pitot-tube was used.This measurement was arduous, since the flow was a two-phase flow and bleeding had to be performed regularly.The velocity was measured on both sides and in the middle at 0.6 meters of depth.Finally, averages of the left, right, and middle sections were averaged and set as the average velocity of that step.
The piezometer was used to statically measure pressure oscillation.These data points were used to enter pressures as initial values in a numerical simulation.In this chute, only four steps (39-42) were inclined and had sill.Depth and velocity values were measured on these steps, as well as on the 38 th step.To determine energy loss, data for the 38 th and 42 nd steps are sufficient.Velocity and depth values were recorded on other steps for future investigations.The entrance flow from the reservoir into the chute flows through the transferring pipes from the pumps and falls into the reservoir.At times that the reservoir completely fills, the vent opens and flow falls into the chute.The utilized pump has a capacity of 220 litters per second and swing check valves were used to open and close the pipes.The following equations were employed to determine upstream and downstream energy for the amended steps Where H is flow energy, y is flow depth, V is flow velocity in the proper sections, g is gravity acceleration and Z is the height above the baseline (The bottom of the last modified step has been assumed to be the baseline, and Z=0).

Numerical Approach
The numerical approach has been applied to simulate flow and estimate energy loss over stepped spillway equipped with an inclined step and an end sill, together.FLUENT 6.3 has been selected for this purpose.In this section, applied grid characteristics, the grid adaption method, boundary conditions, the applied turbulent model and the solver method details are discussed.

Geometry
The Finite Volume Method is an integral-based method.This means that initial differential equations should be integrated into a physical space, and then should be solved using numerical methods.Therefore, point's grid should be developed directly in a physical space.Due to this, there is a non-zero probability that points in the grid don't relate to lines in the grid.These kinds of grids are known as unstructured grids.A two-dimensional unstructured grid applies triangular elements, whereas a three-dimensional one uses tetrahedral elements.Both triangular elements and tetrahedral elements are compatible with complex geometry.A structured grid is easy to develop in simple geometry [18].It needs less time and memory to resolve.However, in complex geometry, like the situation in this research, it is better to use unstructured grids with triangular elements.The time needed to develop an unstructured grid in this study (complex geometry) is less than the time needed to develop a structured grid.Although an unstructured grid needs more time and memory to solve the problem, it is still worth it to use this kind of grid in the study.The unstructured grid for this research was developed in Gambit software by triangular elements [19,20].

Meshing
The area around the end sill should be considered carefully.Due to this fact, the Region Adaption Technique has been applied.This makes triangular elements smaller around end sills, to achieve greater accuracy.The grid has been divided into small triangles and considered 0.001 up to 0.01 meters, depending on the required precision (Figure 2).Also, to increase grid quality, smoothing and swapping have been applied.Smoothing rearranges cells to improve grid quality.In FLUENT software, there are two methods for smoothing.The first one is the Laplace Method, and the second one is smoothing based on skewness.It is recommended to use smoothing based on skewness when triangular elements have been used (FLUENT 6.3 user's guide).Therefore, smoothing based on skewness has been applied, in this study, to decrease the grid's maximum skewness.This method tries to decrease the skewness in cells with the skewness higher than the maximum limit, which is 0.4 in this study, by rearranging the cells.The maximum skewness decreased, but the average skewness increased.This happens because when this method tries to reduce the skewness of critical cells, with a skewness of more than 0.4, the skewness of other cells will increase.Moreover, face swapping has been applied to the grid to increase quality.In triangle cells, the Delaunay Method has been applied to improve grid quality.In computational geometry, the Delaunay triangulation for a set of P points in a plane is triangulation DT (P), such that no point in P is inside the circumcircle of any triangle in DT (P).The Delaunay triangulation maximizes the minimum angle of all the angles of the triangles in the triangulation.Thus, it tends to avoid skinny triangles [21].

Boundary and Initial Condition
There are two sections in the inlet, as can be seen in Figure 3, which are separated by the water surface.For the first section, the velocity inlet boundary type has been applied, because this part is full of water with a certain velocity.The height of this section is selected as the flow depth in each test.Also, the velocity in each test is entered into the software by magnitude and direction.The x-direction is selected as a velocity direction.The second section is above flow, and the pressure in this section is equal to atmospheric pressure.Therefore, the pressure inlet boundary type is selected for this section.A similar situation to the second section is the area above the steps.So, the pressure inlet boundary condition is selected for this area.In the outlet, there are two sections again.The first section incudes flow and the second section includes air.Because pressure inlet is selected in the inlet, the outflow boundary condition cannot be used and the pressure outlet boundary condition should be used as boundary condition in these sections.In the end, the wall boundary condition has been applied to the steps.

Turbulence Modeling
The stepped spillway is a self-aerated hydraulic structure, and a two-phase flow (water and air) can be seen over the chute.To simulate flow over the stepped spillways, both the water and air phases should be defined in the software.The volume fraction is defined as "1" when the cell volume is full of water.The volume fraction is defined as "0" when the cell volume is full of air.The software will calculate the volume fraction of each cell based on the water and air percentage in that particular cell.When the Reynolds number is large, or the geometry is complicated, it is impossible to solve time-dependent Navier-Stokes Equations completely.There are two general methods to the Navier-Stokes Equations in the way that turbulent flow fluctuations don't directly enter into the equations.

Governing Equations
The first method is the Reynolds Average Navier-Stokes Equations (RANS).This time average-based method uses the average values of all turbulent flow fluctuations in transitional quantities.Therefore, it needs less time and memory to solve the problem.Large Eddy Simulation (LES) is the second method.In the LES Method, large eddies and turbulent flow fluctuations will be calculated time-dependently and directly by filtered Navier-Stokes Equations.It means that, the LES is using space filtering to filter small eddies and just considers and calculates large eddies carefully with details.This method needs parallel processing to calculate eddies and is expensive to use.In both methods, new variables and extra quantities are introduced in governing equations.These new variables and quantities should be calculated by the turbulent models [22].In this study, the RANS Method is used, because it can solve flows with a strong vortex and naturally unsteady flows with reasonable cost.Among turbulent models, the Spalart-Allmaras, k-ε , k-ω, and Reynolds Stress Model (RSM) use the RANS Method and all of them are available in FLUENT 6.3.In the RANS Method, continuity and momentum equations are as follows [23]: Where  is density,  is molecular viscosity,  is time,   is the coordinate component,   is the velocity component, and  is pressure.Reynolds stresses should be calculated to solve the problem with the RANS Method.One of the methods that can be used to calculate Reynolds stresses is the Boussinesq Hypothesis (Equation 14), which has a reasonable calculation cost.

𝜌𝑢 𝑖
Where   is the turbulence viscosity.  = 1 when = and   = 0 when .Some turbulent models like Spalart-Allmaras, k-, and k-use the Boussinesq Hypothesis.The RSM Turbulent Model uses another method to calculate Reynolds stresses.In this method, five equations should be solved in two dimensional flows.So, it is a very accurate, but expensive method.In this study, methods which use the Boussinesq Hypothesis (i.e Spalart-Allmaras, k-, and k-) are selected to be considered.Among these turbulent models, Spalart-Allmaras solves one equation to calculate turbulent viscosity.The k-style models (standard, RNG, and Realizable) and the k-turbulent models solve two equations to calculate turbulent viscosity.The number of equations denotes the number of additional Partial Differential Equations (PDEs) that should be solved.The k-style models and the k-are selected for more consideration, in this research, because they are using two equations to calculate turbulent viscosity.
Among the k-style models (standard, RNG, and Realizable) and the k-turbulent models, the k-standard is selected to be considered in this study as the first turbulent model attempt, due to its simplicity and reasonable cost.Moreover, it has been previously used to successfully simulate flow over stepped spillways [3].This turbulent model solves "k" and " "equations as additional PDEs to find turbulent viscosity." ()" is an instantaneous kinetic energy and is the sum the mean kinetic energy (K) and the turbulent kinetic energy (k).and (simplified) are as follows: There are five terms in both equations.From left to right, the terms are the rate of increase, convective transport, diffusive transport, the rate of production, and the rate of destruction, respectively.Where   is calculated by the following equation using and (turbulence energy dissipation rate): Where ρ is the fluid density,   is the experimental constant and equals to 0.09.Prandtl's turbulence numbers for and consist of =1.0, =1.3, 1 and which are the constants of the relation [24].Whenever the VOF method is selected, the k-turbulent model works as a two-phase flow model.Therefore, in this research (incompressible and turbulent flow), there is no need to use modified equations, like the continuity equation for twophase flow [3].
The Second Order Implicit pressure based solver has been used in this research.Also, the segregated algorithm has been applied, because it needs less memory than the coupled algorithm.This algorithm solves equations separately.Also, the second order method has been used for discretization to reach better results.This method is selected because the flow is not in the same direction with the triangular elements.The Standard Method and the First Order Upwind Method have been applied in pressure interpolation and density interpolation, respectively.Moreover, the First Order Upwind Method has been applied in turbulent kinetic energy and turbulent dissipation rate equations.The PISO method has been used as a pressure velocity coupling method.The PISO will give good results when the flow is unsteady.The Neighbour Correction and the Skewness Correction also have been used in PISO with an underrelaxation factor equal to one for all equations.The Green-Gauss Node Based Method has been applied to calculate derivations.This method is accurate in unstructured grids with triangular elements.This method uses the average values of the nodes in each cell to calculate derivations.The velocity inlet boundary condition and the velocity value in x-direction in this boundary condition region in each test have been used for initializing the solution across the entire flow field.Also, patching values, volume fraction equal to zero, have been applied in selected cells to reduce calculations.These selected cells are located above the chute, and it is obvious that there is no water in these cells.

Results and Discussion
The test was conducted on the horizontal steps, before applying the changes on the steps; to calculate the effects of the changes on the energy dissipation rate.The energy loss obtained from the horizontal steps was 0.5128 and 0.5112, respectively, derived from the Chamani and Rajaratnam equations (Equations 2 to 5) and from the test (Equations 10 and 11), which shows acceptable agreement.The changes were then applied to the steps, and the energy dissipation was measured for different slopes and end sills.In this section, the parameter m= (p+w) is used, where p is the height of the end sill, and w is the height of the step inclination (Figure 4).In Figure 5, an overall conclusion of the data regarding all slopes has been presented for 5 and 10 millimeter thicknesses.Additionally, the fitted curve between different points has been presented.As it can be seen, the curve type is quadratic.It is obvious that if the "m/h" ratio increases, up to 0.7, the dissipation energy rate increases.After that, it decreases.This graph suggests that the best ratio for "m/h" is approximately 0.7, and an excessive increase negatively affects the energy dissipation rate.
In practical design, to receive an acceptable energy loss, it is better to eliminate "m/h" with amounts greater than 0.7.When the "m/h" ratio goes past 0.7, the flow jumps from one or more steps and the energy loss rate decrease (Figures 6a and 6b).This step plays practically no role in energy dissipation.As a consequence, the energy dissipation rate decreases.The results for the inclined steps and end sills for the "m/h" equal to or greater than 0.7 (effective end sills) have been presented in Table 2.The obtained results show that the end sill height and the thickness, as well as the step height incline (w), affect the energy dissipation rate.Table 2, also, indicates that "the maximum energy loss increase" is 16.36%, and "the average energy loss increase" is 14.52%.This amount of increase in the energy loss is considerable and proves that step modification (inclined step and end sill) is an effective technique to increase the stepped spillways' efficiency in terms of energy dissipation.Numerical analysis results are indicated as velocity vectors and water surface profile (Figures 7c and 7d).As can be seen in Figure 7a, after the flow encounters the end sill in the first modified step, it jumps from the second modified step.During the jump, some parts of the flow can't reach the third step and fall into the second step.Because the step is reverse inclined, this flow will stay in the second step and make a pool.The flow direction in the pool is from the end sill to the wall (from right to left) due to the step slope.Other parts of the flow jump to the third step, and during the jump combine with air, making highly aerated flow over the chute.Figure 7c shows the flow profile simulation.This flow profile is the same as the flow profile obtained in the laboratory.In Figure 7c, the flow jumps from the second step combined with the air (highly aerated flow) and the pool can be seen in the second step.Also, flow direction is illustrated in Figure 7d as velocity vectors.The flow direction is also similar to experimental results (Figure 7b).The velocity magnitude is larger in the jump flow and is smaller in the pool (Figure 7d).Experimental results confirm these values.Some parts of the jumped flow, also, jump from the second step, but other parts of the jumped flow come back to the chute and impact the third modified step.Flow is highly aerated in the third step due to this impact.This can be seen in Figure 7c.Moreover, the flow direction in the third step pool is similar to the second step pool, because of the step slope (Figure 7d).This is confirmed by experimental results.a b  3, there are eight tests which have been employed to approach numerical calibration to simulate energy loss in the Nappe Flow regime.Moreover, the minimum error for calibration data is 4.58%, and the maximum error is 6.58%.Also, the average error for calibration data is 5.62%.4. shows numerical simulation results in which validation data have been applied.The minimum error for validation data is 4.59%, and the maximum error is 6.72%.Furthermore, the average error for validation data is 5.48%.Figures 8a and 8b indicate the energy loss rate obtained from the measurement and simulation versus the "m/h" dimensionless parameter for the calibration and validation data sets.Although there is a difference between simulation results and measurement, the simulation results are deemed agreeable.As illustrated in Figures 8a and 8b, the energy loss versus the "m/h" follows the same trend in measurement data and simulation results in both the calibration and validation data sets.

Conclusion
In this study, four steps (39 to 42) of a sixty-step physically stepped spillway model have been equipped with inclined steps and an end sill, together.End sills with various heights and thickness have been examined in three reverse stepped slopes (7, 10 and 12 degrees, respectively).The dimensionless parameter, "m/h," has been used to determine the energy dissipation rate.Results indicate that in the discharge 30 L/s ( = 0.0225  ), which represents the Nappe Flow regime, inclined steps and an end sill with an "m/h" less than 0.7 effectively increases the energy dissipation rate.The average of this incensement is obtained 14.52% in effective end sills.Moreover, FLUENT 6.3 software was used in two-dimensional mode in this research to simulate the flow over a stepped spillway equipped with inclined steps with end sills.A numerical approach using the finite volume method, the k-ε standard turbulent model, and triangular elements and an unstructured grid has been successfully applied to simulate flow pattern, velocity vectors, and flow direction over the complicated geometry of stepped spillways equipped with inclined steps and end sill together.Also, the energy loss has been estimated using the numerical model.As a first step, eight tests were selected for calibration; afterwards, four other tests were conducted for validation.Numerical outcomes for depth and velocity were used to calculate energy loss.The results illustrate some discrepancies between the numerical and experimental outcomes.Average errors, generated by comparing the energy loss results between the experimental and numerical methods, were 5.62% and 5.48% for the calibration and validation data, respectively.The magnitude of these errors indicates that instead of time-consuming procedures and high-priced laboratory work, numerical models can be used as a low cost and expeditious method to estimate the energy dissipation rate in a stepped spillway equipped with inclined steps and end sills together.

Acknowledgment
The authors wish to acknowledge the support received in the form of a graduate assistantship from the Civil and Environmental Engineering Department, Florida International University.Also, the authors wish to acknowledge the Water Research Institute, Iran for providing the model and instruments to conduct experiments in the hydraulic laboratory.

Notation
The following symbols are used in this paper:

Figure 7 .
Figure 7. (a) Flow profile -Steep 10 o -End sill 6-5 mm, (b) Flow profile -Steep 10 o -End sill 6-5 mmdetails, (c) Flow profile -Steep 10 o -End sill 6-5 mm, (d) Velocity vectors -Steep 10 o -End sill 6-5 mm Twelve tests have been performed.The first eight tests are used for calibration, and four other tests are used for validation.In Table3, there are eight tests which have been employed to approach numerical calibration to simulate energy loss in the Nappe Flow regime.Moreover, the minimum error for calibration data is 4.58%, and the maximum error is 6.58%.Also, the average error for calibration data is 5.62%.

Figure 8 .
Figure 8.(a) Measurement and simulation energy loss in percentage terms -Calibration data, (b) Measurement and simulation energy loss in percentage terms -Validation data

1 1 ) 1 ) 1 1 )
Constant in the k-turbulent model ∆ Dissipated energy in the length of the chute Constant in the k-turbulent model a Variable in the Chamani and Rajaratnam Equation Constant in the k-turbulent model b Variable in the Chamani and Rajaratnam Equation c Height of slope + end sill (mm) 2) Basic variables (Pi Buckingham Theorem) Total energy (H dam +3/2h c ) (m) n Number of variables in the Pi Buckingham Theorem Residual height (m) (Fratino and Piccini) Number of steps h s Step height (cm) End sill height (mm) s Horizontal step length (cm, m) Unit discharge (m 3 /m.s) Velocity component in the x-direction (m/s) r Basic (main) parameters in dimension analysis Velocity component in the y-direction (m/s) Chute slope Velocity component in the z-direction (m/s) T Time in dimension analysis (second) , ′ Velocity component in the i-or j-direction (turbulent) (m/s) End sill thickness (mm) 2) time (s) Velocity in section 1 (m/s) Height caused by the slope in inclined step (mm) Velocity in section 2 (m/s) Direction Depth in section 2 (cm) Distance to the baseline (m) Water Density (kg/m 3 ) All parameters in dimension analysis Constant in the k-turbulent model Energy loss coefficient Constant in the k-turbulent model Turbulence energy dissipation rate  1 Flow depth on step (cm) 2) Depth in section 1 (cm)  Dynamic viscosity (kg.m/s)   Turbulent viscosity In Buckingham Method

Table 1 . End sill specifications Figure 1. Schematic of the end sill
)