## Abstract

Water distribution systems are basically designed to convey pressurized flow; however, in some situations such as the intermittent operation of the system, the network may experience a transition between free-surface and pressurized flow. On the other hand, combined sewer systems, designed basically for free-surface flow, may undergo pressurization due to extreme rainfalls. During transient flow, free-surface flow changes into a pressurized flow (and vice versa) which could be accompanied by intensive transient pressures causing structural damages to the system. In a pipe filling process, the existing air can be entrapped due to improper ventilation, intensifying the transient pressures in some cases. In this paper, the two-component pressure approach (TPA) and a Harten–Lax–van Leer Riemann solution are applied to model transient flow. The model is validated by comparing its results with the analytical solution of three simple examples, and then the model with an air chamber as the downstream boundary is used to simulate a literature experimental setup that includes air pocket entrapment. As the volume of the air pocket decreases, the errors of the model increase due to the inherent deficiency of the one-dimensional model. Furthermore, it is recommended to limit the Courant number to 0.5 for high acoustic wave speeds.

## INTRODUCTION

Urban water systems may undergo extreme events due to high pressures caused by a transition between free-surface flow and pressurized flow (Zhou *et al.* 2002; Bousso *et al.* 2012). When the ventilation of the system is not appropriate, air entrapment can intensify these pressures (Liu *et al.* 2018). In a combined sewer system, designed as a free-surface flow network, severe storms may fill up the pipes and change the flow into a pressurized flow. The resulting transient pressures can cause structural damages and geysering events (Zhou *et al.* 2002; Wright *et al.* 2017). In 1995 in Edmonton, Alberta, Canada, a big storm led to flooding and infrastructure damage, because of rapidly filling up the combined sewer system. On the other hand, a water distribution system, designed as a pressurized network, is under intermittent operation, and the frequent transition between free-surface flow and pressurized flow can cause damages, water loss, and contamination (Christodoulou & Agathokleous 2012; Kumpel & Nelson 2016). Christodoulou & Agathokleous (2012) in their study showed that during a 2-year intermittent operation of the urban water distribution system in Cyprus, the breaks increase 30–70 percent per year.

Free-surface flow and pressurized flow have different features leading to some challenges in numerical modeling. Acoustic wave speed in pressurized flow is two or three times higher than in free-surface flow (Vasconcelos *et al.* 2006). Moreover, the possible negative pressure in pressurized flow does not exist in the free-surface flow (Vasconcelos *et al.* 2006). In addition, free-surface flow can have different regimes, which are difficult to be included in the model, especially in boundary conditions (Vasconcelos *et al.* 2006). There are two different governing equations in each pressurized and free-surface flow. Accordingly, there are two different approaches to deal with the transient flow: shock-fitting approach and shock-capturing approach.

In the shock-fitting approach, or in other words interface tracking method, in each time step, the location of the pressurization is detected, and in every side of the interface, the flow is modeled by its own governing equations (Fuamba 2002). With respect to the numerical solution, this approach has three main categories: method of characteristics, the finite volume method, and the finite difference method. The ability to estimate the location of the pressurization surge, the ability to consider a wide range of acoustic wave speed, and the ability to model subatmospheric pressure are the advantages of this approach (Bousso *et al.* 2012). The most prevalent method with this approach is the rigid column method (Martin 1976; Zhou *et al.* 2002; Vasconcelos & Leite 2012; Hatcher *et al.* 2015; Fuertes-Miquel *et al.* 2019). However, some references have separated this method from the shock-fitting approach (Vasconcelos *et al.* 2006).

In the shock-capturing approach, there is no need to track the pressurization front by using a set of governing equations for both pressurized and free-surface flow. This approach is fast and simple to implement for complex systems despite its drawbacks and can be used in commercial software if its defects are handled (Bousso *et al.* 2012). Two major methods with this approach are the Preissmann slot method (PSM) and the two-component pressure approach (TPA).

