2D Modelling and energy analysis of entrapped air-pocket propagation and spring-like geysering in the drainage pipeline system

Abstract

Transient (highly unsteady) air–water two-phase flows and spring-like geysers have been one of the critical concerns in drainage pipeline systems, which may cause or exacerbate drainage flooding problems and associated damage consequences. In this paper, the flow dynamics and energy evolution mechanism of the induced spring-like geysers are innovatively investigated through a two-phase full-2D numerical model developed in this study. After full validation by laboratory experimental tests conducted in this study, the proposed 2D model is systematically applied to simulate transient air–water flows in drainage pipelines. The results have shown acceptable accuracy of this full-2D model to capture the complex flow interactions between the air and water phases, and indicated that the velocity and pressure distribution patterns are highly relevant to the air–water interface deformation and energy exchange. The in-depth energy analysis demonstrates that the intermittent eruption of geysers could be attributed to the conservation and release of different energy forms during the transient air–water two-phase flow process. Besides, the numerical applications for the systems with different boundaries and initial conditions indicate that the different ventilation conditions and initially entrapped air volume may significantly affect the velocity distribution of the air phase, thereby playing an essential role to provide effective measures to mitigate unexpected geyser events and pressure oscillations in the system. The results and findings of this paper could provide insights to improve the theory and practice of transient air–water two-phase flows in drainage pipeline systems.

KEYWORDS:

• Energy evolution
• transient air–water flow
• air–water interaction
• spring-like geyser
• full-2D modelling

1. Introduction

With the influences of global climate change and rapid urbanization, extreme rainstorm events are increasing in intensity and frequency worldwide (Feng et al., 2022), causing significant economic losses as well as social and environmental problems (Sun et al., 2014). Recently, flooding risks and drainage disasters have become more severe in the urban stormwater drainage system (USDS) (Li et al., 2015). Actually, previous studies and field observations have evidenced that such urban flooding accidents could be caused and exacerbated by the system drainage capacity reduction due to the occurrence of transient air–water flows within the systems under sudden and fast-changing drainage conditions (Bousso et al., 2013). For instance, Hamam and McCorquodale (1982) reported that several cities in Canada suffered surcharging problems with manhole covers popping off by air–water eruptions during severe storms. Similar phenomena could be found in the previous literatures such as Wang et al. (2003), Vasconcelos and Wright (2005), Lewis (2011), and Qian et al. (2020). Meanwhile, the consecutive spring-like eruptions (defined as ‘geyser’ by Chan et al., 2018) of the air–water mixture from manholes in drainage systems have received substantial attention in recent years (Cong et al. 2017; Liu et al. 2020), which may pose destructive damage on USDS (Leon, 2019). From this perspective, the study of transient air–water flows is important for the design and operation of USDS, which will ultimately lead to the increased resilience of these critical infrastructure networks to natural disasters. Noted that some of the hydraulic problems in drainage pipelines can be solved by theoretical methods (Chen et al., 2021; Li et al., 2020), but the present study only focuses on the simulation methods to solve more complex problems about the entrapped air-pocket in drainage pipelines.

Many researchers have performed considerable efforts in this field to observe and understand the physics of the transient air–water flows in drainage pipelines with different system configurations (e.g. Cong et al., 2017; Guo & Song, 1990; Lewis, 2011; Li & McCorquodale, 1999; Shao, 2013; Wright et al., 2011; Zhou et al., 2004, 2013). Specifically, previous experimental examinations can well reproduce the field-observed geyser phenomenon with the main features of consecutive eruptions and similar pressure oscillation intensity (Leon, 2019; Leon et al., 2019), and the strength of the geyser was found to be highly relevant to the diameter ratio of vertical riser and horizontal pipe, the initial condition such as the upstream pressure head and air-pocket size in specific experimental configuration (Lewis, 2011; Vasconcelos & Wright, 2011). However, the sewer geyser duration in realistic USDS may last for the bulk of the time during the air-phase rising period (e.g. a couple of minutes with several wave propagation cycles), which is inconsistent with many previous laboratory observations (usually within several seconds to reproduce the peak flow/pressure only). Wright (2021) deduced that different sewer configurations could incorporate diverse features and no one set of circumstances will control the formation of all sewer geysers. Consequently, a more universal method/ theory may be valuable to understand the intriguing geyser phenomenon.

Accordingly, different types of numerical schemes have been developed in previous studies for the simulation and analysis of transient air–water interactions and complex geyser processes in drainage pipeline systems. In the early stage, the existing single-phase models, such as the rigid column model (Hamam & McCorquodale, 1982), and the two-phase coupling TPA (Two-component Pressure Analysis) model developed by Vasconcelos and Wright (2007), were widely applied for simple and low-speed mixed-flow simulation with strengths of less time and memory consuming, but the accuracy of the simulation results was inevitably less satisfactory due to the simplification or ignorance of the internal air phase effect. As a result, the full two-phase flow models tend to be more promising to reproduce the complicated air–water interactions. Recently, the 3D computational fluid dynamics (CFD) based numerical models have been applied to better reproduce the unsteady air–water interactions (including capturing the interface reconstruction dynamics during air-pocket formation and release process) with more accurate simulation of pressure transients (Chan et al., 2018; Chegini & Leon, 2020; Choi, 2018; Leon, 2019; Leon et al., 2019). For the flow pattern investigation, He et al. (2022) discussed the velocity distribution along the radial direction of pipe cross-section at the midpoint of an isolated slug flow driven by pressurized air, while Wang et al. (2023) revealed that the given Taylor bubble shape is closely related to the relative bubble velocity. These observed mechanisms of specific air–water interactions in horizontal and vertical directions are referable for further sewer geyser investigation, and detailed velocity distribution and evolution of the air phase should also be taken into consideration. However, it is noted that these 3D models are too time-consuming for practical cases, which may be more favourable for meticulously capturing the main characteristics of the violent erupted geyser for specific research purposes. Actually, Molina and Ortiz (2020) simulated several cases on the propagation of air cavities and concluded that the axisymmetric model for circular ducts preserves the main features of the 3D model, so that the application of 2D numerical models solving the two-phase equations can be the cost-effective choice to reproduce the transient air–water flows with acceptable accuracy and efficiency. For instance, Shao (2013) developed an incompressible 2D model coupled with the volume of fluid (VOF) method (introduced by Hirt and Nichols (1981)) to simulate the air pocket rising process; subsequently, Shao and Yost (2018) modified this scheme for further inspection on the influential factors of spring-like geyser surge. Both proposed models can well capture the distinct air–water interaction profile of the interface before the explicit eruption of the geyser, but the movement of an entrapped air-pocket propagating in the horizontal sewer and the reconstruction of the bubble in the vertical riser was not fully simulated. With these achievements and findings in the literature, and considering the scope of this research within the relatively moderate air–water interaction such as entrapped air-pocket propagation and spring-like geyser, a full-2D model would be adequate to capture more dynamics of air phase with efficient numerical simulation.

