## Abstract

In this paper, the augmented version of finite volume Harten–Lax–van Leer (HLL) solver, including source terms, is extended to free-surface and pressurized mixed pipe flows over complex and frictional topography. This augmented HLL Riemann solver is employed for the flux approximation at the cell interface, where source terms are split into two parts based on the wave propagation speed. The friction term is treated using a splitting implicit method to maintain stability over dry beds. The Preissmann slot method is adopted to reproduce pressurized flow in the conduit. The performances of the numerical model are investigated by several numerical tests and compared with existing methods showing clear improvements.

## HIGHLIGHTS

An augmented HLL Riemann solver is extended to mixed pipe flows over complex topography.

Superiority simulating steady flows over complex topography with different hydraulic flow regimes.

Capability simulating transient dambreak pipe flow, including wetting or drying front and the calculation of the friction term.

## INTRODUCTION

In the context of climate change and rapid urbanization over the last decades, risks of urban flooding have greatly increased worldwide with increased runoff at local catchment, that devastating consequences occurred on both infrastructures and human safety (Zheng *et al.* 2015). To improve risk assessment and evacuation management, flood inundation models are crucial tools to predict flood hydrodynamics in cities (Mignot *et al.* 2006; Hunter *et al.* 2007; Arrault *et al.* 2016; Teng *et al.* 2017; Chen *et al.* 2018). As an integrated component of urban flooding processes, drainage networks composed of junctions and pipes are commonly modeled by different numerical approaches (Martins *et al.* 2017; Li *et al.* 2020). Pipe flows, which are the basic element of drainage network, are often modeled by the 1D shallow water equations (Kerger *et al.* 2011a, 2011b).

During urban flooding, transitions in pipes between free-surface and pressurized flow conditions can take place repeatedly. This phenomenon has been studied by many researchers (Bousso *et al.* 2012) and can be modeled using a method like the shock-fitting method (Wiggert 1972) or the shock-capturing method (Abbott & Cunge 1982; Vasconcelos *et al.* 2006; Vasconcelos *et al.* 2009). The extension of free-surface shallow water equations to the pressurized flow to capture the interface has made the shock-capturing method more attractive than the cumbersome shock-fitting method. Among them, the Preissmann slot method (PSM) has been widely used for its simplicity in model implementation. Lately, Kerger *et al.* (2011a, 2011b) extended the classical PSM by developing a negative slot to reproduce subatmospheric pressurized flow, and Maranzoni *et al.* (2015) extended 1D PSM to 2D accounting for 2D pressurized flow under the bridge. Except the PSM, Vasconcelos *et al.* (2006) and Vasconcelos *et al.* (2009) developed a two-component pressure approach with the assumption that the pipe walls are elastic. In this way, various unsteady flows can be simulated effectively, including free-surface flow, pressurized flow, mixed flow and sub-atmospheric pressurized flow. PSM is used in the current study for mixed pressurized and free-surface pipe flows.

1D shallow water equations may be solved numerically by using a Godunov-type finite volume scheme (Toro & Garcia-Navarro 2007). When tackling complex topography in the framework of Godunov-type scheme, the slope source term needs to be treated properly to reproduce flow hydrodynamics. A number of numerical techniques dealing with the source term have been reviewed in the literature (Guinot 2012) to construct well-balanced schemes. For instance, Liang & Marche (2009) proposed a pre-balanced shallow water equation that no special treatment is needed for the source term considering pressure balancing. For circular cross-sections like pipe flow, Capart *et al.* (2003) reconstructed the momentum flux to account for slope and non-prismaticity by considering the balance of hydrostatic pressure with the approximated water surface level to achieve hydrostatic equilibrium, which was also employed for mixed pipe flow modeling by Sanders & Bradford (2010) and Aureli *et al.* (2015) over complex topography. The balanced property of the above-mentioned schemes is achieved by forcing equilibrium among fluxes and source terms in cases of quiescent water. However, errors on computed discharge will accumulate around the bottom variation, which can be observed in the results published in Liang & Borthwick (2009). To overcome this problem, Murillo & García-Navarro (2012) presented augmented versions of the HLL Riemann solvers for shallow water flows, by incorporating an extra wave that accounts for the presence of the source term and improved results have been obtained. In the work of Murillo & García-Navarro (2010), LeVeque (2011) and Murillo & García-Navarro (2012), only rectangular cross sections in 1D open channel flows are considered.

