## Abstract

The accurate computer simulation of pipe flow is of great importance in the design of urban drainage. The Preissmann box scheme is usually used to model a wide range of subcritical and supercritical flows. However, care must be taken over the modelling of transcritical flows since, unless the correct internal boundary conditions are imposed, the scheme becomes unstable. In this paper, using the scheme in conjunction with the reduced momentum equation and applying boundary condition structure inherent to subcritical flow to all regimes, is an approach that enables efficient numerical simulation of transcritical flows in pipe networks. The approach includes three steps. First, a unified mathematical model which is based on the Preissmann slot model is derived. Second, the Preissmann box scheme is used to solve the set of equations, by analyzing and discussing the origin of the invalidity of applying the scheme, and a numerical model suitable for transcritical flow is proposed by the method of changing the convection acceleration term. Third, the numerical model is assessed by comparison with analytical, experimental and numerical results. The proposed models verified that this method can make the Preissmann box scheme applicable to the computation of transcritical flow in pipes.

## INTRODUCTION

Prevention of flooding in urban areas caused by insufficient sewer systems has become an important issue. With increased property values of buildings and other facilities, potential damage from prolonged flooding can easily extend into the millions of euros (Schmitt *et al.* 2005). Therefore, numerical simulation of unsteady flows in urban storm sewer systems is a very crucial task (Casulli & Stelling 2013).

The water flow is generally unpressurized, but the pipe has a limited storage capacity. When an extreme rainfall occurs, surface run-off from an urban basin can be relatively intense, exceeding the capacity of sewer pipes. In such a situation the pipes may become completely filled with water, which results in a transition from free-surface to pressurized flow (Szydłowski 2014). During such a transition, high transient phenomena appear and may cause structural damage to the systems, and particularly at the beginning of the transition bore, water interactions may raise. So it is necessary to overcome dissimilarity among the sets of equations describing pressurized and unpressurized situations and the transient between them. For numerical simulation of the hydraulic transient, it is difficult to deal with the mixed interface of transitions between the free-surface and pressurized flows, and some mathematical approaches have been developed to date, among which the shock-capturing approach and shock-tracking approach are the two most frequently used ones. The shock capturing approach computes free-surface and pressurized flows by the Preissmann slot method, that is, by adding a narrow slit at the top of the pipe, and the pressurized flow is regarded as the free-surface flow. The main advantage of this approach is that there is no need to track the interface explicitly. It has been applied to urban drainage systems (Ji 1998; Szydłowski 2014; An *et al.* 2018). However, this method exposed a major shortcoming: the mathematical model cannot be used for subatmospheric pressures (Politano *et al.* 2007). The shock tracking approach consists in the description of obtaining separately free-surface and pressurized flows through different sets of equations, and bubbles and negative pressure in the pipe network can be treated (Fuamba 2002). Since the interface between free-surface and pressurized flow is a moving interface, it assumes the introduction of additional equations (compatibility conditions and characteristics), which requires a very complicated algorithmic structure to track the existence and current position of discontinuities and accordingly to switch from one solving procedure to another. Trajkovic *et al.* (1999) and Hatcher & Vasconcelos (2014) achieved satisfactory results by using the shock tracking approach. However, the water depth of the flow cross-section in free-surface could affect the calculation of the free surface, and it may unduly limit the strength of the prediction interface.

