Discharge of a siphon spillway under submerged exit condition

The background of this study is that theoretically, the profile and size of a siphon spillway can be designed for one head only, much like an overflow spillway (Chow 1959). Although the spillways could also be capable of operating under other heads slightly different from the design head, substantial changes to the originally planned operating conditions have occurred or will soon occur in future hydrologic events (Ferdowsi et al. 2020). The changes will possibly cause most existing spillways to operate not only under heads much larger than they were designed for, but also under submerged exit condition (Roushangar et al. 2019). A good understanding of siphon flow characteristics under the new conditions of increased head and submerged exit is essential to further evaluations of cavitation aggravation and a need for retrofits or upgrades. The new conditions mean rising upstream and downstream water levels, which is a gap from previous research. The result is that the siphon carries blackwater flow, which differs from aerated weir flow or sub-atmospheric weir flow (Ali & Pateman 1980). There have been cases of climate change increasing runoffs in certain seasons (Blöschl et al. 2019). In the case the increase produces submerged inlet and submerged outlet conditions, it is physically meaningful to investigate the blackwater flow. When the blackwater flow exits from a siphon at an increased flow rate, the exit flow can possibly cause flash floods at downstream.

In the study of the siphon flow problem, one major challenge facing researchers is how to deal with turbulence, separation of the primary flow from the conduit walls and secondary flow (SF). To overcome the challenge, this paper adopts robust numerical techniques, physically sound boundary-layer treatments and progressive assessments of the techniques and treatments.

This paper aims to explore reliable numerical methods for predicting the flow characteristics. The intent also includes laboratory experiments of submerged siphon flow for measurements of the discharge coefficient C_d, useful for a validation of predicted C_d values. Numerical modelling is cost-effective and efficient and is convenient to be applied to different sites of siphon operation. Previously, Savage & Johnson (2001) used a computational fluid dynamics (CFD) model to successfully predict time-averaged flow over an overflow spillway. Taking advantage of the recent rapid growth of computing power and advancement of CFD techniques, this paper predicts detailed flow characteristics.

Siphon operation undergoes five distinct flow stages: free weir flow, deflected nappe, depressed nappe, air-partialized flow and blackwater flow (Head 1971). The geometric features (siphon crest level, entrance shape and level, exit lip level and deflector shape and position) influence the head–discharge relationship (Head 1975). Ervine & Oliver (1980) experimentally investigated the relationship for air-regulated siphons. They reported underestimates of discharge from experiments of scaled model siphons, when compared to field data from prototype siphons, and attributed the error to vortices at the model siphon entrance. This indicates the importance of turbulence, which is a research gap and is treated in this study. Babaeyan-Koopaei et al. (2002) studied the performance of the Brent Reservoir siphon spillway, aiming at improving discharge capacity and stability. Ali & Pateman (1980) analysed the relationship using several simplifications. Some of them (the flow is irrotational, free of turbulence and vortices; the flow is critical at the crest and has no loss of energy) are crude. Tadayon & Ramamurthy (2013) numerically predicted the C_d values under free flow condition. The influence of submerged exit condition on the head–discharge relationship and other subtle characteristics of flow remains to be discovered.

In the following, experiments of turbulent flow passing through a model siphon, conducted in the Hydraulics Laboratory at Concordia University, are described. Then, numerical methods for three-dimensional (3D) flow predictions are presented. Next, the results are discussed, along with data comparison. Lastly, conclusions are drawn.

We carried out experiments of flow passing through a model siphon under submerged exit condition in order to determine C_d. The siphon has a well-rounded entrance or rounded lips (Figure 1). The conduit has a constant rectangular cross section, with a width b = 25.1 cm and a throat depth d = 5.71 cm at the crest section, giving a cross-sectional area A = 143.32 cm². The crest arc radius is R_i = 1.91 cm. The crown arc radius is R_o = 7.4 cm. Downstream of the arcs, there is a straight, tangential (to the arcs) section. The siphon conduit has a rectangular exit at the downstream end.