The PSM is the first and most common shock-capturing method, presented in 1961. In this method, in pressurized flow, a virtual slot and a hypothetical free-surface are considered on the crown of the conduit (Trajkovic *et al.* 1999; Vasconcelos *et al.* 2006; Lieb *et al.* 2016). Simple concept and the feasibility to distinguish transient flow in a location where there is no bore front are the advantages of this method (Vasconcelos & Wright 2007). However, it has some drawbacks such as the inability to model low pressures and some numerical oscillations at the pressurization front (Vasconcelos *et al.* 2006; Vasconcelos & Wright 2007; León *et al.* 2009). Many studies have tried to overcome these defects (Arai & Yamamoto 2003; León *et al.* 2009; Kerger *et al.* 2011; Malekpour & Karney 2014a).

Vasconcelos *et al.* (2006) have presented the novel TPA to model transient flow, based on which the pressure head term in the Saint-Venant equations is divided into a hydrostatic and a pressurized component. The fundamental assumptions of this approach to calculate the pressurized component are the incompressibility of the flow and flexibility of the pipe. Therefore, the celerity in free-surface flow and acoustic wave speed in the pressurized flow are equal, whereby changes in cross-sectional area can be converted to a pressurized component. This method has been developed in recent studies (Vasconcelos & Wright 2007; Sanders & Bradford 2010; Vasconcelos & Marwell 2011; Trindade & Vasconcelos 2013; Malekpour & Karney 2014b).

Although the shock-capturing approach is simpler and more efficient to model transient flow, spurious numerical oscillations, particularly around the pressurization front, have a significant impact on the stability and accuracy of the model, which can cause failure in the model. Thus, considerable efforts have been made to diminish these oscillations (León *et al.* 2009; Vasconcelos *et al.* 2009; Kerger *et al.* 2011; Malekpour & Karney 2014a, 2014b).

Although such a complex phenomenon is three-dimensional, when it comes to a pipeline or network scale a one-dimensional model is more useful. However, the basic assumptions of the one-dimensional models limit the ability of the model to predict the first peak pressure that is the highest one. These assumptions are: (1) the interface between the air and water is vertical; (2) each section is occupied with water or air; (3) the steady flow friction factor is used; (4) the pressure in the air pocket is uniform; (5) the dissolution and decomposition of the air pocket is discarded (Zhou *et al.* 2002).

In this paper, the TPA applied with the Harten–Lax–van Leer (HLL) solver is used to model transient flow. To treat the numerical oscillations of the TPA model, the numerical viscosity around the filling bore is increased using the proposed method in Malekpour & Karney (2014b). The aim is to investigate the capability of these methods in more complex problems with air pocket entrapment. The model is validated by solving three scenarios for a simple frictionless pipe system and comparing the results with the analytical solution based on the principles of the TPA. Then the model is used to simulate a series of literature experiments (Vasconcelos & Leite 2012) wherein the initial flow rate was not zero, the system contained both pressurized flow and free-surface flow, and pressurization occurred due to a backward flow caused by the downstream valve closure. This experimental study can represent the actual filling process more properly in comparison to previous studies such as the Zhou *et al.* (2002) experiment where the initial flow was zero and the air could release through a downstream orifice. Therefore, the model should be able to simulate pressurized flow at upstream, free-surface flow at the downstream, backward pressurization bore, and air entrapment with no ventilation. An air chamber is considered as a boundary condition to reflect the effect of air entrapment during the filling process. Moreover, the effect of considering different Courant numbers as stability criterion has been studied.

## METHODS

### Governing equations

*A*is the cross-sectional area of the flow,

*Q*is the flow discharge,

*h*is the pressurized water head,

_{s}*h*is the distance between free-surface and the centroid of the cross-section of the flow,

_{c}*g*is the acceleration of gravity,

*S*

_{0}is the pipe slope, and

*S*is the energy grade line slope. Moreover, according to the flow regime

_{f}*h*and

_{c}*h*are calculated based on Equations (3) and (4) (Trindade 2012):where

_{s}*D*is the pipe diameter,

*θ*is the angle between the surface of the flow and the centerline of the pipe,

