Numerical Investigation of Different Stepped Spillway Geometries over a Mild Slope for Safe Operation Using Multi-Phase Model

Abstract

The appropriate design and operation of spillways are critical for dam safety. To enhance design practices and gain insights into flow hydraulics, both experimental and numerical modeling are commonly employed. In this study, we conducted a numerical investigation of flow over a mildly sloping (1V:3H) stepped spillway with various step geometries using a multi-phase mixture model with dispersed interface tracking in ANSYS Fluent. The model was validated against experimental data from Utah State University, focusing on water surface profiles over the crest, velocities, and air concentrations. The validated numerical model was used to simulate flow over different step geometries (i.e., 0.2 m H uniform Step, 0.1 m H uniform step, non-uniform steps, adverse slope steps, and stepped pool) for a range of discharges from 0.285 m³/s/m to 1.265 m³/s/m. While flow depths over the crest and velocities in the chute compared well with experimental results, air concentrations exhibited some deviation, indicating numerical limitations of the solver. The shift in the location of the inception point was found to be mainly influenced by a higher flow rate than the different design configurations over an identical mild slope. The downstream non-linear flow velocity curve with different flow rates indicated less effectiveness of the step roughness over a high flow rate as a result of the reduction in relative roughness. The theoretical velocity ratio indicated the least reduction in downstream velocity with the stepped pooled spillway due to the formation of a “stagnant pool”. A higher negative-pressure region due to flow separation at the vertical face of the steps was obtained by adverse slope steps, which shows that the risk of cavitation is higher over the adverse slope step spillway. Turbulent kinetic energy (TKE) was found to be higher for uniform 0.2 m H steps due to the strong mixing of flow over the steps. The least TKE was found at the steps of the stepped pool spillway due to the formation of a “stagnant pool”. Uniform 0.2 m H steps achieved the maximum energy dissipation efficiency, whereas the stepped pool spillway obtained the least energy dissipation efficiency, introducing higher flow velocity at the stilling basin with a higher residual head. The adverse slope and non-uniform steps were found to be more effective than the uniform 0.1 m H steps and stepped pool spillway. The application of uniform steps of higher drop height and length could achieve higher TKE over the steps, reducing the directional flow velocity, which reduces the risk of potential damage.

1. Introduction

Dam safety issues are often associated with spilling facilities during extreme flood events. The safe release of excess water from the dam is essential to ensure the safety of the spillway. The spillway exhibits higher flow velocity due to the transformation of potential head to kinetic energy. The higher flow velocity zone creates a low-pressure area, leading to cavitation [1], whereas water-borne particles like rocks, sediment, and debris can deteriorate the spillway surface due to abrasion. The incident of Oroville Dam [2], Toddbrook Dam [3], Ruskin Dam [4], Glen Canyon Dam [5,6] and Nagarjuna Sagar Dam [7] spillways are some well-known examples of the potential for overall structural damage. Aeration in spilling facilities helps to reduce the chances of cavitation damage; therefore, artificial aeration or self-aeration structures are advantageous for the safe spillway design [8,9]. In recent decades, different types of spillways have been implemented for releasing the excess water from the reservoir. Most hydraulic issues in the spillway surface and stilling basin have been associated with smooth spillways as compared to rough ones due to higher chances of cavitation and abrasion as a result of higher flow velocity and granular flow. Several studies [10,11,12,13] have been carried out to dissipate the energy with the application of stepped spillways, which helps to reduce energy and minimize the chances of cavitation at the spillway surface as well as in the stilling basin. The stepped spillway has been implemented in several dams due to its high energy dissipation efficiency as it helps to reduce the length of the stilling basin [11].

Chanson [14] and Boes [15] made a detailed experimental investigation of flow regimes (i.e., nappe, transition, and skimming) and aerated flow properties, including energy dissipation over stepped spillways. It was highlighted that stepped spillways cause an early occurrence of the inception point of aeration and facilitate a significant increase in energy dissipation compared to smooth chutes. Felder et al. [16] tested the application of non-uniform step height to evaluate the efficiency of energy dissipation. The rate of energy dissipation was found to be more or less similar and suggested some observed insatiability for lower flow rates due to the application of non-uniform step height. Nina et al. [17] carried out a detailed experimental investigation for air–water flow velocity measurements using dual-tip phase detection probes, and optical flow-data-based techniques using high-speed video on self-areated regions. Nina et al. [17] also identified the reliable results of a dual-tip phase detection probe velocity measurement compared to the optical flow data-based measurement, Whereas the side wall optical flow presented good flow behaviors of cavity recirculation over a highly turbulent regime.

With recent advancements in computation power, computational fluid dynamic (CFD) modeling is also being used in addition to experimental investigations to study flows in engineering applications. Several self-aerated flows such as plunging jets [18,19,20] and hydraulic jumps [21,22] have been investigated using CFD techniques. Furthermore, numerical modeling has also been applied to study complex multi-phase flow over stepped spillways. Salmasi [23] compared the effect of stepped and smooth spillways on Zirdan Dam for the skimming flow regime using a multi-phase mixture model. A significant reduction in the boundary layer length was observed over a stepped spillway with quick self-aeration. Ghaderi [24] studied the influence of different stepped pool geometry on energy dissipation with RNG k-$ϵ$ turbulence model using FLOW-3D. The efficiency of energy dissipation was improved by approximately 5.8% with the notch pool than the other arrangement. Raza [25] investigated the impact of a stepped spillway slope in air entertainment and the location of the inception point using the volume of fluid (VOF) multi-phase model. The non-aerated length for the mild slope spillway was found to be larger than the steep slope stepped spillway. Ma [26] studied the interval pooled stepped spillway to evaluate the energy dissipation using a VOF multi-phase model. The effect of increasing pool height was found to be lower on energy dissipation. Morovati [27] carried out a detailed study of the design effects on the opening size of the pool step as well as step length using a VOF multi-phase model with RNG k-$ϵ$ turbulence model. A larger opening in the pool step was observed with less flow resistance. The flow resistance was enhanced with higher pool depth.

