Numerical investigation of flow characteristics over stepped spillways

Dam structures are the most important projects around the world to store water or to transport water because protection of water is the key to living. And the spillway is classified as one of the most important parts of a dam. A spillway is constructed to protect a dam from destruction or damage by flood. Dam building and flood control can be considered a very important issue across the world given the importance of hydroelectric power generation, navigation, recreation and fishing. There are many types of spillway, but the most common types are: ogee spillways, free over-fall spillways, siphon spillways, chute spillways, side channel spillways, tunnel spillways, shaft spillways and stepped spillways. And every spillway consists of four necessary components: an inlet channel, control structure, discharge carrier and outlet channel. A large number of stepped spillways were constructed through the recent decades, as especially associated with the technique of roller compacted concrete (RCC) dam construction and construction of stepped spillways classified as an easier, quicker and cheaper technique of construction (Chanson 2002; Felder & Chanson 2011). A stepped spillway structure increases the rate of energy dissipation which decreases cavitation risk (Boes & Hager 2003b). And stepped spillways have advantages which make them more attractive under various conditions.

The flow behaviour in stepped spillways is generally classified by three different regimes: nappe, transition and skimming flow regimes (Chanson 2002). When the flow rate is low the nappe flow regime occurs and it is characterized as a sequence of free-falling nappes, while in a skimming flow regime, water flows over steps as a coherent stream on the pseudo-bottom above the outer step edges. It is also clear that under the main flow three-dimensional recirculating vortices occur (e.g., Chanson 2002; Gonzalez & Chanson 2008). At the pseudo-bottom near the step edge the direction of the flow is virtually aligned with the pseudo-bottom. According to Takahashi & Ohtsu (2012), for a given flow rate in a skimming flow regime the flow impacts the horizontal step face near the step edge and, with a decreasing chute slope, the area of the impact region will increase. Transition flow regime occurs between the nappe flow and skimming flow regime. In design of stepped spillways, the skimming flow regime should be considered (e.g., Chanson 1994; Matos 2000; Chanson 2002; Boes & Hager 2003a).

Computational fluid dynamics (CFD), namely numerical models of hydraulic engineering, generally reduce the amount of total cost and time that will be spent on physical models. So numerical models are classified as faster and cheaper than experimental models and can also be used for more than one purpose at the same time. There are many CFD software packages available but the most widely used one is Flow-3D. In this study, Flow 3D software is used to simulate the air concentration, velocity distribution and dynamic pressure distribution on a stepped spillway for two different models with different flow rates.

Roshan et al. (2010) studied the investigation of flow regimes and energy dissipation over two physical models of stepped spillways with different numbers of steps and discharges. The slope of the experimental models was 19.2% and the number of steps 12 and 23, respectively. The results illustrated that the observed flow regimes in the 23-step physical model was considered more acceptable than the 12-step model. However, the energy dissipation on the 12-step model was more than on the 23-step model. And the experiments observed that in the skimming flow regime the energy dissipation in the 23-step model was less than in the 12-step model by about 12%.

Ghaderi et al. (2020a) carried out experimental studies of stepped spillways to investigate the influence of scouring parameters with different step sizes and flow rates. The results showed that the flow regime affected the scour-hole dimensions, such as the minimum scouring depth which happened under the nappe flow regime. Moreover, the tailwater depth and step size are actual parameters for maximum scouring depth. By increasing the depth of tailwater from 6.31 cm to 8.54 and 11.82 cm this increased the scouring depth by 18.56% and 11.42% respectively. Also, this increasing tailwater depth decreases the scouring length by 31.43% and 16.55% respectively. In addition, the Froude number increases by increasing the flow rate, and the increased momentum of the flow promotes scouring. Also, the results indicated that the scouring in the middle is less than at the sidewalls of the cross-section. An empirical formula was suggested to predict the maximum depth of scouring downstream of stepped spillways and then compared with the experimental results. And the comparison illustrated that the suggested formula can predict the depth of scouring within 3.86% and 9.31% relative and maximum errors, respectively.

