## Abstract

Flow through a siphon is difficult to predict due to inherent turbulence, separation and secondary circulation. This paper overcomes the difficulty by using advanced numerical techniques and rigorously assessing their suitability. The aim of this paper is to explore reliable numerical methods for predicting submerged siphon characteristics. Using the Reynolds-averaged Navier–Stokes equations, we predicted three-dimensional velocity, pressure and turbulence quantities. We also conducted laboratory experiments for measurements of the submerged discharge coefficient *C*_{d}. The mean value of *C*_{d} predicted matches the measured mean value. The numerical results show flow separation in the siphon upper leg, causing secondary flow (SF) and increasing velocity above the crest. The SF shows complicated patterns and multiple turbulent eddies and reaches a maximum relative strength as large as 16%. The relative pressure has negative values in the crest region. The profiles of predicted longitudinal velocity in the crest region resemble the theoretical solution. The numerical methods and computation strategies from this paper are useful for investigating the performance of submerged siphons of various dimensions and/or geometric configurations under a wide range of hydraulic conditions. The RNG *k*-*ε* model is more suitable than the standard *k*-*ε* model and the Realizable *k*-*ε* model for turbulence closure.

## HIGHLIGHTS

Submerged siphon discharge coefficient from laboratory experiments is reported.

Reliable numerical methods for predicting submerged siphon characteristics are developed.

Complex primary and secondary flows, flow separation and turbulent vortices are predicted.

### Graphical Abstract

## INTRODUCTION

*et al.*(2020) experimentally studied the energy dissipation.

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.

## METHODS

### Laboratory experiments

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^{2}. 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^{3}) 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.

*Q*. This tank is 304.8 cm long, 60.96 cm wide and 91.44 cm high. The V-notch angle is

*θ*= 60°. A honeycomb plate was installed inside this tank to damp turbulence and thus improves the measurement accuracy. The static head

*H*in the approach water over the V-notch weir was measured using point gauges with an accuracy of ±0.1 mm. The discharge over the weir is given bywhere

*C*is the V-notch discharge coefficient, equal to 0.6 (Henderson 1966);

*g*is the gravity. The accuracy of

*Q*measurements is better than 97%. This V-notch setup has consistently produced accurate measurements of

*Q*in previous studies (e.g., Ramamurthy

*et al.*2017; Thapa

*et al.*2018).

*Q*for the flow is given by Equation (1). The discharge coefficient

*C*

_{d}was determined as the ratio of the actual to the theoretical discharge

Case . | Q
. | Δh/d
. | h_{1}/d. | h_{2}/d. | Re . | Fr . | C_{d}. | |||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

(m^{3}/s)
. | Exp^{a}
. | CFD^{b}
. | Exp . | CFD . | Exp . | CFD . | Exp . | CFD . | Exp . | CFD . | ||

1 | 0.024 | 6.340 | 9.019 | 13.012 | 2.680 | 6.673 | 247,943 | 247,943 | 0.08 | 0.05 | 0.628 | 0.628 |

2 | 0.025 | 5.762 | 8.879 | 12.872 | 3.117 | 7.128 | 236,372 | 236,012 | 0.09 | 0.05 | 0.687 | 0.688 |

3 | 0.024 | 6.340 | 8.196 | 12.207 | 1.856 | 5.867 | 247,943 | 247,943 | 0.10 | 0.05 | 0.628 | 0.628 |

4 | 0.022 | 5.412 | 9.089 | 13.100 | 3.678 | 7.671 | 229,075 | 229,445 | 0.07 | 0.04 | 0.623 | 0.622 |

5 | 0.024 | 5.902 | 9.842 | 13.853 | 3.940 | 7.933 | 239,228 | 239,583 | 0.07 | 0.05 | 0.651 | 0.650 |

6 | 0.024 | 6.165 | 9.229 | 13.240 | 3.065 | 7.075 | 244,494 | 244,494 | 0.08 | 0.05 | 0.637 | 0.637 |

Case . | Q
. | Δh/d
. | h_{1}/d. | h_{2}/d. | Re . | Fr . | C_{d}. | |||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

(m^{3}/s)
. | Exp^{a}
. | CFD^{b}
. | Exp . | CFD . | Exp . | CFD . | Exp . | CFD . | Exp . | CFD . | ||

1 | 0.024 | 6.340 | 9.019 | 13.012 | 2.680 | 6.673 | 247,943 | 247,943 | 0.08 | 0.05 | 0.628 | 0.628 |

2 | 0.025 | 5.762 | 8.879 | 12.872 | 3.117 | 7.128 | 236,372 | 236,012 | 0.09 | 0.05 | 0.687 | 0.688 |