Li [28] applied a two-phase mixture model to study the effect of rounded and trapezoidal steps; higher air concentration was observed in rounded steps. The air–water velocity was found to be higher for trapezoidal steps than for rounded steps. The chamfering edges on steps obtained a slight energy dissipation enhancement and were considered less effective. Saqib [29] studied the influence of curved-tread steps regarding energy dissipation and pressure fluctuation using the volume of fluid (VOF) multi-phase method through FLOW-3D. The enhancement of energy dissipation was found up to 5–7% for low flows, whereas an unnoticeable difference was observed with high flow rates. Also, the curved treads enhanced the negative pressure. Overall, the effect of the curve treads was found to be an unfit design regarding the engineering aspect. Chen [30] carried out numerical simulation using the Realizable-k-$ϵ$ model, the tested influence of converging sectional width over the steps resulted in higher downstream flow velocity, inducing less energy dissipation; lower air concentration, causing a potential risk of cavitation; and higher water level, indicating the requirement of the higher side wall to avoid over-topping, which was shown to be an inappropriate engineering design.

Jahad [31] investigated the velocity and pressure variation due to the effect of an end sill on the steps of a spillway for several flow conditions using the VOF multi-phase model. Negative pressure was observed at the nappe flow regime for flat steps but not in curved steps. Velocity and pressure distribution were obtained with $\pm 6\%$ tolerance with an experimental data set. Kaouachi [32] studied the influence of six different step spillway widths for different flow regimes, using the multi-phase VOF along with the SST k-$\omega$ turbulence model. The wider steps were observed with alternating skimming flow due to the formation of asymmetrical vorticity patches in the step cavity.

The utilization of mildly sloped stepped spillways has increased due to their suitability for embankment dams. Nevertheless, enhancing the energy dissipation rate over such spillways with shorter chute lengths, particularly for small embankment dams, remains a significant challenge. Numerical investigations employing the ANSYS Fluent 2023 R1 solver have been conducted on identical mild slope spillways, exploring various spillway geometries to assess energy dissipation under extreme flow conditions.

The main objectives of the study of extreme flow events over different geometry of spillways are as follows.

• To assess the suitability of multiphase models to simulate flow over the stepped spillway;
• To evaluate the general flow behavior, velocity, and pressure distribution;
• To examine the distribution of turbulent kinetic energy (TKE) over the steps;
• To analyze the influence of step geometries on energy dissipation efficiency.

2. Methodology

For a better understanding and visualization of the complex flow behavior, as well as to capture the detailed effects, a combination of physical and numerical investigation is necessary. Both models have their strengths and weaknesses.

2.1. Multi-Phase Models

The application of the multi-phase model considers phenomena such as inter-phase energy and momentum transfer, phase distribution, phase transformation, and dynamics of the phase interface. Accurately capturing these types of effects is necessary to safely design and optimize the hydraulic structure. Different multi-phase models were tested to validate the numerical results with experimental data sets. The Volume of Fluid (VOF), Mixture, and Eulerian multi-phase flow models were adopted for model comparison during the simulation.

• VOF model

The VOF model shares a single set of momentum equations among the fluids. In each computational cell, the fluid volume of the fraction is tracked throughout the computational domain. It is primarily developed to track the free surface between the immiscible fluid, the flow behavior with free-surface flow, stratified flow, sloshing, and large bubble motion in the fluid. The volume fraction equation for tracking the interface phases is provided by the solution to the continuity equation, as represented below for qth phase [33]:

$$\frac{1}{{\rho }_q}\left[\frac{𝜕}{𝜕t}\left({\alpha }_q{\rho }_q\right)+\nabla ·\left({\alpha }_q{\rho }_q{\overrightarrow{v}}_q\right)=S_{{\alpha }_q}+\sum _{p=1}^n\left({\dot{m}}_{pq}-{\dot{m}}_{qp}\right)\right]$$

(1)

where ${\dot{m}}_{qp}$ and ${\dot{m}}_{qp}$ define phase mass transfer from qth to pth and vice versa, $S_{{\alpha }_q}$ is the additional mass source term, and ${\rho }_q$ is the density of qth phase,

$$\sum _{q=1}^n{\alpha }_q=1$$

where ${\alpha }_q$ is the volume of the fraction. The total volume of the fraction becomes 1 when the cell is full with the qth fluid.

• Mixture Model

The mixture model is a simplified Eulerian approach based on the assumption of the small Stokes number (St $<<$ 1). The mass average mixture’s velocity is solved through one continuity and one momentum mixture equation. The volume of the fraction transport is solved for each secondary phase. The mixture model is suitable for flow behaviors exhibiting dispersed flow, such as bubbly flow, droplets, and slurry flow. The following equation represents the continuity equation, volume transport fraction equation of each secondary phase, and momentum equation [33]:

Continuity equation:

$$\frac{𝜕}{𝜕t}\left({\rho }_m\right)+\nabla ·\left({\rho }_m{\overrightarrow{v}}_m\right)=0$$

(2)

Volume fraction transport equation of each secondary phase (pth phase):

