Numerical studies regarding the influence of entrapped air on the hydraulic performance of gullies are nonexistent. This is due to the lack of a model that simulates the air-entrainment phenomena and consequently the entrapped air. In this work, we used experimental data to validate an air-entrainment model that uses a Volume-of-Fluid based method to detect the interface and the Shear-stress transport k-ω turbulence model. The air is detected in a sub-grid scale, generated by a source term and transported using a slip velocity formulation. Results are shown in terms of free-surface elevation, velocity profiles, turbulent kinetic energy and discharge coefficients. The air-entrainment model allied to the turbulence model showed a good accuracy in the prediction of the zones of the gully where the air is more concentrated.

INTRODUCTION

Towards efficient performance of urban drainage systems, the study of the hydraulic behaviour of their elements is of crucial importance. Gullies, also known as storm inlets, are the most common linking element to convey the water from surface runoff to the pipe networks, with multiple designs and applications (Figure 1(a)1(c)). In the United Kingdom, drop manholes can be used as substitutes for gullies (Figure 1(d)). Several studies regarding the flow characteristics of pipe networks and the surface runoff can be found in the literature; however, there are few works on linking elements due to the numerous costs of constructing such experimental facilities and the difficulty of performing calibration/validation of numerical models.
Figure 1

Some examples of linking elements: (a) grated inlet close to a side walk (Coimbra, Portugal), (b) curbside inlet (Coimbra, Portugal), (c) continuous grated inlets placed transversely to the street (Coimbra, Portugal) and (d) grated drop manhole in a park (Sheffield, United Kingdom).

Figure 1

Some examples of linking elements: (a) grated inlet close to a side walk (Coimbra, Portugal), (b) curbside inlet (Coimbra, Portugal), (c) continuous grated inlets placed transversely to the street (Coimbra, Portugal) and (d) grated drop manhole in a park (Sheffield, United Kingdom).

Djordjević et al. (2013) used the 3D Computational Fluid Dynamics (CFD) OpenFOAM® code to model a United Kingdom (UK)-type gully to investigate the interactions between surface flood flow and surcharged pipe flow. Carvalho et al. (2011) used a 2D-V Volume-Of-Fluid (VOF) in-house numerical model to simulate a typical Portuguese gully (Figure 1(a)) under drainage conditions, work that was further extended in Carvalho et al. (2012) to cover both drainage and surcharge conditions using the OpenFOAM® model. An experimental investigation of turbulence characteristics and mean flow properties was made by Romagnoli et al. (2013) for a gully in surcharge conditions. More detailed work on both the drainage and surcharge performance of a gully was made by Leandro et al. (2014b) by means of the characterisation of a 2D middle plane of the gully and a qualitative interpretation of the entrapped air. Martins et al. (2014) conveyed numerical and experimental investigation of a gully under drainage for a wide range of inflows in order to understand the hydraulic performance of this device and investigate the drainage coefficients. Similar work was extended by Lopes et al. (2015) for surcharge conditions. Experimental works on drainage efficiency also exist for grated inlets (Russo & Gómez 2011; Sabtu et al. 2016) and continuous transverse gullies (Gómez & Russo 2009; Russo et al. 2013; Lopes et al. 2016a). An example of a continuous transverse gully can be seen in Figure 1(c).

The admission of air in the drainage systems is due to multiple sources as hydraulic jumps, turbulent free-surfaces, or dropshaft elements (gullies and manholes), or it is admitted through air vacuum valves. Its presence may cause damage and considerable operational disruption, but also contributes to some beneficial aspects. The downsides are essentially related to the appearance of air pockets in pluvial drainage systems, which can increase the roughness of the conduits through corrosion of metallic components, reduce pump and turbine efficiency with associated growth of electrical consumption, and reduce the drainage capacity (Pothof & Clemens 2010), which can occur both in the drainage elements and in the conduits. The movement of air in pressurized conduits and its sudden release can cause high-pressure spikes (Ramezani et al. 2016) and consequent disruption of the systems. On the other hand, deficient ventilation of the drainage elements leads to negative pressures and an increase in the pool depths of drop manholes (Figure 1(d)) (Granata et al. 2014), which might happen in gullies on a smaller scale; whereas scarce aeration will diminish the capacity of the sewer to auto-treat the residuals due to the low levels of dissolved oxygen. The latter is even more significant in the case of a combined sewer overflow. Furthermore equilibrium of air quantity is needed.

On the basis of the previous results and investigations, this work aims to verify the accuracy of an air-entrainment model and recognize the influence of the entrapped air on the hydraulic performance of a gully under a specific drainage condition. A new numerical solver implemented in the CFD toolbox OpenFOAM® (named airInterFoam) is used to predict the air-entrainment mechanism at the free surface and to detect the amount of air that is trapped inside the gully box. The free surface is calculated using the VOF technique and the turbulence statistics with a Shear-stress transport (SST) k-ω turbulence model. The next section describes the experimental apparatus of a real scale gully box located at the University of Coimbra. The following section describes the numerical model used amongst the mathematical formulations. The methodology used in this work is presented in the subsequent section, then the next section describes and discusses the results. The final section summarizes the work.

EXPERIMENTAL APPARATUS

