Numerical investigation of alternating skimming flow over a stepped spillway Open Access

Stepped spillways (SP) are a famous technique used for their simplicity of shape to transfer high-energy water and to decrease flow velocities (Chanson 2015). However, some physics are not yet understood like the formation of an alternating skimming flow (Lopes et al. 2017). In this paper, the limit width, from which the formation of the alternating skimming flow occurs, is investigated, as well as the reasons and what is behind this formation.

The flow over SP occurs in three different regimes: nappe flow, transition, and skimming flow (Wan et al. 2019). The SP is typically designed to operate in a skimming flow regime (SK), which can be one of the two sub-regimes: SK1 and SK2 (Chanson 2002; Ohtsu et al. 2004). The SK regime is characterized by a main water flow over a pseudo-bottom, which connects all the steps where the development of a boundary layer along the length appears, reaching the free surface and allowing air-entrainment. This is usually called the inception point. Toro et al. (2017) and Zabaleta et al. (2020) presented details of the mechanisms leading to air entrainment in the skimming flow, as it is likely that the phenomenon occurs because of vorticity patches that were generated in between the steps, which then become released to the flow and reach positions close to the free surface. In the SK1 sub-regime, the mixing layer does not reach the step end, an undulant free surface appears, and in extreme cases, this effect is observed partially parallel to the step tread. The reason for this behaviour is the wake zone formed in each step downstream and the recirculation vortex underneath. In the SK2 sub-regime, there is an interaction between the wake and the recirculation vortex of one step and the wake formation of the subsequent one. Moreover, in this sub-regime (SK2), the time-averaged free-surface is roughly parallel to the pseudo-bottom being independent of this sub-regime's effect. New research in the field of SP, based either on experimental or numerical modelling, has focused on the SK regime (Wan et al. 2017; Bayon et al. 2018; Nóbrega et al. 2020).

Nowadays, several research methods based on physics and computational fluid dynamics (CFD) have greatly promoted the development of studying the characteristics of complex fluid flow in the entire hydraulic structure (Ghaderi et al. 2020a, 2020b; Zhou et al. 2020). Numerical simulation contributes to reduce laboratory tests and to conserve time/money (Morovati & Eghbalzadeh 2018; Valero & Bung 2018; Ghaderi et al. 2020a, 2020b; Güven & Mahmood 2021).

By using this approach, Lopes et al. (2017) presented different turbulence closures via steady-state and time-dependent approaches in three dimensions (3D) and verified a zig-zag behaviour of the flow close to the pseudo-bottom, which they named an alternating skimming flow regime. The phenomenon occurs when a combination of the flow nature or the dimensional characteristics of the SP convene. They studied two different SP with two different widths, 0.30 m and 0.50 m, and they found an alternating skimming flow for the largest 0.50 m-wide SP (SSW0.50), which was not visible in the 0.30 m-wide SP (SSW0.30).

In this paper, our focus is to continue the previous work of Lopes et al. (2017), and to understand the reason behind the formation of the alternating skimming flow and whether it is related to the vorticity patches generated between the steps discovered by Zabaleta et al. (2020); also it is to identify the critical width for this alternating skimming flow.

To achieve this purpose, we used the experimental data of Lopes et al. (2017), to validate the 3D numerical models established in Ansys Fluent^®.

Then, concentrating on the parameter of the width, we constructed similar models with different widths (w) and the same step height (s), expressing a non-dimensional ratio w/s to identify the widths for which the alternating skimming flow occurs. We also investigated the influence of the spillway width over the alternating skimming flow regime, and for doing this, the spillway was changed for different stages, starting from 0.1 m and increasing to 0.6 m. As a result, the alternating skimming flow regime in several steps can be defined by the alternated presence of the SK1 and SK2 sub-regimes.

