Numerical Research of Flows into Gullies with Different Outlet Locations

Gullies are sewer inlets placed in pavements usually covered by bar grates. They are the most common linking-element used to drain a wide range of flows from surface runoff into the buried drainage system. Their hydraulic behavior and their overall hydraulic performance is dependent on the flow conditions, the gully dimension, geometry, and location of the outlet device. Herein a numerical research based on Volume Of Fluid ( V O F ) to detect the interface, and on the Shear Stress Transport S S T k - ω turbulence model was conducted to study the importance of the outlet location and characterize flows through them in drainage conditions. Results provided detailed information about flow features, discharge coefficients, and efficiencies for different outlet locations. The authors identified three different regimes, R 1 , R 2 , and R 3 , and concluded that the outlet location influences the velocity field along the gully, the discharge coefficient, and the drainage efficiency. This allows for the estimation of uncertainty and its variation for different outlet positions.


Introduction
In the past decade, there has been an increasing use of dual-drainage models, linking storm water pipe flow (the buried system or minor system) to the overland flow (the major system). Models of the minor system increasingly allow for detailed modelling of the network through representing all manholes, combined sewer overflows, outlet devices, control facilities, among others, usually by defining discharge coefficients [1,2]. On the other hand, some particular structures and devices, such as gullies and manholes, have been analyzed using either experimental investigation [3,4] or computational fluid dynamics (CFD) [2] in order to understand flows and search for parameters that characterize such flows. The understanding of the underlying physics and the capability to model these physical phenomena allowing to adequately predict the flow field is of paramount importance to perform numerical environmental studies [5]. This detail has not been extended effectively to the interface devices between the overland drainage system and the buried drainage infrastructure, although it is well known that flooding can be caused or aggravated by insufficient flow capacity of the inlet devices [6][7][8]. The hydraulic capacity of gullies has been related with factors such as the flow depth on the surface, surface flow conditions, grate type and clogging conditions, inlet area and dimensions, geometry and slope of streets and gullies [8][9][10][11]. Factors such as the local geometry and the location of the gully outlet are also known to influence the flow. CFD is becoming increasingly used as part of simulation schemes for hydraulic structures in urban drainage systems [2,[11][12][13].
In particular, Beg et al. [2] tested different turbulence models and compared them with experimental measurements, concluding that Shear Stress Transport (SST) k-ω and Renormalization Group (RNG) kturbulence models reached best accuracy. Also Lopes et al. [13] presented simulations of a specific gully using SSTk-ω with a high accuracy.
Different gully geometries may be found throughout the world. It seems that their geometry and details follow some regional traditional criteria and that sometimes they are constrained by the street and sewer system location. Figure 1 illustrates two gullies, where the connections to the buried system are vertical and located in two different positions along the gully longitudinal axis. The inlet efficiency of gullies is not well known, and their discharge coefficients depend mainly on the geometry, the hydraulic head, and the velocity field [14]. In that work, the authors found that the initial conditions have a great influence in the time to attain steady conditions. They concluded that the outlet position can influence discharge coefficients. This work aims to study the hydraulic behavior of gullies in usual drainage conditions featuring different locations of the vertical outlet pipe that connects the gully to the buried drainage system, taking into consideration different outlet dimensions, as well as a large range of discharge flows, and 2D and 3D analyses. The 3D model constructed using OpenFOAM R code was validated in [10,11,15,16], and [13] by comparing with experimental measurements. Different gully models for different outlet locations and dimensions were now constructed and those were used to simulate a wide range of discharge flows. We present the qualitative description of the gully flow in steady and unsteady conditions as well as the quantitative analysis of the following parameters along the gully: water depths, streamlines, velocity, pressure fields, and the inlet discharge capacity.

Numerical Model
The numerical model is based on the Navier-Stokes equations/Reynolds-Averaged Navier-Stokes (RANS) equations governing the motion of the 3D incompressible and isothermal flows in which the free surface is described using a Volume-Of-Fluid method (VOF). According to this description [17,18], the VOF-function, F = F(x, y, z, t), ranging from 0 to 1 and corresponding respectively to cells without water and full occupied by water, is included in the mass and momentum conservation equations. It is also updated using an advection equation for F. Some improvements of VOF models were developed including surface tension and the interface curvature, as well as the artificial compression of the interface, to improve accuracy of the interface [19]. The VOF method used in the interFoam solver implemented in the OpenFOAM R (Equations (1)-(3)) has two particularities: a volumetric surface force, explicitly estimated by the Continuum Surface Force (CSF) function of the surface tension, and the interface curvature, which are included in the momentum equation [20]; the compression of the interface is achieved by introducing an extra, artificial compression term in the advection equation [19,21]. ∇ ·ū = 0 (1) whereū is the mean velocity vector, p * is the modified pressure adapted by removing the hydrostatic pressure from the total pressure, α is the VOF F function, t is the time, ρ is the fluid density, g is the acceleration due to gravity, τ is the shear stress tensor, f is the volumetric surface tension force (where CSF and interface curvature are included) andū c is the compression velocity. Further in the mass and momentum conservation equations, VOF-function is included through physical properties such as density and viscosity, which are defined by a weighting of the values for air and water. Lopes et al. [13] developed an air-entrainment model. The new solver, airInterFoam, considers the air entrainment triggered by an additional advection equation for dispersed gas phase. Lopes et al. [13] thus simulated a gully with air and water flow, led to conclusion on whether to consider the air-entrainment or not. Although important air-entrainment occurs, it revealed small influence on the hydraulic performance of the gully. They also used a turbulence model, where turbulence variables were calculated with the SSTk-ω turbulence model. This is known for the best combination of two Reynolds-Averaged Simulation (RAS) formulations, using high-Reynolds-number formulation of kmodel for the free-stream region and taking advantage of the accuracy and robustness of k-ω model in the near-wall zone. This study follows the Lopes et al. [13] methodology. However, the air-entrainment model was not used as it is not needed for the detailed requirements, saving computer time.
Total Variation Diminishing (TVD) limited form of central-differencing is used for convective terms in momentum equation. The Van-Leer scheme is used for the convective term in VOF-advection equation and 'Interface Compression' scheme is used in order to bound the solution of the compressive term between 0 and 1. To ensure boundedness of the phase fraction and avoid interface smearing, the solution of the VOF equation is done with the Multidimensional Universal Limiter for Explicit Solutions (MULES). The Pressure-Implicit with Splitting of Operators (PISO) procedure proposed by [22] is used for pressure velocity coupling in transient calculations with 3 loops. We used the same discretization schemes employed in previous works on gullies tested with positive outcome [10,13] as well as we used turbulence models with the best accuracy [2].
The computational domain of the flow, aiming to represent a length (L) × width × height = 0. The gully was placed between 1.2 m and 1.8 m of the street (1.2 m < x < 1.8 m). The mesh analysis was done and presented in [13]. We tested three meshes. The two finer meshes showed equivalent vortex details, therefore we choose an intermediate mesh size, which retains the main features of the vortices, while also keeping the calculation time acceptable. We constructed the 3D geometry using blockMesh utility from OpenFOAM R adapting the work of [10,13]. Figure 2 presents the longitudinal cross sections of the eight gully configurations studied indicating the outlet vertical pipe position center (P op ). We consider five configurations for an outlet pipe of diameter D op = 80 mm (Configurations 1 to 5, which corresponds to five different locations from upstream to downstream) and three configurations for an outlet pipe of diameter D op = 200 mm (Configurations 6 to 8, which corresponds to three different locations from upstream to downstream). Table 1 summarizes the outlet pipe characteristics for each configuration.       To simplify the mesh construction, we considered xx-axis parallel to the channel bottom instead of it being horizontal. Thus, the acceleration due to gravity, instead of being considered vertical, was set to g = (0.00017; 0; −9.81200) to consider a 1% slope. The different configurations of the gully outlet were simulated for different flow rates (3 discharges for Config. 1 to 5: Q10, Q20, and Q50) and 4 discharges for Config.6 to 8: Q20, Q50, Q100, and Q200).
The upstream inflow boundary condition (BC) defined as Dirichlet-BC, at x = 0, was set according to supercritical uniform conditions in a hypothetical 0.5 m wide channel and 1% slope street. Table 2 presents the inflow conditions for the different discharge flows in supercritical uniform conditions, supercritical uniform height, h in , and uniform velocity, U in , which was considered constant along the depth, as well as Froude and Reynolds numbers ( with µ being the viscosity). The inlet and outlet at the right boundary (downstream channel) were defined by a specific gradient for dynamic part of pressure (p * = p total − p hydrostatic ). In the outlet at the bottom (outlet pipe) a hydrostatic pressure was assumed. The top boundary considered hydrostatic pressure and specific gradient for VOF and velocity. The remain boundaries were considered as walls, imposing zero velocity in the vicinity of the face using Dirichlet-BC. Table 2 were run to understand the flow behavior during surcharge transient conditions. It was found that 15 s is enough to achieve the convergence of the main flow properties in steady conditions (gully and downstream water depth and outlet flows, as well as the total water volume in the domain). However, for low flows, the flow continues to decrease progressively, making it no longer interesting to continue with the simulations.