In addition, for subcritical flows, the Preissmann box scheme is one of the most widely used in industrial models (Freitag & Morton 2006). The implicit finite-difference method is unconditionally stable and extremely robust, and it is also a valid method for fully supercritical flow. But if transcritical flow is taken into account under unsteady conditions, the scheme is invalid, for flows in which two regimes coexist (Sart *et al.* 2010; Kane *et al.* 2017). A common procedure allows the treatment of transcritical flows with a numerical method created for subcritical flows. This method can be applied to a degenerate system of equations, reducing or totally dropping the influence of the convective momentum term as the flow becomes supercritical (Kutija 1993). To reduce artificial diffusion and ensure numerical stability, the convective acceleration term was reduced locally by multiplying it by a factor (where *F* is the Froude number of the flow) (Djordjević *et al.* 2004). In urban flood plains, Abebe *et al.* (2016) took the simplified models which ignore convective acceleration terms in the momentum equations to simulate one-dimensional flows for supercritical and transcritical flows, and the research showed that it may be effectively used to simulate urban flood plains without a significant loss of accuracy. Many scholars have done a lot of work to solve the set of equations under the coupling of free-surface and pressurized ﬂows. Some of them are based on the finite volume method, and the other part is based on the Preissmann slot method. For the former, using the finite volume and kinetic framework, Bourdarias *et al.* (2011) proposed an approximation of the source terms following the principle of interfacial upwind with a kinetic interpretation; based on the ﬁrst-order Roe scheme in the frame of ﬁnite volume methods, Fernández-Pato & García-Navarro (2014) also built a numerical model; Abebe *et al.* (2016) used the MIKE 11 flow model, Kutija method, and the Roe scheme for modelling one-dimensional flow. For the latter, Ji (1998) used the ‘superlink’ algorithm, staggered grid and implicit scheme to improve the stability and speed of the computation; Jhai *et al.* (2000) proposed the flux-difference splitting scheme of Roe to simulate flows in a closed conduit; Kerger *et al.* (2011) used a first-order explicit finite volume Godunov-type scheme to solve the set of equations; An *et al.* (2018) proposed a new hybrid numerical ﬂux solver, in which the upwind ﬂux solver was combined with the centered ﬂux solver. In a word, their main goal was to use the same algorithm designed for pure subcritical flow for supercritical and transcritical flows, and these algorithms were based on imposing one boundary condition at each end of the computational domain. However, very little attention was devoted to discussing or analyzing the origin of the difficulties associated with the application of the Preissmann box scheme to transcritical flow.

## GOVERNING EQUATIONS AND DISCRETIZATION SCHEME

### Governing equations

*R*is the hydraulic radius. The mass equation reflects the water balance in the pipeline, that is, the change rate of storage (the first term) should be equal to the change rate of flow along the pipeline (the second term). The momentum equation reflects that the combined action of gravity and pressure makes the water flow overcome the energy loss caused by the inertia force and friction and obtain acceleration. Among them, the first term reflects the local acceleration at a fixed point, and the second term reflects the convective acceleration caused by the uneven spatial distribution of velocity. The above two terms are called inertia terms. The third term reflects the influence of water depth, which is called the pressure term. The fourth term is the friction loss in flow inside the boundary. In this paper, a horizontal pipeline is studied, so the gravity term is neglected.

In the numerical simulation of an urban pipe network, the governing equations are a good description for the partially full flowing pipe. For the completely full flowing pipe, the Preissmann slot concept was introduced: extending along the entire length of a closed pipe is a hypothetical slot extending infinitely upward. The use of the slot allows for the application of a single form of 1D shallow water to be applied to describe three distinct flow situations: partially full, completely full or transient between both flows. The width of the slot is very narrow, so that the slot does not introduce appreciable error in the volume of water (see Supplementary Material 1 for detailed demonstration, available with the online version of this paper). Therefore, the free-surface/pressurized flows of the pipe network are governed by Equation (1). Meanwhile, the water height () of an open channel is equivalent to the water level of a pipe flowing partially full and pressure head of a pipe flowing completely full; the width of the water surface () of the open channel is equivalent to the width of the water surface of a pipe flowing partially full and the slot width of a pipe flowing completely full.

### Discretization scheme

*f*is any function,, which are temporal weighting coefficients. The grid points are variably spaced with for all . are spatial weighting coefficients. For the scheme is fully explicit and for it is fully implicit. The Preissmann box scheme is unconditionally stable for , .

In order to implement the Preissmann box scheme, this study considers a finite domain and uses () grid points [] which are variably spaced with distance . Since there are two unknowns and at each computational node , a total number of () are unknown. Since every reach connecting two grid points in the computational domain can get two discretized equations (Equation (3); Equation (4)), () grid points can get discretized equations. Thus, two more equations are needed, which are obtained by boundary conditions.

## METHODS

*Q*is usually prescribed upstream and the height

*h*downstream. Both