*a*is the acoustic wave speed, and

*A*is the initial cross-sectional area of the pipe.

_{pipe}### Finite volume solver

*n*refers to the previous time step,

*n*

*+*

*1*refers to the current time step,

*i*

*−*1/2 and i + 1/2 are related to the left and right boundaries of the cell

*i*respectively.

In the finite volume method, it is required to calculate the flux at the border of two adjacent cells with appropriate accuracy. To this purpose, the first-order Godunov numerical scheme is utilized based on which the fluxes are computed by means of the Riemann solution. Because of some numerical difficulties in implementing the precise solution, many approximate methods are introduced to solve the Riemann problem (Toro & Garcia-Navarro 2007; Malekpour & Karney 2014b), one of which is the HLL solution presented by Harten *et al.* (1983) for the first time.

### Harten-Lax-van Leer Riemann solution

*F*and

_{L}*F*;

_{R}*S*and

_{L}*S*;

_{R}*U*and

_{L}*U*are representatively fluxes; wave speeds; and conservative variable values in the left and the right cells, respectively. In the current model, the wave speeds are calculated with the presented equations in León

_{R}*et al.*(2009) (Equations (9) and (10)):where

*V*and

_{L}*V*are flow velocities in the left and the right cells respectively.

_{R}*P*is pressure,

*A*is the cross-sectional area,

*c*is celerity, and index

*G*represents the variables of a point with a depth of

*Y*. Malekpour & Karney (2014b) estimated this depth in a way to decrease the oscillations by increasing numerical viscosity around the pressurization front locally (Equation (11)):where

_{G}*K*is a constant coefficient increasing viscosity, which is recommended to be taken as 1.4 when the water depth is more than the pipe diameter and 1.001 for other situations.

_{a}*NS*determines the search span; based on numerical experiments, Malekpour & Karney (2014a) suggest that this span should be two or three times longer than the pipe diameter and the value of

*NS*should be at least three. Therefore, in this paper, when the ceiling value of three pipe diameter over two computational cell length is lower than three, a minimum search span is assumed as seven computational cells.

### Friction source term

*Q*. First,

_{i}^{n+1}*Q**in each cell is computed by neglecting the friction term; then

*Q*can be derived from Equation (13):where

_{i}^{n+1}*c*is the resistance parameter calculated using Equation (14):where

_{D}*n*is Manning's coefficient and

_{m}*R*is the hydraulic radius (Sanders & Bradford 2010).

_{h}### Controlling the wet/dry status of a cell

*A*) (Trindade & Vasconcelos 2013). After calculating the new cross-sectional area of the flow (

_{min}*A*) and comparing it with the minimum cross-sectional area, the wet/dry status of the cell can be distinguished according to Equation (15), in which is a Boolean variable demonstrating the wet/dry status of each cell.

_{i}^{n+1}### Stability criteria

*et al.*2006). Therefore, the time step size is not constant during each simulation. At first, an initial time step duration has been assumed; then in subsequent iterations, having the celerity and velocity in the previous time step, the time interval is derived from Equation (17):where

*Cr*is the Courant number,

*dt*is the time step duration,

*dx*is the cell length,

*u*is the velocity, and

*c*is celerity.

### Air entrapment

*et al.*(2015) (Figure 1) and the boundary water head and discharge are derived from Equations (18)–(21). Since the air pocket forms at downstream with no orifice, this is a rational choice.where

*Q*is the discharge at the downstream boundary,

_{nc}*Q*is the discharge of the air chamber,

_{c}*H*is the water head at downstream boundary,

_{nc}*K*is polytropic coefficient (assumed equal to 1.2 according to the reference paper),

*H*is the air pressure,

_{air}*H*is the atmospheric pressure,

_{atm}*V*is the air pocket volume,

_{air}*dt*is the time step length,

*a*is acoustic wave speed,

*g*is acceleration of gravity,

*A*is the cross-sectional area at the boundary cell,

_{nc}*H*is the water head at left cell of the boundary,

