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.
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.
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.
General flow equations
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
Mesh, initial and boundary conditions
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.
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).
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)|
|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)|
|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
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 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
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.
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.
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.
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.