*Q*and

*h*are prescribed at the upstream boundary for supercritical flow. The model descriptions below, in order to use the one-point boundary condition at each end for supercritical and transcritical, focus on the momentum equation and discuss how to deal with supercritical or transcritical flow conditions by the convective acceleration term (Longxi

*et al.*2015). In order to deal with the convection acceleration term more clearly, this study introduces two coefficients and in front of the expansion term of the convective acceleration term, respectively: where , . Because , this yields:

The eigenvalues of the equation are obtained: . Thus, the slopes of the two characteristic lines are given by .

The boundary conditions of the governing Equation (1) are provided by the method of characteristics, and the number of characteristics entering the domain in the forward time direction at a point on the boundary gives the number of boundary conditions required at that same point. Firstly, in order to determine the physical boundary types of the transcritical flow, the slopes of the two characteristic lines should be identified. Secondly, in order to quickly identify the slopes of the two characteristic lines, the model is based on the previously remarked range of coefficient and coefficient , and the limit values 0 and 1 are used respectively. Finally, these limit values of the coefficients have the following four combinations:

- (1)
When and , the complete convective acceleration term is preserved. The slopes of the two characteristic lines are given by . For subcritical flow, , where exceeds

*u*, the slopes of the two characteristic lines have opposite signs and therefore enter the domain from two different boundaries, which correspond to the use of a one-point boundary condition at each end. For the supercritical flow condition, , where*u*exceeds , the slopes of the two characteristic lines are two positive signs and therefore enter the domain from upstream only, which correspond to the use of two boundary conditions upstream. - (2)
When and , the complete convective acceleration term is deleted. The slopes of the two characteristic lines are given by . It is apparent that these slopes are not independent of the Froude number, and moreover, these always have opposite signs, which correspond to the use of a one-point boundary condition at each end.

- (3)
When and , the term is preserved. The slopes of the two characteristic lines are given by . From this expression for the slopes of the two characteristic lines, it is seen that they do not depend at all on the Froude number. The slopes of the characteristic lines always have opposite signs, which correspond to the use of a one-point boundary condition at each end.

- (4)
When and , the term is preserved. The coefficient matrix

*M*has no real eigenvalue. Thus, the slope of the characteristic line does not exist.

## RESULTS AND TEST CASES

In all the results based on the equations with changing convective acceleration term, it is known that combination (1) is suitable for a single flow regime such as pure subcritical and pure supercritical. For transcritical flow, the slope of the characteristic line is affected by the Froude number. If the flow is from supercritical to subcritical, that is, the physical boundary type is 2 + 1, the number of governing equations is (), which is different from the number of unknown variables (), and a unique solution cannot be obtained. The analysis presented shows that the Preissmann scheme cannot be used to simulate transcritical flow using the combination (1) method, and vice versa. In the case of combination (2) and combination (3), the slopes of the two characteristic lines are not affected by the Froude number, and these always have opposite signs, which corresponds to the use of a one-point boundary condition upstream and downstream, and sets up a closed system in the computational domain. For a domain of () computational points, there are two unknowns at each node. Making sure that () unknowns correspond to () discretized equations, a unique solution can be obtained. The last combination (4) is not suitable for any flow regime. To sum up, combination (1) and combination (4) should be abandoned, combination (2) and combination (3) are suitable for any flow regime, including transcritical. Next, in the verification section, for the numerical simulation of transcritical flow, the model based on combination (2) and combination (3) is validated by comparison with analytical solutions and experimental data for various cases (Wiggert 1972; MacDonald *et al.* 1995), and the influence of the value of the coefficient on the model is discussed.

### Numerical simulation of transcritical flow in a channel

For transcritical flow in a channel, the four study cases are used to verify the proposed model. These cases are four combinations of the inflow/outflow boundary conditions, including subcritical inflow to supercritical outflow, supercritical inflow to subcritical outflow, subcritical inflow to subcritical outflow and supercritical inflow to supercritical outflow. The channel is 1,000 m long, and the bottom width is 10 m. To describe the rough channel bottom, Manning's roughness coefficient is set to be 0.02 s/m^{1/3}, and an upstream inflow water discharge of 20 m^{3}/s is given herein.