Ghaderi et al. (2020b) made a numerical investigation of trapezoidal labyrinth shaped (TLS) steps. The results observed that these types of spillways have better performance because they increase the magnification ratio L_T/W_t (L_T is the total edge length, W_t is the width of spillway). Also, the trapezoidal labyrinth shaped stepped spillway has a larger friction factor and a lower residual head. The friction factor differs from 0.79 to 1.33 for various magnification ratios, while for the flat stepped spillway it is approximately equal to 0.66. Also, the ratio of the residual head (H_res/dc) is approximately 2.89 in a TLS stepped spillway, while it is approximately equal to 4.32 for a flat stepped spillway.

Shahheydari et al. (2015) investigated the flow over a stepped spillway in skimming flow regime by using Flow-3D software, the RNG k-ε model and Volume of Fluid (VOF) method to study the profiles of free surface flow such as discharge coefficient and energy dissipation and compared them with the expermental results. The results showed that the relationship between the energy dissipation rate and discharge coefficient rate was inverse, and were in good agreement with the results from the experimental model.

Mohammad Rezapour Tabari & Tavakoli (2016) investigated the effect of various parameters such as step height (h), step length (L), number of steps (Ns) and the discharge in unit width (q) on the energy dissipation in a stepped spillway. They used Flow-3D software in the analysis to evaluate the relationship between the energy loss and the critical depth of flow in a stepped spillway. Moreover, the finite volume method was applied to solve the equations and the standard k-ɛ model used for flow turbulence. According to the results, the energy dissipation decreases when the number of steps increases and the flow discharge increases. The gained results were compared with other studies, empirical and mathematical investigations were conducted, and eventually passable results were acquired.

Turbulence modelling is one of three key elements in CFD (Gunal 1996). There are many types of turbulence models, but the most common are Zero-equation models, One-equation models, Two-equation models, Reynolds Stress/Flux models and Algebraic Stress/Flux models. In FLOW-3D software, five turbulence models are available. The formulation used in the FLOW-3D software differs slightly from other formulations that includes the influence of the fractional areas/volumes of the FAVOR^TM method and generalizes the turbulence production (or decay) associated with buoyancy forces. The latter generalization, for example, includes buoyancy effects associated with non-inertial accelerations.

The available turbulence models in Flow-3D software are the Prandtl Mixing Length Model, the One-Equation Turbulent Energy Model, the Two-Equation Standard Model, the Two-Equation Renormalization-Group (RNG) Model and large Eddy Simulation Model (Flow Science Inc. 2012).

VOF and FAVOR are classifications of volume-fraction methods. In these two methods, firstly the area should be subdivided into a control volume grid or a small element. Each flow parameter like velocity, temperature and pressure values within the element are computed for each element containing liquids. Generally, these values represent the volumetric average of values in the elements.

The FAVOR method is a different method and uses another volume fraction technique, which is only used to define the geometry, such as the volume of liquid in each cell used to determine the position of fluid surfaces. Another fractional volume can be used to define the solid surface. Then, this information is used to determine the boundary conditions of the wall that the flow should be adapted for.

In this study, the experimental results of Ostad Mirza (2016) was simulated. In a channel composed of two 4 m long modules, with a transparent sidewall of height 0.6 m and 0.5 m width. The upstream chute slope (i.e. pseudo-bottom angle) Ɵ1 = 50°, the downstream chute slope Ɵ2 = 30° or 18.6°, the step heights h = 0.06 m, the total number of steps along the 50° chute 41 steps, the total number of steps along the 30° chute 34 steps and the total number of steps along the 18.6° chute 20 steps.

The flume inflow tool contained a jetbox with a maximum opening set to 0.12 meters, designed for passing the maximum unit discharge of 0.48 m²/s. The measurements of the flow properties (i.e. air concentration and velocity) were computed perpendicular to the pseudo-bottom as shown in Figure 1 at the centre of twenty stream-wise cross-sections, along the stepped chute, (i.e. in five steps up on the slope change and fifteen steps down on the slope change, namely from step number −09 to +23 on 50°–30° slope change, or from −09 to +15 on 50°–18.6° slope change, respectively).

Pressure sensors were arranged with the x/l values for different slope change as shown in Table 1, where x is the distance from the step edge, along the horizontal step face, and l is the length of the horizontal step face. The location of pressure sensors is shown in Table 1.