Experimental tests were carried out in a full scale gully box placed in a flume inside a multipurpose channel at the Hydraulic, Water Resources and Environment Laboratory of the University of Coimbra. Since the experimental results have been published in Leandro et al. 2014b, only a brief summary is provided here. The water is initially pumped from an underground tank to the multipurpose channel and returned in a closed circuit. The flow is controlled by four butterfly valves, associated to each of the four pumps, and measured by an electromagnetic flow meter. Valves are electronically controlled by a SCADA (Supervision, Control and Data Acquisition) system. Inside the multipurpose channel is placed another flume that contains our gully box. This flume has a rectangular shape with 1% bottom slope, and is 8 m long and 0.5 m wide. The gully box has dimensions of 0.6 m length, 0.3 m depth and 0.3 m width (Figure 2) and is placed in the middle of the flume. The dimensions of the gully box chosen are the most representative among the different types of gullies that can be found in Portugal (an example can be seen in Figure 1(a)). At the bottom of the box, a circular hole with a diameter of 0.08 m and an additional 0.05 m pipe, allows drainage of the gully. This pipe diameter is smaller than the one specified in the Portuguese by-laws (0.2 m), however 0.08 m is chosen as it can simulate the worst case scenario of an obstructed outlet pipe during a flood event. Figure 2 is a schematic representation of the experimental facility. The Cut A-A′ represents the centre plane of the gully in which the results are taken.
Figure 2

Sketch of the experimental gully installation. Inside the gully, the spatial grid with uniform 0.03 m spacing is used to measure the velocities. Cut A-A′ represents the centre-plane of the gully. Units in metres.

Figure 2

Sketch of the experimental gully installation. Inside the gully, the spatial grid with uniform 0.03 m spacing is used to measure the velocities. Cut A-A′ represents the centre-plane of the gully. Units in metres.

The grate inlet, normally placed on top of the gully, was removed in the present work. Although this component is known for decreasing the drainage capacity of the gully, the removal was done because our main goal is to verify the accuracy of the air-entrainment model and further compare, against other works (also without a grate), the influence of air on the drainage capacity. Also, grate removal (in gullies and manholes) is a normal procedure during flood events to increase the drainage efficiency.

Flow photographs are captured using a digital Sony Alpha DSLR-A350 camera. The camera is fixed on a tripod, about 2 m in distance outside the channel's side wall, as such a camera lens can be adjusted to 50 mm focal distance to avoid a ‘fisheye’ effect. The perspective distortion in the image is neglected in this work as images are just used to qualitatively define the position of the bubbles inside of the gully. To provide flow illumination, a reflector with 1,000 W light is positioned above the gully. A video of the flow is taken using a Panasonic DMC-FS16 in a frame rate of 24 fps. Instantaneous velocities are measured with an Acoustic Doppler Velocimeter (ADV) probe. The frequency and sampling period are set to 25 Hz and 180 s, respectively.

NUMERICAL MODEL

General flow equations

The flow in the gully is a complex problem to simulate, due to the highly complex three-dimensional rotational flow and the turbulent surface. As such, the solution cannot be achieved using shallow water approaches, but requires a full description of the flow properties just possible with the solution of the Navier-Stokes equations. For an incompressible and isothermal fluid, the mean flow properties can be achieved by the solution of the incompressible Reynolds-Averaged Navier-Stokes equations: 
formula
1
 
formula
2
where is the fluid density, the mean velocity vector, the gravitational acceleration, the shear stress tensor, p the pressure, the surface forces and t the time. The decomposition of the viscous stress term is given by constitutive relation: 
formula
3
where is the dynamic viscosity given by the sum of the mean and turbulent dynamic viscosity () admitting Boussinesq hypothesis. In this model, modified pressure () is used, 
formula
4

Free-surface position

The free-surface position is calculated using the interface capturing VOF method. The original VOF (Nichols & Hirt 1975; Hirt & Nichols 1981) is based on the introduction of the advection equation (Equation (5)). This equation introduces an indicator function that records the volume fraction of each fluid, 
formula
5
Further, a distinct advection technique is used to control the numerical diffusion of the interface while ensuring boundedness and conservation of the phase fraction (Rider & Kothe 1997; Carvalho et al. 2008). Many alternatives to the initial donor-acceptor scheme (Nichols & Hirt 1975; Hirt & Nichols 1981) have been proposed as the geo-reconstruction techniques: SLIC (Noh & Woodward 1976), PLIC (Youngs 1984) and isoAdvector (Roenby et al. 2016) and high-resolution compressive schemes: flux-corrected transport, Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) (Ubbink & Issa 1999) and High Resolution Interface Capturing (HRIC) (Muzaferija et al. 1998). Gopala & van Wachem (2008) and Waclawczyk & Koronowicz (2008) provide a useful comparison between different advection schemes.
The VOF method used in this work, and implemented in the OpenFOAM® as interFoam solver, has slight modifications to the original VOF. The free-surface boundary condition (BC), which was implemented in the original VOF in order to solve just one phase, is removed and the two phases are now solved together. This is possible due to the implementation of a volumetric surface force function (), explicitly estimated by the continuum surface force (Brackbill et al. 1991), which is added to the momentum equation (Equation (2)), 
formula
6
where is the surface tension and the interface curvature. The physical properties () density and viscosity need to be further defined by a weighting of the values for air and water, 
formula
7
Finally, in OpenFOAM®, the compression of the interface is achieved by introducing an extra, artificial compression term (Rusche 2002; Weller 2008) into the VOF equation (Equation (5)). This term is given by: 
formula
8
The sub-term ensures that interface compression is calculated just where the mixture happens and is known as compressive velocity, acting perpendicular to the interface. The compression term appears as an artificial contribution to the convection of the phase fraction, 
formula
9
The term is the interface unit normal vector, which yields the direction of the compressive velocity. The coefficient is an adjustable coefficient (standard ), which determines the magnitude of the compression.