Free-Surface and Velocity Field
Detailed free-surface and 2D flow velocity pattern in central plane over time is illustrated in Tables 3-5, showing the following three hydraulic behaviors, respectively: (1) the water volume starts decreasing immediately (Config. 1 to 5 for Q10, Config. 6 to 8 for Q20 and Q50); (2) the water volume remains relatively stable since the beginning (Config. 1 to 5 for Q20 and Config. 6 to 8 for Q100); (3) the water volume increases up to a point that the downstream water depth smoothly tends to stabilize; (Config. 1 to 5 for Q50 and Config. 6 to 8 for Q200).
Tables 6 and 7 illustrate 3D velocity field for all simulations under steady conditions for all configurations and discharges. We calculate average results from 15 s ≤ t ≤ 20 s. It is clear that while for some simulations the flow entering the gully is mostly discharged by the bottom outlet, in other situations the water flows over the gully suggesting that the gully efficiency depends not only on the dimensions but on the location of the outlet pipe. For Q200, Config. 8 (D op = 200 mm) presents lower water depth downstream of the gully and for Q100, Config. 7 and 8 (D op = 200 mm) present more stable condition than Config. 6.     It was found that the natural direction of the jet, as well as the downstream channel water depth can change over time for a constant flow discharge. For Behavior 1, the gully initially (t = 0.5 s) shows to be drowned especially under Config. 7 and 8, and a considerable amount of inflow goes downstream causing large variation of water depth. For Behavior 2, in configurations where the outlet pipe is placed in an upstream or central position (Table 4), it is possible to distinguish 3 phases: 1. at the beginning, the inflow causes a complex velocity field in the gully and the jet tends to reach the gully downstream of the outlet pipe (Table 4, t = 0.5 s); 2. the flow through the outlet pipe drives the jet into the gully bottom outlet controlling the vortices placed upstream and downstream of the outlet pipe. As the main flow is deflected to the outlet pipe, vortices change dimensions (Table 4, t ≤ 5 s); 3. further, in almost configurations a hydraulic jump forms at surface above the outlet pipe, also pushing the jet into the outlet pipe. Thus, the discharge capacity increases and the velocity field in the gully downstream stabilizes ( Table 4, t ≥ 15 s). When the hydraulic jump disappears, it forms a step wave which is responsible for a higher water depth downstream.
While this happens, in Config. 3 and 7, the jet tends to move slightly in the outlet pipe which induce the variation of discharge capacity and the water depth in the gully until equilibrium is attained with oscillating characteristics (Table 4). However, this behavior with oscillatory characteristics cannot be seen in Config. 6 and 8 nor in Config. 1, 2, and 4 to 5. The hydraulic jump occurs for D op = 80 mm and Q20 in Config. 1 to 4 and for D op = 200 mm and Q100 in all configurations (Config. 6 to 8). However, for Config. 1 to 5 a stable hydraulic jump only occurs for Q20 and Config. 1 to 3. For Behavior 3 (Table 5, simulations for D op = 80 mm and Q50, as well as for D op = 200 mm and Q200) when the outlet pipe is placed at the center or upstream, the main flow is not deflected to the outlet pipe. In the first seconds, a wave is formed, which increases the water depth downstream. A substantial part of the flow is drained to channel downstream and there is not a hydraulic jump formation. Therefore, the flow through the outlet pipe always comes from the gully downstream area. Similar behavior is observed for simulations from Config. 1 to 5 for Q20 and Q50 and from Config. 6 to 8 for Q100 and Q200.

