Introduction

Energy dissipator structures are used to mitigate the damage that could potentially happen downstream from energy of the supercritical flow and extreme kinetic energy in such a flow. Different types of these structures can be used in the ocean to provide calm areas or downstream of dams in rivers. These structures can control the hydraulic jump by creating a jump in a designated location. Among energy dissipators, the stilling basin downstream from the dam’s chute is one of the most common structures.

So far, extensive experimental and numerical research have been carried out on hydraulic jumps and energy dissipation in stilling basins. Wu and Rajaratnam (1996) examined the flow transitions from the hydraulic jump to fully developed flow in open channels. The velocity profile, at this distance, changes from a tangential jet profile to a semi-logarithmic profile. Chanson and Brattberg (2000) conducted new experiments on developing a shear layer in a hydraulic jump with a semi-developed upstream flow. In these experiments, the air concentration distribution, the average velocity of water-to-air flow, and air bubble frequencies were measured. Wang and Liu (2000) compared different numerical methods for two-dimensional shallow water, a 2D-dam break, and a hydraulic jump. They used four different finite volume methods on unstructured triangular networks called Roe-MUSCL, the Roe-Upwind, the Hll-MUSCl, and the hybrid methods. The results of the numerical solutions, the computational speed rate, and the stability of the methods were compared together. Ead and Rajaratnam (2002) obtained the hydraulic jump parameters considering bed roughness. They reported that the required tailwater depth for a hydraulic jump on the conical floor is less than for a flat floor. Zhao and Misra (2004) simulated the turbulent hydraulic jump by numerical methods. They used the VOF method to investigate jump turbulence. Torabi et al. (2019) used VOF method to simulate flow to consider erosion and sediment in channel junction. Ni and Liu (2005) conducted a study on the depth ratio, energy loss ratio, and jump length. Madsen et al. (2005) developed the implicit formulation of nonlinear finite-difference in shallow water for application in tidal wave behavior and hydraulic jumps. Yan and Zhou (2006) carried out a statistical analysis of the pressure fluctuations under spatial hydraulic jumps. Shafai Bejestani and Neisi (2009) mentioned that by choosing a specific bed roughness, the length of the rectangular stilling basin could be reduced 40%.

Despite numerous numerical studies on energy dissipators, this study can determine the effect of grid size on the simulated flow profile for different step heights. In this research, a numerical model for the stilling basin of the Upper Siah Bisheh Dam is created, then calibrated with results obtained from the physical model. In this model, the effect of the step’s height modification at the end of the chute on the flow profile and the flow velocity, in the stilling basin, is examined. The results indicate that increasing step height raises flow depth and decreases flow velocity. Decreasing step height generates the opposite situation. In the next section, a numerical approach is introduced to determine the optimum grid size by performing sensitivity analysis. The results show a significant amount of time can be saved by optimizing grid size.

Methodology and governing equations

Physical model

The numerical model has been created based on the stilling basin of the real Siah Bisheh Dam information. The physical model of this dam is located in the Hydraulic Laboratory of the Iranian Water Research Institute. Siah Bisheh is an earthen dam located 10 km north of the Kandovan Tunnel (in the Chalus Road) on the Chalus River; it has a height of 85 m. At the end of the chute, there are some steps that are 0.5 to 0.7 m in height. The stilling basin terminated at a trapezoidal channel. The length of the stilling basin is 22.90 m in the floor and 24.90 m, including the steps, at the end of the stilling basin. To verify, optimize, and ensure the operation of the flood discharge system, a hydraulic model with a 1:15 scale has been designed and tested. This model is made of Plexiglas. (Water Resources Research Institute 2005).

In the real dam, three steps with a height of 0.6 m and a length of 1 m were installed to dissipate energy at the end of the stilling basin that leads to the channel. The rest of the basin features and the longitudinal section of the real dam’s stilling basin is shown in Fig. 1 (Water Resources Research Institute 2005). In this research, three heights (0.5, 0.6 and, 0.7 m) of the last step of the chute were modeled numerically to determine the effect of step height on the flow profile in the stilling basin. The original step height (last step of the chute) in the dam is 0.5 m.

Fig. 1
figure 1

Longitudinal section of the studied basin

Full size image

Numerical model