3 | 0.024 | 6.340 | 8.196 | 12.207 | 1.856 | 5.867 | 247,943 | 247,943 | 0.10 | 0.05 | 0.628 | 0.628 |

4 | 0.022 | 5.412 | 9.089 | 13.100 | 3.678 | 7.671 | 229,075 | 229,445 | 0.07 | 0.04 | 0.623 | 0.622 |

5 | 0.024 | 5.902 | 9.842 | 13.853 | 3.940 | 7.933 | 239,228 | 239,583 | 0.07 | 0.05 | 0.651 | 0.650 |

6 | 0.024 | 6.165 | 9.229 | 13.240 | 3.065 | 7.075 | 244,494 | 244,494 | 0.08 | 0.05 | 0.637 | 0.637 |

Measurements of flow depth and siphon conduit width were true to mm. Discharge measurement errors were limited to 3%.

^{a}For the laboratory experiments.

^{b}For the numerical simulations using the CFD model.

### Governing equations for numerical computations

*i*= 1, 2, 3);

*t*is the time;

*ρ*and

*μ*are, respectively, the density and dynamic viscosity of the fluid mixture;

*p*is the Reynolds-averaged relative pressure;

*τ*is the specific Reynolds stress tensor;

_{ij}*g*is the gravitational acceleration in

_{i}*x*direction. The flow velocity vector iswhere is a unit vector in

_{i}*x*direction. For the fluid mixture,

_{i}*ρ*and

*μ*are calculated as the volume-weighted averages:

*ρ*=

*α*

_{1}

*ρ*

_{1}+

*α*

_{2}

*ρ*

_{2}and

*μ*=

*α*

_{1}

*μ*

_{1}+

*α*

_{2}

*μ*

_{2}, where the subscripts 1 and 2 refer to air and water, respectively, and

*α*is the volume fraction. For any computational cell,

*α*=

*α*

_{1}+

*α*

_{2}= 1.

The stress tensor is , where *u _{i}*’ and

*u*’ are the fluctuating velocity components in

_{j}*x*and

_{i}*x*directions, respectively. Based on Boussinesq eddy-viscosity approximation, the stress tensor is calculated as the product of an eddy viscosity,

_{j}*ν*, and the mean strain-rate tensor,

_{t}*S*= (∂

_{ij}*u*/∂

_{i}*x*+ ∂

_{j}*u*/∂

_{j}*x*)/2. The specific turbulence kinetic energy (TKE) is defined as . For turbulence closure, this paper used schemes which relate

_{i}*ν*to

_{t}*k*and its dissipation rate

*ε*. Values of these two turbulence quantities were computed using two-equation models of turbulence.

*k*-

*ε*model (Jones & Launder 1972), renormalization group (RNG)

*k*-

*ε*model (Yakhot & Orszag 1986) and Realizable

*k*-

*ε*model (Shih

*et al.*1995), were compared to assess their suitability for the problem of siphon flow. The standard

*k*-

*ε*model assumes an algebraic formula for

*ν*and uses the modelled equations for

_{t}*k*, and TKE dissipation rate

*ε*of the following form:where the dissipation rate is defined as . Equation (6) is based on dimensional arguments (Buckingham's

*Π*-theorem). Equations (7) and (8) describe the changes of

*k*and

*ε*due to advection, production, dissipation and diffusion. The closure constants (Launder & Sharma 1974) are

*C*

_{1}

*= 1.44;*

_{ε}*C*

_{2}

*= 1.92;*

_{ε}*C*= 0.09;

_{μ}*σ*= 1.0 (Prandtl number for

_{k}*k*);

*σ*= 1.3 (Prandtl number for

_{ε}*ε*). In Equation (7),

*k*is in the denominators of the production and dissipation terms, creating a singularity in flow regions of low

*k*values or of locally weak turbulence at low values for the Reynolds number Re. Nevertheless, the standard

*k*-

*ε*model is widely used.

*k*-

*ε*model uses Equations (7) and (8) like the standard

*k*-

*ε*model but different closure constants:

*C*

_{1}

*= 1.42; C*

_{ε}_{2}

*= , where = 1.68,*

_{ε}*η*= (2

*S*)

_{ij}S_{ij}^{1/2}

*k*/

*ε*,

*C*= 0.085,

_{μ}*η*

_{0}= 4.38 and

*β*= 0.012. The model replaces the diffusion coefficient (

*μ*+

*μ*/

_{t}*σ*) in Equation (7) with

_{k}*α*

_{k}μ_{eff}and replaces

*μ*+

*μ*/

_{t}*σ*in Equation (8) with

_{ε}*α*

_{ε}μ_{eff}, where

*α*=