The siphon conduit was fitted into a recirculating flume. The flume setup has a large head tank whose storage capacity (about 1.8 m³) is much larger than the volume of the conduit. This tank was attached to the siphon on the entrance side. This mimics an upstream reservoir that continuously supplies water to a siphon on a dam. The inside of the head tank had a ramp and a horizontal contraction. During the experiments, these devices forced flow streamlines to contract vertically and horizontally and thus damped background turbulence in the reservoir. Further, the flow passed through a honeycomb plate installed inside the tank, creating smooth streamlines approaching the entrance of the siphon. The head tank has an overflow opening and becomes a constant head tank if the water surface reaches the overflow level. In the channel section downstream of the exit, a tailgate controlled the water level and created target submergence. A floating polystyrene plate was placed at the water surface at a short distance from the exit, to attenuate surface waves and thus reduce water-level fluctuations during an experiment.

The stress tensor is ⁠, where u_i’ and u_j’ are the fluctuating velocity components in x_i and x_j directions, respectively. Based on Boussinesq eddy-viscosity approximation, the stress tensor is calculated as the product of an eddy viscosity, ν_t, and the mean strain-rate tensor, S_ij = (∂u_i/∂x_j + ∂u_j/∂x_i)/2. The specific turbulence kinetic energy (TKE) is defined as ⁠. For turbulence closure, this paper used schemes which relate ν_t to k and its dissipation rate ε. Values of these two turbulence quantities were computed using two-equation models of turbulence.

Water flows from an upstream reservoir through a siphon spillway into a downstream reservoir (Figure 1). The position of the air–water interface or the free surface is computed from a continuity equation: ⁠. For any computational cell, α₂ = 1 – α₁; α₂ = 1 means that the cell is occupied completely by water, whereas α₂ = 0 means that the cell is occupied completely by air. The free surface was discerned using the criterion α₂ = 0.5. The volume of fluid method (Hirt & Nichols 1981) is efficient and convenient for modelling two-phase flow.

The model domain consisted of two reservoirs, connected by a siphon section (Figure 1(b)). This section matched the design and geometry of the laboratory model siphon. The domain was bounded by solid walls and seven open boundaries, labelled as boundaries 1–7 in Figure 1(b). At the solid walls, u₁ = u₂ = u₃ = 0. Conditions imposed at the other boundaries were: (1) at boundary 1 (x₂ ≤ η₁), u₁ = U₁, u₂ = u₃ = 0 and α₂ = 1; p was relative hydrostatic pressure. The upstream reservoir water level η₁ was specified (η₁ = 0.964 m). The value U₁ was specified such that it matched the discharge (0.025 m³/s) used in the experiments; (2) at boundary 2 (x₂ ≤ η₂), u₂ = u₃ = 0, α₂ = 1 and p = 0; u₁ was extrapolated from u₁ values for the adjacent interior cells. The downstream reservoir water level η₂ was specified (η₂ = 0.612 m); (3) at boundaries 3 and 5, u₁ = u₃ = 0 and p = 0; u₂ was extrapolated from u₂ values for the adjacent interior cells; (4) at boundary 4 (x₂ > η₂), u₂ = u₃ = 0, α₁ = 1 and p = 0; u₁ was extrapolated from u₁ values for the adjacent interior cells; (5) at boundary 6, u₁ = u₃ = 0, α₂ = 0 and p = 0; u₂ was extrapolated from u₂ values for the adjacent interior cells, and was limited to u₂ ≤ 0; 7) at boundary 7, u₂ = u₃ = 0, α₂ = 0 and p = 0; u₁ adjusted itself to achieve a total of zero mass flow rate.

Finite volume mesh (hexahedra) was created to cover the model domain (Figure 1(b)) for numerical solutions to the governing equations. The spatial derivatives in these equations were approximated using second-order upwind schemes, and the time derivatives using first-order schemes. A pressure-correction equation was derived from the continuity and momentum equations. The velocity and pressure fields were solved using the pressure-based coupled algorithm. The momentum and pressure-correction equations were solved consecutively.

When creating the mesh, refinement was applied to near-wall regions, giving the dimensionless wall distance, y⁺, in the range 2.5 ≤ y⁺ ≤ 35 for the siphon conduit. Mesh independence of the solutions was tested using mesh sizes Δx = 2, 3 and 4 mm. Individual runs of the numerical model commenced from a state of rest and proceeded until a target period of time was reached and solutions converged. The time step was Δt = 0.001 s, satisfying the Courant–Friedrichs–Lewy criterion (Courant et al. 1967) for numerical stability, as estimated using Δx and the cross-sectional average velocity (v_s = 1.75 m/s) of the siphon flow. The convergence criterion per Δt was 1.00 × 10⁻⁶.

