Application of two-component pressure approach and Harten–Lax–van Leer (HLL) solver to model transient flow with regard to air entrapment

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.

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.

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.

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.

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.

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.

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.