In this research, the computational fluid dynamics (CFD) package is used to simulate flow on a stilling basin. The applied numerical method is based on solving the Reynolds-averaged Navier–Stokes (RANS) equations with the finite-volume method. This method can solve flows with a strong vortex and naturally unsteady flows with reasonable cost. Detailed information on the derivation of the RANS equations that were used in this study is obtained from Daneshfaraz and Ghaderi (2016), Daneshfaraz et al. (2017) and Zeidi and Mahdi (2014,2015). The continuity equation (principle of conservation of mass) and the momentum equations are presented in Eqs. (1) and (2), respectively:

$$\frac{\partial {U}_{i}}{\partial {x}_{i}}=0$$
(1)
$$\frac{\partial {U}_{j}}{\partial t}+{u}_{i}\frac{\partial {U}_{j}}{\partial {x}_{i}}=-\frac{1}{\rho }\frac{\partial P}{\partial {x}_{j}}-\frac{1}{\rho }\frac{\partial {\tau }_{ij}}{\partial {x}_{i}}+{g}_{j}$$
(2)

where \({u}_{i}\) is the momentum component in the direction of xi, \(\nu\): is kinematic viscosity \(\left(\frac{\mu }{\rho }\right)\), \(\rho\): is the volume mass of fluid, gxi: is the gravity acceleration components in the direction of xi, P: is the pressure at each point of the fluid, i: varies from 1 to 3 and represents the direction of the axis x, y, z, respectively.

To solve the Navier–Stokes equations, the Reynolds averaged method is used. Therefore, turbulent flow fluctuations are indirectly inserted into equations. In this case, the modified forms of the Reynolds-averaged Navier–Stokes equations (RANS) are presented in Eqs. (3) and (4), respectively.

$$\frac{\partial }{\partial x}\left(u{A}_{x}\right)+\frac{\partial }{\partial y}\left(\nu {A}_{y}\right)+\frac{\partial }{\partial z}\left(w{A}_{z}\right)=0$$
(3)
$$\frac{{\partial u}_{i}}{\partial T}+\frac{1}{{V}_{F}}\left({U}_{j}{A}_{j}\frac{{\partial u}_{i}}{\partial {x}_{j}}\right)=\frac{1}{\rho }\frac{\partial p}{\partial {x}_{j}}+{g}_{i}+{f}_{i}$$
(4)

where Ax, Ay, and Az are the fractional areas open to the flow in the x, y, z directions; u, v, and w are the velocities in the directions of x, y, and z, respectively; T is time and VF is fractional volume open to the flow; Uj and Aj are velocity and the face area of the cell, respectively; ρ is density of fluid; p is pressure; gi is gravitational force in i-direction; and fi represents the Reynolds stress to close the turbulence model (Sangsefidi et al. 2019).

Reynolds stress models or RSM

In the Reynolds stress model (RSM) model, a transition equation is used to calculate each of the Reynolds stresses. The Reynolds stress tensor is a symmetrical tensor, which means, that in the two-dimensional model (Eq. 5) only three transformation equations and in the three-dimensional state (Eq. 6), only six equations are required to determine the distribution Reynolds stress (Fluent Inc 2006).

$$\left[ {\begin{array}{*{20}c} {\overline{{u^{{\prime^{2} }} }} } & {\overline{{u^{\prime } v^{\prime } }} } \\ {\overline{{v^{\prime } u^{\prime } }} } & {\overline{{v^{{\prime^{2} }} }} } \\ \end{array} } \right]$$
(5)
$$\left[ {\begin{array}{*{20}l} {\overline{{u^{{\prime^{2} }} }} } \hfill & {\overline{{u^{\prime } v^{\prime } }} } \hfill & {\overline{{u^{\prime } w^{\prime } }} } \hfill \\ {\overline{{v^{\prime } u^{\prime } }} } \hfill & {\overline{{v^{{\prime^{2} }} }} } \hfill & {\overline{{v^{\prime } w^{\prime } }} } \hfill \\ {\overline{{w^{\prime } u^{\prime } }} } \hfill & {\overline{{w^{\prime } v^{\prime } }} } \hfill & {\overline{{w^{{\prime^{2} }} }} } \hfill \\ \end{array} } \right]$$
(6)

Two-dimensional mode