_{nc−1}*V*is the velocity at left cell of the boundary,

_{nc−1}*f*is the Darcy–Weisbach friction factor, and

*D*is the pipe diameter.

## RESULTS AND DISCUSSION

### Numerical model validation

To validate the numerical model, a simple system with a frictionless horizontal pipe and two reservoirs, upstream and downstream, are considered. The pipe length is 400 m and its diameter is 2 m. Different scenarios are investigated by changing the water depths in the pipe and reservoirs. First, the basic ability of the model is tested by simulating pressurized flow and free-surface separately, and then the model is implemented for the pipe filling process with both positive and negative hydraulic bore by modeling pressurization of a pipe containing an initial water depth with a positive hydraulic bore from upstream to downstream and with a negative hydraulic bore from downstream to upstream, resembling the pipe filling process when a huge amount of water enters a sewer pipeline. A horizontal frictionless pipe has been chosen so that the results could be validated with an analytical solution. The results are compared with an analytical solution of the energy equation, the continuity equation, and the momentum equation using the TPA concept (variable cross-sectional area along with assuming incompressibility of water and elasticity of the pipe).

In the first scenario, the initial water depth in the pipe and reservoirs is 1.5 m, and then the water depth of the upstream reservoir abruptly is increased creating a positive hydraulic bore filling and pressurizing the pipe. The acoustic wave speed is assumed as 10 m/s. In this study, different Courant numbers were considered (0.3, 0.4, 0.5, 0.7, 0.8, 0.9, 1). The resulting trajectories for two different Courant numbers (0.5, 0.8) were compared with the analytical solution (Figure 2). When the maximum Courant number is assumed as 0.8, some oscillation can be observed around the hydraulic bore and the average difference between analytical and numerical results is 2 percent. By decreasing the maximum Courant number to 0.5, this average has declined to 0.7 percent but the difference between the average velocities of the hydraulic bore in the analytical solution and model is increased from 0.08 percent to 2 percent. These errors have been used to show the effect of changing the Courant number, in order to reduce numerical instability, on the accuracy of the results. According to Figure 3, lowering the Courant number decreases the speed of propagation of the filling bore in the model. On the other hand, in lower Courant numbers the time step duration became very small, increasing the computational time with no specific effect on the stability of the solution.

In the following results, the Courant number is limited to 0.5. In Figure 4 the velocity and pressure changes of the model are compared with the analytical solution over time in different locations. The average difference for velocity is significant (12 percent) and for pressure is 1.1 percent.

In the case of higher acoustic wave speeds, the model undergoes high numerical oscillations because the TPA is based on the assumption that the celerity of the free-surface flow is equal to acoustic wave speed in the pressurized flow and the additional numerical viscosity cannot eliminate the oscillation at the first seconds of the simulation. Figure 5 demonstrates the comparison of the trajectories between numerical and analytical solutions over time assuming a higher acoustic wave speed of 50 m/s. There are some oscillations around the pressurization front, declining over time. The 15 percent difference in maximum water height between the two solutions in the fifth second is decreased to 1 percent in the twentieth second; however, the difference in the average hydraulic bore velocity has increased 1 percent. The same pattern has been observed in the second scenario when we have a negative hydraulic bore from downstream toward the upstream.

### Model application

The model is applied to simulate the experiments conducted by Vasconcelos & Leite (2012). In this study, the experimental apparatus consists of two reservoirs, a pipe, a knife gate valve, and a pump system providing the water circulation. In the initial steady situation, flow is pressurized in the upstream part of the pipe and it is free-surface flow at the downstream part of the pipe. By closing the downstream valve, a backward filling bore is generated causing air pocket entrapment.

Different input discharges resulted in the formation of air pockets with different sizes. Three experiments, with a horizontal pipe, of these series, are used to simulate in this study (Table 1). According to Vasconcelos & Leite (2012), the values of the pressure head, time, flow discharge, and volume of the air pocket are normalized. The major distinction between these experiments is the initial volume of the air pocket.

Experiment . | Velocity (m/s) . | . | . |
---|---|---|---|