$$\frac{𝜕}{𝜕t}\left({\alpha }_p{\rho }_p\right)+\nabla ·\left({\alpha }_p{\rho }_p{\overrightarrow{v}}_m\right)=-\nabla ·\left({\alpha }_p{\rho }_p{\overrightarrow{v}}_{dr,p}\right)+\sum _{q=1}^n\left({\dot{m}}_{qp}-{\dot{m}}_{pq}\right)$$

(3)

Momentum equation:

$$\frac{𝜕}{𝜕t}\left({\rho }_m{\overrightarrow{v}}_m\right)+\nabla ·\left({\rho }_m{\overrightarrow{v}}_m{\overrightarrow{v}}_m\right)=-\nabla p+\nabla ·\left[{\mu }_m\left(\nabla {\overrightarrow{v}}_m+\nabla {\overrightarrow{v}}_m^T\right)\right]+{\rho }_m\overrightarrow{g}+\overrightarrow{F}-\nabla ·\left(\sum _{k=1}^n{\alpha }_k{\rho }_k{\overrightarrow{v}}_{dr,k}{\overrightarrow{v}}_{dr,k}\right)$$

(4)

where n is the number of phases, ${\mu }_m$ is the viscosity of the mixture, $\overrightarrow{F}$ is the body force, ${\overrightarrow{v}}_{dr,k}$ is the drift velocity for the secondary phase (k) obtained from the volume fraction transport equation, and ${\overrightarrow{v}}_m$ is the mass average velocity.

• Eulerian Model

In the Eulerian multi-phase model, each phase is modeled separately with the interaction. The conservation equations for each phase contain a single phase term. The interfacial terms (i.e., lift, drag, mass transfer, etc.) are applied in the set of equations. For a model with complex flow behavior, the solution might face convergence issues due to the complexity of the model. The model is suitable for flow behavior such as bubbly flow, slurry flow, droplet flow, and particle-laden flow. The following equation represents the continuity and momentum equation of the Eulerian approach for the multi-phase flow model for qth phase [33]:

Continuity equation:

$$\frac{𝜕}{𝜕t}\left({\alpha }_q{\rho }_q\right)+\nabla ·\left({\alpha }_q{\rho }_q{\overrightarrow{v}}_q\right)=\sum _{p=1}^n\left({\dot{m}}_{pq}-{\dot{m}}_{qp}\right)+S_q$$

(5)

Momentum equation:

$$\frac{𝜕}{𝜕t}\left({\alpha }_q{\rho }_q{\overrightarrow{v}}_q\right)+\nabla ·\left({\alpha }_q{\rho }_q{\overrightarrow{v}}_q{\overrightarrow{v}}_q\right)=-{\alpha }_q\nabla p+\nabla ·{\overline{\overline{\tau }}}_q+{\alpha }_q{\rho }_q\overrightarrow{g}+\sum _{p=1}^n\left({\overrightarrow{R}}_{pq}+{\dot{m}}_{pq}{\overrightarrow{v}}_{pq}-{\dot{m}}_{qp}{\overrightarrow{v}}_{qp}\right)+\left({\overrightarrow{F}}_q+{\overrightarrow{F}}_{lift,q}+{\overrightarrow{F}}_{wl,q}+{\overrightarrow{F}}_{vm,q}+{\overrightarrow{F}}_{td,q}\right)$$

(6)

where ${\overline{\overline{\tau }}}_q$ is the strain–stress tensor, ${\overrightarrow{F}}_{lift,q}$ is the lift force, ${\overrightarrow{F}}_{wl,q}$ is the wall lubrication force, ${\overrightarrow{F}}_{vm,q}$ is the virtual mass force, ${\overrightarrow{F}}_{td,q}$ is the turbulent dissipation force for turbulent modeling, ${\overrightarrow{R}}_{pq}$ is the phase interaction force, and p is the pressure sheared by the phases.

2.2. Interface Tracking

The interface surface tracking methods (i.e., the sharp and dispersed interface) are essential to understand before implementing the multi-phase flow problem. The selection of interface techniques depends upon the type of flow behaviors. The Dispersed interface is likely to be suitable for highly disturbed flow, whereas the sharp interface is applied to the distinct free surface multi-phase flow [33], as shown in Figure 1. The sharp interface tracking technique is primarily developed for the application of the VOF method, where the volume fraction equation, representing the portion of the fluid on the particular cell, is solved for each cell, whereas the dispersed tracking technique is suitable for the Mixture and Eulerian multi-phase model.

2.3. Experimental Setup

The results of the numerical model were validated against the experimental results from the investigation of a large sectional model of a stepped spillway performed at the Utah Water Research Laboratory at Utah State University, USA. The experimental facility featured a stepped spillway with a series of 8 uniform steps, designed with 0.2 m step height, as represented in Figure 2a. The sectional model of the stepped spillway consists of a 1 m width and a mild overall slope of 1V:3H, which is mostly a typical slope of the embankment dam. Air–water flow data were recorded using a dual-tip phase detection conductivity probe attached to an automatic moving arm. The flow velocity and the concentration of air were estimated at each aerated step edge along the vertical line perpendicular to the direction of flow. The sampling frequency (F) of 200 kHz and a sampling duration (T) of 90 s were adopted using a LabVIEW interface. Data acquisition was performed at an interval of 5 mm along the vertical profile, considering the center-line axis of the stepped chute for two discharges of 0.425 m³/s/m and 0.565 m³/s/m. Air water flow properties such as void fraction (C), instantaneous velocities (V_inst), flow depths (y), and turbulence intensities (T_u) were then obtained by analyzing the conductivity probe data.

The conductivity probe consisted of a leading and a trailing tip with a longitudinal separation, $\Delta x$ = 5.24 mm, as shown in Figure 2b.