Although significant achievement has been accomplished based on experimental and numerical investigations, research on the geyser phenomenon is still in the primary stage for ascertaining dominating factors during the geyser formation process and its profound influence on realistic drainage structures. Actually, the further relationship of energy transformation has not been taken into consideration for the complicated geyser system. As suggested by Karney (1990), the energy approach could support an integrated insight for the transient dynamic propagation in closed-conduit flow, which is physically consistent with the realistic pipeline system and could well interpret the sophisticated transition process in a simple, efficient way. Some previous literature research was carried out for investigating the energy constitution and transformation pattern in the field of transient pipe flows (e.g. Duan et al., 2010, 2017; Jelev, 1989), and recently Rokhzadi and Fuamba (2022) discussed a particular variable of the maximum pressure of the entrapped air-pocket on rapid filling issues in a partially pressurized system. Still, the detailed analysis of the influence of air phase involvement on the energy dissipation mode during transient air–water transition has not been sufficiently investigated to fully understand the overall geyser revolution of the process, which is explicable as the incapability of calculation for previous research with insufficient collected data.

This paper aims to investigate transient air–water interactions to derive applicable evolution mechanisms to improve the understanding of the realistic sewer geyser by examining the detailed pressure and velocity distributions so as to discuss the relationship between the flow pattern and the interface deformation, and perform energy evolution analysis for both air and water phase with the aid of the applied full-2D numerical model. To this end, governing equations to simulate the hydraulic characteristics as well as the coupled air–water interface tracing scheme are first introduced in Sections 2.1 and 2.2, followed by the method for energy revolution assessment in section 2.3. Based on the validated 2D models, a more visible and quantitative inspection of the essential stage during the dynamic process both in the horizontal pipeline and vertical riser could be studied in detail in section 3. Besides, the different forms of energy (i.e. the potential and kinetic energy of water and air phase), as well as flow dynamics such as the pressure oscillation pattern would be investigated to figure out the underlying mechanisms of the transient air–water mixed flows. Finally, the practical implications of the results and the main achievement of this study are summarized in the conclusion part. The numerical analysis and findings achieved in this study are expected to improve and enhance the theory and practice of the transient air–water two-phase flows in drainage pipe systems.

2. Methodology and models

To capture the complex flow dynamics and evolution processes of transient air–water interactions, a Coupled Level Set and Volume of Fluid (denoted as CLSVOF) model (Sussman & Puckett, 2000; Son & Hur, 2002) and the corresponding numerical scheme are proposed for the numerical simulations by solving the full 2D Navier-Stokes equation. With this full-2D model and CLSVOF scheme, it is expected that the variation of surface curvature and surface tension during the transient air–water flow process can be captured accurately by the Level Set (LS) scheme (Kan et al., 2023; Osher & Sethian, 1988; Sussman et al., 1994), and meanwhile, the mass loss occurred in LS function can be corrected precisely by incorporating the VOF approach (Chang & Wei, 2022; Zhu et al., 2022). The complete model and numerical scheme are provided in detail as follows.

2.1 Governing equations