Those previous works are extended in the present contribution to circular cross-sections and mixed pipe flows over complex topography. Based on Murillo & García-Navarro (2010) and Murillo & García-Navarro (2012), a source term adapted to circular sections has been added to the continuity equation that along with the numerical scheme, and upgraded celerity calculations, leads to improved simulations of free-surface and mixed circular pipe flows. Following Murillo & García-Navarro (2010) and Murillo & García-Navarro (2012), the source term is discretized approximately as a succession of discontinuities separating regions of constant state and split into two parts based on the wave propagation speed. Governing equations are solved using a Godunov-type scheme with an augmented HLL approximate Riemann solve developed for mixed flow (León *et al.* 2006). The friction term is included and discretized by a splitting implicit method (Liang & Marche 2009) to maintain stability for flows over dry beds.

The paper is organized as follows: 1D shallow water equations with conserved variables *A* (flow area) and *Q* (total discharge) are presented in the ‘Governing equations’ section with the introduction of PSM for mixed flow modeling. In the ‘Numerical scheme’ section, a first-order Godunov-type scheme with augmented HLL Riemann solver, including source terms, is extended to the studied governing system for mixed pipe flows. The adaptations of the numerical method presented for simulating free-surface and mixed circular pipe flows are tested against the method from Sanders & Bradford (2010) in the ‘Numerical experiments’ section. Conclusions are drawn in the ‘Conclusions’ section.

## GOVERNING EQUATIONS

*A*and

*Q*may be expressed in a matrix form as:where

*t*denotes time and

*x*is the distance. , and are the vectors of conserved variables, fluxes and source terms, respectively, and are given by:where is the flow area, is the total discharge, is the gravitational acceleration, is the bottom slope term ( is the bottom elevation) and

*I*is the hydrostatic pressure force term, which is defined as:where

*h*represents the water depth, and is the cross-sectional width at elevation above the bottom (Guinot 2012). is the friction slope and may evaluated as:where is Manning's roughness coefficient, and is the hydraulic radius.

*a*is the conduit acoustic wave velocity. For circular pipe with diameter , the flow area

*A*and pressure term

*I*for free-surface and pressurized flow based on PSM are given, respectively, as follows (Sturm 2010):where is the wetted angle.

## NUMERICAL SCHEME

*i*is the cell index, and the superscript

*n*is the time index. and are the time step and cell size, respectively. and are the fluxes calculated at the interface and , respectively, and is the source term.

### Numerical flux calculation

*et al.*1983; León

*et al.*2006; Sanders & Bradford 2010). Taking the cell interface as an example, the fluxes can be evaluated by solving a local Riemann problem, i.e., , in the context of a Godunov-type scheme. and are the left and right Riemann states, which is used to define the Riemann problem and hence calculate the fluxes. The solution of the Riemann problem is approximated by an intermediate region of constant state separated from the left and right states and by two shocks. Therefore, the numerical flux is (Toro 2001):where and are the left and right characteristic wave speeds calculated by (Toro 1992; Fraccarollo & Toro 1995):in which

*u*is the averaged flow velocity defined as and is calculated by:

*c*calculated based on the slot width .

### Bottom slope source term calculation

*et al.*(2003) for natural channel and further developed by Sanders & Bradford (2010) for mixed flow modeling. As illustrated in Figure 1, by assuming the slope of the water surface is mild and the streamwise component of the sidewall inclination is small, momentum balance is approximated using an average elevation of water surface within the control volume instead of the local free-surface elevation . The overall pressure thrust acting on control volume is approximated as:

*et al.*(2003), when the channel is rectangular in shape and has a constant width (regardless of the slope of the water surface) and when the water surface is horizontal (regardless of channel geometry), momentum balance is achieved exactly. The bottom slope source term is treated explicitly and is calculated as:where .

In the following, the above-described schemes with Equation (17) are referred to as the ‘Sanders & Bradford (2010)’ scheme.

### Augmented HLL Riemann solver including slope source term

In this section, the numerical scheme on source term discretization presented by Murillo & García-Navarro (2010) and Murillo & García-Navarro (2012) is extended from rectangular shape to the currently studied pipe flow system, which is based on the projection of the source term onto the Jacobian's eigenvector basis.

