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

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

*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.

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

*p*the pressure, the surface forces and

*t*the time. The decomposition of the viscous stress term is given by constitutive relation: 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,

### Free-surface position

*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.

^{®}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)), 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, 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: 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, 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

*et al.*(1978) model, 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), 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, 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).

### Mesh, initial and boundary conditions

*V*= 64 cm

_{cell}^{3}) whereas the smallest cell volume is generated inside the gully box (

*V*= 1 cm

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

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).

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

*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).

### Free-surface position

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

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 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 m^{2}/s^{2}. 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

*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.

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

*et al.*2014), following the equations: 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: 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).

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 R

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