Three runs (RN2, RN3 and RN4) were carried out to investigate mesh independence. These runs used the same conditions and the same turbulence closure (RNG k-ε model), except the mesh size, which was 4 mm for RN4, refined to 3 mm for RN3 and further refined to 2 mm for RN2.

The computed velocities and relative pressures at a series of other locations were consistent among RN4, RN3 and RN2, confirming that the mesh size Δx = 4 mm (used in RN4) was sufficiently fine, and further refinement was not necessary. Thus, subsequent runs of the CFD model used the 4-mm mesh. On the one hand, it incurred lower computing costs than the 2- or 3-mm mesh. On the other hand, the model results were of independence of the computational mesh used.

Detailed knowledge of turbulence, eddies and SF in a siphon is of relevance to design improvement. To accurately compute the foregoing flow details, along with boundary-layer separation, how complex a turbulence model must be? This paper took the practical approach that compared the computed flows among three well-referenced models for turbulence closure.

The crest region is the most critical to siphon design, and relative pressure drop at the crest can lead to cavitation. Thus, the Realizable k-ε model was less desirable. Compared to the standard k-ε model, the RNG k-ε model is more sophisticated and more suitable for use in order to handle complicated boundary-layer phenomena. Discussion of the flow characteristics for RN4 follows.

The complexity of flow is further demonstrated by vertical profiles of u₃ (Figure 6(b)), which is a component of SF. The solid line showed a zero value of u₃ in the POS. In the middle width, the dashed curve showed the circulation of fluids in the direction from the sidewall toward the POS (u₃ > 0) over the lower throat depth and in the opposite direction (u₃ < 0) over the upper throat depth. Near the sidewall, the dotted-dashed curve showed changing directions of u₃ from the upper to the central and from the central to the lower portion of the throat depth. In cs2, u₃ had a maximum magnitude of as high as 0.36v_s.

The relative pressure decreased with decreasing distance from the crest and had negative values close to the crest (Figure 6(c)). The lowest value of p at the crest was proportional to the maximum speed shown in Figure 6(a). The TKE k was very strong within a very thin boundary layer near the crest (Figure 6(d)). The maximum speed of u₁ in the vicinity of the crest gave rise to strong velocity shear and hence the production of turbulence. The profiles of k showed strong turbulence in the upper throat depth, where 3D eddy motions and reversed flow took place, as indicated by changing directions of u₃ (Figure 6(b)). There was little turbulence in the portion of the throat depth where the profiles of u₁ resembled the potential flow solution.

The flow near the crest gained kinetic energy at the expense of local pressure energy. At a given Δh, the most severe drop of p occurs in the POS, compared to other planes parallel to it. For the scaled model siphon in this study, the negative relative pressures in the crest region are not low enough to create cavitation, but they imply a potential problem of cavitation at the prototype scale.

At most cross sections, the SF showed complicated structures with multiple eddies. At some cross sections, the SF reached a maximum strength of 16% of the 3D flow velocity. Given that SF reduces the discharge efficiency, it would be beneficial to retrofit existing siphons and adjust the design of new siphons to reduce SF.

One main finding from the aforementioned results is that the primary flow highly varies in siphon velocity (Figures 2, 4(a) and 5(a)), not to mention the SF (Figure 8), from one point to another at a given cross section and from one cross section to another along the flow path. The reasons for the spatial variations are twofold – boundary-layer flow separation, coupled with spiral/eddy motion. The variations are inherent features of siphon flow but have been ignored in the traditional approach to the siphon flow problem. The traditional approach relies on the energy relationship and resorts to the two key assumptions: (1) The siphon velocity is constant without spatial variations; (2) the energy losses caused by friction and by the siphon entrance and exit can be expressed in terms of siphon velocity head, which are to make the problem mathematically trackable. However, the inherent features of siphon flow revealed in this paper imply that the assumptions are very crude approximations and thus the traditional approach is questionable. In comparison, the new approach discussed in this paper leads to much-needed improvement over the traditional approach.