The work is organized as follows: ‘Experimental facility’ describes the physical experiment in the constructed facility to measure the flow depths and velocities. ‘Numerical model’ presents the equations used in the numerical simulation research, the computational models (mesh generation and its analysis and the definition of constants and procedures), as well as its numerical validation using experimental data, and introduces the methodology followed to study the free surface position. ‘Results’ shows the numerical simulation results for the different widths focusing on their comparison. The transition and the characteristics of the alternating skimming flow regime are discussed in ‘Discussion’. The final section summarizes the main findings of this work.

Experimental tests were previously conducted at the Hydraulic Engineering Section, FH Aachen, University of Applied Sciences (Germany) by Lopes et al. (2017). The setup consisted of an SP model (scale 1:10) with 0.5 m width. This model was composed of 28 identical steps, the step height was 0.06 m and the length 0.12 m, placed in a channel with a slope of φ = 26.6° (Figure 1). The SP model was assembled in a closed loop facility. The spillway crest featured a filter for minimizing the flow turbulence for discharges q_w = 0.07 m²s⁻¹ (measured by inductive flow meter). Further details regarding the experimental facility can be found in Lopes et al. (2017).

The bubble image velocimetry (BIV) measurement technique and an ultrasonic sensor (micUS) were used for measuring the velocity profiles at step niches and water depth at step edges, respectively. The operating range of the used ultrasonic sensors was 200–1,300 mm, their resolution was 0.18 mm and they were moved on to the top of the step edge to measure the water elevation. The data was sampled for 180 s at 200 Hz. The setup also used the bubbles as tracers, illuminating them with halogen spotlights. The technique helped the high-speed camera to capture the flow velocity field near the sidewall through the BIV technique. The bubbles can leave focused planes and distort the resulting velocity field. Also, following Felder & Chanson (2014a), sub-sample durations between 1 and 45 s were selected so the maximum cross-correlation coefficient does not change significantly.

In this research, the SP model studied experimentally is simulated using the commercial code Ansys Fluent^®, where special attention is given to alternating skimming flow behaviour and the numerical results are compared with the experimental data. Then, the model is used to studied different stepped spillway widths. The following sections describe the overall characteristics of the proposed numerical model regarding its definition and implementation.

The SP numerical model uses a set of 3D Reynolds-Averaged Navier–Stokes equations, together with the turbulence model closure: SST k-ω (Menter 1994), and volume of fluid (VOF) method, proposed by Hirt & Nichols (1981) and used in Carvalho & Amador (2009). The SST k-ω is a hybrid model that merges the best of two RANS formulations, k-ω (Menter 1994) and k-ɛ (Launder & Sharma 1974). The model works as a k-ω model in the near-wall region to glean the vortices on the spillway steps, and as k-ɛ in the free stream region (above the pseudo-bottom) where the high Reynolds number formulation of the model is super-effective.

The mechanical properties of the fluid used in this paper are: ρ_w = 998.78 kg m⁻³ and μ_w = 1.1094 × 10⁻⁶ m²s⁻¹ as density and kinematic viscosity of water, and ρ_a = 1.225 kg m⁻³ and μ_a = 1.4657 × 10⁻⁵ m²s⁻¹ for air. Regarding surface tension, it was defined as 0.072 kgs⁻². The discretization scheme for momentum, turbulent kinetic energy, and dissipation rate second-order methods were used as they are relevant when representing swirl flows. Finally, we used a SIMPLE coupled phase algorithm, to couple the pressure velocity. Moving to the boundary conditions, we considered the hydrostatic pressure for inlet on the left, with input water and air velocity, and outlet on the right. Also, the atmosphere was on the top and a no-slipping and stationary wall at the spillway. At the bottom, standard wall functions were applied to avoid the use of excessively fine meshes and to save computational resources while modelling the boundary layer. The initial conditions considered at t₀ are water just at the inlet with 0.1 m water depth (h) and velocity of 0.7 ms⁻¹ (u_x).

