Finite volume simulation of unsteady water pipe flow

Abstract. The most commonly used hydraulic network models used in the drinking water community exclusively consider fully filled pipes. However, water flow numerical simulation in urban pipe systems may require to model transitions between surface flow and pressurized flow in steady and transient situations. The governing equations for both flow types are different and this must be taken into account in order to get a complete numerical model for solving dynamically transients. In this work, a numerical simulation tool is developed, capable of simulating pipe networks mainly unpressurized, with isolated points of pressurization. For this purpose, the mathematical model is reformulated by means of the Preissmann slot method. This technique provides a reasonable estimation of the water pressure in cases of pressurization. The numerical model is based on the first order Roe’s scheme, in the frame of finite volume methods. The novelty of the method is that it is adapted to abrupt transient situations, with subcritical and supercritical flows. The validation has been done by means of several cases with analytic solutions or empirical laboratory data. It has also been applied to some more complex and realistic cases, like junctions or pipe networks.


Introduction
The numerical simulation of water pipe networks is characterized by the difficulty of simulating both transient and steady states.Most of the time, the water flow is unpressurized but, under exceptional conditions, the limited storage capacity of the pipes could cause a momentary pressure peak if the water level raises quickly.Under these conditions, a pressurized flow should be taken into account, implying a change in the mathematical model in order to solve the problem accurately.
The most commonly used hydraulic network models used in the drinking water community exclusively consider fully filled pipes.However, a complete pipe system model should be able to solve steady and transient flows under pressurized and unpressurized situations and the transition between both flows (mixed flows).Many of the models developed to study the propagation of hydraulic transients solve both system of equations, e.g.Gray (1954), Wiggert and Sundquist (1977).For this purpose, numerical schemes derived from the Method of Characteristics (MOC) are employed.This method transforms the continuity and momentum partial differential equations into a set of ordinary differential equa-tions, much easier to solve.MOC schemes can handle complex boundary conditions.However, the main issue with these kind of methods is the need for interpolation in some cases, so that they are not strictly conservative.This results in a diffusion of the wave front that implies a wrong arrival time of the waves to the boundaries.Hence, a new efficient numerical scheme is justified in order to solve accurately the wave front propagation.
Two main points of view have been considered by the researchers to study mixed flow problems along the last decades.The first one consists of solving separately pressurized and shallow flows, treating the transitions between both flow types as internal boundary conditions.These models are capable of simulating sub-atmospheric pressures in the pipe, but the need for solving both mathematical models implies a high level of complexity.Other authors, e.g.García-Navarro et al. (1994), León et al. (2009), Kerger et al. (2010) developed simulations applying the shallow water equations in a slim slot over the pipe (Preissmann slot method).An estimation of the pipe pressure can be obtained from the water level in the slot.The great advantage of this model is the use of a single equation system for solving the complete prob-Published by Copernicus Publications on behalf of the Delft University of Technology.