A 3D numerical model of hydraulic phenomena was simulated based on an experimental study by Ostad Mirza (2016). The water surcharge and flow pressure over the stepped spillway was computed for two models of a stepped spillway with different discharge for each model. In this study, the package was used to simulate the flow parameters such as air entrainment, velocity distribution and dynamic pressures. The solver uses the finite volume technique to discretize the computational domain. In every test run, one incompressible fluid flow with a free surface flow selected at 20̊ was used for this simulation model. Table 2 shows the variables used in test runs.

For stepped spillway simulation, several parameters should be specified to get accurate simulations, which is the scope of this research. Viscosity and turbulent, gravity and non-inertial reference frame, air entrainment, density evaluation and drift-flux should be activated for these simulations. There are five different choices in the ‘viscosity and turbulent’ option, in the viscosity flow and Renormalized Group (RNG) model. Then a dynamical model is selected as the second option, the ‘gravity and non-inertial reference frame’. Only the z-component was inputted as a negative 9.81 m/s² and this value represents gravitational acceleration but in the same option the x and y components will be zero. Air entrainment is selected. Finally, in the drift-flux model, the density of phase one is input as (water) 1,000 kg/m³ and the density of phase two (air) as 1.225 kg/m³. Minimum volume fraction of phase one is input equal to 0.1 and maximum volume fraction of phase two to 1 to allow air concentration to reach 90%, then the option allowing gas to escape at free surface is selected, to obtain closer simulation.

The flow domain is divided into small regions relatively by the mesh in Flow-3D numerical model. Cells are the smallest part of the mesh, in which flow characteristics such as air concentration, velocity and dynamic pressure are calculated. The accuracy of the results and simulation time depends directly on the mesh block size so the cell size is very important. Orthogonal mesh was used in cartesian coordinate systems. A smaller cell size provides more accuracy for results, so we reduced the number of cells whilst including enough accuracy. In this study, the size of cells in x, y and z directions was selected as 0.015 m after several trials.

Figure 3 shows the 3D computational domain model 50–18.6 slope change, that is 6.0 m length, 0.50 m width and 4.23 m height. The 3D model of the computational domain model 50–30 slope changes this to 6.0 m length, 0.50 m width and 5.068 m height and the size of meshes in x, y, and z directions are 0.015 m. For the 50–18.6 slope change model: both total number of active and passive cells = 4,009,952, total number of active cells = 3,352,307, include real cells (used for solving the flow equations) = 3,316,269, open real cells = 3,316,269, fully blocked real cells equal to zero, external boundary cells were 36,038, inter-block boundary cells = 0 (Flow-3D report). For 50–30 slope change model: both total number of active and passive cells = 4,760,002, total number of active cells equal to 4,272,109, including real cells (used for solving the flow equations) were 3,990,878, open real cells = 3,990,878 fully blocked real cells = zero, external boundary cells were 281,231, inter-block boundary cells = 0 (Flow-3D report).

When solving the Navier-Stokes equation and continuous equations, boundary conditions should be applied. The most important work of boundary conditions is to create flow conditions similar to physical status. The Flow-3D software has many types of boundary condition; each type can be used for the specific condition of the models. The boundary conditions in Flow-3D are symmetry, continuative, specific pressure, grid overlay, wave, wall, periodic, specific velocity, outflow, and volume flow rate.

There are two options to input finite flow rate in the Flow-3D software either for inlet discharge of the system or for the outlet discharge of the domain: specified velocity and volume flow rate. In this research, the X-minimum boundary condition, volume flow rate, has been chosen. For X-maximum boundary condition, outflow was selected because there is nothing to be calculated at the end of the flume. The volume flow rate and the elevation of surface water was set for Q = 0.1 and 0.235 m³/s respectively (Figure 2).

The bottom (Z-min) is prepared as a wall boundary condition and the top (Z-max) is computed as a pressure boundary condition, and for both (Y-min) and (Y-max) as symmetry.

The air concentration distribution profiles in two models of stepped spillway were obtained at an acquisition time equal to 25 seconds in skimming flow for both upstream and downstream of a slope change 50°–18.6° and 50°–30° for different discharge as in Table 2, and as shown in Figure 4 for 50°–18.6° slope change and Figure 5 for 50°–30° slope change configuration for dc/h = 4.6. The simulation results of the air concentration are very close to the experimental results in all curves and fairly close to that predicted by the advection-diffusion model for the air bubbles suggested by Chanson (1997) on a constant sloping chute.