Following these boundary conditions, we looked for a difference that is less than 1% between the inlet and outlet water flows. Also, we established as conditions the values k, ω, volume fraction, and mass staying constant (allowing a tolerance of 10⁻⁶ for the residual). With these conditions and by using a time step variable to find the compatibility with the Courant number, we found that 22 seconds is the value for the steady-state solution (Kaouachi et al. 2019).

To optimize the process and reduce numerical problems, the best choice was to work using structured mesh algorithms. They also offer more regular memory access and reduce significantly its latency (Keyes et al. 2000). To use this method, we considered a static structured rectangular hexahedral mesh. To be sure that the mesh was not affecting the study, three different meshes – N1 (0.01 m × 0.01 m × 0.01 m) with 797,600 cells, N2 (0.012 m × 0.012 m × 0.012 m) with 481,488 cells, N3 (0.015 m × 0.015 m × 0.015 m) divided in 256,734 cells – were tested. These meshes were created by using Meshing Modeler tools from Ansys Fluent^®. Once the meshes were defined, Richardson extrapolation was performed, following Celik et al. (2008). The results obtained for the velocity profile at step niche 12 (Figure 2) showed reliability on N1 (0.01 m × 0.01 m × 0.01 m) as the grid convergence indexes were GCI₂₁ 8.22%, and GCI₃₂ 36.73%.

To validate the numerical model, we did two tests: (1) validation of the numerical time-averaged free surface against flow depths measured with the ultrasonic sensors along the centreline of the stepped spillway chute for q_w = 0.07 m²s⁻¹ (Figure 3); (2) comparison of the velocity profiles (taken by BIV) at steps niches 8 and 13 with the numerical results. The use of ultrasonic sensors to measure flow depths is a quite reliable method even in aerated flows (Felder & Chanson 2014b). It should be noted that validation of the numerical model may also include other features of skimming flows on stepped spillways, e.g. the air–water flow characteristics. However, air–water properties are not the focus of this work. Furthermore, the VOF model cannot be used to simulate entrapped air.

The results of the turbulence models SST k-ω are compatible with the experimental measurements obtained by the ultrasonic sensor; wherever the flow depth decreases gradually from the crest it equals 48 mm until becoming constant and roughly equals 32 mm at step 5 and 6. We realized that at step 5 the water depth rises due to the air-entrainment. The SST k-ω model is a hybrid model working as a k-ω at the wall and as a standard k-ɛ in the free stream. This is done to improve the performance in flows under adverse pressure gradients (such as those with separation as on a stepped spillway in skimming flow).

Figure 4 compares BIV velocity profiles against numerical results at step niches 8 and 13, showing very similar profiles. As an example, looking at the step 13 profiles the mean velocity's RMSE was 4.4

A slight difference is found in the step 8 (Figure 4(a)) profiles (BIV measurements and numerical results) of approximately 5%, which is very good.

Table 2 presents the different widths of all stepped spillways studied using Ansys Fluent^® with the SST k-ω model and VOF method. Results analysis is presented in the following sections.

Figure 5 presents the cross-section free-surface profile in all stepped spillway models corresponding to different widths. There is a gradual decrease in water depth along the steps, which is a safety factor and considered the main advantage of the stepped spillway. Regarding variations along the width, the free surface is quasi-stable in SSW0.10 (Figure 5(a)) and SSW0.20 (Figure 5(b)) on all the steps, while SSW0.35 (Figure 5(d)) and SSW0.45 (Figure 5(f)) present a nearly horizontal free-surface, and SSW0.30 (Figure 5(c)) and SSW0.40 (Figure 5(e)) present an oscillatory free-surface. In wider stepped spillways, SSW0.60 (Figure 5(h)), there is an increase of the number of free-surface oscillations and their amplitude along the width, which is just smooth in the first and second steps. In SSW0.40, SSW0.50 (Figure 5(g)) and SSW0-60 after the third step the amplitudes of free-surface oscillations are significant.