Shallow water equations
The unsteady open channel water flow are usually modelled by means of the 1D shallow water or St. Venant equations.These equations represent mass and momentum conservation along the main direction of the flow and become a good description for the behaviour of the most of the pipe-flow-kind problems (Kundu et al., 2012;Abbott, 1979).The conservative form of this systems of equations is presented in Eqs. ( 1) and (2): where A is the wetted cross section, Q is the discharge, g is the acceleration due to gravity, I 1 represents the hydrostatic pressure force term and I 2 accounts for the pressure forces due to channel/pipe width changes: in which η and h(x, t) represent vertical coordinate and the free surface depth, respectively.The remaining terms, S 0 and S f , represent the bed slope and the energy grade line (defined in terms of the Manning roughness coefficient), respectively: where R = A/P , P being the wetted perimeter.The coordinate system used for the formulation of the shallow water equations is shown in Fig. 1.A vectorial form of the Eqs.( 1) and ( 2) is widely used: (2) e, g is ostatic forces (3) (4) (5) in which η and h(x, t) represent vertical coordinate and the free surface depth, respectively.The remaining terms, S 0 and S f , represent the bed slope and the energy grade line (defined in terms of the Manning roughness coefficient), respectively: where R = A/P, P being the wetted perimeter.The coordinate system used for the formulation of the shallow water equations is shown in Fig. 1.A vectorial form of the equations ( 1) and ( 2) is widely used: where In those cases in which F = F(U), when I 2 = 0, it is possible to rewrite the conservative system by means of the Jacobian matrix of the system 11: where c is the wave speed (analogous to the speed of sound in gases), defined as follows: The set of real eigenvalues (λ 1,2 ) and eigenvectors (e 1,2 ) of the system matrix can be used for diagonalizing it.These meagnitudes represent the speed of propagation of the information: It is very common to characterize the flow type by means of the Froude number Fr = u/c (analogous to the Mach number in gases).It allows the classification of the flux into three main regimes: subcritical Fr < 1, supercritical Fr > 1 and critical Fr = 1.
www.drink-water-eng-sci.net where In those cases in which F = F(U), when I 2 = 0, it is possible to rewrite the conservative system by means of the Jacobian matrix of the system Eq.( 11): where c is the wave speed (analogous to the speed of sound in gases), defined as follows: The set of real eigenvalues (λ 1,2 ) and eigenvectors (e 1,2 ) of the system matrix can be used for diagonalizing it.These magnitudes represent the speed of propagation of the information: It is very common to characterize the flow type by means of the Froude number F r = u/c (analogous to the Mach number in gases).It allows the classification of the flux into three main regimes: subcritical F r < 1, supercritical F r > 1 and critical F r = 1.

Water-hammer equations
The instant response for changes in a pipe flow is closely related to the elastic compresibility of both the fluid and the pipe wall material.Unsteady flow in pipes is commonly described by the cross-section integrated mass and momentum equations (Chaudhry and Mays, 1994):

Preissmann slot model
The Preissmann slot approach consist in assuming that the top of the pipe/closed channel is connected to a fictitious narrow slot, which is open to the atmosphere.This fact allows to consider a free surface flow situation, even when the pipe is pressurized, by including the slot in the shallow water equations (see Fig. 2).The slot width b s is ideally chosen equaling the speed of gravity waves in the slot to the water hammer wavespeed.Therefore, the water level in the slot corresponds to the pressure head level.This model is based on the previously remarked similarity between the wave equations which describe free surface and pressurized flows.The water hammer flow comes from the capacity of the pipe system to change the area and fluid density, so forcing the equivalence between both models requires that the slot stores as much fluid as the pipe would by means of a change in area and fluid density.The pressure term and the wavespeed for A ≤ bH are: and for A > bH: The ideal choice for the slot width results in: in which A f stands for the full pipe cross-section.boundary conditions (BC).The number of boundary conditions depends on the flow regime (subcritical or supercritical).Hence, there are four posibilities for a 1D numerical problem (see Table 2).
In the cases that consider pipe junctions, additional internal boundary conditions are necessary in order to model this feature.Fig. 3 shows an example of a 3 pipe junction.The water level equality condition is imposed at the junction for all the pipes: Discharge continuity condition is formulated depending on the flow regime.In the cases presented, water flows in a subcritical regime: Storage wells are another common junction type.For these cases, the boundary conditions are modified as follows: where A well and H well are the well top area and depth, respectively.
The presented method use an explicit time discretization and it requires a control on the time step in order to avoid being e the pipe thickness and E, K the elastic modulus of the pipe material and fluid, respectively.By neglecting convective terms, it is possible to reach a linear hyperbolic equation system: These equations conform a hyperbolic system, analogous to the shallow water equations, considering the pressure and the velocity as the conserved variables.

Wiggert test case
The experimental case designed by Wiggert (Wiggert (1972)) and widely numerically reproduced (e.g.Kerger et al. (2010), Bourdarias et al. (2007)) consists in a horizontal 30 m long and 0.51 m wide flume.A 10 m roof is placed in the middle section, setting up a closed rectangular pipe 0.148 m in height (see Fig. 6).A 0.01m −1/3 s Manning roughness coefficient is assumed.As initial conditions, a water level of 0.128 m and zero discharge are considered.Then a wave coming from the left causes the pressurization of the pipe.The imposed downstream boundary conditions are the same values measured by Wiggert (Fig. 7).Fig. 8 shows the numerical results for the four gauges.The experimental comparison is done only for the second one, due to data availability issues.An overall good agreement is observed with only small numerical oscillations presented.