Turbulence closure

Turbulent variables were calculated with a modified version of the SST k-ω model (Menter 1993) implemented in OpenFOAM. It combines the best of two Reynolds-Average Simulation (RAS) formulations (Blazek 2001; Versteeg & Malalasekera 2007): it takes advantage of the accuracy and robustness of the k-ω model in the near-wall zone, whereas the free-stream region, i.e. the outer part of the boundary layer, is simulated with a high Reynolds number formulation of the k-ɛ model. By operating with the ω-equation next to the wall, this model becomes substantially more accurate in the near wall zone (Menter et al. 2003), does not require wall-damping functions in low-Reynolds number flows (Versteeg & Malalasekera 2007), and improves the accuracy of prediction of flows with strong adverse pressure gradients (Blazek 2001).

Sub-grid air-entrainment model

The air-entrainment model is triggered by solving an additional advection equation for the dispersed gas phase (), 
formula
10
where is a switching function that returns the position of the free-surface and is the bubble advection velocity, which can be calculated by: 
formula
11
in which is the bubble-relative velocity, calculated according to the Clift et al. (1978) model, 
formula
12
where is the averaged bubble radius. The source calculates the rate of dispersed gas that penetrates through the free surface. In this study we adopt the formulation described in the work of Ma et al. (2011), 
formula
13
where is the normal velocity component to the free-surface, a is the amplitude of the surface disturbances formed at the free-surface, is the interface thickness, given by , where L is the characteristic linear dimension (equal to the pipe diameter for pipe flows or the water depth in free-surface flows). In this work, L is equal to the water depth upstream of the gully box () (Figure 2). The symbols 〈 〉 are used to turn the inner parcel to zero if its value is negative. As mentioned in Lopes et al. 2016b, the amplitude of the surface disturbances () is considered as having the same order of magnitude as the radius of the turbulent eddies at the free-surface, calculated according to the formulation for the SST k-ω turbulence model, 
formula
14
where . The transport of the gaseous-single-phase is completed when bubbles present in the water phase are able to break up along the free-surface area. The dispersed gas phase is turned to continuous when the volume of the water at each cell is residual. In this work, we assumed this residual as = 0.1, as this is commonly assigned to be the water volume ratio at the hypothetical free surface (Wood 1991).
The air-entrainment model used in this work is also able to influence the dynamics of the continuous water phase () (i.e. two-way coupling model). This is done by imposing a new conservation law to the volume fractions, guaranteeing the mass conservation for variables and during the calculation of the exceeded fraction of the continuous air phase (), 
formula
15
Finally, the physical properties () of the two-fluid mixture, such as density or viscosity, are weighted in the volume of each fluid, 
formula
16

Mesh, initial and boundary conditions

A regular non-uniform mesh with variable grid sizing space from 0.01 m to 0.04 m is used in this work (Figure 3). The largest cell is placed on top of the domain (Vcell = 64 cm3) whereas the smallest cell volume is generated inside the gully box (Vcell = 1 cm3). The mesh is adapted from the work of Martins et al. (2014), created with the blockMesh utility from OpenFOAM®. The initial condition comprehends a column of water with the same level and velocity as the inlet.
Figure 3

Numerical mesh and BC definition.

Figure 3

Numerical mesh and BC definition.

Five types of BC are applied to this case (Figure 3). The ‘inlet surface’ boundary allows the flow to enter the domain by setting Dirichlet-BC, whereas the ‘outlet surface’ allows the flow to leave the domain by fixing the pressure. The boundary ‘outlet pipe’ has the value of pressure imposed using Dirichlet-BC. The ‘atmosphere’ boundary allows exchanges of air. The ‘wall’ boundary imposes zero velocity using Dirichlet-BC.

Model discretisation and solution methods

As the study of the effect of discretisation schemes in the solution is not the objective of this work, we used the same discretisation schemes as the one employed in previous works on gullies (Lopes et al. 2015; Martins et al. 2014) and tested, with positive outcomes, in many other works (Zhao & Wan 2015; Zhang et al. 2016). The governing equations were discretised in time using an explicit Euler scheme. The gradient terms are discretised using linear interpolation from cell centres to face centres. Different schemes are used to discretise the convective terms. The term of the momentum equation uses the TVD-limited form of central differencing in the ‘V’ (vector-field) version. The ‘V-schemes’ limiter is calculated based on the direction of the most rapidly changing gradient, and is applied to all components of the vectors, resulting in a less accurate but stable solution (OpenFOAM Foundation Ltd 2016). The Van Leer scheme is used for and the upwind scheme is used in . For the second convective term of the VOF equation (Equation (8)), the so-called ‘Interface Compression’ scheme is used in order to bound the solution of the compressive term between 0 and 1.

In order to ensure boundedness of the phase fraction and avoid interface smearing, the solution of the VOF equation is achieved with the Multidimensional Universal Limiter for Explicit Solutions (MULES). In this study, the PISO procedure proposed by Issa (1985) is used for pressure-velocity coupling in transient calculations. The PISO was calculated three times in each time step.

METHODOLOGY

Experimental methodology