Figures 6 and 7 present respectively the V/Vc distribution (where V is the average velocity and Vc is the critical velocity), and average turbulent kinetic energy (k) at two different longitudinal sections from step 1 to step 19. For all the studied models, we present two sections named W/4 and W/2, where W is the total width. For example, in the SSW0.10 model, W = 0.10 m and we present Z (W/4) = 0.075 m, Z (W/2) = 0.05 m, considering different conditions of oscillations in all widths (see Figure 5).

Analysing features for the increasing stepped spillway widths, we observed that while turbulent kinetic energy near the free-surface starts to differ at Z (W/4) and Z (W/2) from 0.2 m, the velocity becomes different at Z (W/4) and Z (W/2) from 0.35 m at the last steps. The difference at Z (W/4) and Z (W/2) increases for wider stepped spillways and suggests an existing variation of velocity and kinetic energy along the width, which is in accordance with the observation of an oscillating free-surface referred to in the above section.

Figure 8 illustrates the pressure distribution on the step surfaces for the different stepped spillway widths along the length. In SSW0.10 (Figure 8(a)), the pressure on different horizontal step surfaces is uniform, but it seems very clear that the pressure increases longitudinally along the x-axis. In SSW0.20 (Figure 8(b)) no significant change related to SSW0.10 can be observed. Even the locations where the minimum and maximum pressure occurred are similar.

The pressure distributions in SSW0.30 (Figure 8(c)) and SSW0.35 (Figure 8(d)) are also similar on all steps, with a slight increase/decrease in the middle of the last steps.

For stepped spillway SSW0.40 (Figure 8(e)) the oscillation is clear, and it is obviously indicating that there are two types of pressure distribution, from step 1 to step 7, with an organized band distribution of pressure, and an alternating skimming flow after step 7 until the end.

For wider stepped spillways SSW0.45 and SSW0.50 as illustrated in Figure 8(f) and 8(g), two types of the pressure distribution can be observed: a distribution along the first steps with an irregular pattern appearing on the steps, and a distribution appearing from the final steps, where we can observe an alternating skimming flow pattern.

In the particular cases SSW0.50 and SSW0.60, as shown in Figure 8(g) and 8(h), after step 12, the pressure distribution presents alternating oscillations with two or three peaks along the step width until step 19, leading to an increase of pressure intensity, reaching the maximum at step 19.

Figures 9–14 present the velocity field in transversal cross-section at different steps, for all the studied models. All these models show regular velocity until step 6/8, then an irregularity is observed, in particular in the models SSW0.50 and SSW0.60

In all stepped spillways, the velocity at the free surface increases along the longitudinal direction, growing from the crest until reaching the maximum at step 17. These results are in accordance with previous research works (Leandro et al. 2014; Lopes et al. 2018). According to Meireles et al. (2014), the location of the inception point occurs when the k at the free surface reaches 0.22 m²/s². Following this, we estimated the inception point at step 10 for w < 0.3 (SSW0.10, SSW0.20, SSW0.30); at step 9 for 0.35 < w < 0.4 (SSW0.35, SSW0.40); at step 7 for w = 0.5 (SSW0.50); and at step 6 for w = 0.6 (SSW0.60), which agrees with the previous results. Moreover, Figures 6 and 7 show a global perspective of the steps, where it is noticeable that the flow skims over all the steps, from 1 to 10 faced in both plans, indicating the presence of an SK2 regime. Nevertheless, this is not seen at steps 11 to 19, where the flow covers all the even steps faced in W/4, which is a typical description of an SK1 regime.

Velocity results are fully consistent with the pressure contour and kinetic energy dissipation. Where there is a distribution showing an alternation of high- and low-pressure values along the width and when elevated pressure appears in the centre of odd steps and near walls in even steps, the velocity is larger.

For the stepped spillway SSW0.10, the free surface shape is steady at all the steps (Figure 5(a)). The pressure field shows a uniform distribution along each section (Figure 8(a)). The velocity field (Figure 6(a)) shows perfect similarity in both profiles' W/4 and W/2, which leads to identical values of kinetic energy turbulence at the free surface and in step cavities (Figure 7(a)).