A | 0.511 | 0.4021 | 3.4 |

B | 0.605 | 0.476 | 1.078 |

C | 0.644 | 0.507 | 0.38 |

Experiment . | Velocity (m/s) . | . | . |
---|---|---|---|

A | 0.511 | 0.4021 | 3.4 |

B | 0.605 | 0.476 | 1.078 |

C | 0.644 | 0.507 | 0.38 |

Figure 6 shows the pressure hydrograph for the middle point and endpoint of the pipe in experiment A. The percentage error of the first peak pressure head between numerical results and experimental measurements at the middle point of the pipe is 6.5 percent and 8.5 percent at the endpoint of the pipe. However, this difference for the next peak pressures is significant, because, after the first period, the air pocket is compressed and deformed making the assumptions of the one-dimensional questionable. Besides, when air pocket is considered as an ideal gas regardless of its energy absorption and dissolution in water, it causes such significant error. Furthermore, at downstream a pressure peak delay is observable. This delay might be because of assuming an air chamber at the endpoint instead of the knife gate valve and ignoring the water impact on the valve.

The pressure hydrograph of experiment B of the middle and endpoint of the pipe is depicted in Figure 7. The peak pressure is overestimated about 12.7 percent at the middle point. It is underestimated about 30.8 percent by the model. The error is more significant in comparison to experiment A, since in this experiment initial air volume is smaller, absorbing less energy and intensifying the effect of the collision between water and valve.

As Figure 8 demonstrates, in experiment C the results of the model are not acceptable due to the very small initial size of the air pocket (57.6 percentage error for the middle point and 72.2 percentage error for the endpoint). Apart from the aforementioned reasons, this meaningful error could have been because of the difficulty in measuring the volume of the small air pocket.

As a recommendation for future work, combining one of the shock-capturing methods with a multilayer Saint-Venant system (Audusse & Bristeau 2007) could be compelling. Using a multilayer Saint-Venant system in the TPA might be able to surpass this assumption. Since we have different layers each one has vertical hydraulic bore but there is no need for the total bore to be vertical. In addition, assuming a thin upper layer near the crown of the pipe and locally increasing the acoustic wave speed of the layer might decrease the numerical oscillations.

## CONCLUSIONS

The TPA has an acceptable ability to capture the first peak pressure in the middle of the pipe as a one-dimensional model with its inherent deficiencies; however, at the end of the pipe, the error is considerable, which might have occurred because of the simple air chamber boundary conditions at downstream. Moreover, limiting the Courant number to 1 cannot guarantee the stability of the solution especially for the higher acoustic wave speed. It is recommended to limit the Courant number to 0.5, although it causes a delay in the propagation of the hydraulic bore. Besides, in higher acoustic wave speeds, the Courant number criterion might result in selecting infinitesimal time steps, increasing the duration of the computation dramatically, without considerable effect on the stability and accuracy of the results.

As mentioned before, air entrapment is basically a three-dimensional phenomenon and the basic assumption of the one-dimensional approach (vertical air–water interface) limits the ability of the model to anticipate the highest peak pressure, which can be used in designing and operating urban water systems. While three-dimensional models might be able to model transient two-phase flow, they are not efficient to model this flow in a real water or wastewater networks. To model air pocket entrapment in detail, one should consider air–water interaction, air incorporation through a hydraulic jump, and the air pocket split and motion. These can hardly be incorporated in a one-dimensional model. These sorts of studies are tiny steps toward finding a powerful and efficient tool to model transient flow to design, manage, and maintain intermittent water systems and combined sewer systems. As climate change causes more extreme events such as floods and droughts, well-designed and well-operated stormwater management systems and intermittent water supplies, as a demand management strategy, become noteworthy to provide safe and equitable public services.

## ACKNOWLEDGEMENTS

The authors gratefully acknowledge Dr Jose G. Vasconcelos for providing valuable experimental data. The authors also wish to thank Dr Ahmad Malekpour for his helpful advice and assistance during the course of this research.