The discharge tested corresponds to = 22 l/s and = 0.037 m. Just one flow rate is chosen, as the aim is to verify the accuracy of the numerical model. For this combination, the mean velocity is = 1.18 m/s and Froude is = 1.94 at the free-surface. Unlike the work of Martins et al. (2014), the bottom outlet of the gully is slightly under surcharge to simulate the real conditions of pressurised systems. This surcharge condition corresponds to a rough hydrostatic head of about 0.02 m of water in the bottom outlet pipe, as experienced in Páscoa et al. 2013 and Leandro et al. 2014b (see also Figure 2). This value was calculated based on the free-surface level at the bottom reservoir, and considering the distance to the end of the outlet pipe in the captured video.

The captured video is used to trace the free surface position using a Computational Vision Model (Roque 2011). Video records a total of 720 images. Each image is firstly divided in several horizontal sub-domains. Then, image treatment and segmentation techniques are used to trace the two largest lines found in each sub domain. While one of these lines represents the channel bottom, the other detects the free-surface position, and the distances between the lines' midpoints correspond to the flow water depths for each sub-domain.

A spatial grid defined with uniform 0.03 m spacing was chosen to measure the velocities (Figure 2). The bottom of the grid is situated 0.05 m above the base of the gully box, which is well above the minimum (of 0.02 m) required to avoid erroneous ADV signals that we found when measuring near the acrylic wall. A total of 171 points were measured. This resolution has been shown to be sufficient to characterize the large vortices structures inside the gully box (Leandro et al. 2014b). The ADV signal is post-processed in this work by imposing minimum thresholds on two data quality indicators normally associated with acoustic Doppler data: the correlation (COR) and the signal to noise ratio (SNR). The COR measures the level of similarity between two consecutive pulses of the ADV. The SNR is calculated using the signal amplitude and background noise level. In the present study, we used SNR > 15 and COR > 70, as such values are typically recommended by many ADV manufacturers. According to Wahl (2000), samples with COR < 70 can still be used when the SNR is high and the flow is relatively turbulent. The presence of spikes in the water velocity time series were detected and deleted by means of the phase-space thresholding method proposed by Goring & Nikora (2002).

Numerical methodology

In total, three different simulations are set up to study the 3D gully under drainage conditions (Table 1). Simulation S1 corresponds to the simulation of the gully using the VOF technique to calculate the free-surface calculation and the RANS equation to model the flow characterization. Simulation S2 includes the SST k-ω turbulence model with the previous simulation in order to calculate the influence of the turbulence on the mean flow. Simulation S3 adds the air entrainment model to S2. In order to overcome some numerical instabilities of the air-entrainment model when the water drops into the gully box, S3 was solely activated on top of S2 when this model has simulated 10 s of real time (see also ‘Steady-state’ section).

Table 1

Summary of simulation cases tested in this work

Simulation ID Model equations
 
RANS + VOF model (sec.3.1 + sec.3.2) Includes turbulence model (sec.3.3) Includes air-entrainment model (sec.3.4) 
S1 Yes No No 
S2 Yes Yes No 
S3 Yes Yes Yes 
CPU time (h) 103.94 79.19 102.42 
Simulation ID Model equations
 
RANS + VOF model (sec.3.1 + sec.3.2) Includes turbulence model (sec.3.3) Includes air-entrainment model (sec.3.4) 
S1 Yes No No 
S2 Yes Yes No 
S3 Yes Yes Yes 
CPU time (h) 103.94 79.19 102.42 

As referred to before, since the bottom outlet of the gully is slightly under surcharge, the pressure of this boundary is set to a rough hydrostatic considering a head of 0.02 m of water.

RESULTS AND DISCUSSION

Steady-state

Achieving a steady-state condition in numerical simulations is sensible for initial and boundary conditions, model properties and geometry. In the authors’ previous works, it was shown that running the simulation for 15 s is enough to achieve convergence of the main flow properties (Carvalho et al. 2011; Leandro et al. 2014b; Martins et al. 2014; Lopes et al. 2015). Figure 4 shows that herein, all simulations also converged to the steady state within 15 s of real time. In S2 and S3, the volume of water converged to a much more constant and stable value than in S1. The detected oscillations in simulation S1 are related to some unsteadiness observed in the outlet pipe. However, for the last 5 s, the relative errors in the mean water volume fraction remain smaller than 0.8% for all simulations. The inclusion of the turbulence model preserved more 0.8% of water in the domain when compared with S1 (the difference between S1 and S2 is 0.8%). The small collapse of water volume, observed at 10 s in simulation S3, is due to the activation of the air-entrainment model. This suddenly imposed a higher quantity of air inside the domain, which recovered after 3 s (13 s into the run time).
Figure 4

Steady-state convergence check for numerical simulations S1 (RANS + VOF), S2 (RANS + VOF + turbulence) and S3 (RANS + VOF + turbulence + air).

Figure 4

Steady-state convergence check for numerical simulations S1 (RANS + VOF), S2 (RANS + VOF + turbulence) and S3 (RANS + VOF + turbulence + air).

Free-surface position

Figure 5 shows the comparison between the numerical and experimental free-surface position. The experimental free-surface is detected using a Computational Vision Model (Roque 2011). The mean water level is marked with dots in Figure 5, whereas error bars are used to bound the maximum and minimum values of the free-surface elevation. The numerical free surface is given by isolines of α. Three isolines of = 0.1, 0.5 and 0.9 are plotted, representing respectively the free-surface levels in which the fraction of water is 10, 50 and 90%.
Figure 5

Free-surface position detected by α = 0.1, 0.5 and 0.9 at y = 0 m for simulations: (a) S1 (RANS + VOF), (b) S2 (RANS + VOF + turbulence) and (c) S3 (RANS + VOF + turbulence + air). Dots mark the experimental free-surface position detected by the Computational Vision Model.