But as is shown in all above mentioned figures it is clear that at the pseudo-bottom the CFD results of air concentration are less than experimental ones until the depth of water reaches a quarter of the total depth of water. Also the direction of the curves are parallel to each other when going up towards the surface water and are incorporated approximately near the surface water. For all curves, the cross-section is separate between upstream and downstream steps. Therefore the (-) sign for steps represents a step upstream of the slope change cross-section and the (+) sign represents a step downstream of the slope change cross-section.

The dimensionless velocity distribution (V/V₉₀) profile was acquired at an acquisition time equal to 25 seconds in skimming flow of the upstream and downstream slope change for both 50°–18.6° and 50°–30° slope change. The simulation results are compared with the experimental ones showing that for all curves there is close similarity for each point between the observed and experimental results. The curves increase parallel to each other and they merge near at the surface water as shown in Figure 6 for slope change 50°–18.6° configuration and Figure 7 for slope change 50°–30° configuration. However, at step numbers +1 and +5 in Figure 7 there are few differences between the simulated and observed results, namely the simulation curves ascend regularly meaning the velocity increases regularly from the pseudo-bottom up to the surface water.

Figure 8 (50°–18.6° slope change) and Figure 9 (50°–30° slope change) compare the simulation results and the experimental results for the presented dimensionless dynamic pressure distribution for different points on the stepped spillway. The results show a good agreement with the experimental and numerical simulations in all curves. For some points, few discrepancies can be noted in pressure magnitudes between the simulated and the observed ones, but they are in the acceptable range. Although the experimental data do not completely agree with the simulated results, there is an overall agreement.

The pressure profiles were acquired at an acquisition time equal to 70 seconds in skimming flow on 50°–18.6°, where p is the measured dynamic pressure, h is step height and ϒ is water specific weight. A negative sign for steps represents a step upstream of the slope change cross-section and a positive sign represents a step downstream of the slope change cross-section.

Figure 10 shows the experimental streamwise development of dimensionless pressure on the 50°–18.6° slope change for dc/h = 4.6, x/l = 0.35 on 50° sloping chute and x/l = 0.3 on 18.6° sloping chute compared with the numerical simulation. It is obvious from Figure 10 that the streamwise development of dimensionless pressure before slope change (steps number −1, −2 and −3) both of the experimental and simulated results are close to each other. However, it is clear that there is a little difference between the results of the streamwise development of dimensionless pressure at step numbers +1, +2 and +3. Moreover, from step number +3 to the end, the curves get close to each other.

Figure 11 compares the experimental and the numerical results for the streamwise development of the dimensionless pressure on the 50°–30° slope change, for dc/h = 4.6, and x/l = 0.35 on 50° sloping chute and x/l = 0.17 on 30° sloping chute. It is apparent that the outcomes of the experimental work are close to the numerical results, however, the results of the simulation are above the experimental ones before the slope change, but the results of the simulation descend below the experimental ones after the slope change till the end.

In this research, numerical modelling was attempted to investigate the effect of abrupt slope change on the flow properties (air entrainment, velocity distribution and dynamic pressure) over a stepped spillway with two different models and various flow rates in a skimming flow regime by using the CFD technique. The numerical model was verified and compared with the experimental results of Ostad Mirza (2016). The same domain of the numerical model was inputted as in experimental models to reduce errors as much as possible.

Flow-3D is a well modelled tool that deals with particles. In this research, the model deals well with air entrainment particles by observing their results with experimental results. And the reason for the small difference between the numerical and the experimental results is that the program deals with particles more accurately than the laboratory. In general, both numerical and experimental results showed that near to the slope change the flow bulking, air entrainment, velocity distribution and dynamic pressure are greatly affected by abrupt slope change on the steps. Although the extent of the slope change was relatively small, the influence of the slope change was major on flow characteristics.

The Renormalized Group (RNG) model was selected as a turbulence solver. For 3D modelling, orthogonal mesh was used as a computational domain and the mesh grid size used for X, Y, and Z direction was equal to 0.015 m. In CFD modelling, air concentration and velocity distribution were recorded for a period of 25 seconds, but dynamic pressure was recorded for a period of 70 seconds. The results showed that there is a good agreement between the numerical and the physical models. So, it can be concluded that the proposed CFD model is very suitable for use in simulating and analysing the design of hydraulic structures.

All relevant data are included in the paper or its Supplementary Information.