Transient mixed flow
The aim of the next test case is to prove the model in the simulation of large-scale strong transients.It considers a uniform slope (0.1%) rectangular pipe connected to a downstream valve (León (2007) and León et al. (2009)).The lenght, width and height are 10 km, 10 m and 9.5 m, respectively, and the chosen Manning roughness is 0.015.
As initial condition, a uniform water depth (8.57m) and discharge (Q = 240m 3 /s) are assumed.From this state, the downstream valve is closed, generating a shock wave moving upstream and pressurizing the pipe.Fig. 9 shows the numerical results obtained for the pressure head at three different times.For this simulation, a 500 cells mesh and a timestep given by CFL=0.5 have beed used.It is clearly observed that the pipe pressure raises gradually with the upwards shockwave movement.

Pressurized flow transients
A full-pressurized (León (2007) and León et al. (2009)) 10 km flat, frictionless pipe of rectangular section (10 m x 7.853 m) is connected upstream to a water reservoir that guarantees a constant pressure head of 200 m.A closure valve is placed downstream.When the valve is suddenly closed, a waterhammer wave is generated moving upwards.Fig. 9 shows the numerical results with CFL=0.8 and a 500 cells mesh.
This case simulates a completely pressurized pipe, so it is necessary to be very precise in the slot width selection in order to ensure the numerical solution accuracy.A waterhammer speed of 1000 m/s and a initial flow speed of 2.0 m/s are assumed.Following (Eq.( 25)) condition: 6 Pipe networks

Transient flow in a 3-pipe junction
A case proposed in Wixcey (1990) is presented first.A 5 km long prismatic main channel (1) dividing into two secondary branches (2 and 3) also 5 km long with the same 1 m wide rectangular geometry have been considered, as presented in fig. 3. Manning roughness coefficient for all the pipes is 0.01m −1/3 s.The main branch has a slope of 0.002 and the secondary ones 0.001.

Preissmann slot model
The Preissmann slot approach consist in assuming that the top of the pipe/closed channel is connected to a fictitious narrow slot, which is open to the atmosphere.This fact allows to consider a free surface flow situation, even when the pipe is pressurized, by including the slot in the shallow water equations (see Fig. 2).The slot width b s is ideally chosen equaling the speed of gravity waves in the slot to the water hammer wavespeed.Therefore, the water level in the slot corresponds to the pressure head level.This model is based on the previously remarked similarity between the wave equations which describe free surface and pressurized flows.The water hammer flow comes from the capacity of the pipe system to change the area and fluid density, so forcing the equivalence between both models requires that the slot stores as much fluid as the pipe would by means of a change in area and fluid density.The pressure term and the wavespeed for A ≤ bH are: www.drink-water-eng-sci.net Drink.Wate and for A > bH : The ideal choice for the slot width results in: in which A f stands for the full pipe cross-section.It is necessary to build an approximate Jacobian matrix J whose eigenvalues and eigenvectors satisfy: The wave strenght coefficients αk represent the variable variation coordinates in the Jacobian matrix basis.The eigenvalues and eigenvectors are expressed in terms of the average flux velocity and wave speed: where Following (García-Navarro and Vázquez-Cendón ( 2000)), the source terms of the equation system (Eq.7) are also discretizated using an upwind scheme: in which the source streghts βk take the role of αk coefficients for the fluxes.Then, the complete discretization of the system becomes:

Boundary conditions and stability conditions
In order to solve a numerical problem, it is necessary to characterize the domain limits by imposing some physical boundary conditions (BC).The number of boundary conditions depends on the flow regime (subcritical or supercritical).Hence, there are four posibilities for a 1-D numerical problem (see   Drink.Water Eng.Sci.www.drink-water-eng-sci.net In the cases that consider pipe junctions, additional internal boundary conditions are necessary in order to model this feature.Figure 3 shows an example of a 3 pipe junction.
The water level equality condition is imposed at the junction for all the pipes: Discharge continuity condition is formulated depending on the flow regime.In the cases presented, water flows in a subcritical regime: Storage wells are another common junction type.For these cases, the boundary conditions are modified as follows: where A well and H well are the well top area and depth, respectively.
The presented method uses an explicit time discretization and it requires a control on the time step in order to avoid instabilities.Hence the Courant-Friedrichs-Lewy (CFL) condition is applied on every branch (j ) of the network: Then, the minimum of all the computed time steps is chosen:

Test cases
The model is applied to several analytic or experimental test cases in order to evaluate its accuracy.

Steady state over a bump
Following Murillo and García-Navarro (2012) a frictionless rectangular 25 m × 1 m prismatic channel is considered.The variable bed level is given by: and the initial conditions for the water depth and discharge: The different boundary conditions for the test cases simulated are summarized in Table 3.
Figure 4a shows the generation of hydraulic jump that connects both subcritical and supercritical regimes.In the second test case, shown in Fig. 4b the connection is made without any shock wave and it can be mathematically proved (and numerically checked) that this transition takes place in the highest part of the bump.ckwave propagation in transient mixed flow at t=100 s (red), t=200 s (green) and t=300 s (blue).ckwave propagation in transient mixed flow at t=100 s (red), t=200 s (green) and t=300 s (blue).Pressure head (m) x (m) aterhammer simulation.Pressured head at t=3 s (red), t=6 s (green) and t=9 s (blue).

Dam-break
Dam-break is a classical example of non-linear flow with shocks.It has been widely used to test the accuracy and conservation of the numerical scheme, by comparing with its analytical solution.In the case presented, an initial discontinuity of 1 m : 0.5 m ratio with no friction was considered.
Figure 5 shows the good agreement between numerical and analytical solution at the given time.

Wiggert test case
The experimental case designed by Wiggert (1972) and widely numerically reproduced (e.g.Kerger et al., 2010;Bourdarias and Gerbi, 2007) consists in a horizontal 30 m long and 0.51 m wide flume.A 10 m roof is placed in the middle section, setting up a closed rectangular pipe 0.148 m in height (see Fig. 6).A 0.01 s m −1/3 Manning roughness coefficient is assumed.As initial conditions, a water level of 0.128 m and zero discharge are considered.Then a wave coming from the left causes the pressurization of the pipe.The imposed downstream boundary conditions are the same values measured by Wiggert (1972) (Fig. 7). Figure 8 shows the numerical results for the four gauges.The experimental comparison is done only for the second one, due to data availability issues.An overall good agreement is observed with only small numerical oscillations presented.Drink.Water Eng.Sci.www.dri

Transient mixed flow
The aim of the next test case is to prove the model in the simulation of large-scale strong transients.It considers a uniform slope (0.1 %) rectangular pipe connected to a downstream valve (León, 2006;León et al., 2009).The lenght, width and height are 10 km, 10 m and 9.5 m, respectively, and the chosen Manning roughness is 0.015.As initial condition, a uniform water depth (8.57m) and discharge (Q = 240 m 3 s −1 ) are assumed.From this state, the downstream valve is closed, generating a shock wave moving upstream and pressurizing the pipe.Figure 9 shows the numerical results obtained for the pressure head at three different times.For this simulation, a 500 cells mesh and a timestep given by CFL = 0.5 have beed used.It is clearly observed that the pipe pressure raises gradually with the upwards shockwave movement.

Pressurized flow transients
A full-pressurized (León, 2006;León et al., 2009) 10 km flat, frictionless pipe of rectangular section (10 m × 7.853 m) is connected upstream to a water reservoir that guarantees a constant pressure head of 200 m.A closure valve is placed downstream.When the valve is suddenly closed, a waterhammer wave is generated moving upwards. Figure 9 shows the numerical results with CFL = 0.8 and a 500 cells mesh.
This case simulates a completely pressurized pipe, so it is necessary to be very precise in the slot width selection in order to ensure the numerical solution accuracy.A waterhammer speed of 1000 m s −1 and a initial flow speed of 2.0 m s −1 are assumed.Following (Eq.25) condition: 6 Pipe networks

Transient flow in a 3-pipe junction
A case proposed in Wixcey (1990)     Fig. 3. Manning roughness coefficient for all the pipes is 0.01 s m −1/3 .The main branch has a slope of 0.002 and the secondary ones 0.001.First, a steady state has been calculated from the initial conditions: using as boundary condition the inlet discharge: and free outflow boundary condition at the outlet of pipes 2 and 3.
Using this steady state as initial condition for the transient calculation, a triangular function of peak discharge Q MAX = 3.2 m 3 s −1 and a period of 600 s (Fig. 11) is imposed at the beginning of the pipe 1. Figure 12 shows the numerical results using a 100 cells mesh and a CFL = 0.9 condition.
Junctions J 1 and J 2 are treated as normal confluences (Eq.34) and a storage well of 5 m 2 top surface is assumed in junctions W 1 and W 2 (Eq.36).The previously computed steady state is now used as an initial condition for a second calculation.A triangular function of peak discharge Q MAX and a period of 600 s (Fig. 11) is imposed at the beginning of the pipe 1.The minimum discharge is 0.1 m 3 s −1 .
Two different situations were studied.In the first one the peak discharge is fixed to Q MAX = 2.0 m 3 s −1 and the flow remains unpressurized all over the pipe system.In the second case, the maximum discharge is increased to Q MAX = 3.0 m 3 s −1 , so pressurization occurs in some points of pipe 1. Figure 14 shows the results for the pressurized case.The results in pipes 3 and 6 have been omitted because of symmetry reasons.The discharge at the center of the pipe 4 is constantly zero and equal in magnitude but opposite in sign in the symmetric points, which verify the simmetric character of the pipe network.

Conclusions
The present work leads to the conclusion of the reasonably good applicability of the Preissmann slot model for an estimation of the pressure values in unsteady situations between both shallow and pressurized flows.In this second case, a small deviation of the ideal slot width induces a mass and momentum conservation error.On the other hand, it is clear that wider slots increases the stability of the numerical model.Hence, in benefit of the results accuracy, it is remarkably important to find an agreement between both facts.
The model takes advantages of the similarity of both equation systems (shallow water and pressurized flow) and provides a simple way to simulate occasionally pressurized pipes, treating the system as an open channel in the Preissmann slot.This results in an easier implementation of the model because it avoids the managing of two separate systems of equations (shallow water and water hammer) in order to model separately the pressurized and the free surface flows.
The use of an explicit scheme implies a limitation on the computational time step, which increases the calculation time, specially in the in the pressurized cases, due to the small width of the slot.
Finally, the authors want to remark the possibility of adapting the method to more realistic systems, like pipe networks which can be punctually pressurized.Some additional internal boundary conditions are necessary in these cases, in order to represent pipe junctions and storage wells.

1 .
Coordinate system for shallow water equations.

Figure 1 .
Figure 1.Coordinate system for shallow water equations.

Figure 8 .Figure 9 .
Figure 8. Comparision between numerical results (blue), experimental data (green dots) and numerical results from Kerger et al. (2010) (red dots) for the four gauges of the Wiggert test case setup.

Figure 8 .
Figure 8. Comparision between numerical results (blue), experimental data (green dots) and numerical results from Kerger et al. (2010) (red dots) for the four gauges of the Wiggert test case setup.
results (blue), experimental data (green dots) and numerical results fromKerger et al. (2010) (red our gauges of the Wiggert test case setup.
results (blue), experimental data (green dots) and numerical results fromKerger et al. (2010) (red our gauges of the Wiggert test case setup.

Figure 13 .
Figure 13.(Scheme of the looped pipe network.

Figure 14 .
Figure 14.Time histories at grid points i=N/2 (red) and i=N (blue) for punctually pressurized flow: (a) Water depth; (b) Discharge.The discharge for pipe 4 is recorded at locations i=1 (red), i=N/2 (grey) and i=N (blue).

Figure 14 .
Figure 14.Time histories at grid points i = N/2 (red) and i = N (blue) for punctually pressurized flow: (a) Water depth; (b) Discharge.The discharge for pipe 4 is recorded at locations i = 1 (red), i = N/2 (grey) and i = N (blue).

Table 1 .
List of symbols and units.

Table 3 .
Steady state over a bump.