The implemented best practice time-averaged void fraction, the Krammer [34] relation, is

$$C=\frac{1}{T}{\int }_{t=0}^T{c}_1\left(t\right)dt$$

(7)

where C is the time average void fraction, T is the sampling duration, and $c_1\left(t\right)$ is the instantaneous void fraction.

The velocity was estimated according to the ratio of the distance between tips ($\Delta x$), with the interfacial travel time obtained from the cross-correlation of the leading tip and trailing tip signals using the following relations, Krammer [34]:

$$R_{12,i}\left(\tau \right)=\frac{{\sum }_{t=t_{s,i}}^{t_{s,i}+{\mathcal{W}}_{T,i}}\left(\left[S_1\left(t\right)-〈S_{1,i}〉\right]·\left[S_2(t+\tau )-〈S_{2,i}〉\right]\right)}{\sqrt{{\sum }_{t=t_{s,i}}^{t_{s,i}+{\mathcal{W}}_{T,i}}{\left(S_1\left(t\right)-〈S_{1,i}〉\right)}^2}·\sqrt{{\sum }_{t=t_{s,i}}^{t_{s,i}+{\mathcal{W}}_{T,i}}{\left(S_2(t+\tau )-〈S_{2,i}〉\right)}^2}}$$

(8)

$${\left[\langle V\rangle \right]}_{t_{s,i}}^{t_{s,i}+{\mathcal{W}}_{T,i}}\approx V_i=\frac{\Delta x}{{\mathcal{T}}_i}$$

(9)

where $U_i$ is the instantaneous velocity (m/s), ${\mathcal{T}}_i$ is the interfacial travel time indicating the peak of $R_{12},i$ (${\mathcal{T}}_i\left(=argmax\left(R_{12,i}\right)\right)$), $S_1$ and $S_2$ are the tip signal of the leading and tailing (V), $W_{T,i}$ is the time window duration (s), $\tau$ is the time lag (s), and $\Delta x$ is the distance between the tip (m).

2.4. Numerical Schemes

The numerical schemes of the spatial and temporal discretization adopted for the model are listed in tabular form in

Table 1. The model’s implemented numerical schemes.

Multi-phase model	Mixture-Dispersed Interface
Turbulence model	Realizable-k-$ϵ$
Pressure–velocity coupling	Coupled
Gradient	Least squares cell based
Pressure	Body force weighted
Momentum	First order upwind
Volume fraction	First order upwind
Turbulent kinetic energy	First order upwind
Turbulent dissipation rate	First order upwind
Transient formulation	First order implicit

. The ”Coupled” scheme was implemented for pressure–velocity coupling, which solves the continuity and momentum equation simultaneously in a coupled way. The coupled schemes also improve the convergence by solving pressure–velocity fields together, reducing the number of iterations. The ”Least squares cell-based” scheme was used during the simulation to evaluate the cell center gradient of a variable to minimize the sum of squares of the errors. It helps to increase the accuracy of the gradient for a complex geometry mesh. For the discretization of the pressure gradient term, a ”Body force weighted” scheme was used as flow behavior over the stepped spillway, which has a significant impact on gravity flow. To solve the momentum equation, the volume fraction equation, the turbulent kinetic energy equation, and the turbulent dissipation rate, the ”First-order upwind” schemes were used to reduce the computational cost and ensure the stability of the model, as higher-order numerical schemes often result in model instability issues. Also, the ”First-order implicit” scheme was adopted for time discretization in the transient formulation [33].

The realizable k-$ϵ$ turbulence model was implemented to accurately capture the effects of the flow recirculation, massive separation, and boundary layer separation that occurs in the steps of the spillway. Details about the selection of the appropriate model, depending upon the flow behaviors, are discussed in Shaheed [33,35].

Mesh Quality and Grid-Independent Test

The accuracy of the results is often associated with the quality of the mesh and grid size; a low-quality mesh might distort the elements affecting the solution. On the other hand, grid size might also impact the solution concerning the water surface elevation, the velocity distribution, and the air concentration. The mesh quality was improved by evaluating the skewness, orthogonality, and aspect ratio, as presented in

Table 2. Mesh quality assessment.

Qualitiy	Model	Acceptable Threshold [33]
Skewness	0.01	<0.8
Orthogonality	1 (90°)	>0.6
Aspect ratio	1.26	<3

. The structured mesh was implemented using ANSYS ICEM-CFD to avoid these effects and improve the mesh quality. A typical mesh arrangement of a stepped spillway is presented in Figure 3 with inflation layers.

The grid-independent test was carried out by considering different mesh sizes of 0.004 m, 0.005 m, 0.008 m, and 0.01 m with a fixed time step of 0.01 s for the flow rate of 0.425 m³/s/m. The flow velocity at spillway step 7 and the water surface elevation were evaluated with experimental results through different mesh sizes. The best practice phase technique of ISO-surface generation through volume fraction with a fractional value of 0.5 was implemented to calibrate the model so as to obtain the expected water surface elevation. The water surface elevations measured at the consistent difference of 0.2 m over a crest were found to be more or less similar among all the meshes, as presented in Figure 4.

The Grid Convergence Index (GCI) concept suggested by Roache [36] is often considered an essential tool in CFD for identifying convergence of the numerical solution between the fine, medium, and coarse grid sizes. Therefore, a fine grid size of 0.005 m, a medium grid size of 0.008 m, and a coarse grid size of 0.01 m were considered for the grid convergence test. The simulated depth average velocity listed in

Table 3. Estimated velocity for different grid sizes.