Figure 5

Free-surface position detected by α = 0.1, 0.5 and 0.9 at y = 0 m for simulations: (a) S1 (RANS + VOF), (b) S2 (RANS + VOF + turbulence) and (c) S3 (RANS + VOF + turbulence + air). Dots mark the experimental free-surface position detected by the Computational Vision Model.

A good agreement is found between the numerical and the experimental free surface in a zone where the free surface is smooth, i.e. x < 0.3 m. The position of the hydraulic jump is also well detected in all simulations.

Velocity and turbulence

Figure 6 presents the time-averaged velocity profiles at centre-channel measured with ADV (Figure 6(a)) and simulated by S1 (Figure 6(b)), S2 (Figure 6(c)) and S3 (Figure 6(d)). In the experimental profile (Figure 6(a)), the points measured near the water surface were discarded due to low COR found in the ADV signal. Pearson's-R2 COR coefficients are calculated between numerical and experimental, horizontal () and vertical () components of the velocity.
Figure 6

Time-averaged velocity fields at the centre plane of the gully: (a) measured with ADV, and simulations (b) S1 (RANS + VOF), (c) S2 (RANS + VOF + turbulence) and (d) S3 (RANS + VOF + turbulence + air). The blue area on top of the ADV measurements indicate a zone where the signal COR was low and no measurements were made. The full colour version of this figure is available in the online version of this paper, at http://dx.doi.org/10.2166/wst.2017.071.

Figure 6

Time-averaged velocity fields at the centre plane of the gully: (a) measured with ADV, and simulations (b) S1 (RANS + VOF), (c) S2 (RANS + VOF + turbulence) and (d) S3 (RANS + VOF + turbulence + air). The blue area on top of the ADV measurements indicate a zone where the signal COR was low and no measurements were made. The full colour version of this figure is available in the online version of this paper, at http://dx.doi.org/10.2166/wst.2017.071.

A small vortex was measured by the ADV next to the left side of the pipe outlet, but not simulated by the numerical model. Two reasons may justify this behaviour: (1) since we are not solving DNS, some small turbulent structures are likely not being solved; (2) the ADV is an intrusive method, which in this case occupies a large portion of the gully box, as such it must influence the flow field and be influenced by near wall effects. In any case, the large vortex found inside the gully is reproduced by both the numerical model and ADV measurements.

Simulations S2 and S3 present some slight improvements regarding the accuracy of the velocity patterns. A higher difference is verified in terms of the velocity along the z-axis between S1 and S2/S3, as is , whereas and are and , respectively. The difference is mostly found on the left side of the gully box, where S2 and S3 detect a vertical velocity of ∼0.25 m/s, which is more in accordance with the experimental velocity. The horizontal velocity decreases in its accuracy from S2 to S3. However, if we observe the flow pattern, it looks quite similar. The Pearson's-R2 coefficients between S2 and S3 agree with the previous assumption as and .

The turbulent kinetic energy (TKE) yields the value of kinetic energy generated by the fluctuation velocity field (), 
formula
Figure 7 shows the distribution of TKE at the centre channel measured with ADV (Figure 7(a)) and simulated using S3 (Figure 7(b)), which is very similar to S2. In the experimental profile (Figure 7(a)), the points measured near the water surface were discarded due to low COR found in the ADV signal. Figure 7(c) presents the absolute deviation between values of experimental TKE and numerical TKE (from S3).
Figure 7

Time-averaged TKE profile at centre-plane of the gully: (a) measured with ADV, (b) simulation S3 (similar results to S2) (RANS + VOF + turbulence + air) and (c) absolute deviation.

Figure 7

Time-averaged TKE profile at centre-plane of the gully: (a) measured with ADV, (b) simulation S3 (similar results to S2) (RANS + VOF + turbulence + air) and (c) absolute deviation.

The overall magnitudes of TKE are in accordance with experimental results despite some deviations in the top-left corner and the right side of the gully (Figure 7(c)). On the top-left corner, below the stream that feeds the gully box, the numerical model returns higher values of TKE than in the experiments. However, as TKE is mostly generated due to fluid shear, friction or buoyancy, the numerical results are making more sense than the experimental ones in this specific zone. Moreover, due to the presence of bubbles, low COR signals from the ADV were found in this zone, which could lead to misinterpreted values of experimental TKE. On the right side of the gully, the experimental profile has values of TKE in the order of 0.14 m2/s2. These high values are comprehensive as the flow collides with the top-right corner of the gully and turns down, contributing substantially to the generation of a large eddy inside the structure.

Time-averaging of air-concentration profiles

Figure 8 presents the time-averaged air-concentration profiles in the time interval 45–50 s using simulation S3. Profiles are taken at y = −0.1 (front panel), −0.04, 0, 0.04 and 0.1 m (back panel). Figure 8(f) presents the time-averaged air-concentration profile, integrated over the y-axis for every 5 s of simulation. For this integration we used a discrete air-profile every 0.02 m in the y direction.
Figure 8

Time-averaged air-concentration profiles (45–50 s) using simulation S3, at positions (a) y = 0.1 (back panel), (b) y = 0.04, (c) y = 0, (d) y = −0.04 and (e) y = −0.1 m (front panel). (f) Time-averaged air-concentration profile, integrated over y-axis using simulation S3.

Figure 8