The governing mass and momentum equations of fluid motions (air and water) in transient air–water flows, by considering incompressible air content (Shao, 2013), are as follows: (1) $\begin{array}{rl}& \mathrm{\nabla }\cdot U=0\end{array}$(1) (2) $\begin{array}{rl}& U_t+U\cdot \mathrm{\nabla }U=-\frac{\mathrm{\nabla }P}{\rho \left(\phi \right)}+\frac{\mathrm{\nabla }\cdot \left(\mu \left(\phi \right)\left(\left(\mathrm{\nabla }U\right)+{\left(\mathrm{\nabla }U\right)}^T\right)\right)}{\rho \left(\phi \right)}\\ & \phantom{U_t+U\cdot \mathrm{\nabla }U=}-\frac{\gamma \kappa \left(\phi \right)\mathrm{\nabla }H\left(\phi \right)}{\rho \left(\phi \right)}+F\end{array}$(2) (3) $\begin{array}{rl}& {\phi }_t+U\cdot \mathrm{\nabla }\phi =0\end{array}$(3) (4) $\begin{array}{rl}& f_t+U\cdot \mathrm{\nabla }f=0\end{array}$(4) where U = (u, v) is the velocity field (u is the horizontal velocity, and v is the vertical velocity), P is the pressure, $\gamma$ is the surface tension coefficient, $\kappa \left(\phi \right)$ is the local mean curvature, F is the body force due to gravity F = (0, g), $\phi$ is the LS function defined as a signed distance from the interface (i.e. the zero-level-set), $f$ is the volume fraction applied in VOF method. In this case, the LS function $\phi$ is positive in the water and negative in the air; the density $\rho \left(\phi \right)$, viscosity $\mu \left(\phi \right)$ and curvature $\kappa \left(\phi \right)$ are dependent on the LS function ($\phi$), and $H\left(\phi \right)$ is the Heaviside function as defined below: (5) $H\left(\phi \right)=\{\begin{array}{ll}1& \text{if }\phi >0\\ 0& \text{otherwise}\end{array}$(5)

2.2 Numerical methods

To achieve an efficient solution process, a staggered numerical method as shown in Figure 1 is adopted in this study for solving the 2D NS equation. Specifically, velocity U and some scalar variables such as pressure, volume fraction, body force, density, and viscosity are defined at the centre of a computational cell; meanwhile, the edge-based velocity field u_{i ± 1/2,j} and v_{i,j ± 1/2} is located at the centre of cell edges, and the pressure difference term $q=p^{n+1}-p^n$ is defined at the cell corners.

Figure 1. Computational grid and location of different discrete variables.

The main algorithms and numerical schemes for solving the proposed 2D model are illustrated in Figure 2, and the application procedures are as follows:

1. discretizing the momentum equation for an explicit prediction of the edge-based velocity field $U^{ADV}$ and advection term $(\left(U\cdot \mathrm{\nabla }\right)U{)}^{n+1/2}$ at half time step t^n + 1/2 in the 2D model;

2. with the updated $U^{ADV}$, coupling the LS and VOF functions (${\phi }^{n+1}$ and $f^{n+1}$) which is then corrected synchronously by the flux split advection schemes, so that the density and viscosity at half time step can be updated with ${\phi }^{n+1}$;

3. based on the above half-time-step updated results, solving the NS momentum equation (except for the pressure term) to obtain the intermediate velocity field $U^{\ast }$;

4. updating the pressure field $p^{n+1}$ by implicitly solving the Poisson equation with the multigrid projection method, followed by the correction and updating of the divergence-free velocity field $U^{n+1}$ by the obtained pressure gradient.

Figure 2. The algorithm for the overall computation procedure.

In this study, the proposed full-2D model and the CLSVOF numerical scheme are implemented by the in-house codes in the MATLAB platform, which are to be tested and validated firstly by the experimental tests and then applied to further examine the underlying mechanism and dynamics of transient air–water flows later in this paper.

2.3 Energy constitution and evolution analysis

Based on the available data of velocity and pressure distributions from the proposed full-2D model, different forms of energy could be obtained for drainage pipeline systems in different variation stages. However, compared to the kinetical and potential energy forms which have pervasive applications in engineering projects, other forms of energy (such as internal energy) may not be so convenient to quantify for engineers in system design. To this end, only the variation patterns of kinetical and gravitational potential energy between air and water phases are retrieved from the total energy variations and analysed in detail in the following study. For simplicity, the kinetic and gravitational potential energy forms of the water phase are denoted as ‘E_kw’ and ‘E_pw’, and those of the air component are written as ‘E_ka’ and ‘E_pa’ respectively. The mathematical forms as given as, (6) $\begin{array}{rl}& E_{pw}=\sum _{i,j}{\rho }_w\cdot g\cdot \mathrm{\Delta }x\cdot \mathrm{\Delta }y\cdot y_{i,j}\cdot f_{i,j}\end{array}$(6) (7) $\begin{array}{rl}& E_{kw}=\sum _{i,j}\frac{1}{2}{\rho }_w\cdot \mathrm{\Delta }x\cdot \mathrm{\Delta }y\cdot U_{i,j}^2\cdot f_{i,j}\end{array}$(7) (8) $\begin{array}{rl}& E_{pa}=\sum _{i,j}{\rho }_a\cdot g\cdot \mathrm{\Delta }x\cdot \mathrm{\Delta }y\cdot y_{i,j}\cdot \left(1-f_{i,j}\right)\end{array}$(8) (9) $\begin{array}{rl}& E_{ka}=\sum _{i,j}\frac{1}{2}{\rho }_a\cdot \mathrm{\Delta }x\cdot \mathrm{\Delta }y\cdot U_{i,j}^2\cdot \left(1-f_{i,j}\right)\end{array}$(9)

Note that the above energy expressions are approximately averaged for each computational cell, which ignores the accurate air–water interface profile for estimating the actual mass of either water or air phase within each cell.

3. Experimental validations and numerical applications

Prior to the applications of the proposed full-2D model, the presented numerical model and scheme are first verified and evaluated through different experimental tests from previous studies and the well-designed laboratory system in this study. Specifically, the overall schematic setup of the experimental system used in this study is shown in Figure 3. This system is composed of an upstream supply tank, horizontal and vertical pipelines, a downstream drainage tank, and control devices (valves and recycling facilities). The upstream tank is 0.6 m long, 0.6 m wide, and 0.8 m height, with its bottom connected with the horizontal pipe with an outer diameter = 0.1 m and inner diameter = 0.094 m by a hosepipe with the same diameter. The vertical riser with also inner diameter = 0.094 is installed at a distance of 2.255 m from the upstream tank. The butterfly valve ‘B’ is located at a distance of 3.018 m away from the centre of the riser to hold the pressurized water, while the distance of each valve ‘A’, ‘B’, ‘C’, and ‘D’ is 1.006 m.

Figure 3. Schematic of laboratory experimental test system for this study.

During the tests, the sections of air and water are initially isolated by the butterfly valves, and an air–water interface could be formed once after the fast opening of the valve. For example, the uses of valve A–D can provide different sizes of entrapped air-pocket to study the air–water interaction in rapid filling process, which will be further investigated in section 3.1. For more complicated geyser cases analysed in section 3.2, if the propagated air-pocket arrived at the vertical riser, a complex process with drastic air–water interaction and air–water mixture release through the vertical riser will be induced, which is recorded by the high-speed camera and measured by the high-frequency pressure transducers (1kHz) installed at the bottom of the vertical riser in the system.

3.1 Air-pocket propagation in the horizontal pipe

To investigate the air-pocket advancing effect during rapid drainage process, the experimental test settings of the air cavity in Figure 4, which are similar to the previous experimental study in Atrabi et al. (2015) and Cong (2016), are applied here for validating the proposed method and model. The detailed parameters including the initial water length L_w, air length L_a, the height of initial water column h, and applied grid sizes for test cases are listed in Table 1.

Figure 4. Experimental system settings for air-pocket advancing process.

Table 1. Settings of test cases for rapid filling process in horizontal pipes.

Test case	Top boundary (cover wall)	Lw (m)	La (m)	h (m)	Grid size (length × height) of numerical model
Case 1	Closed wall	1.000	1.000	0.100	0.001 m × 0.001 m;0.0025 m × 0025 m;0.005 m × 0.005 m
Case 2	Closed wall	1.250	0.450	0.050	0.005 m × 0.001 m
Case 3	Closed wall	2.000	1.006	0.100	0.005 m × 0.005 m
Case 4	Open outlet	1.000	1.000	0.100	0.005 m × 0.005 m

For the implementation of the proposed model, the computational domain was discretized using a structured quad grid element, and the grid density was enhanced particularly in regions of high flow gradients or of smaller size (such as regions near the wall and within the vertical riser for geyser simulation), so as to capture the features of the flow field and the complex interface deformation. Before carrying out all the simulation cases, a grid independence verification was fully performed through coarse, medium, and fine grid sizes. As shown in Table 1, the grid density with acceptable accuracy and computational time was applied in all simulation configurations.

The simulated results of the air–water interface profile with coarse grids in case 1 are plotted in Figure 5(a)–(d) together with the experimental data for comparison, which generally demonstrate the consistent trends of numerical results with the experimental data obtained from Atrabi et al. (2015). In particular, although the flow levels fluctuate with time during the overall air pocket propagation process, the distinct interface profiles simulated by different grid sizes were found highly overlapped with each other especially before the interface tail arrived at the end, and the boxplots of the water level difference simulated by different grid sizes are illustrated in Figure 5(e). As the propagation velocity of the air pocket front is the more commonly used parameter for model validation, the difference in this parameter is also plotted in Figure 5(f) for grid-size comparison. Both these results show that the differences between coarse and fine grids are less than 10%, and most of the average difference is within 5%. Therefore, the coarse grid can be adopted for the investigation of test cases in this study.

Figure 5. Results of air-water interface profiles from experimental measurement and numerical modelling in case 1 with different grid sizes: (a) interface profile at t = 0.6 s; (b) interface profile at t = 1.2 s; (c) interface profile at t = 1.6 s; (d) interface profile at t = 2.2 s; (e) interface profile difference with different grids; (f) difference of propagation velocity with different grids.

Moreover, the temporal variations of air-pocket front advancing obtained from the model are plotted in Figure 6 for comparison with experimental data of cases 1, 2, and 3 (Cong, 2016). The results comparison shows clearly a very good accuracy of the model prediction with the measured upper front positions of entrapped air-pocket. Consequently, these comparisons confirm again the suitability and reliability of the developed 2D model and numerical method for the study of the transient air–water interaction process in drainage pipeline systems. The detailed hydraulic characteristics with different boundary conditions of case 1 and case 4 are discussed later in section 4.1.

Figure 6. Temporal variation of the air-pocket front under different initial water length conditions in the closed system.

3.2 Dynamics of air-pocket in the vertical riser filled with water

To further study more complex air–water interactions with the vertical riser as shown in Figure 3, different parameters and settings of experimental tests for the geyser case are tested and listed in Table 2, and the detailed dynamics of air–water interaction in the vertical riser in case 5 is simulated and validated with the measured data in this section.

Table 2. Parameter settings for different geyser test cases.

Test case	The boundary condition at the top of the riser	Height of upstream water Hw (m)	Length of air pocket La (m)
Case 5	Open (pressure outlet)	0.98	1.006
Case 6	Closed (wall)	0.98	1.006
Case 7	Open (pressure outlet)	0.98	2.012
Case 8	Open (pressure outlet)	0.50	1.006

The no-slip condition is imposed for the wall boundary condition, and the constant pressure head condition is applied for the inlet of the upstream tank and the open end of the vertical riser.

The snapshots of experimental tests during different evolution stages of the air-pocket rising process along with the numerical simulation results are compared in Figure 7. Once the air-pocket arrives at the riser, the pressure at the centre bottom of the riser decreases accordingly due to the involvement of the relatively low-pressure air phase. As shown in Figure 7(a), the bottom of the riser is filled with air–water mixture during the initial air bubble rising process, with an upper air–water interface formed distinctly (as highlighted by the red circle region in Figure 7(a)). Meanwhile, the top of the free-surface water level in the riser is slightly increased with the elevation promotion of the air-pocket. Afterwards, the powerful buoyancy force instantly drives the air-pocketing rising and leaving away from the bottom of the riser, as shown in the stage of Figure 7(b), which further elevates the level of free-surface top and results in the instantaneous escalation of the bottom pressure. During this process, the air–water interface (circled in the red region in Figure 7(b)) also accelerates to catch up with the free-surface level, which makes the distance even more closed. At the later stage of the release process, the air–water interface is very close to the free-surface level, so that the air-pocket will merge quickly with the outer atmosphere and thus break down into broken expansive bubbles (as shown in Figure 7(c)). After this moment, the dispersed bubble will fall down to the bottom of the riser and form new air-pocket(s) to reproduce the similar process of the air–water interaction and release, which actually contributes to the pressure oscillation in the vertical riser.

Figure 7. Numerical and experimental results of air-pocket rising process in the vertical riser at: (a) t = 8.1 s; (b) t = 9.0 s; (c) t = 11.7 s.

For validating the model, more results of numerical simulation and experimental measurements are plotted in Figure 8 during the first bubble-rising process (7.5∼9.7s) in this test. Overall, the simulated results are consistent with the measurement data, demonstrating the bubble’s rising speed exceeds the ascending velocity of the free-surface water lever especially at the later rising stage, which complies with the similar experimental and numerical tests in Chan et al. (2018). Consequently, it is promising to conclude from the comparison of the results that the developed model can well reproduce the air-pocket rising process, which achieves remarkable performance in tracing the complex deformation of transient air–water evolution, including clear reconstruction of the distinct air–water interface and free-surface level in such complicated dynamic processes.

Figure 8. Comparison of simulated and measured water levels and bubble interface heights at the centreline of the riser.

Furthermore, the pressure variations at the bottom of the vertical riser with different models and experimental results during the transient air–water evolution process are provided in Figure 9 for inspection, and the oscillations of the pressure traces are corresponding to the observed phenomena in Figure 7 regarding the air–water mixture motion and air-pocket release process. The overall results of the proposed full-2D incompressible model (in black solid line) can offer acceptable agreement with experimental tests, in both maximum amplitudes and oscillation phases of the pressure trace. Specifically, the simulated maximum pressure of the first air–water eruption at t = 9.7s was about only 1.79% lower than the measured data (shown in Table 3), which performs much better than the simplified analytical model proposed by Liu et al. (2020), and comparable to the accuracy of their improved model (with a difference of 2.3%) in Liu et al. (2022).

Figure 9. Comparison of pressure oscillations at the bottom of the vertical riser by different models with experimental data.

Table 3. A detailed comparison between measured pressure and simulated data.

		Experimental measurement of pressure P0 (Pa)	2D incompressible model P1 (Pa)	2D compressible model P2 (Pa)	(P1 − P0)/P0 × 100% (%)	(P2 − P0)/P0 × 100% (%)
t < 7s	Average	9584.60	9875.18	9468.19	3.03	−1.21
	Maximum	9767.81	9581.44	10,665.37	1.91	9.19
	Minimum	9259.11	9567.46	7837.25	3.33	−15.36
t = 9.7 s	Average	10,678.33	10,487.50	10,975.91	−1.79	2.79
t = 13.6 s	Average	10,191.01	10,097.40	10,404.07	−0.92	2.09

However, it is also noticeable for the difference in the pressure results during the first 7s for different models. In particular, compared with the measured data, a commercial CFD model of ANSYS FLUENT based on a 2D compressible model under similar settings applied in Chan et al. (2018), has overestimated the pressure oscillations, while the proposed full-2D incompressible model in this study has underestimated the results. To explain the discrepancies in Figure 9, on one hand, the pressure oscillation simulated by the compressible model may be dominated by the water hammer effect (featured with unsustainable peak pressure only occurred at the first 1∼2 periods and then weakened instantly), which is quite different from the measured pressure variation pattern. On the other hand, the influence of the discs of the used butterfly valves has been ignored in the numerical model for simplification, which may affect the air-pocket propagation although valve ‘C’ and valve ‘D’ are fully open during the test process.

After the initial setup and formation of the air-pocket, the proposed full-2D model behaves much better than the CFD-based 2D compressible model in terms of both oscillatory amplitude and phase. Based on these above results and analysis, the developed numerical model can perform superior advancement in efficient air–water simulation for capturing the detailed air–water interaction and evolution, which is very useful to in-depth understand the physical mechanism and process of such complex phenomena in practical engineering systems.

4. Results and discussions

4.1 Investigation of physical dynamics of the rapid filling process

To further investigate the detailed interaction process, the pressure and velocity distributions at an intermediate moment before the wave phase arrives at the downstream end (e.g. t = 0.6 s) for case 1 of the closed system are visualized in Figure 10. In this figure, a brief configuration of grid mesh about the model structure is plotted in Figure 10(a). Figure 10(b) shows the pressure distribution contour, and the air-pocket approaching upstream is driven by the pressure difference at each side of the stagnant point, as the top of the upstream end is filled with a negative pressure plotted with dark blue colour, while the other side of the stagnant point is close to the standard atmospheric pressure. The detailed pressure variation in Figure 10(b) clearly shows that there exists a noticeable pressure gap across the stagnant point, which is the driving force of the upper air phase propagation. For the air phase on the upper-side, this propelling force makes the contour lines mainly distributed in the horizontal direction (shown in Figure 10(c)), which is distinctly different from the contour lines of the water phase distributed in the vertical direction (mainly driven by hydraulic pressure). Besides, the bottom tail of the water wave perceives a faster moving speed compared with the front of the air-pocket advancing, as the tail of water wave is reaching 1.6 m on the bottom while the front of the air-pocket only approaches 0.7 m towards the upstream. Moreover, the velocity distribution is highly correlated to the air–water interface profile, both the motions of air and water phases are regularly changed along the interface profile. As shown in Figure 10(d), for the region near the downstream tail of the air–water interface, the velocity vector tends to be tangent to the interface profile, and the direction of the vector reverses at the upper side of the interface. As the air–water interface is highly related to the sectional area of the pipeline in the drainage system, and the regular-distributed velocity can be further investigated to estimate the flow rate variation along the air–water interface, which is promising to support the quantification of the flow attenuation effect of the entrapped air-pocket, so as to assess the reduction of drainage capacity in realistic USDS.

Figure 10. (a) brief configuration of grid mesh about the model structure; (b) contours of pressure (Pa) for case 1 at t = 0.6 s; (c) contours of velocity magnitude (m/s) for case 1 at t = 0.6 s; (d) the velocity vector field for case 1 at t = 0.6 s.

Afterwards, as shown in Figure 11 for case 1 with the closed end at t = 1.2 s (after the water tail is reflected by the downstream end), the pressure and velocity distributions still comply with the similar patterns of this configuration, and the air–water interface profile tends to be the medium of exchange for the complex variations of different dynamic properties (e.g. pressure and velocity) in the system, as illustrated in Figure 11. However, the drastic deformation of the tail of the air–water interface greatly affects the exchange process, and the high-speed area (in yellow and red colour in Figure 11) decreased quickly in this period.

Figure 11. The contour of dynamic variables for case 1 of the closed system after the water arrives at the downstream end: (a) contours of pressure (Pa) for case 1 at t = 1.2 s; (b) contours of velocity magnitude (m/s) at t = 1.2 s; (c) velocity vector field at t = 1.2 s.

By contrast, the result of pressure and velocity distributions with the open outlet condition of case 4 at this same time moment (t = 0.6 s) is illustrated in Figure 12, which shows the distinct difference from the stagnant front in a closed system in case 1. Actually, the propagating front in this case 4 can fleetly reach the upstream end, which leads to a very quick collapse of the water column. During this process, the air phase appears to be less influenced by this interaction with uniform atmospheric pressure. As the driven force on the bottom tail of water is similar to the condition in case 1, the contour lines distribution pattern and the velocity vector field in Figure 12(b,c) are quite similar with the region in Figure 10(c,d); while the distribution pattern of air phase is totally different. Meanwhile, the significant velocity variation region is confined to the downstream tail of the air–water interface profile, and the region of the high-speed air phase (e.g. velocity > 0.5 m/s) is much less than the situation of the closed system in case 1.

Figure 12. The contour of dynamic variables for case 4 of the open system at t = 0.6 s (before the water arrives at the downstream end): (a) pressure contour; (b) velocity distribution; (c) velocity vector field.

4.2 Energy transformation during the rapid filling process

Based on the validated full-2D model, the dynamic evolutions of different energy forms during the rapid filling process of case 1 and case 4 are systematically investigated in this section. As stated in section 2.3, the total mechanical energy including the kinetical and potential energy of both water and air components over the entire computational domain is denoted as ‘E_t’ (i.e. $E_t=E_{pw}+E_{kw}+E_{pa}+E_{ka}$). The variations of ‘E_t’ for open and closed systems are plotted and compared in Figure 13(a), which reveals that the boundary with an open outlet (case 4) experiences more decrease of energy ‘E_t’, and this descending trend performs even obvious in the later stage of the rapid filling process. To further examine the energy variations, the energy forms of ‘E_kw’ and ‘E_pw’ of the water phase, along with ‘E_ka’ and ‘E_pa’ of the air phase, are retrieved from the model results and discussed as shown in Figure 13(b,c), respectively.

Figure 13. Variations of different energy forms for both the closed (case 1) and open (case 4) systems: (a) loss of total mechanical energy ‘E_t’; (b) kinetical and potential energy variations of water phase; (c) kinetical and potential energy variations of air phase.

Through comparisons of closed and open system results in Figure 13, the energy change in the water phase occurs relatively slowly and gently for both kinetical and potential energies (‘E_kw’ and ‘E_pw’) in the closed system, which is in accordance with the trend of the total energy ‘E_t’. On the contrary, the magnitude of energy constitution in the air phase is rather negligible compared to that of the water phase, with a value below 0.1 J (about only 0.2% of the value in the water phase). Nevertheless, the entrapped air-pocket in the closed system still introduces its typical features: the accumulation of potential energy (‘E_pa’) grows slower and the kinetical energy of air escalates higher than those in open outlet condition, which is in accordance with the findings made by Rokhzadi and Fuamba (2022) that the pressurized flow tends to store more kinetic energy and free-surface flow may contribute less in the energy exchange process. Moreover, the temporal variation of kinetical energy also agrees with the velocity distribution discussed in Section 4.1, which is highly related to the extent of interface deformation. This correlation implies that less kinetical energy can be stored in the later stage as the distorted interfaces are usually accompanied with more energy dissipation.

4.3 Energy evolution analysis for transient air–water geysering process with different riser ventilation conditions

It can be observed from the previous results in Figure 8 that the whole air–water flow evolution process of that experimental case 5 could be divided into four key stages for their different formulation and evolution mechanisms, which are shown in Figure 14 and elaborated as follows.

Figure 14. Evolution of different energy forms with time in case 5: (a) potential energy of water ‘E_pw’; (b) kinetic energy of water ‘E_kw’; (c) potential energy of air ‘E_pa’; (d) kinetic energy of air ‘E_ka’.

During the initial stage (e.g. t < 7 s, the light blue region marked as ‘stage 1’ in each subplot of Figure 14), the air-pocket propagates gradually towards the upstream (close to the riser) under the pressurized condition, which is quite similar with the variation pattern as clarified in Figure 13. During this period, the potential energy of water ‘E_pw’ (see Figure 14(a)) transfers gradually to other energy forms: the kinetic energy of water ‘E_kw’, the potential energy of air ‘E_pa’, and the kinetic energy of air ‘E_ka’, shown in Figure 14(b,c and d), respectively, with the continuously decreasing of energy ‘E_pw’, increasing of energy ‘E_pa’, and the synchronized rise-fall variation mode of kinetic energy for both ‘E_kw’ and ‘E_ka’.

When the pressurized air-pocket arrives at the vertical riser (e.g. t ∼ 7 s), the air bubble is formed promptly in the vertical riser to replace the corresponding volume due to the relatively larger pressure of the arrived air-pocket, which thereafter starts to rise up because of the difference of pressure and density between air and water phases. Since this moment, the transient air–water interaction enters into the air bubble rising process as studied in the former Section 3.2. As a result, this transition of air–water interaction stages may have induced an instant decreasing of potential energy ‘E_pw’ and immediate soaring of kinetical energy ‘E_kw’ until t = 8.7 s in the results of Figure 8 (i.e. the pink region marked as ‘stage 2’ in each subplot of Figure 14).

Afterwards, the decelerated rising bubble is approaching the peak-height location of the vertical riser with drastic decreasing kinetical energy ‘E_kw’ from t = 8.7 s to t = 9.6 s (i.e. the brown region marked as ‘stage 3’ in each subplot of the figure). During this period, since part of the water column in the vertical riser is lifted up by the rising bubble, the potential energy ‘E_pw’ is changed with a slightly escalating trend with a higher position of water in the vertical riser.

Finally, it comes to stage 4 with the marked yellow region in each subplot of Figure 14, where the air bubble becomes fractured at the peak location in the vertical riser due to the pressure and velocity oscillations in the system. During this stage, a spring-like geysering phenomenon of air–water mixing eruption could be formed at the top of the riser, if only the formed motion (velocity and pressure) is large enough and at the same time the riser is open to allow the mixed fluid (air and water) to be released/erupted (otherwise the mixed fluid is further compressed due to covered top end of the riser). Immediately after the air–water mixture eruption/geysering, the splitting air–water mixture (or part of its total volume due to the ‘loss’ at the top outlet of the riser) will fall down with decreasing ‘E_pw’. Meanwhile, the next round of fluctuation of energy transformation (e.g. E_kw) will repeat until all the air is released from the riser (as revealed by the oscillations after stage 4 in Figure 14).

Through the comparison of the energy transformation trend during the whole air–water interaction process in Figure 14, it is clearly shown that the variation of potential energy turns out to keep always in the opposite trends for air and water phases, as shown in Figure 14(a,c). However, for the kinetic energy, a similar (consistent) variation trend is kept for both air and water phases during each evolution process in spite of their different variation amplitudes, as shown in Figure 14(b,d). Moreover, for the specific energy form such as potential energy, the variation period between different phases is also exactly the same due to the synchronous interactions of air and water components, e.g. the same variation frequency between ‘E_pw’ in Figure 14(a) and ‘E_pa’ in Figure 14(c). Compared with potential energy transformation, the oscillation of kinetic energy becomes more frequent for both phases. By estimation, it turns out that the variation period of kinetic energy is approximately half of that of potential energy (or the oscillation frequency of kinetic energy is about as twice as that of potential energy). This is again because of the locally drastic interactions between the air and water phases during the air bubble moving along the pipeline and rising upward in the riser.

To study the influence of the ventilation condition of the top/cover of the vertical riser, the temporal variations of different energy forms with open (case 5) and sealed risers (case 6) are calculated and plotted in Figure 15 for comparison. During the initial stage (‘stage 1’) before the air-pocket reaches the vertical riser (t < 7 s), the variation pattern is almost the same since the formulation and evolution of the rapid filling process are actually identical for both cases. However, the air–water interaction in the vertical riser becomes quite different once the air-pocket arrived at the riser. Specifically, under sealed-riser conditions, the upper air of the vertical riser is continuously compressed by the pressure oscillation during the air bubble rising process, so that the air–water dynamic oscillation can be greatly suppressed in both the pipeline and the riser. Actually, the air-pocket entrapped in the rapid filling process mixes quickly with and passes through the water phase due to the buoyance force during the bubble rising and releasing period (e.g. 7 ∼ 12 s in Figure 15) with drastic distortion of the air–water interface. Meanwhile, the potential and kinetic energy forms of the water phase, and the kinetic energy of the air phase in Figure 15(a,b, and d), respectively all decrease quickly in this period. Afterwards, the air–water mixture becomes relatively stable and no further oscillation of energy is observed. Unsurprisingly, no drastic air–water mixture eruption/geysering was observed in this system due to the fully closed top of the vertical riser. In fact, similar experimental tests with different ventilation sizes of manhole covers were conducted by Zhang et al. (2022), which confirmed that decreasing cover ventilation reduces the number of geysers. Additionally, Zhou et al. (2023) revealed that air–water interaction within a vertical dead end may evolve more energy dissipation (with convective heat transfer and transient wall shear stresses) after the first two pressure oscillations, which justified that less residual potential and kinetic energy can be reserved for geyser eruptions.

Figure 15. Temporal variations of different energy forms with open and sealed risers: (a) potential energy of water ‘E_pw’; (b) kinetic energy of water ‘E_kw’; (c) potential energy of air ‘E_pa’; (d) kinetic energy of air ‘E_ka’.

By comparison, the variations of the dynamic pressure at the centre of the riser bottom with open and sealed risers are shown in Figure 16. With the sealed riser (in blue dotted line), the entrance of the air bubble at the bottom of the riser induces a slight decline of the dynamic pressure, which is thereafter recovered rapidly during the air–water mixture rising stage. Unlike the result with an open riser (in black solid line), no clear pressure surge oscillations were observed during the air–water mixing and rising process. Accordingly, it is found that the boundary control method (e.g. riser top cover) could be a suitable way to reduce pressure failure risks.

Figure 16. Variations of static and dynamic pressures at the centre of the riser bottom for systems with open and sealed risers.

Furthermore, for the dynamic bubble rising process with an open riser (e.g. t > 7 s), the water level converted hydrostatic pressure (in red dashed line in Figure 16) shows a significant difference with the dynamic pressure extracted from the model due to the transient air–water mixing motion. Specifically, the dynamic pressure is less than the equivalent static pressure (converted by water level) due to the existence of an air bubble and its interaction with water in the riser. This amplitude difference becomes particularly profound for the first 1∼2 rounds of the transient air–water mixture rising and geysering process. It is also noted that the oscillation phases of these two parameters almost coincide with each other. This finding is actually consistent with the former results of energy analysis. In fact, the rising bubble in the vertical riser could significantly ascend the upper water and thus could greatly amplify the potential energy, so that the gained energy is then conveyed by the connecting air bubble to influence the bottom pressure. Once the rising bubble breaks down and the water on the top falls towards the bottom, the bottom pressure would decrease accordingly with the lower water level, and therefore, the oscillation periods/phases of these two parameters are consistent during the whole air–water rising and geysering process. Consequently, it is reasonable to conclude that the motion of air-pocket is of vital importance in the overall air–water evolution, which concurs with the point raised by many other researchers, such as Vasconcelos and Wright (2005), who asserted that the hydrodynamics in the vertical riser is greatly involved with the presence of air-pockets.

4.4 Energy evolution analysis for transient air–water geysering process with different initial conditions

To evaluate the significance of initially trapped air volume at the downstream end, a different length of air column La is applied in former case 7 from other cases in Table 2. Moreover, in order to explore the influence of the imposed upstream pressure (e.g. stormwater collection/drainage condition), a different height of the upstream water Hw is employed in former case 8. For further demonstration and clear comparison, the results of energy evolutions for different cases listed in Table 2 and the corresponding pressure oscillation are plotted in Figure 17, and the following conclusions can be drawn:

1. When a larger air pocket is involved in the air–water interaction (e.g. case 7), the intensified air–water dynamic evolution process (illustrated by the greatly increasing magnitude of kinetic energy in Figure 17(b)) would result in the subsequent decreasing trend of the potential energy of water ‘E_pw’ (shown in Figure 17(a)). During this process, a relatively large-size air bubble may not afford continuous rising to lift the upper water column as the surrounding hydrodynamics is highly unstable. Under this condition, such a large bubble becomes breaking down quickly, and the pressure will thus suffer a larger reduction in the amplitude than in other cases with normal air bubble sizes (Figure 17(c,d)).

2. For the influence of the initial upstream water level (drainage pressure condition), the results indicate that a lower upstream pressure head (in case 8) can significantly affect the decreasing ratio of the potential energy of water (Figure 17(a)). Besides, the drainage pressure condition may have a direct influence on the oscillation magnitude of transient air–water flows in the system (Figure 17(c)). Actually, the slightly increased kinetic energy of air (in the black-dash line in Figure 17(b)) implies that the high-speed air bubble in the vertical riser can promote a higher pressure head (Figure 17(d)).

Figure 17. Energy evolution and pressure oscillation for different cases: (a) ratio of the potential energy (E_pw) of water to the initial total energy (E_t0); (b) kinetic energy of air; (c) pressure variations at the centre of the riser bottom; (d) ratio of pressure head (pressure/ρg) to the initial water level (H_w).

In summary, the energy evolution analysis is considered a universal method applicable to various system configurations, which can either provide practical methods to mitigate the intensity of the formed transient air–water flows or improve the previous experimental configurations with further extended applications to better reproduce the field-observed sewer geysers. For example, it is quite possible that the pumping pressurized air cavity propagating towards the surcharged riser with the transformable boundary on the riser top could make a difference in the dynamically varied air–water flow evolution to sustain consecutive energy transformation of potential and kinetic energy between air and water phases.

5. Summary and conclusions

This paper studies the transient air–water interactions and the induced severe spring-like geysering consequences in urban drainage pipeline systems. To this end, a full-2D model and corresponding CLSVOF scheme are proposed and validated for the investigation, which has made significant breakthroughs to simulate complicated geyser eruption process capturing the interface reconstruction dynamics with acceptable accuracy and less time-consuming. Based on the validated full-2D model, the main contributions and significances of this study include:

1. The detailed velocity and pressure distribution patterns are retrieved and found to be highly relevant to the air–water interface deformation and energy-exchanging process. The velocity vector field of the air phase was firstly taken into consideration, which indicates that the velocity distribution of the air phase with different driven forces under different ventilation conditions would perform quite distinct patterns. Both the magnitude and flow direction of the velocity vary regularly along the interface profile. This finding can be further extended to quantify the reduction of drainage capacity due to air cavities in realistic USDS.

2. Energy evolution analysis for both water and air phases during the overall process including the air pocket propagation and bubble oscillation is conducted for the first time in this study, which reveals that the intermittent eruption of geyser could be considered as the conservation and release of energy including the transformation of potential and kinetical energy of both air and water. Accordingly, this universal finding is helpful to provide practical methods to mitigate drastic sewer geysers for better design and management of stormwater collection and drainage conditions in realistic USDS.

3. For the frequently rising surge in the vertical riser, the corresponding pressure oscillation is highly correlated to the intensity of the geyser. Based on the exploration of comparative pressure variation at the bottom of the vertical riser from the extensive simulations, the energy analysis is found to be the preferable method to clarify the geyser’s intrinsic properties. Besides, the initially entrapped air volume and riser ventilation conditions can greatly affect the oscillation pattern and peak pressure, while the upstream water height is highly related to the oscillation magnitude.

Finally, the preliminary results and findings of this study confirm the feasibility and applicability of the proposed full-2D model and numerical method for efficient and accurate simulations of transient air–water flow in drainage pipeline systems. Nevertheless, it is also worthy of noting that more investigations on transient air–water flows are required to fully consider and include more influence factors (such as the tunnel relative length, ventilation size, and location conditions) in realistic drainage pipelines, so as to further complete the theoretical development and improve the practical applications.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Natural Science Foundation of China [grant number 52000142 and 51808396]; Shanghai Pujiang Program [grant number 20PJ1417500]; Research Institute for Sustainable Urban Development, Hong Kong Polytechnic University [grant number 1-ZVWM and 4-ZZNF].