Velocity and Pressure Profiles
To map the flow and quantify the water depths, velocity, and pressure, values were evaluated at the gully boundaries and top, as well as at the central profile. Figure 3 illustrates velocity and pressure profiles at the gully bottom and top sections (z = 0.3 and z = 0.6 m) for the different configurations and discharge flows. Figure 3 demonstrates that for all the outlet locations with diameter 80 mm, independently of the flow rate, the bottom pressure (Figure 3a), the bottom velocity profile (Figure 3c) and the pressure at the top (Figure 3e) show the same trend. However, while for discharge Q10, the flow at the outlet pipe is drained through upstream side of the bottom outlet pipe, for Q20 and Q50, the flow is drained from downstream side to the outlet pipe (maximum velocity at the right).
At the top entrance (Figure 3g), different patterns can be detected for the velocity longitudinal profile: (i) outlet pipe located upstream with lower discharges presenting velocity into the gully in upstream area and smaller velocity variation along the gully; (ii) outlet pipe located in center or downstream with larger discharges presenting velocity into the gully in upstream area and velocity out of the gully in the gully center; and (iii) outlet pipe with larger discharges presenting small velocity upstream, as well as the highest values for velocity downstream. Considering all the results, it is clear that the 200 mm outlet pipe in the different locations (Figure 3b,d,f,h) present larger variations along the gully than the 80 mm outlet pipe in the different locations. Pressure values at the gully bottom are larger downstream the outlet for the largest discharges. Different discharges induce different maximum velocity locations into the outlet pipe and different velocity profiles for each location and discharge. The flow in higher discharges for each outlet pipe diameter comes from downstream, similarly in D op = 80 mm and D op = 200 mm outlets. Figure 3f shows different pressure profiles at the gully top entrance for different locations and discharges. Pressure values on the gully center are higher for larger discharges and upstream outlet pipe location.
Higher pressure values at the right wall (large water depth) are presented for larger discharges and location at the center and downstream. On the left (at the gully upstream wall) a separation zone is verified as negative pressure occurs for almost simulation except for higher discharges and outlet pipe location in center or downstream. Figure 3h (velocity profiles at the gully top) shows a very smooth concavity of free surface for the Q200 whereas in the Q100 there is an abrupt rise, which allows the identification of a hydraulic jump occurrence.