The SSW0.20 model shows that the free surface along a cross-sectional profile is also quasi-stable at all the steps (Figure 5(b)). The pressure on different horizontal step surfaces shows uniform distribution (Figure 8(b)), as well as velocity longitudinal distribution (Figure 6(b)). However, looking to the k values, we can see a larger rise in this model than in SSW0.10 simulations, especially for step 7, where it reached its higher value (Figure 7(b)).

The SSW0.30 model also shows regular velocity until step 6/8, then a change in the velocity values near the wall can be observed (Figure 9), which corresponds to a change in the pressure distribution and in the kinetic energy near the surface from step 6 (Figures 8(c) and 7(c)). A bell-shaped distribution of the free surface starts from step 8 across the model width (Figure 5(c)). However, a perfect similarity of the main flow velocity shows up in the longitudinal profiles W/4 and W/2 (Figure 6(c)). An identical k value from the crest to step 5 is shown in the profiles W/4 and W/2 (Figure 7(c)). There is an increase of 20% of k in W/2 compared with W/4 for steps 6 to 28. Focusing on the step cavity, a slight drop of k is found in odd steps for W/2.

The SSW0.35 model shows regular velocity distribution until step 12 (Figure 10). After that, a velocity irregularity starts to appear on both sides in all the steps, which may be noted in step 13. In the pressure contour (Figure 8(d)) no irregularity in the pressure is detected.

However, the longitudinal sections W/4 and W/2 present a significant similarity of the main flow velocity as well as in the step cavity (Figure 6(d)). Regarding k profiles for W/4 and W/2 (Figure 7(d)), the k values at W/2 are higher than the k values at W/4 by 13% for the steps between 4 and 8. After that, the k values increase up to 20% for steps 9 to 19. From SSW0.40, it is clear for every two steps there appears one peak of free-surface velocity at the centre, while in alternating steps two appear on both sides showing a similar location (Figure 11). Figure 5(e) shows a bell-shaped oscillation in the cross-section free-surface profile, noted with a high peak starting from step 8. However, in SSW0.45, this is not so evident. Small differences are noticed in both velocity and k profiles (Figure 6(f)). Only the k value at the free surface rises above step 7 in the W/2 profile (Figure 7(f)). The cross-section free-surface profile is quasi-stable on the first five steps (Figure 5(f)), then a slight oscillation appears at the centre and near the walls. This also suggests a dependence of the ratio between width and length of oscillations, which will not be discussed in this paper.

Figure 15 illustrates pressure contours and streamlines for (a) SSW0.35, (b) SSW0.50 and (c) SSW0.60, where in SSW0.35 the uniform scale of the vortices over the width can be observed (Figure 15(a)), which can explain a significant similarity of the main flow velocity as well as the step cavity (Figure 6(d)). Figure 15(b) (SSW0.50) clearly shows two alternating vortices across the step face 12 and 13. This is consistent with the cross-section free-surface profile, where we observe intensive oscillation (Figure 5(g)).

The longitudinal section of velocity W/4 (Figure 6(g)) demonstrates that the flow does not skim until the end of each step (SK1). This effect is not seen for W/2 where, from step 10 and on, the flow reaches the end of the odd steps (SK2). In even steps, the mixing layer does not reach the end of the step (SK1). On the horizontal step faced, there is a large pressure differential (Figure 8(g)). This is a phenomenon to be taken into consideration during the design of masonry stepped spillways as described by Winter et al. (2010), as a pressure variation in short distances may cause the masonry block to rotate out of position. The SSW0.60 model is similar to SSW0.50 in terms of velocity distribution. However, in step 12 three vortices are formed which turn into only two in step 12, and after that they return to only two vortices again in step 13 (Figure 15(c)). From step 8 to step 20 (Figure 8(h)) three vortices are formed in even steps and just two in odd steps.