Contours of the TKE k are plotted in Figure 9 for cs1–cs6 (Figure 3). Upstream of the crest, the occurrence of turbulence was limited to the upper part of cs1 and cs2 (Figure 9(a) and 9(b)). At these two cross sections, had a maximum value of 0.876 and 0.947 m/s (or 0.50v_s and 0.54v_s), respectively. At cs3 downstream of the crest, the maximum k value was slightly lower than that at cs1, but turbulence also occurred in the corner between the sidewall and the siphon's lower surface (Figure 9(c)). Turbulence intensified at cs4; reached a maximum value of 0.67v_s. Further downstream, dropped to 0.55v_s and 0.42v_s at cs5 and cs6, respectively.

In the upper leg at cs1 and cs2, the regions of strong turbulence (Figure 9) and the flow separation zone coincided, whereas in the lower leg at cs3–cs6, the regions corresponded with the cores of SF. The turbulence activities were strongly driven by the large eddies associated with either boundary-layer separation or SF or both. The large eddies effectively interacted with the mean flow strain S_ij, feeding the turbulence. This is the key mechanism of turbulence production (Wilcox 2006), captured very well by the CFD model. The region of the strongest production was the elongated zone of flow separation immediately downstream of the crest. At any cross section in the lower leg, turbulence activities occupied a large percentage of the flow area. No variations of 3D velocity or turbulence are considered in classical analyses of siphon.

Another finding from the computational results is the coexistence of different scales of motion in siphon flow (see Figures 2 and 8). A key to the success of turbulence computations in this paper lies in accounting for the contributions of all scales of motion to the turbulent diffusion, including the effects of relatively small scales of motion. This has been achieved by mathematical techniques in the RNG k-ε approach, through changes to the production term in the modelled equation for ε. For the problem of siphon flow computations, it is not sufficient to determine the eddy viscosity (Equation (6)) from a single turbulence length scale, as is specified in the standard k-ε approach. Recently, Aminoroayaie Yamini et al. (2021) successfully applied the RNG k-ε approach to the problem of flow passing a bottom outlet of a dam.

Arguably, determining the coefficient C_d is practically the most important task in solving siphon flow problems. The C_d values from the experiments and numerical simulations in this study are presented in Table 1. The experimental values had an arithmetic mean of 0.64 and a standard deviation of 0.024 under a range of dimensionless heads ⁠. It would be sufficient for practising engineers to use C_d = 0.64 in the head–discharge relationship to estimate the siphon discharge. The numerical simulations covered the same dimensionless heads as the experiments and predicted C_d values in good agreement with the experimental values (Table 1). In this regard, the numerical methods and computation strategies discussed in this paper are reliable.

It is difficult to make measurements of flow velocity from the siphon conduit. In the absence of such measurements for validating the CFD model, this paper has taken rigorous steps to ensure mesh quality, turbulence closure suitability and numerical convergence and stability. Moreover, the predicted profiles of above-crest velocity (Figure 5(a)) are supported by the theoretical solution of Dressler (1978).

The design of siphon used in this study was a modification of the model of Head (1975), and siphons with different geometric features (e.g., the shape of entrance, deflector size and position and crest radius) may lead to somewhat different values of C_d. The numerical methods and computation strategies discussed in this paper offer cost-effective, efficient, robust tools for the exploration of the influence of different geometric features.

The methods and strategies would be applicable to prototype siphons; there is little scaling problem related to the Reynolds number Re or Froude number Fr (Table 1), which are defined as where D is the hydraulic diameter of the siphon conduit, and is the viscosity of water, and as Fr = ⁠, where is a characteristic velocity in the upstream reservoir. The values of Re in the test of the robustness of the methods and strategies were sufficiently high (Table 1), ensuring negligible viscous effects. In the upstream reservoir, Fr was very small (Table 1), and thus, no significant waves entered the siphon conduit from the reservoir. The test realistically reflects the field condition of siphon discharge from a quiescent reservoir.