Experimental Velocity (m/s)		4.524
Grid Size (m)	CFD-Velocity (m/s)	Remarks
0.005	4.513	Fine
0.008	4.573	Medium
0.01	4.691	Coarse

for the fine, medium, and coarse grid sizes was tested against the experimental results for step 7 of the spillway for a flow rate of 0.425 m³/s/m using the following relation [36]:

$$\begin{array}{ccc}\hfill & GC{I}_{\mathrm{Fine}}=\frac{1.25\left|\frac{{\alpha }_F-{\alpha }_M}{{\alpha }_F}\right|}{{\left(\frac{\Delta x_M}{\Delta x_F}\right)}^p-1}\hfill & \hfill p=\frac{ln\left[\frac{{\alpha }_C-{\alpha }_M}{{\alpha }_M-{\alpha }_F}\right]}{ln\left[\frac{\Delta x_M}{\Delta x_F}\right]}\end{array}$$

(10)

where, $\Delta x_C$, $\Delta x_M$, and $\Delta x_F$ are the grid sizes of coarse, medium, and fine. ${\alpha }_C$, ${\alpha }_M$, and ${\alpha }_F$ are the parameter values to be tested for coarse, medium, and fine. In this case, depth average velocity was considered at step 7.

Asymptotic range of convergence relation:

$$\mathrm{Asymptotic}\mathrm{range}=\frac{GC{I}_{\mathrm{coarse}}}{GC{I}_{\mathrm{Fine}}{\left(\frac{\Delta x_M}{\Delta x_F}\right)}^{p_{\mathrm{Fine}}}}$$

(11)

The GCI value was calculated as 1.7%, whereas the absolute relative error was calculated as 0.24%, considering the fine mesh and expected value of the experimental data set. Also, the asymptotic range of convergence was found to be 0.986, nearly 1, which defines the model results as acceptable. The plot of the solution convergence using different mesh sizes has been illustrated in Figure 5. The result shows a converged solution below the mesh size 0.005 m.

2.5. Multi-Phase Model Comparison

The model comparison between the physical laboratory test and different multi-phase models is presented in Figure 6, for a flow rate of 0.425 m³/s/m. Skimming flow was observed in both the physical and numerical models. The length of the inception point was more or less adequately captured by Mixture and Eularien–Dispersed interface than the VOF–sharp interface model, with the experimentall flow result represented in Figure 6a,b. The length was estimated based on the air–volume fraction over the steps of the spillway. To maintain consistency, void fraction values ranging from 0.09 to 0.18 at the step edges were deemed indicative of the inception point for the results of the numerical simulations. This approach was also supported by the observations of the TKE, which are expected to increase with the increase in air and subsequent mixing. For instance, the initiation of full-flow excitation due to the rapid increase in TKE at step 4 indicates the beginning of an aerated zone, as presented in Figure 6b. The dispersed interface model shows the suitability of strong mixing behaviors of air and water over the steps of the spillway. However, the numerical model was unable to produce the dispersed bubbly flow as it is in laboratory results, which can be considered a limitation of the volume fraction transport model.

Additionally, the air concentration along the flow depth was tested against the different multi-phase models at step 6 for a flow rate of 0.425 m³/s/m. Air concentrations were plotted against normalized flow depths (y/$y_{90}$), as shown in Figure 7, where y is the flow depth and $y_{90}$ is the depth at 90% air concentration. The concentration of air using the VOF–sharp interface highly deviated from the experimental air concentration due to the absence of air concentration until the water depth ratio ($y/y_{90}$) around 0.65, as a result of distinct water flow, which is clearly illustrated by reference Figure 6b. On the other hand, the Mixture and Eulerian models were able to obtain the air concentration on the water column. However, the model still struggles to adequately capture the concentration of the air in the water column compared to the experimental result. Several studies have been carried out to attain closer results to the experiment through different CFD volume fraction techniques, and replicating the results with the volume fraction method was found to be a challenging task, according to Valero [37], Almeland [38], and Li [39]. A simplified Mixture multi-phase model was adopted based on the optimum results that were closer to the Eulerian results and the experiment results to save on computational cost, as the full Eulerian equation increases the computational cost of the model for multiple test run case studies.

Model Validation

To enhance the accuracy and robustness of the numerical model, proper model validation for different flow parameters is necessary before implementing several test cases. Therefore, the function of dimensionless velocity (V/$V_{90}$) versus depth ($y/y_{90}$) and water surface elevation was evaluated at different steps of the spillway, especially in the aerated regions at the downstream end of the spillway, as shown in Figure 8. The two different independent flow rates of 0.425 m³/s/m and 0.565 m³/s/m were considered with the adopted mesh size of 0.005 m.

• Water Surface Profile

The water surface profiles for a flow rate of 0.425 m³/s/m and 0.565 m³/s/m were implemented to validate the model with the experimental results (water surface elevations measured at a consistent difference of 0.2 m over a crest, as referenced in Figure 4), as represented in Figure 9. The model adequately captures the profile of the water surface with the adopted mesh size of 0.005 m for different flow rates. The empirical relation of brink depth ($y_b$) at the end of the spillway crest also demonstrates the model’s efficiency, which was all found with negligible differences, as shown in

Table 4. Brink depth ($y_b$) comparison between the models and the empirical relation.

Brink Depth (${\mathit{y}}_{\mathit{b}}$)
Sp. Discharge (m3/s/m)	Experimental	CFD	Emperical
0.425	0.19	0.178	0.188
0.565	0.22	0.216	0.228

.

The empirical relation of brink depth ($y_b$) (Rouse, 1936) [40] for the rectangular channel is precisely estimate by