Regarding the TKE profiles in SSW0.40 for W/4 and W/2 (Figure 7(e)) starting from step 6, the k value is higher in W/2 than the k value at W/4 by 20%, and the k value of odd steps at W/4 is equal to the k of even steps at W/2. In SSW0.50, k profiles at W/4 and W/2 (Figure 7(g)), the value of k increases in direction to the free surface. For the W/4 section, the k in even steps is higher than in odd steps. At SSW0.60 the k value increases compared with the previous cited models, but a high similarity shows up in the profile's W/4 and W/2 (Figure 7(h)).

In summary, when the relative stepped spillway width is low, no differences of pressure and velocity across the channel width are noted. With increasing width, some differences across the width appear, first in kinetic energy near the free surface, pressure and velocity field. From relative stepped spillway width around 6 (w/s = 0.35/0.06 = 6, SSW0.35), velocity differences across the width appear near the wall on both sides, leading to a higher free-surface velocity in the centre and some oscillation near the bottom. For wider stepped spillways (w > 0.4 m, SSW0.40) clear differences are observed in both pressure and velocity distributions, either along the stepped spillway or along its width starting in the non-aerated region. For the models SSW0.50, SSW0.60, generally the differences occur earlier (always in the non-aerated region with the exception of SSW0.60, where the inception point and the beginning of differences occur at step 6) and the number of oscillations increases as well as the amplitude.

Differences observed for 0.3 < w < 0.5 suggest that the number of oscillations and their amplitude are dependent on the relative width per oscillation length (Figures 5 and 9–13). For w > 0.5 a new mechanism concerning the formation and distribution of vortices occurs in the aerated region, which has not been seen in the previous studied models. The dominant vorticity inside the step cavities described in this research was mentioned in the work of Amador et al. (2006) and observed also in numerical results executed by Toro et al. (2016), and Nóbrega et al. (2020) for identical conditions to those of Amador (2005). These vorticity patches generated in between the steps were recently mention by Zabaleta et al. (2020).

Towards the understanding of the formation of alternating skimming flow, noted by Lopes et al. (2017), as well as to achieve a width limit for its occurrence, we did numerical research on the skimming flow regime in different stepped spillways with various widths. using CFD with fully three-dimensional turbulence model closure SST k-ω and employing the VOF method for free-surface detection. Then we used laboratory data from a specific stepped spillway width to validate the model. The following conclusions can be retrieved from this work:

• CFD based on VOF and RANS is a valid methodology to simulate the flow over a stepped spillway in the laboratory with an appropriate computational mesh, and the SST k-ω gives good results.

• While when the stepped spillway width is less than 0.35/0.06 (w/s = 6), no differences of pressure and velocity across the stepped spillway width are relevant, from that limit, a clear alternating pattern can be identified in the distribution of velocity, kinetic energy and pressure.

• The relative stepped spillway width of w/s = 0.35/0.06 = 6 (w between 0.35 and 0.4), can be identified as a transition (critical limit) since it leads to the formation of vortices along the width which are still of uniform size and still keep the similarity of the main flow velocity.

• When increasing stepped spillway width, an alternating skimming flow and inception point appear early.

• Asymmetrical distribution of vortices in terms of shape and size are present.

• A dependence from the ratio of width and oscillation length is suggested, however more simulations are needed.

Future work on alternating skimming flow will focus on research of more simulations using different widths and conjugate with more flow rates. Research on the influence of air concentration profiles and bubble behaviour relation with the alternating patterns is also an interesting topic.

The first author would like to acknowledge the support of MESRS (Ministry of Higher Education and Scientific Research, Algeria) through the scholarship PNE. The second and third authors acknowledge FCT (Portuguese Foundation for Science and Technology), through the Project UIDB/04292/2020 – MARE, which was financed by MEC (Portuguese Ministry of Education and Science) and the FSE (European Social Fund), under the programs POPH/QREN (Human Potential Operational Programme from National Strategic Reference Framework) and POCH (Human Capital Operational Programme) from Portugal 2020.

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