The total source term in Equation (23) is split into two parts:

one part that corresponding to negative wave speeds, , is assigned to the cell on the left-hand side of the initial discontinuity.

another part that corresponding to positive wave speeds, , is assigned to the cell on the right-hand side of the initial discontinuity.

In subcritical flow, the wave propagates to both the upstream and the downstream. Therefore, and .

Three flow conditions related to dry bed are illustrated in Figure 3:

, and

, and

and

### Friction calculation

The value of will be replaced by the critical one when it is computed beyond the limit.

## NUMERICAL EXPERIMENTS

In this section, the performance of the proposed numerical model is tested on a number of numerical examples with a circular cross-section. Pipe flow is simulated in steady or transient state over complex topography. The numerical results obtained by using the ‘Sanders & Bradford (2010)’ scheme on a very fine grid are adopted as references for model comparisons. The stability of the numerical scheme is governed by the Courant–Friedrichs–Lewy criterion as the scheme is explicit. For all the test cases considered in this work, fixed timestep is used and the Courant number is less than 0.9 with and .

### ‘Lake at rest’ steady states

The initial conditions are defined as , over the whole computational domain so that the bump is completely submerged. The computed hydraulic variables ( and ) at are shown in Figure 4. Results show that both methods can preserve ‘lake at rest’ steady states for water surface elevation, but the proposed model outperforms the ‘Sanders & Bradford (2010)’ model considering the total discharge computed.

### Free-surface pipe flow over a bump

The reference solution for the test case is obtained on a very fine grid using the ‘Sanders & Bradford (2010)’ model. The reference for test case T2a (subcritical flow) is established using and . For the hydraulic jump (T2c), the discharge and . For the transcritical case (T2b), the discharge , and no water depth is necessary at the boundary since the flow is controlled by a critical depth (such as ) that appears at (middle of the bump).

For three hydraulic regimes, the domain length is with a mesh of uniform cell . The initial velocity of this test is set to be zero. An initial water surface elevation ζ_{ini} = 0.4 m is set for test case T2a and T2c, and ζ_{ini} = 0.2 m is given for test case T2b. The computational parameters are given in Table 1. The simulation is carried out until steady flow conditions are reached over the computational domain (after barely ).

Symbol . | Meaning . | Value . |
---|---|---|

Pipe diameter | ||

Domain length | ||

Bump length | ||

Middle of bump | ||

Minimum bottom elevation | ||

Maximum bottom elevation | ||

Cell size | ||

Initial water surface elevation (T2a) | ||

Initial water surface elevation (T2b) | ||

Initial water surface elevation (T2c) | ||

Prescribed total discharge at the upstream boundary | ||

For subcritical flow (T2a) | ||

For transcritical flow (T2b) | ||

For transcritical flow with a jump (T2c) | ||

Prescribed water depth at the downstream boundary | ||

For subcritical flow (T2a) | ||

For transcritical flow (T2b) | Free-outlet | |

For transcritical flow with a jump (T2c) |

Symbol . | Meaning . | Value . |
---|---|---|

Pipe diameter | ||

Domain length | ||

Bump length | ||

Middle of bump | ||

Minimum bottom elevation | ||

Maximum bottom elevation | ||

Cell size | ||

Initial water surface elevation (T2a) | ||

Initial water surface elevation (T2b) | ||

Initial water surface elevation (T2c) | ||

Prescribed total discharge at the upstream boundary | ||

For subcritical flow (T2a) | ||

For transcritical flow (T2b) | ||

For transcritical flow with a jump (T2c) | ||

Prescribed water depth at the downstream boundary | ||

For subcritical flow (T2a) | ||

For transcritical flow (T2b) | Free-outlet | |

For transcritical flow with a jump (T2c) |

The simulation results of three flow regimes are presented in Figure 5. Comparing with the reference, numerical results show that the proposed model provides more accurate predictions of the computed hydraulic variables ( and ) than the ‘Sanders & Bradford (2010)’ model, whatever the hydraulic flow regime is. For subcritical flow, the proposed method provides better estimations of the water surface elevation above and upstream of the bump, which can be clearly observed in the zoomed comparison in rectangular. For transcritical and hydraulic jump, the proposed method has better estimation of water surface elevation upstream the bump and more accurate shock position for hydraulic jump. Regarding the water discharge, the proposed method provides the exact value for sub- and transcritical flow, except one cell for hydraulic jump, which is located downstream the top of the bump. For the ‘Sanders & Bradford (2010)’ model, the errors accumulated upstream and downstream of the top point are due to the HLL solver adopted, which can also be observed in the results presented by Liang & Borthwick (2009).