$$y_{\mathrm{b}}=0.715{\left(\frac{q^2}{g}\right)}^{\frac{1}{3}}$$

(12)

• Velocity Distribution

The velocity profiles obtained from the dual-tip conductivity probe at spillway steps 6 and 7 in the highly aerated water column were evaluated with the numerical model results for the flow rate of 0.425 m³/s/m and 0.565 m³/s/m. The dimensionless function of velocity (V/$V_{90}$) and water depth (y/$y_{90}$) were evaluated regarding the velocity plots, as illustrated in Figure 10 and Figure 11.

The model was able to reproduce the velocity distribution on the water column with adequate results in comparison to the experimental data set by following a similar distribution, as illustrated in Figure 10 and Figure 11. Furthermore, the depth average velocity was also evaluated for the adopted flow rate of 0.425 and 0.565 m³/s/m at steps 6 and 7 using the following empirical relation (Chanson [41]):

$$\begin{array}{c}{y}_{clear}={\int }_{y=0}^{Y_{90}}(1-C)\times dy\end{array}$$

(13)

$$\begin{array}{c}V=q/{\int }_0^{Y_{90}}(1-C)\times dy\end{array}$$

(14)

where $y_{clear}$ is the clear water depth $\left(\mathrm{m}\right)$, q is the specific discharge $({\mathrm{m}}^3/\mathrm{s}/\mathrm{m})$, V is the depth average velocity $(\mathrm{m}/\mathrm{s})$, and C is the air concentration.

The evaluated depth average velocity using Equations (13) and (14) is presented in

Table 5. Depth average velocity at steps 6 and 7.

Depth Average Velocity (m/s)
Sp. Discharge (q) (m3/s/m)	Spillway Steps	CFD	Experiment
0.425	6	4.11	4.14
	7	4.51	4.52
0.565	6	4.53	4.55
	7	4.59	4.61

. The model’s results were found to be more or less similar to the experimental values for flow rates of 0.425 and 0.565 m³/s/m over steps 6 and 7.

2.6. Run Cases

The different step geometries were adopted over a mild slope spillway, with a horizontal distance of 4.8 m and a vertical distance of 1.8 m from the crest. The horizontal lengths for small steps uniform and non-uniform steps were made up to 5.1 m to maintain the overall slope of the spillway, as shown in Figure 12. The uniform steps with different step heights, non-uniform steps, 8° adverse slope steps, and stepped pool were tested against the eight different flow rates of 0.285, 0.425, 0.565, 0.705, 0.845, 0.985, 1.125, 1.265 m³/s/m. The low flow rate of 0.285 m³/s/m was increased by almost 4.5 times with a discharge until 1.265 m³/s/m, concerning the risk of high flow effects over a spillway. The energy dissipation efficiency of the different step geometry was tested against low and high flows for the safe operation of the spillway over a mild slope.

3. Results and Discussion

The flow behaviors in terms of velocity, pressure distribution, cavitation index, turbulent kinetic energy, and energy loss due to different step geometries are fundamental in characterizing the flow in depth. Therefore, a detailed study was carried out to identify the efficiency of the different step geometry over a mild slope for safe spillway operation during low and high flows.

3.1. Flow Behavior

The sudden change in the water surface streamline–flow profile along the spillway due to steps causes highly disturbed flow. The low-pressure zones are created behind the steps as a result of changes in flow direction and velocity due to flow separation, which leads to the formation of recirculation zones. The recirculating eddies were observed during the numerical simulation, as expected, which can be visualized in Figure 13b. The recirculation plays a key role in energy dissipation to reduce the flow directional velocity.

The flow behaviors over different step geometries of the spillway were studied to identify the location of the inception point with incremental discharges. The inception point separates the regions between non-aerated and aerated flow, as illustrated in Figure 8. The elongated length of the inception point increases the risk of cavitation damage due to a less aerated flow region over a spillway. The impact of the increased discharge was evaluated for 0.425 and 0.565 m³/s/m over different geometries, as shown in Figure 14.

The flows exhibit a skimming flow type for all the geometries, as expected for critical depth; step height ratio ($\frac{y_c}{h}\ge 0.8$) and step height; and length ratio ($0.4\le \frac{h}{l}\le 0.9$), Rajaratnam [42]. The location of the inception point was found to be shifted downwards as discharge increased from 0.425 to 0.565 m³/s/m for all the geometries, which can be visualized in Figure 14. This indicates the elongation of the non-aerated region and increases the risk of cavitation during high flows over stepped spillways. The influence of different geometry over the length of the inception point was found to be lower, considering an equal overall slope geometry (reference Figure 12), as the length of the inception point was found to be more or less similar, which can be observed in Figure 14.

3.2. Velocity and Pressure Distribution

The average depth velocity was evaluated at the downstream end of the spillway. The eight different flow rates were tested against the adopted geometry. The impact of the geometry roughness was analyzed against the downstream depth average velocity, as the reduction of velocity helps to reduce the design length of the stilling basin, increasing the energy dissipation rate. The flow velocity against the sp. discharge was plotted for a typical section of the 0.2 m H uniform steps spillway, as shown in Figure 15, considering low and high flows. Lower effectiveness of the geometry step roughness height (k, reference Figure 13a) was observed for high flows due to higher flow depth, as the flow–velocity curve behaves non-linearly. This implies that the relative roughness height $(k/y)$ becomes lower as flow depth $\left(y\right)$ increases due to higher flow rates, ultimately decreasing the energy dissipation efficiency of the spillway.