Time-averaged air-concentration profiles (45–50 s) using simulation S3, at positions (a) y = 0.1 (back panel), (b) y = 0.04, (c) y = 0, (d) y = −0.04 and (e) y = −0.1 m (front panel). (f) Time-averaged air-concentration profile, integrated over y-axis using simulation S3.

The air-profiles at y = −0.1 and y = 0.1 m present quite similar features. These profiles are characterised by a large ring of air that occupies almost all the gully and a peak of air-concentration in the centre (maximum of 0.1%). In the positions y = −0.04 and y = 0.04 m, the air-profiles do not show such similarity. Although the profiles are characterised by a central region occupied by air, their shape is different. Nevertheless, the asymmetry can be the result of the calculation of the turbulent statistics of the flow. In the middle profile (y = 0 m), the air is concentrated at the point (0.62, 0.55), reaching a peak of 5%. This was also the maximum value of air concentration found in the gully. Although no experimental data is present to compare these values, the given order of magnitude is plausible since the ADV instrument would not give good correlations (COR > 70%) when the air-concentration is above 10%.

Regarding Figure 8(f), the first aspect that should be highlighted is that this integrated profile is much more understandable in terms of dispersion of air concentration when compared to the remaining profiles. For this flow rate, the hydraulic jump occurring on top of the gully is the only source of air, a result that is in line with previous assumptions (Leandro et al. 2014a). This air is further carried down to the gully following the downstream wall. While some bubbles are trapped in the vortex, some bubbles are allowed to escape through the bottom orifice or the free surface. The latter is likely to happen when the outer part of the vortex reaches the interface. The swirling regime of the vortex is characterised by zero velocity on the centre (Figure 6) and high tangential velocities in the outer part. This is the reason why a large number of bubbles are trapped in its centre, whereas a small percentage flow across the entire gully.

Discharge coefficients

Weir and orifice coefficients can be obtained experimentally and numerically as in (Martins et al. 2014), following the equations: 
formula
17
 
formula
18
where the subscripts w and o stand for weir and orifice respectively, Q is the discharge flow, A the area, b the width, g the acceleration due to gravity and h the uniform water depth upstream of the weir or orifice. Figure 9 presents through lines, the drainage coefficients and according to the power law: 
formula
19
where and are coefficients given in Martins et al. 2014. and are the drainage coefficients for the weir and orifice formulated based on the experimental (E) and numerical (N) data. In the experiments and simulations of Martins et al. (2014), the bottom outlet was under atmospheric pressure as the gully had free outflow. In the present work, the pipe outlet is under surcharge; therefore, the drainage coefficients previously formulated cannot be directly compared to our case. Simulations S1*, S2* and S3* use the same numerical methodology as S1, S2 and S3, except that the bottom outlet pipe BC is now a free outlet boundary. This means that S1* is a copy of the simulations performed by Martins et al. (2014).
Figure 9

(a) Orifice () and (b) weir () discharge coefficients for simulations S1, S2 and S3. Lines are the coefficients found by Martins et al. (2014). Asterisks (*) are simulations where the outlet BC characterises the experiments of Martins et al. (2014).

Figure 9

(a) Orifice () and (b) weir () discharge coefficients for simulations S1, S2 and S3. Lines are the coefficients found by Martins et al. (2014). Asterisks (*) are simulations where the outlet BC characterises the experiments of Martins et al. (2014).

The drainage coefficient decreases for both weir and orifice formulations when the turbulence and air entrainment model was activated. The introduction of an air entrainment formulation still underpredicts when compared to the experimental coefficient . As such, the reason for the underestimated coefficients given by the simulations of Martins et al. (2014) cannot be associated with the non-existence of an air-entrainment model. In the simulations where the gully has a free-outlet (simulations marked with asterisk), the drainage coefficient does not change too much with the inclusion of turbulence or the air-entrainment model. In contrast, for the pressurized outlet (simulations without asterisk), it can be noticed that there is a large difference between S1 and S2/S3. As we have seen in the ‘Velocity and turbulence’ section, the velocity fields for S2 and S3 are more accurate than S1, therefore such a difference is associated with the inclusion of the turbulence model.

Another relevant aspect we can discuss is the influence of air on the drainage. For this particular case, the amount of air found in the gully box almost does not influence the drainage coefficient, and consequently, the drainage efficiency. The inclusion of air should decrease the drainage efficiency, as the volume occupied by the air would at some point interrupt the quantity of water that passes through the outlet (Pothof & Clemens 2010; Ramezani et al. 2016). Nevertheless, for S3 and S3* the inclusion of air has increased compared to S2 and S2*. It could be that, while considering the contribution of air, an increase of momentum may occur as the air within the boundary layer is known to be a reducer of the shear stress (Ackers & Priesley 1985; Chanson 1993). However, more tests on other discharges are needed to verify this finding.

By imposing pressure on the pipe outlet, the simulations return the coefficients without asterisks. As expected, the drainage coefficients were reduced. Such reductions are in the order of 22.4% for S1, 64.4% for S2 and 63.8% for S3 in comparison to the coefficients for the simulations with a free outlet.

CONCLUSIONS

In this work, a 3D gully structure is simulated using a solver that combines the interface tracking VOF model for free-surface detection with a sub-grid air-entrainment model, to simulate the self-aeration process of the flow and the RAS SST k-ω model turbulence model to predict the turbulent statistics.