_{k}*α*= 1.393 (inverse effective Prandtl numbers), and

_{ε}*μ*

_{eff}is the effective eddy viscosity. Instead of Equation (6), the eddy viscosity is solved from a differential equation:where ;

*C*

_{ν}≈ 100. This equation describes effective turbulent transport varying with effective Re (eddy scale), and thus the model better handles low Re and near-wall flows. At high Re, Equation (9) with

*C*= 0.0845 gives slightly lower values of

_{μ}*ν*, compared to Equation (6).

_{t}*k*-

*ε*model, the transport equation for

*k*is the same as in the standard

*k*-

*ε*model. The transport equation for

*ε*is different, given bywhere

*C*

_{1}= max[0.43,

*η*/(

*η*+ 5)];

*η*is defined the same way as in the RNG

*k*-

*ε*model. The eddy viscosity is computed using the same formula (Equation (6)) as in the standard

*k*-

*ε*model. However, instead of a constant value for

*C*, the coefficient is computed from

_{μ}*C*= 1/(

_{μ}*A*

_{0}+

*A*

_{s}

*U**

*kε*), where , and

*Ω*= (∂

_{ij}*u*/∂

_{i}*x*− ∂

_{j}*u*/∂

_{j}*x*)/2 is the mean rate-of-rotation tensor. The closure coefficients for the realizable

_{i}*k*-

*ε*model are

*A*

_{0}= 4.04,

*A*

_{s}= , , . Compared to the standard

*k*-

*ε*model, the Realizable

*k*-

*ε*model is known to improve predictions of the spreading rate of planar and round jets, boundary layers under strong adverse pressure gradients, flow separation and recirculation. For more details about the rationale on the choice of the closure constants, parameters and their optimized values, refer to Jones & Launder (1972), Yakhot & Orszag (1986) and Shih

*et al.*(1995).

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, *α*_{2} = 1 – *α*_{1}; *α*_{2} = 1 means that the cell is occupied completely by water, whereas *α*_{2} = 0 means that the cell is occupied completely by air. The free surface was discerned using the criterion *α*_{2} = 0.5. The volume of fluid method (Hirt & Nichols 1981) is efficient and convenient for modelling two-phase flow.

### Model domain and boundary conditions for numerical computations

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*_{1} = *u*_{2} = *u*_{3} = 0. Conditions imposed at the other boundaries were: (1) at boundary 1 (*x*_{2} ≤ *η*_{1}), *u*_{1} = *U*_{1}, *u*_{2} = *u*_{3} = 0 and *α*_{2} = 1; *p* was relative hydrostatic pressure. The upstream reservoir water level *η*_{1} was specified (*η*_{1} = 0.964 m). The value *U*_{1} was specified such that it matched the discharge (0.025 m^{3}/s) used in the experiments; (2) at boundary 2 (*x*_{2} ≤ *η*_{2}), *u*_{2} = *u*_{3} = 0, *α*_{2} = 1 and *p* = 0; *u*_{1} was extrapolated from *u*_{1} values for the adjacent interior cells. The downstream reservoir water level *η*_{2} was specified (*η*_{2} = 0.612 m); (3) at boundaries 3 and 5, *u*_{1} = *u*_{3} = 0 and *p* = 0; *u*_{2} was extrapolated from *u*_{2} values for the adjacent interior cells; (4) at boundary 4 (*x*_{2} > *η*_{2}), *u*_{2} = *u*_{3} = 0, *α*_{1} = 1 and *p* = 0; *u*_{1} was extrapolated from *u*_{1} values for the adjacent interior cells; (5) at boundary 6, *u*_{1} = *u*_{3} = 0, *α*_{2} = 0 and *p* = 0; *u*_{2} was extrapolated from *u*_{2} values for the adjacent interior cells, and was limited to *u*_{2} ≤ 0; 7) at boundary 7, *u*_{2} = *u*_{3} = 0, *α*_{2} = 0 and *p* = 0; *u*_{1} adjusted itself to achieve a total of zero mass flow rate.

### Procedures and accuracy of numerical computations

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^{−6}.

## RESULTS AND DISCUSSION

### Mesh independence of numerical results

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.

*h*

_{1}= 2.25 m and

*h*

_{2}= 0.429 m are shown in Figure 2, where

**v**(Equation (3)) vectors are projected onto the plane

*x*

_{3}=

*b*/2, along with contours of the projected vectors. The flow converged from the upstream reservoir toward the siphon entrance, separated from the upper surface of the siphon upper leg and from the lower surface downstream of the crest, and had a maximum speed of 7.44 m/s in the vicinity of the crest (Figure 2(a)). This maximum speed was more than four times the cross-sectionally average flow speed of