### Mixed pipe flow over a bump

The reference solution for this test case is obtained on a very fine grid using the ‘Sanders & Bradford (2010)’ model. The reference for test case T3a (subcritical flow) is established using and . For the hydraulic jump (T3b), the discharge and . For the pressurized flow, , and free-outlet is used at the right side of the boundary. Same numerical parameters as free-surface pipe flow are used for this test case, which are listed in Table 2. The simulation is carried out until steady flow conditions are reached over the computational domain (after barely ).

Symbol . | Meaning . | Value . |
---|---|---|

Pipe diameter | ||

Preissmann slot width | ||

Domain length | ||

Bump length | ||

Middle of bump | ||

Minimum bottom elevation | ||

Maximum bottom elevation | ||

Cell size | ||

Initial water surface elevation | ||

Prescribed total discharge at the upstream boundary | ||

For subcritical flow (T3a) | ||

For hydraulic jump (T3b) | ||

For pressurized flow (T3c) | ||

Prescribed water depth at the downstream boundary | ||

For subcritical flow (T3a) | ||

For hydraulic jump (T3b) | ||

For pressurized flow (T3c) | Free-outlet |

Symbol . | Meaning . | Value . |
---|---|---|

Pipe diameter | ||

Preissmann slot width | ||

Domain length | ||

Bump length | ||

Middle of bump | ||

Minimum bottom elevation | ||

Maximum bottom elevation | ||

Cell size | ||

Initial water surface elevation | ||

Prescribed total discharge at the upstream boundary | ||

For subcritical flow (T3a) | ||

For hydraulic jump (T3b) | ||

For pressurized flow (T3c) | ||

Prescribed water depth at the downstream boundary | ||

For subcritical flow (T3a) | ||

For hydraulic jump (T3b) | ||

For pressurized flow (T3c) | Free-outlet |

The simulation results of mixed pipe flow are presented in Figures 6 and 7 for pressurized flow. Numerical results show that the proposed model provides better estimations of the computed water surface elevation than the ‘Sanders & Bradford (2010)’ model, including upstream water surface elevation, shock position and transition between free-surface and pressurized flow. Regarding the water discharge, transition between free-surface and pressurized flow has caused additional oscillations for both models. For the ‘Sanders & Bradford (2010)’ model, the errors accumulated around the bump can still be observed for mixed or pressurized flow due to the adoption of HLL solver.

### Dambreak with bottom discontinuity

This test case aims at validating the convergence ability of the model on a transient dambreak problem over a discontinuous bottom with no friction. Two cases are tested with different bottom topography: step-up case (T4a) and step-down case (T4b). The domain length is with a mesh of uniform cell . The bottom step (discontinuity) is positioned at . For the step-up case (T4a), the bottom elevations at left and right sides of the discontinuity are and , respectively. For the step-down case (T4b), the bottom elevations at left and right sides of the discontinuity are and , respectively. The dam is located at , and water velocity is zero everywhere at . For the test case T4a, the water surface levels at the left and right sides of the dam are and , respectively. For the test case T4b, the water surface levels at the left and right sides of the dam are and , respectively. The dam-break flows instantaneously start at . The computational results with different mesh sizes are compared at when the wave fronts are not reached the domain boundaries. The simulation parameters of this test are summarized in Table 3. The pipe diameter is .

Symbol . | Meaning . | Value . |
---|---|---|

Pipe diameter | ||

Domain length | ||

Cell size | ||

Bottom elevation on the left-hand side of the dam | ||

Bottom elevation on the right-hand side of the dam | ||

For test case T4a | ||

For test case T4b | ||

Initial water depth on the left-hand side of the dam | ||

Initial water depth on the right-hand side of the dam | ||

Computation time |

Symbol . | Meaning . | Value . |
---|---|---|

Pipe diameter | ||

Domain length | ||

Cell size | ||

Bottom elevation on the left-hand side of the dam | ||