The following conclusions can be derived from this work:

  • The air-concentration profiles are in qualitative good agreement when compared to the photographs.

  • The numerical accuracy of the velocity patterns were more influenced by the turbulence closure than the air-entrainment model as given by the R2 coefficient of fit.

  • Numerical values of air-concentration do not surpass 6% in the middle profile. Given such a low value, it is not surprising that both the velocity profiles and the drainage coefficients do not change with the simulation of air in the gully.

  • The increment of pressure on the gully's outlet diminishes the drainage coefficient and efficiency. For this specific case, imposing a hydrostatic pressure of 0.02 m in the outlet pipe, we found decrements of about 60% on the drainage.

  • In the case of a gully with a surcharged outlet, it was verified a slight increment of drainage efficiency was gained by activation of the air model. Such a conclusion might be connected to the reduction of shear stress and the increase of momentum inside the gully.

ACKNOWLEDGEMENTS

This study had the support of FCT (Portuguese Foundation for Science and Technology) through the Project UID/MAR/04292/2013 and PhD Grant SFRH/BD/85783/2012 (owned by Pedro Lopes). All the numerical results here showed were performed on the Centaurus Cluster of the Laboratory for Advanced Computing of the University of Coimbra, Portugal.

REFERENCES

REFERENCES
Ackers
P.
Priesley
S. J.
1985
Sef-aerated flow down a chute spillway
. In:
Proceedings of 2nd International Conference on the Hydraulics of Floods and Flood Control, British Hydraulics Research Association, Fluid Engineering
,
Cambridge, UK
, pp.
1
16
.
Blazek
J.
2001
Computational Fluid Dynamics: Principles and Applications
.
Elsevier Science B.V.
,
Switzerland
.
Brackbill
J. U.
Kothe
D. B.
Zemach
C.
1991
A continuum method for modeling surface tension
.
Journal of Computational Physics
100
(
2
),
335
354
.
Carvalho
R. F.
Lemos
C. M.
Ramos
C. M.
2008
Numerical computation of the flow in hydraulic jump stilling basins
.
Journal of Hydraulic Research
46
(
6
),
739
752
.
Carvalho
R. F.
Leandro
J.
Martins
R.
Abreu
J. M.
de Lima
J. M. L. P.
2011
2DV numerical modelling of different flows occurring in gullies
. In:
12th International Conference on Urban Drainage (12ICUD)
,
IWA Publishing
,
London, UK
,
September 11–16, 2011
,
Porto Alegre/RS, Brazil
, pp.
1
8
.
Carvalho
R. F.
Leandro
J.
Martins
R.
Lopes
P.
2012
Numerical study of the flow behaviour in a gully
. In:
4th IAHR International Symposium on Hydraulic Structures
,
APRH Publishing
,
Lisbon, Portugal
,
February 9–11, 2012
,
Porto, Portugal
.
Chanson
H.
1993
Self-aerated flows on chutes and spillways
.
Journal of Hydraulic Engineering
119
(
2
),
220
243
.
Clift
R.
Grace
J. R.
Weber
M. E.
1978
Bubbles, Drops and Particles
.
Academic Press Inc.
,
New York
,
USA
.
Djordjević
S.
Saul
A. J.
Tabor
G. R.
Blanksby
J.
Galambos
I.
Sabtu
N.
Sailor
G.
2013
Experimental and numerical investigation of interactions between above and below ground drainage systems
.
Water Science and Technology
67
(
3
),
535
542
.
Gómez
M.
Russo
B.
2009
Hydraulic efficiency of continuous transverse grates for paved areas
.
Journal of Irrigation and Drainage Engineering
135
(
2
),
225
230
.
Gopala
V. R.
van Wachem
B. G. M.
2008
Volume of fluid methods for immiscible-fluid and free-surface flows
.
Chemical Engineering Journal
141
(
1–3
),
204
221
.
Goring
D. G.
Nikora
V. I.
2002
Despiking acoustic doppler velocimeter data
.
Journal of Hydraulic Engineering
128
(
1
),
117
126
.
Granata
F.
de Marinis
G.
Gargano
R.
2014
Flow-improving elements in circular drop manholes
.
Journal of Hydraulic Research
52
(
3
),
347
355
.
Hirt
C. W.
Nichols
B. D.
1981
Volume of fluid (VOF) method for the dynamics of free boundaries
.
Journal of Computational Physics
39
(
1
),
201
225
.
Lopes
P.
Leandro
J.
Carvalho
R. F.
Páscoa
P.
Martins
R.
2015
Numerical and experimental investigation of a gully under surcharge conditions
.
Urban Water Journal
12
(
6
),
468
476
.
Lopes
P.
Leandro
J.
Carvalho
R. F.
Russo
B.
Gómez
M.
2016a
Assessment of a VOF model ability to reproduce the efficiency of a continuous transverse gully with grate
.
Journal of Irrigation and Drainage Engineering
142
(
10
).
Lopes
P.
Tabor
G.
Carvalho
R. F.
Leandro
J.
2016b
Explicit calculation of natural aeration using a volume-of-fluid model
.
Applied Mathematical Modelling
40
(
17–18
),
7504
7515
.
Ma
J.
Oberai
A. A.
Drew
D. A.
Lahey
R. T.
Hyman
M. C.
2011
A comprehensive sub-grid air entrainment model for RaNS modeling of free-surface bubbly flows
.
The Journal of Computational Multiphase Flows
3
(
1
),
41
56
.
Martins
R.
Leandro
J.
de Carvalho
R. F.
2014
Characterization of the hydraulic performance of a gully under drainage conditions
.
Water Science and Technology
69
(
12
),
2423
2430
.
Menter
F. R.
1993
Zonal two-equation k-ω turbulence model for aerodynamic flows
. In:
AIAA 24th Fluid Dynamics Conference
,
Orlando, Florida
.
Menter
F. R.
Kuntz
M.
Langtry
R.
2003
Ten years of industrial experience with the SST turbulence model
.
Turbulence, Heat and Mass Transfer
4
,
625
632
.
Muzaferija
S.
Peric
M.
Sames
P.
Schelin
T.
1998
A two-fluid Navier-Stokes solver to simulate water entry
. In:
Proceedings Twenty-Second Symposium on Naval Hydrodynamics
,
Washington, DC
.
Nichols
B. D.
Hirt
C. W.
1975
Methods for calculating multidimensional, transient free surface flows past bodies
. In:
1st Intern. Conf. on Numerical Ship Hydrodynamics
,
Gaithersburg, MD
,
October 20, 1975
.
Noh
W. F.
Woodward
P.
1976
SLIC (simple line interface calculation)
. In:
Proceedings of the Fifth International Conference on Numerical Methods in Fluid Dynamics
(
Vooren
A. I.
Zandbergen
P. J.
, eds),
June 28 – July 2, Lecture Notes in Physics
,
Springer
,
Berlin Heidelberg
,
Berlin, Heidelberg
.
OpenFOAM Foundation Ltd
2016
OpenFOAM User Guide v.4.0
.
Páscoa
P.
Leandro
J.
Carvalho
R. F.
2013
Characterization of the flow in a gully: average velocity and air entrainment
. In:
International Workshop on Hydraulic Design of Low-Head Structures (IWLHS 2013), Bundesanstalt für Wasserbau (BAW)
,
Karlsruhe, Germany
,
February 20–22, 2013
,
Aachen, Germany
, pp.
141
148
.
Pothof
I.
Clemens
F.
2010
On elongated air pockets in downward sloping pipes
.
Journal of Hydraulic Research
48
(
4
),
499
503
.
Ramezani
L.
Karney
B.
Malekpour
A.
2016
Encouraging effective air management in water pipelines : a critical review
.
Journal of Water Resources Planning and Management
142
(
12
),
1
11
.
Rider
W. J.
Kothe
D. B.
1997
Reconstructing volume tracking
.
Journal of Computational Physics
141
(
2
),
112
152
.
Roenby
J.
Bredmose
H.
Jasak
H.
2016
A Computational Method for Sharp Interface Advection
.
Royal Society
,
London, UK
.
Romagnoli
M.
Carvalho
R. F.
Leandro
J.
2013
Turbulence characterization in a gully with reverse flow
.
Journal of Hydraulic Engineering-ASCE
139
(
7
),
736
744
.
Roque
J. M.
2011
Medição de alturas de água usando visão computacional num modelo Simulink®
.
MSc Thesis
,
University of Coimbra
,
Portugal
[in Portuguese]
.
Rusche
H.
2002
Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions
.
PhD Thesis
,
Imperial College of Science, Technology & Medicine Department of Mechanical Engineering
.
Russo
B.
Gómez
M.
2011
Methodology to estimate hydraulic efficiency of drain inlets
.
Proceedings of the ICE – Water Management
164
(
2
),
81
90
.
Russo
B.
Gómez
M.
Tellez
J.
2013
Methodology to estimate the hydraulic efficiency of nontested continuous transverse grates
.
Journal of Irrigation and Drainage Engineering, American Society of Civil Engineers
139
(
10
),
864
871
.
Sabtu
N.
Saul
A. J.
Sailor
G.
2016
Hydraulic interaction of a gully system
.
American Scientific Research Journal for Engineering, Technology, and Sciences
21
(
1
),
202
209
.
Ubbink
O.
Issa
R. I.
1999
A method for capturing sharp fluid interfaces on arbitrary meshes
.
Journal of Computational Physics
153
(
1
),
26
50
.
Versteeg
H. K.
Malalasekera
W.
2007
An Introduction to Computational Fluid Dynamics – The Finite Volume Method
.
Pearson Education
,
Harlow, UK
.
Waclawczyk
T.
Koronowicz
T.
2008
Comparison of CICSAM and HRIC high-resolution schemes for interface capturing
.
Journal of Theoretical and Applied Mechanics
46
(
2
),
325
345
.
Wahl
T. L.
2000
Analyzing ADV data using WinADV
. In:
Joint Conference on Water Tesources Engineering and Water Resources Planning & Management
,
Minneapolis, Minnesota
, pp.
1
10
.
Weller
H.
2008
A New Approach to VOF-based Interface Capturing Methods for Incompressible and Compressible Flows
.
Report TR/HGW/04
,
OpenCFD Ltd
.
Wood
I. R.
1991
Air Entrainment in Free-Surface Flows. IAHR Hydraulic Structures Design Manual Nr.4, Hydraulic Designs Considerations
.
A. A. Balkema
,
Rotterdam
,
The Netherlands
.
Youngs
D. L.
1984
An Interface Tracking Method for a 3D Eulerian Hydrodynamics Code
.
Zhang
T.-T.
Yang
W.-J.
Lin
Y.-F.
Cao
Y.
Wang
M.
Wang
Q.
Wei
Y.-X.
2016
Numerical study on flow rate limitation of open capillary channel flow through a wedge
.
Advances in Mechanical Engineering
8
(
4
),
1
11
.
Zhao
W.
Wan
D.
2015
Numerical study of interactions between phase II of OC4 wind turbine and its semi-submersible floating support system
.
Journal of Ocean and Wind Energy
2
(
1
),
45
53
.