The ratio of the theoretical velocity $\left(V_{max}\right)$ assumed with smooth surface and simulated velocity $\left(V_1\right)$ due to step roughness at the downstream end, in relation to upstream total head $\left(H_o\right)$ and brink depth $\left(y_b\right)$ ratio, was plotted for different flow rates, as presented in Figure 16. The implemented energy relation for the plot is

$$H_o=y_b+\frac{v^2}{2g}+\Delta Z$$

(15)

$$V_{\max }=\sqrt{2\times 9.81\times (H_o-y\times cos\left(\theta \right))}$$

(16)

A noticeable reduction in the velocity ratio $(V_1/V_{max})$ was observed with a 0.2 m H stepped spillway, indicating higher velocity reduction over a mild slope among the different geometrical steps. The stepped pool spillway configuration was found to be less effective in terms of velocity reduction as a result of an increased velocity ratio due to higher flow velocity at the downstream end of the spillway. The stepped pool spillway acts as a “stagnant pool”, as shown in Figure 17a, which results in less reduction in the flow velocity, thereby reducing the efficiency. The sequential difference in velocity ratio $(V_1/V_{max})$ is reduced with increasing flow rates for all the adopted geometries, due to less effectiveness of the roughness, for higher flows in the spillway. Increasing the uniform roughness over the spillway span with 0.1 m H steps was found to be less effective than with 0.2 m H steps, as illustrated in Figure 16 concerning the velocity ratio $(V_1/V_{max})$. Furthermore, adverse slope steps and non-uniform steps were found to be more effective than the 0.1 m H step and the stepped pool configuration. However, the result also shows a more or less similar reduction in the velocity ratio during the low flow for $(H_0/y_b)\ge 10$ as a result of less difference in the downstream flow velocity. This implies the overall roughness impact over the steps during the low flow depth results in a more or less similar reduction in the downstream flow velocity for adopted geometries, except for stepped pool configuration.

The variation in the pressure distribution over the steps of the spillway for different geometries with two different flows was evaluated. The flow incremental low-pressure zones identified in the case of 0.285 m³/s/m and 0.425 m³/s/mare presented in Figure 18. Negative pressure was observed at the vertical steps of the spillway for all the configurations, indicating the risk of cavitation at the vertical face of the spillway. The low-pressure zones in front of the steps were created due to the flow separation as a result of recirculating secondary flow. The static pressure contour plot shows the incremental pressure with the increase in flow rate for all the geometries, as expected. The minimum values of the negative pressure were observed for the low flow (0.285 m³/s/m) for all the configurations. This implies a risk of cavitation for low flow over the steps spillway due to the localized low-pressure zones with the presence of cavities. A higher negative value of pressure (−2106.54 Pa) was observed for the adverse slope steps among the other configurations, which revealed that the design configuration might be inappropriate in the case of cavitation risk. A minimum negative pressure value (−802.24 Pa) was observed for a uniform 0.1 m H; 0.3 mL step spillway, with a lower impact concerning the risk for cavitation.

Furthermore, the cavitation index was evaluated for the safety assessment of the spillway, as presented in

Table 6. Cavitation index for different geometry configurations for 0.285 m³/s/m.

Geometry Type	Static Pressure (Pa)	Velocity (m/s)	Cavitation Index
Adverse slope steps	−2106.54	4.51	9.51
Non-uniform steps	−2028.65	4.44	9.82
Stepped pool	−1955.6	4.62	9.08
0.2 m H steps	−1301.05	4.39	10.12
0.1 m H steps	−802.24	4.4	10.13

. The permissible safe cavitation damage index is greater than 0.2 for concrete, as per Falvey (1990) [6]. The cavitation index was estimated based on the 2.33 kPa value of the saturated vapor pressure of water at 20 °C, the density of 998.2 kg/m³, and the minimum static pressure obtained during the simulation through 0.285 m³/s/m. The cavitation index for different geometries with associated flow velocity was calculated using the following relation [6]:

$$\mathrm{Cavitation}\mathrm{index}=\frac{P_o-P_v}{\frac{1}{2}{\rho }_{water}{V}^2}$$

(17)

$$P_o=P_a+P_g$$

(18)

where $P_o$ is the reference pressure (Pa), $P_g$ is the gauge static pressure (Pa), $P_a$ is the atmospheric pressure (Pa), $P_v$ is the vapor pressure (Pa), V is the velocity (m/s), and ${\rho }_{water}$ is the density of water (kg/m³).

Uniform steps of 0.2 m H and 0.1 m H were found with a higher cavitation index value, indicating optimal design for cavitation risk. However, a lower chance of cavitation damage was identified for all the geometries. The cavitation index was much higher than the critical cavitation index value of 0.2 due to the negligible negative pressure values compared to the atmospheric pressure of 101.0 kPa. Nevertheless, the assurance of cavitation damage risk for small physical models might impact the prototype. The scaling up of pressure and velocity values reduces the cavitation index, ultimately increasing the risk of cavitation, as the cavitation index is highly dominated by higher flow velocities. The cavitation damage due to high-pressure values with higher flow velocity in Pedrógão and Dona Francisca Dam’s stepped spillway during the operation is a typical example of potential damage [43].

3.3. Turbulent Kinetic Energy (TKE)

The study of turbulence kinetic energy ensures structural integrity and efficiency of the hydraulic structures through detailed flow impact assessment of higher flow rates over different geometry spillways. The turbulent kinetic energy per unit mass can be expressed in the form of fluctuating mean square velocities as follows [33]:

$$TKE=\frac{1}{2}\left(\overline{u^{\prime 2}}+\overline{v^{\prime 2}}+\overline{w^{\prime 2}}\right)$$

(19)