For transcritical flow in the channel, the bottom elevation, critical level, analytical solution, and the simulated water depth at m, s for combination (2) and combination (3) are exhibited (see Figure 1), and the results of combination (2) and combination (3) are exactly the same. This is because the flow is constant, and so , regardless of whether its coefficient is 0 or 1, and it will not affect the result of the model, so the results of the two new models are exactly the same. However, subcritical flow changes to supercritical flow at (see Figure 1(a)), supercritical flow changes to subcritical flow at (see Figure 1(b)), subcritical flow changes to supercritical flow at , back to subcritical flow at (see Figure 1(c)), and supercritical flow changes to subcritical flow at , back to supercritical flow at (see Figure 1(d)). To compare with the analytical solutions, the flow pattern is slightly changed in advance. This is due to ignoring the term of the convection acceleration term. That is to say, the energy losses of hydraulic jump/hydraulic drop generated during the transformation of the flow pattern are not neglected. Therefore, the water depth is greater than the value of the analytical solution, so that the flow pattern is slightly changed in advance. It can be seen that the overall simulation results are satisfactory, so the proposed models can be applied to get a good prediction of the flow in single transcritical flow (see Supplementary Material 3 for detailed demonstration, available with the online version of this paper).

### Numerical simulation of transcritical flow in the pipe

The experimental apparatus (see Figure 2(a)) used by Wiggert is the fifth study case. The Manning coefficient is 0.01 s/m^{1/3} and the pressure wave speed is 20 m/s. A water level of 0.128 m and no discharge are considered as initial conditions. Then a wave coming from the left side causes the pressurization of the pipe. The water level and pressures in the apparatus are monitored by four gauges within the pipe. To specify boundary conditions for the numerical method, Wiggert measured upstream pressure head and downstream pressure head at the pipe extremities (see Figure 2(b)).

The numerical results of the proposed models and experimental data for pressure heads at the four gauges are shown in Figure 3. On the whole, the results of combination (2) and combination (3) are only slightly different. Due to the unsteady flow, . The difference between the results of the two models is due to these different values of the coefficient of the term. However, the difference is very small, which shows that the influence of the term on the model is almost negligible. When a rise is observed in the upstream water level of the flume, a transition develops upstream and propagates towards the downstream end of the pipe, and when the pressure head reaches the top of the pipe, there is no violent fluctuation. This is because the convective acceleration term tends to expand the fluctuation. However, the proposed new models discard the term, which plays a major role in the convection acceleration term, so only small numerical oscillations are observed. Through comparison, the numerical simulation results are basically consistent with the experimental results.

## CONCLUSIONS

The method of reducing the term of the convective acceleration term and preserving the structure of the boundary conditions inherent to subcritical flow is commonly employed in practice as a means of modelling transcritical flow in hydraulics. According to the numerical simulations, the following conclusions are drawn:

- 1.
In the convective acceleration term, the term plays a decisive role for the numerical model results, while the term has negligible influence on the simulation results.

- 2.
Whether the term of the convective acceleration term is preserved directly affects the occurrence of hydraulic jump/hydraulic drop in the process of flow transition, so as to affect the water depth of flow after the flow regime conversion.

- 3.
When the Preissmann slot method is used to simulate the free-surface and pressurized flow, a numerical oscillation occurs when the flow reaches the top of the pipe, and the convection acceleration term will enlarge the oscillation. However, since the term of the convection acceleration term is omitted from the proposed model, the numerical oscillation at the top of the pipe is reduced, which makes the proposed models more consistent with the real flow regime than previous models.

From the numerical results, this proposed method of reducing the term can be shown to be simple and effective when using the Preissmann scheme to calculate transcritical flow in pipe networks.

## ACKNOWLEDGEMENTS

This research was supported by Major Project of National Natural Science Foundation of China (No. 51739009), the National Natural Science Foundation of China (No. 51509224), the Scientific and Technological Research Program of Henan Province (No. 182102210017), and the Scientific and Technological Research Program of Henan Province (No. 162102310522), for which the authors are grateful.