Bottom elevation on the right-hand side of the dam | ||

For test case T4a | ||

For test case T4b | ||

Initial water depth on the left-hand side of the dam | ||

Initial water depth on the right-hand side of the dam | ||

Computation time |

The simulation results at of the two test cases with bottom stepping up or down, respectively, are presented in Figure 8. It shows that, for both test cases, the solutions converge as mesh is refined. No numerical instability is produced at the position of the bed discontinuity during the dam-break process. Notice that for the step-up case T4a, a spike is observed at , this may be due to that the contact wave is ignored in the HLL Riemann solver adopted in the current study.

### Transient mixed pipe flow of dambreak

In this test case, the capability of the proposed method to reproduce transient flow over complex topography is tested by a mixed pipe flow experiment performed by Aureli *et al.* (2015), involving friction term. The experimental settings are illustrated in Figure 9. The red line represents the initial flow conditions with a pressure head of , and blue lines are the six gauges from G1 to G6 as illustrated in Table 4. The black line represents the pipe's bottom and top. The total length of the pipe is , and the first part of the pipe is with slope . The length of the second pipe is with slope . The inner diameter is , and the wall thickness is . At , the gate from the upstream inlet opens and then the water flows out. The left end of the pipe is partly closed to maintain the initial water height, and the right side is fully open to maintain a ventilated system.

Pressure Gauges . | . |
---|---|

G1 | 1.00 |

G2 | 3.00 |

G3 | 4.50 |

G4 | 6.80 |

G5 | 7.32 |

G6 | 8.52 |

Pressure Gauges . | . |
---|---|

G1 | 1.00 |

G2 | 3.00 |

G3 | 4.50 |

G4 | 6.80 |

G5 | 7.32 |

G6 | 8.52 |

This experiment mainly aims to study simple transient flow behavior in a single sloped Plexiglas pipe as the effect of the boundary conditions has been ignored. The numerical model proposed and the method developed by Sanders & Bradford (2010) is used to reproduce the experimental results. The cell size , and time step . Manning coefficient is set to , and wave celerity is chosen as with . Transmissive and closed boundaries are used in this work (Toro 2001). The numerical results are compared with the measurement data illustrated in Figure 10. The horizontal line in each graph is the pipe crown.

As presented in Figure 10, compared with the experimental data, results show that two methods tackling slope source term can both reproduce transient mixed pipe flow of dambreak, and the curve fitting between them is satisfied. Differences can also be identified that simulated pressure is higher than the experiment in gauges G4 and G5 and lower in G6. This may be due to the minimum water depth used and adopted for free-surface flow calculation. One can observe that both the proposed scheme and Sanders & Bradford (2010) show satisfying agreement with experimental measurements.

## CONCLUSIONS

In this paper, an augmented HLL Riemann solver is extended to circular cross-sections and mixed pipe flows over complex topography. The proposed model provides accurate predictions in reproducing steady flows with different hydraulic flow regimes. The summarized main features of this scheme are presented as follows:

An augmented HLL Riemann solver has been extended from rectangular shape to the currently studied pipe flow system by projecting the source term onto the Jacobian's eigenvector basis; the source term is discretized approximately as a succession of discontinuities separating regions of constant state.

The proposed model performs better than the reference solution obtained with the ‘Sanders & Bradford (2010)’ method when simulating steady flows over a bump, including sub-, trans-critical and hydraulic jump.

More accurate estimations than the reference are generated by the proposed model for mixed and pressurized pipe flow.

Transient dambreak pipe flow has been reproduced by the proposed model, including wetting or drying front and the calculation of the friction term.

The effectiveness and robustness of the proposed numerical scheme are verified by comparing the simulated results with numerical or measurement results. This method can also be applied to other cross-sections like trapezoid with the calculation of Riemann invariant over complex topography.

## ACKNOWLEDGEMENTS

This work is funded by the National Natural Science Foundation of China (Grant No. 51922096) and the Excellent Youth Natural Science Foundation of Zhejiang Province, China (LR19E080003). The authors warmly thank Pierre-André Garambois, INRAE, for his final lecture of the manuscript.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

Engineering Applications of Computational Hydraulics: Hommage to A. Preismann. Vol. II: Numerical Models in Environmental Fluid Mechanics. Pitman Advanced Publishing Program, Boston, MA and London