Flow Behavior and Characterization
Based on the methodology proposed by [3] and also by [23,24] for drop manholes, we propose the description of transitional flow in a gully organized in 3 different regimes, and based on three aspects, related but independent: (i) the natural jet direction and its relationship with the location of the outlet pipe; (ii) the direction of the streamlines around the outlet pipe location and (iii) the relationship between the discharge in the channel and the outlet pipe discharge capacity. The new regimes which are represented in Figure 4 will be described taking into account the unsteady behavior of the flow in the gully, since the regimes for a specific flow changes according to the water depth in the gully. Regime R1 occurs mainly for smaller discharge flows and locations more downstream as the natural jet trajectory tends to reach upstream outlet-The main flow direction hits the box floor upstream of the bottom outlet pipe, and the streamlines are near horizontal, either along gully bottom in the upstream part, R1a, or just near the outlet, R1b, causing small discharge capacity. If the discharge in the channel is lower than the bottom outlet pipe discharge capacity, the water depth and the water volume decrease and streamlines follow the bottom wall, and the regime depends only on the jet trajectory, R1a (e.g., Config. 6 and 8, D op = 200 mm for Q20 and Config. 8, D op = 200 mm for Q50, Tables 3 and 7). If the discharge flow in the channel is similar to the outlet capacity, the main flow is driven to the outlet pushed by the vortices located at the gully box upstream part, and it hits the outlet pipe boundary and a contraction is seen, causing asymmetry in the streamlines entering the pipe, R1b (e.g., Config. 7, D op = 200 mm for Q50, Config. 7, D op = 200 mm for Q100, Tables 4, 6 and 7).
Regime R2 occurs for intermediate discharge flows, when they are similar to the outlet discharge capacity-The main flow direction hits the outlet pipe, the streamlines are almost symmetric and discharge capacity is maximum. If there is water in the gully, the velocity field in it is influenced by the jet into the outlet pipe and upstream, downstream, or both vortices. This is possibly due to the trend of the jet trajectory or the equilibrium in upstream and downstream vortices (e.g., Config. 5, Regime R3 occurs when the natural jet trajectory overshoots the outlet downstream, for larger discharge flows or outlet pipe gullies located upstream-The main flow direction tends to hit the area downstream of the outlet pipe, and the streamlines close to it become non-symmetric, which causes a decrease of the discharge capacity. If the outlet pipe capacity is similar or larger than the upstream discharge flow, the jet streamline is deflected and guided into the outlet pipe. Therefore, a vortex could be created downstream, providing an increase of the discharge capacity, R3a (e.g., Config. 1 to 4, D op = 80 mm for Q20; Config. 6, D op = 200 mm for Q100, Config. 6, 7, D op = 200 mm for Q200 and Config. 8, D op = 200 mm for Q200, and for t ≥ 10 s, Tables 4-7). In this case, during discharge, if the outlet discharge is not enough to discharge the flow from the upstream channel, the water depth and the water volume may initially increase and thus vortices on the gully may be formed deflecting the jet into the outlet pipe, thus, increasing its capacity, R2 could be attained. Alternatively, if the jet tendency is the gully wall to far downstream from the outlet pipe, the vortices upstream are unable to deflect the flow into the outlet pipe, and thus most of the upstream discharge will flow downstream: R3b (e.g., Config. 1 to 4, D op = 80 mm for Q50).
Apart from the regime description, some additional behaviors must be noted from a theoretical standpoint: R1 and R2 should not be physically possible on Config. 1; however, R2 is verified since the flow is forced to curve into the outlet for lower discharges; R3a could be defined for Config. 5, since downstream vertical wall forces the flow into the outlet.

Discharge Coefficients
Discharge coefficients, depending on the bottom outlet pipe diameter and water depth inside the gully, can be evaluated considering either orifice (Equation (4)) or weir (Equation (5)) formulas: where subscripts 'o' and 'w' represent an orifice and a weir variable, respectively; Q o and Q w represent the discharge flow through outlet pipe (either considering orifice or weir), A is the area of the orifice, b w is the channel width, g is the acceleration due to gravity, and, ultimately, H in is the uniform height of the flow upstream the weir/orifice (in the upstream channel). Efficiency is obtained by calculating the ratio between Q o or Q w and Q in , which is the sum of Q outlet_bottom and Q outlet_channel , presented in Tables 8 and 9. Tables 8 and 9, apart discharges, present discharge coefficients considering weir and orifice formulas and bottom outlet efficiency. Figure 5 illustrates the Tables 8 and 9 discharge coefficients and efficiency data. As expected the higher outlet pipe diameter drain higher flows and largest efficiency is found for lower flows. However, for lower flows and Behavior 1, as defined at the beginning of Section 3 and for each outlet pipe diameter, all discharge coefficients and efficiency are similar. Looking at the coefficients and efficiency values of gullies with different outlet locations, Figure 5, the largest efficiency is found for the center and downstream location (95%). In intermediate discharges, those maintaining volume in the gully (Behavior 2), the higher discharge coefficient and efficiency is also found for the central bottom outlet. For higher discharges, flood conditions, it is found for the downstream position, which performs better. Calculating the differences, 26% is found to be the maximum uncertainty range occurring in the D op = 200 mm and Q200 between Config. 6 and 8 related with Config. 7 (c), followed by 25% for D op = 80 mm and Q50 between configurations 3 and 5 related with Config. 3 (c) and 13% for D op = 200 mm and Q100.

Conclusions
An urban drainage gully with 0.6 × 0.3 × 0.3 m 3 (length × width × height) placed in a channel with 1% slope was modelled to evaluate the importance of the bottom outlet pipe location along the gully longitudinal axis. Five different flow rates and two sizes of the outlet pipe were used in the analysis. Different 3D models based on Navier-Stokes/Reynolds equations, VOF method and SSTk-ω were constructed and used in the OpenFOAM R platform, which showed ability to predict complex flow features, air-water interface, and turbulence characteristics. The study was motivated by the increase of applications of dual-drainage models that consider the performance of linking elements between the surface and buried systems, as well as by the need to take uncertainty into account. The different locations of the bottom outlet pipe lead to different limitation of discharge and influence the stability of the free surface and the water depth along the gully and downstream channel. Three regimes were identified based on the velocity field, which is the main responsible for the hydraulic behavior and also determines the symmetry of the streamlines at the outlet and thus the contraction and the discharge capacity. Besides the initial water depth in the gully, the velocity field and, thus, the regime depend on three main factors: (1) the outlet location; (2) the natural jet direction and the relative distance between the area where it hits the gully bottom and the outlet pipe location, which influence streamlines around the bottom outlet pipe; (3) the relationship between the discharge in the channel and the outlet pipe discharge capacity. The largest difference of 26% was found for the same outlet diameter and the same upstream flow with different outlet pipe locations corresponding to largest flow and gully pipe bottom outlet diameter (26% for D op = 200 mm and Q200 related with Config. 7, followed by 25% for D op = 80 mm and Q50 related with Config. 3), showing that unknown outlet pipe location creates a significant range uncertainty in flooding conditions.