Determining the time steps in a fluent software requires numerical modeling convergence. The time step is used to adjust the grid dimensions to achieve sustainability. The discrete differential equation for the non-permanent fluid flow is seen in Eq. (7) (Patankar 1980). Figure 2 shows the control volume for the two-dimensional situation.

$$\rho_{c} \frac{\partial T}{{\partial t}} = \frac{\partial }{\partial x}\left( {k\frac{\partial T}{{\partial x}}} \right) + \frac{\partial }{\partial y}\left( {k\frac{\partial T}{{\partial y}}} \right) + S$$
(7)
Fig. 2
figure 2

Control volume for the two-dimensional situation

Full size image

There are different iterative methods to solve algebraic equations:

$$a_{P} T_{P} = a_{E} T_{E} + a_{W} T_{W} + a_{N} T_{N} + a_{S} T_{S} + a_{T} T_{T} + a_{B} T_{B} + B$$
(8)

The point-to-point Gaussian–Seidel method has been used throughout this work. In this method, variable values for each node are calculated according to the specified network order. The breakdown equation will be written as Eq. (9).

$$a_{P} T_{P} = \sum {a_{{{\text{nb}}}} T_{{{\text{nb}}}} + b}$$
(9)

where nb implies a neighbor's node, then \({T}_{P}\) will be calculated using the following Equation.

$$T_{P} = \frac{{\sum {a_{{{\text{nb}}}} T_{{{\text{nb}}}}^{*} + b} }}{{a_{P} }}$$
(10)

In Eq. (10), \({T}_{\mathrm{n}\mathrm{b}}^{*}\) indicates the values of the neighboring node. For the neighbors that have been considered in the current iteration, the value of \({T}_{\mathrm{n}\mathrm{b}}^{*}\) is for previous iterations. In any case, \({T}_{\mathrm{n}\mathrm{b}}^{*}\) is the last available value for the neighbor's node temperature. When all nodes are considered in this way, an iteration of the point-to-Gauss–Seidel method is completed. The Gauss–Seidel approach will not always converge; in fact, whenever the Scarborough benchmark (1985) is applied, the Gauss–Seidel approach is guaranteed. This criterion is seen in Eq. (11).