The performance of existing geometrically similar small siphons which are subject to submerged conditions can be easily determined by changing downstream boundary conditions at a later time. Indeed, small siphons are used to divert excess water from irrigation channels to drainage. The availability of submerged siphon performance data permits development of validated submerged siphon characteristics. This enables one to develop validated CFD models instead of expensive experimental methods to find siphon characteristics when the downstream is submerged. Because of turbulence, streamline separation, SF vortices present in submerged siphon flow, one should not analyse the flow as a simple pressurized flow within a closed conduit (or a bent pipe).

In summary, the strengths of this study are on two aspects. The first aspect is rigorous procedures implemented for experimental measurements and numerical computations. The experimental part uses suitable devices (large head tank, streamlined vertical ramp, gradual horizontal contraction and honeycomb) for straightening the approach flow to the siphon, and uses instruments with high precision (water-level gauge and V-notch) for measurements. The computational part implements systematic tests about the suitability of mesh configuration, time interval of integration and turbulence closure schemes, and confirms the suitability through progressive data comparisons. The second aspect is the uniqueness of research contributions; a review of the literature shows that there has been little study of submerged siphon flow with detailed measurements and computations.

The exact equation for ε (Wilcox 2006) is complicated, involving multiple unknown double and ripple correlations of fluctuating velocity, pressure and velocity gradients. The Standard, RNG and Realizable k-ε models use closure approximations for the complicated correlations, yielding the modelled equation for ε (e.g., Equation (8)). It is worth noting that the model equation aims to predict the required length scale of the energy-containing, Reynolds stress-bearing eddies, whereas the exact equation represents the much smaller dissipating eddies. Thus, the relation between the modelled equation for ε and the exact equation for ε is weak (Wilcox 2006). This is a limitation of the turbulence models, used in this paper. According to Dimitriadis et al. (2016), there is another limitation, i.e., the models violate basic stochastic asymptotic properties; they do not include such important stochastic parameters as Hurst coefficient and fractal dimension.

In this paper, the CFD model of submerged siphon flow incurs high computing costs. This is another limitation of this study. Computational efficiency is an issue when there is a need to cover a large number of submergence scenarios and/or to accommodate prototype scales. Understandingly, future advances in computing power will increase the computational efficiency. This study is limited to a siphon of a rectangular cross section. Practical applications have used siphons of different shapes. For future research, it is recommended to explore the flow behaviours of prototype scale siphons of different shapes.

This study experimentally and numerically investigates submerged siphon characteristics. The following conclusions have been reached:

• The discharge coefficient C_d has a mean value of 0.64 and a standard deviation of 0.024 from experiments under a good range of dimensionless submergence levels and driving heads. It is sufficient to use the mean value for the purpose of determining submerged discharge (Equation (2)). Small siphons are also used to release excess flow from irrigation channels to drainage area. Even here, if free-flowing siphon exit gets submerged, one can use the submerged flow siphon C_d in place of the free-flowing siphon C_d.

• The numerical predictions give discharge coefficient C_d matching the experimental results and velocity profiles at the siphon crest supported by the theoretical solution (Dressler 1978). The numerical methods and computation strategies from this paper are robust for predicting submerged siphon characteristics. The RNG k-ε model for turbulence closure realistically captures the boundary layers in the siphon conduit and velocity profiles between the boundary layers; it outperforms the standard k-ε and realizable k-ε models.

• The flow velocity displays large spatial variations in not only the primary flow but also the secondary flow. The velocity attains a maximum magnitude in the vicinity of the crest, whereas the relative pressure drops to negative values, implying cavitation risks in a submerged siphon.

• The secondary flow is complex at most cross sections, featuring multiple eddies. It reaches a maximum strength of larger than 16% of the primary flow strength. It would be beneficial to retrofit existing siphons and adjust the design of new siphons to reduce secondary flow.

• The flow separates from the upper surface of the upper leg all the way to the siphon crown and produces elongated eddy motions. The separation is triggered by the lip curvature at the siphon entrance. The resultant shear causes strong turbulence. Due to the presence of eddies, shear and turbulence, one should not treat submerged siphon flow as a simple pressurized flow in a close conduit.

This study received financial support from the National Sciences and Engineering Research Council of Canada through Discovery Grants held by S. S. Li.

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

The authors declare there is no conflict.