where TKE is the turbulent kinetic energy per unit mass $({\mathrm{m}}^2/{\mathrm{s}}^2)$ and $\overline{u^{\prime 2}}$,$\overline{v^{\prime 2}}$, and $\overline{w^{\prime 2}}$ are the mean square velocities of the fluctuating component over x, y, and z directions.

The turbulence kinetic energy indicates the strength of the turbulence and vigorous mixing of the flow, which helps to enhance the energy dissipation and cures the potential damage, especially at the downstream end of the dam in the stilling basin. The turbulence strength was tested considering the flow rate of 0.565 m³/s/m for different geometries through the TKE contour plot, as presented in Figure 19. The higher TKE of 1.42 m²/s² was obtained through the 0.2 m H steps over the other geometrical configurations, which indicates higher chances of energy dissipation rate over the steps of the spillway. Higher TKE values generate vigorous mixing of the flow, which leads to a reduction in the directional flow velocity over every step. The lowest TKE of 1.05 m²/s² was obtained through the stepped pool spillway due to the formation of a “stagnant pool” over the steps. The lower TKE values were observed over the steps of the stepped pool spillway, which indicated less effectiveness concerning the energy dissipation over the steps.

3.4. Energy Dissipation

The dissipation of kinetic energy is essential to avoid excessive erosion at the downstream end of the spillway over the stilling basin. The adopted energy loss relation is presented in Figure 8. The total energy loss over the steps of the spillway was plotted in relation to the loss ratio ($\Delta$H/$H_o$) over different flow rates (q), considering the adopted geometrical configuration, to evaluate the energy efficiency, as presented in Figure 20a. Almost 60% of energy loss was achieved during the low flow of 0.285 m³/s/m by 0.2 m H steps, 0.1 m H steps, and non-uniform steps. Although the adverse step slope and stepped pool configurations achieved less than 50% of the energy, the adverse slope steps and stepped pool steps helped to retain water, causing a smooth flow transition between the steps, reducing the energy-dissipation efficiency. For the higher flow rates, the results clearly indicate that the 0.2 m H step was the most effective one for energy dissipation and is suitable for design consideration.

Furthermore, the residual head ($H_1$) with different flow rates (q) was also plotted for better visualization, as shown in Figure 20b. More or less similar residual heads were observed for low flows. However, a higher reduction in the residual head was observed for the uniform 0.2 m H step spillway than the other geometrical configurations during a high flow. The higher residual head observed for the stepped pool configuration implies that the design configuration is not effective for most of the flows in the spillway. A higher risk of potential damage for the stilling basin was observed with the stepped pool spillway by introducing higher flow velocity at the stilling basin, whereas the least risk was identified by the uniform steps of the 0.2 m H steps spillway as a result of lower residual head. This revealed that a higher drop height (H) and length (L) of uniform steps have better performance over a mild slope during high flows than the other configurations, as it generates a high TKE, which reduces the potential risk of damage in the stilling basin. Furthermore, increasing the number of uniform steps with 0.1 m H and the stepped pool was found to be less effective in comparison to the adverse and non-uniform steps as a result of the higher residual head.

4. Conclusions

The study carried out on the numerical investigation of the safe operation of the mild slope spillway with different geometrical test configurations resulted in the following conclusions:

• The mixture dispersed interface model was found to be suitable for the highly disturbed flow over the steps of the spillway, as it adequately represents the physical behavior of the multi-phase flow. While the VOF method was unable to replicate the self-entrainment of air into the flow, the mixture models successfully demonstrated the multi-phase flow behavior observed in the stepped spillway flow.
• The shift in the location of the inception point is mainly dependent upon the discharge for a more or less identical slope geometry. The inception point was found to be shifted downstream due to a higher flow rate, which increases the non-aerated region. The impact of change in geometry over an identical slope showed less influence in the location of the inception point.
• The downstream flow velocity curve has shown non-linear behavior, as relative roughness was found to be lower with sequential increments in flow rates. This indicates that the impact of roughness will be less for higher flow depth over a stepped spillway, which ultimately decreases the energy dissipation efficiency.
• The lowest velocity ratio $(V_1/V_{max})$ with increasing flow rates was obtained with uniform 0.2 m H steps compared to the rest of the configuration, which shows the higher rate of velocity reduction. A higher velocity ratio $(V_1/V_{max})$ was obtained using the stepped pool design, which illustrates the lower rate of velocity reduction, increasing the downstream velocity as a result of the formation of a “stagnant pool”.
• Higher negative pressure values were observed for adverse slope steps due to the low-pressure zone created near the vertical face of the steps as a result of flow separation. Although negative pressure was observed, the impact of cavitation was found to be lower as a result of low negative values compared to the atmospheric pressure, indicating less risk of cavitation damage for the stepped spillway for the model test.
• Turbulent kinetic energy (TKE) for uniform 0.2 m H steps was found to be higher, indicating a strong mixing of flow over the steps than the other configuration, which helps to reduce the direction flow velocity. Lower TKE was observed for the stepped pool spillway due to the formation of a “stagnant pool”, which reduces the mixing of the flow.
• The maximum efficiency ($\Delta$H/$H_o$) of the energy dissipation was obtained by 0.2 m H steps for lower and higher flow rates over the spillway due to higher TKE over the steps. Almost 60% was achieved during low flows. The lower value of the residual head ($H_1$) also indicated the higher efficiency of the uniform steps 0.2 m H steps. The impact of incremental numbers of uniform 0.1 m H step roughness shows less reduction in residual heads, indicating lower effects for low drops. Furthermore, the stepped pool and uniform 0.1 m H steps were found to be less effective than adverse slope steps and non-uniform step configurations.