$$\frac{{\sum {\left| {a_{{{\text{nb}}}} } \right|} }}{{\left| {a_{P} } \right|}}\left\{ {\begin{array}{*{20}l} { \le 1} \hfill & {{\text{For\,all\,Equations}}} \hfill \\ { < 1} \hfill & {{\text{At\,least\,for\,one\,Equation}}} \hfill \\ \end{array} } \right.$$
(11)

The above convergence relationship is acceptable if time intervals are correctly selected and arranged (Patankar 1980). The selected time interval in this research is 0.0001 s which meets the Seidel-Gauss convergence criterion.

Volume of fraction

Nickels and Hert proposed the volume of fraction model in 1981. To determine the grid level, two fluid phases have been considered in many hydrodynamic problems. In hydraulic phenomena, the free flow is extremely important in the solution of the flow field. In this method, for each cell volume, a differential equation is solved, which ultimately determines the fluid component in each cell. If αw is the volume of water, then the volume of air is equal to Eq. (12) (Fluent Inc 2006).

$$\alpha_{\rm a} = 1 - \alpha_{\rm W}$$
(12)

In this paper, the free flow surface is defined by calculating the cells’ ratio of fluid to void. A 90% ratio is considered to be the surface flow line (this means the cell is 90% water and 10% air).

Boundary conditions

In general, one of the most important steps in the numerical analysis of the flow is determining appropriate boundary conditions and conforming to the problem’s physical condition. Considering that the governing equation for fluid flow is constant, depending on the type of fluid in terms of compressibility, the differences in the computational fluid dynamics arise from differences in the geometry of problems and their boundary conditions. Here, the importance of considering the appropriate boundary conditions becomes more apparent. Hamedi et al. (2016) suggested that velocity inlet and pressure inlet can be applied to inlet water level and air above the inlet water level respectively. Also, they applied pressure outlet for the outlet boundary condition and pressure inlet. The boundary condition of the wall, steps, and floor of the stilling basin, the boundary condition for the input velocity for the flow input, the fixed pressure boundary condition for the output of the model, and the constant boundary condition for the free surface are shown in Fig. 3 (Zahabi et al. 2018).

Fig. 3
figure 3

Boundary conditions of the numerical model of the stilling basin

Full size image

Results and discussion

Sensitivity analysis

In this research, for proper model gridding, an unstructured rectangular grid is used as shown in Fig. 4. Figure 5, the flow profile under this type of grid, is shown at different points in the stilling basin. Also, all calculations in the tables related to the two primary sections 1 and 2 of the desired stilling basin have been obtained from experimental results.

Fig. 4
figure 4

An example of gridding and grid size distribution (0.0085 m)

Full size image
Fig. 5
figure 5

The flow profile on the chute and the stilling basin

Full size image

To access the most suitable grid size for the proper calculation of flow factors on the stilling basin, the results of the numerical method are compared with the experimental results (Eq. 13). According to the relative difference values, the most appropriate option to resemble the grid size in calculating flow parameters on the stilling basin can be selected.

$$\hbox{Relative\,Difference} =\frac{|\hbox{Experimental\,Results} - \hbox{Numerical\,Results}|}{\hbox{Experimental\,Results}}$$
(13)

Accordingly, the numerical and experimental results for the velocity in the different sizes of the grid are shown in Tables 1 and 2. The results obtained for the 0.0085-m grid are the most accurate results with reasonable cost, compared with other cell sizes.

Table 1 Sensitivity analysis on grid dimensions in section 1
Full size table
Table 2 Sensitivity analysis on grid dimensions in section 2
Full size table

Effect of changing the height of the final step by the VOF method

In this section, the effect of height modifications on the last step of the chute on the flow profile is investigated by a numerical evaluation using an irregular rectangular grid with different mesh sizes. The RSM is used for numerical simulations. In Fig. 6a–c, the surface flow profile and pressure values are shown when modifications have been applied to the step height.

Fig. 6
figure 6

Flow profiles on the stilling basin with the height of the final step: a 0.50 m; b 0.60 m; c 0.70 m

Full size image

The depth and flow velocity values were measured by changing the height of the chute’s final step in two Sects. 1 and 2 (based on a grid size of 0.0085 m), on the last step before the basin and the mid-section of the stilling basin, respectively. (Tables 3, 4). Figure 7 indicates the comparison between surface flow profiles in the modes of increasing and decreasing the height of the chute’s last step with the actual flow profile.

Table 3 The effect of changing the height of the step on the flow parameters in section 1
Full size table
Table 4 The effect of changing the height of the step on the flow parameters in section 2
Full size table
Fig. 7
figure 7

Comparison of the water flow profiles in the numerical model with laboratory results

Full size image

The Navier–Stokes equations are used to solve the flow field. The RSM is used to model Reynolds stresses, and volume of fluid is used to determine the flow surface profile. The sensitivity analysis was performed on the numerical results of seven mesh sizes to determine the optimized mesh size, and from the results, it can be inferred that using the grid size 0.0085 m leads to the minimum relative difference of experimental and numerical results. Besides the RSM model for turbulence, the K-Epsilon and K-Omega model were applied. The results show that these models can’t properly capture flow behavior and contour.

Increasing the height of the last step of the chute in this model reduces velocity. Also, when the height of the downstream steps rises, the velocity vector in the x-direction decreases significantly. In step height 0.5 m, the depth and velocity obtained from the numerical solution in Sect. 1 are 0.204 m and 0.713 m/s, respectively, which matches data provided in the report of the Water Research Institute (2005) of 2.05 m and 0.710 m/s respectively. In this model, by increasing the height of the last step of the chute, flow velocity drops.

Conclusion

Based on a comparison of the water flow profiles in a numerical and laboratory model, as the final step of the chute changes, data obtained from the numerical approach are matched with experimental values. Also, by comparing the effect of changing steps’ height on flow parameters, it was found that by increasing the height of the final step of the chute, flow velocity decreases. Moreover, decreasing the height of the final step leads to an increase in flow velocity and a decrease in flow depth inside the basin.

Performing a sensitivity analysis on the grid size showed that a 0.0085 m size can save a considerable amount of time and cost. The relative error of the flow velocity (step height 0.5 m) is equal to 0.42% and 3.07% in sections 1 and 2, respectively. This implies high precision for the numerical method.