*v*

_{s}= 1.75 m/s in the along-siphon direction, and larger than the characteristic velocity

*V*= 5.98 m/s. The two zones of flow separation, delineated by the zero- and negative-value contours (Figure 2(b)), were elongated, where fluids circulated in elongated eddies and reversed direction (

*u*

_{1}< 0). Further downstream, the flow separated downstream of the deflector and discharged as a submerged jet.

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.

### Suitability of turbulence closure scheme

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.

*k*-

*ε*model and Realizable

*k*-

*ε*model, respectively, for turbulence closure, were carried. The results for SK4, RL4 and RN4 were compared. These three runs used the same setup and conditions, except the turbulence closure model. Examples of computed velocity and relative pressure distributions are shown in Figure 5. RN4 and SK4 produced consistent vertical profiles of the horizontal velocity and relative pressure above the crest. The velocity profiles (Figure 5(a)) over the lower half of the siphon throat depth showed strong non-linearity and resembled the potential flow profile of Dressler (1978). Moreover, RN4 and SK4 predicted expected flow separation and reversal near the crown. RL4 gave a profile with a quasi-linear (or weakly non-linear) decrease in velocity from the peak value near the crest to small positive values near the crown; the profile did not show flow reversal in the boundary layer near the crown. RL4 predicted a smaller pressure drop, compared to RN4 and SK4 (Figure 5(b)).

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.

### Flow characteristics above the crest

*u*

_{1},

*u*

_{3},

*p*and

*k*at three

*x*

_{3}locations in cs2 are shown in Figure 6. The head was

*Δh*= 1.821 m. In Figure 6(a), all profiles of

*u*

_{1}showed a shape of potential flow (Dressler 1978) over the lower half of the throat depth, with a decrease in the peak value of

*u*

_{1}close to the solid sidewall at

*x*

_{3}= 0. Over the upper half,

*u*

_{1}was less strong at the middle width (the dashed curve) than both sides (the other two curves); reversed flow occurred near the crown at the middle width, but not on either side. These differences were some aspects of the complicated 3D flow.

The complexity of flow is further demonstrated by vertical profiles of *u*_{3} (Figure 6(b)), which is a component of SF. The solid line showed a zero value of *u*_{3} 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*_{3} > 0) over the lower throat depth and in the opposite direction (*u*_{3} < 0) over the upper throat depth. Near the sidewall, the dotted-dashed curve showed changing directions of *u*_{3} from the upper to the central and from the central to the lower portion of the throat depth. In cs2, *u*_{3} had a maximum magnitude of as high as 0.36*v*_{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*_{1} 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*_{3} (Figure 6(b)). There was little turbulence in the portion of the throat depth where the profiles of *u*_{1} resembled the potential flow solution.

### Pressure distribution

*p*contours in the POS for

*Δh*= 1.821 and 0.579 m are plotted in Figure 7. In both cases,

*p*started to drop as the flow approached the siphon entrance and dropped to the lowest in the vicinity of the crest. The lowest

*p*occurred near the crest due to increased velocity and raised elevation. The increase was more significant in the larger than the smaller

*Δh*case.

*p*was larger at the entrance but dropped to a more severe level near the crest in the former than the latter case.

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.

### Secondary flow

**v**(Equation (5)) onto a cross section gives SF. The SFs in cs1–cs3 (Figure 3) are plotted in Figure 8, where the sidewall is at

*x*

_{3}= 0, and

*n*is the normal distance from the intersection of the cross section with the siphon's lower surface. Like reversed flow in a separation zone, SF gained energy at the expense of the head Δ

*h*, contributed no discharge and thus reduced discharge efficiency. The SF at cs1 had a moderate strength, the strongest vector having a magnitude 0.27

*v*

_{s}(Figure 8(a)). The strongest magnitude increased to 0.70

*v*

_{s}at cs2 (Figure 8(b)), reached the maximum 1.43

*v*

_{s}at cs3 (Figure 8(c)), and decreased to 0.96, 0.67 and 0.31

*v*

_{s}at cs4–cs6, respectively. SF was relatively strong near the siphon's upper surface in the crest region, which is also the case at some distance downstream. Further downstream, strong SF also appeared near the lower surface.

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.

### Turbulence kinetic energy

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.50*v*_{s} and 0.54*v*_{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.67*v*_{s}. Further downstream, dropped to 0.55*v*_{s} and 0.42*v*_{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.

## DATA COMPARISON AND APPLICATIONS

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.

## CONCLUSIONS

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.

## ACKNOWLEDGEMENTS

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

## DATA AVAILABILITY STATEMENT

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

## CONFLICT OF INTEREST

The authors declare there is no conflict.