A Mixed Flow Analysis of Sewer Pipes with Different Shapes Using a Non-Oscillatory Two-Component Pressure Approach (TPA)

Abstract

This paper aimed to justify the performance of a non-oscillatory TPA-based model proposed by the authors for capturing transient mix flow in sewer systems consisting of a variety of pipe shapes. The model utilizes a first-order Godunov Finite volume numerical scheme in which a Harten–Lax–van Leer (HLL) Riemann solver was used for calculating the fluxes at the cells’ boundaries. The spurious numerical solution associated with the transient mix flow analysis is suppressed by enhancing the numerical viscosity of the scheme when the pipe pressurization is imminent. Due to the lack of experimental data for systems with pipe shapes other than circular and rectangular, a hypothetical pipe system for which analytical solutions exist was employed to verify the model performance. The results reveal that for all pipe shapes considered, the model provides oscillation-free solutions even at a high acoustic speed of 1400 m/s. It is also observed that the numerical results are in perfect agreement with the analytical solution. The obtained results conclude that the proposed model can be utilized to capture transient responses of sewer systems with any pipe shape.

1. Introduction

Sewer pipe systems rarely run under steady-state flow conditions as the inflows into such systems change with time. In dry weather conditions, the flow gradually changes with time inside the conduits, and quasi-steady open-channel flow is established across the system. However, when a severe storm occurs, complex transient flows may be onset and the system undergoes flow regime transition from open channel to pressurized flow. The induced pressurization front serves as a moving piston and makes the air expelled out of the system through different components including, drop shafts, manholes, and outfalls. Fortunately, during pressurization, the available storage in the partially filled conduit can accommodate the transient flow energy and does not allow the energy to be stored in the pipe and liquid as strain and compression energy [1]. In such conditions, the elastic feature of the flow does not play an important role even if the pace of the transient is high, and the inertia and mass oscillation mainly govern the transient flow.

However, there are some situations in which the elastic feature of the flow becomes important. Two pressurization hydraulic bores moving in the opposite direction can easily trap a large air pocket in the conduit. As the entrapped air becomes pressurized, it may absorb a good portion of the transient energy and can significantly affect the hydraulics of the system. If the energized air pocket finds its way to the drop shafts or manholes, it can violently leave the system and produce geysering [2,3]. However, the presence of an air vent centered at the right position causes the air pocket to leave the system and remove the risk of geysering. On the other hand, when the last air escapes from the system, the two adjacent hydraulic bores collide, and depending on the velocities of the adjacent water columns, significant water hammer pressures can occur [4]. The generated water hammer pressure spike will propagate and affect the rest of the pipe located between two adjacent drop shafts reflecting the pressure spike. The reflected wave may produce intense negative pressures, followed by column separation in the system.

Available off-the-shelf software such as SWMM [5], InfoWorks [6], etc. can calculate the most dominant transient flow governed by the flow inertia and mass oscillation in the system, but they are crippled in capturing those transient flow phenomena governed by the elastic feature of the flow. To fill this gap, though partially, some in-house computer programs have been emerged for performing surge analysis in sewer pipe systems [7]. These programs use either shock-fitting or shock-capturing numerical strategies to solve the one-dimensional momentum and continuity equations governing the transient flow in sewer systems.

In the shock-fitting strategy, the pressurization front(s) that separate open channel, and pressurized flow is traced with time through applying the discrete form of continuity, momentum, and/or energy on either side of the pressurization front [8,9,10]. Having calculated the location of the front, the open channel and pressurized flow on either side of the interface are calculated using their own set of governing equations. Song et al. [11] employed the shock-fitting method along with the method of characteristic for analyzing transient mixed flow in sewer systems. Although the method of characteristics can replicate the pressurized flow quite accurately, it fails to reasonably capture hydraulic bores in the open-channel flow sections. Abbott and Minns [12] showed that the method of characteristics cannot conserve mass when applied along a hydraulic bore. Leon et al. [9] address this issue by solving the conservative form of the governing equations using a finite volume numerical method. Although the proposed method succeeded in replicating both open channel and pressurized flow accurately, it suffers from the common disadvantage of shock-fitting methods, which is the difficultly in handling several interfaces in the system, as well as in treating the interaction of the interfaces with each other and with the boundaries of the system. This has made the shock-capturing method the center of interest among the researchers.

In the shock-capturing method, both open channel and pressurized flow are treated using the same set of equations utilized for calculating transient open-channel flows. The most popular method of this type is well-known as the Preissmann slot method (PSM) [13,14,15], named after its inventor engineer Preissmann [16]. By a virtual narrow slot on the crown of the pipe, the flow is assumed to always remain in open-channel flow condition. The width of the slot is calculated in a way that open channel waves move as fast as their counterpart elastic waves in pressurized flows; the height of the water in the slot represents the pressure head of the flow. However, this approach cannot replicate negative pressures, and as soon as negative pressures tend to occur, the flow switches to open-channel flow. Kerger et al. [13] resolved this problem by improving the PSM using a virtual negative slot and showed that their method can accurately replicate negative pressures. By splitting the pressure term in the momentum equation, Vasconcelos et al. [17] proposed a novel approach, the two-component pressure approach (TPA), which can capture negative pressures during transient mixed-flow analysis.

Both PSM and TPA, however, generate spurious numerical oscillation when the flow regime changes from open channel to pressurized flow. The numerical oscillation intensifies as the acoustic wave speed considered in the calculation increases, and therefore, beyond a certain level, it makes the solution unstable and worthless. One approach to control the numerical oscillation is to artificially reduce the wave speed, but this approach is acceptable as long as the physics of the problem is not distorted. In a flow condition that is governed by inertia and mass oscillation, transient flow is independent of the wave speed, and this approach works well, though due to a considered wider slot, some artificial storage is added to the system, resulting in extra non-physical pressure attenuation [18]. However, when the elastic feature of the transient flow becomes of significant importance, reducing the wave speed compromises the results by distorting the magnitude and distribution of the resulting pressures.

Vasconcelos et al. [19] proposed numerical filtering and hybrid flux approaches to suppress the spurious numerical oscillations and showed that these approaches can reasonably control the numerical oscillations if the acoustic wave is below a certain level (100 m/s). Malekpour and Karney [14] proposed an approximate HLL solver that can remove the numerical oscillation in the PSM even when the acoustic speed exceeds 1000 m/s; the validity of his method has been independently confirmed by others [20].

Objective and Organization of the Paper

Recently, Khani et al. [21] employed the approximate HLL solver proposed by Malekpour and Karney [14] in conjunction with TPA and showed that the solver can effectively remove the numerical oscillations for as high acoustic wave speed as 1000 m/s, though the validity of the model was tested for pipes with circular sections and rectangular cross sections. However, in real sewer pipe systems, a variety of pipe shapes can be used, and therefore, further investigation is required to find out if this model can be utilized in practice. To fill this gap, this paper aimed to verify whether the performance of the model is retained when applied to other pipe cross sections.

The paper’s organization is as follows: Theoretical background including governing equations, TPA, and the numerical solution is first described, and the hypothetical example and its associated analytical solution utilized for validation of the produced numerical results are then discussed. In the next stage, the numerical results are compared with the analytical solutions and finally, the conclusions are made.

2. Theoretical Backgrounds

2.1. Governing Equations

Unsteady flows in open channels are governed by the following partial differential equations representing the conservative form of the continuity and momentum equations [22]:

$$\frac{\partial \mathit{U}}{\partial t}+\frac{\partial \mathit{F}}{\partial x}=\mathit{S}$$

(1)

where the vectors $\mathit{U}$, F, and S are the flow variables, fluxes, and source terms, respectively, and shown in the following:

$$\mathit{U}=\left[\begin{array}{c}A\\ Q\end{array}\right];\mathit{F}=\left[\begin{array}{c}Q\\ \frac{Q^2}{A}+Agh\end{array}\right]; \mathit{S}=\left[\begin{array}{c}0\\ \left(S_0-S_f\right)gA\end{array}\right]; S_f=\frac{Q^2{n}^2}{A^2{R}^{\frac{4}{3}}}$$

where $A$ = flow cross-sectional area, $Q$ = flow rate, $h$ = distance between the free surface and the centroid of the flow cross-sectional area, $S_0$ = pipe slope, $S_f$ = energy grade line slope, $R$ = flow hydraulic radius, $n$ = Manning coefficient, and $g$ = gravitational acceleration.

2.2. TPA

By a smart trick, the above equation can be also used for calculating pressurized flows. In PSM, a narrow slot above the crown of the pipe allows the flow to remain in an open-channel flow regime even when the pipe is surcharged. The width of the slot is selected narrow enough to make sure that, while the flow is in the slot the open channel wave velocity replicates the acoustic wave speed of the pipe. The following equation relates the slot width and pipe acoustic speed:

$$T_S=\frac{g{A}_f}{a^2}$$

(2)

where $T_S$, $A_f$, and $\mathrm{a}$ are the slot width, the cross-sectional area of the pipe, and the acoustic speed of the pipe, respectively.

In TPA [21], in order to use Equation (1) for pressurized flow, the flow depth in the momentum flux split into two terms, as shown in the following:

$$\mathit{F}=\left[\begin{array}{c}Q\\ \frac{Q^2}{A}+Ag\left(h_c+h_s\right)\end{array}\right]$$

(3)

In this equation, $h_c$ represents the distance between the water surface and the centroid of the flow cross-sectional area, and $h_s$ measures the surcharging head. When the pipe carries open-channel flows, the surcharging head is set to zero, and the flux becomes identical to its original form presented in Equation (1). However, when the pipe is pressurized, $h_c$ remains constant at the maximum height of the pipe, and $h_s$ measures the pressure head of the system. Since the pipe can expand and contract when exposed to positive and negative surcharging pressure heads, the flow cross-sectional area is a function of the pressure head and acoustic speed of the pipe and can be calculated using the following equation:

$$A=A_f\left(1+\frac{g{h}_s}{a^2}\right)$$

(4)

where $a$ = pipe acoustic speed.

Note that when the pipe tends to be depressurized at a specific location with adequate ventilation, the flow is switched back to open-channel flow, and $h_s$ is reset to 0. However, lack of ventilation causes the pipe to receive negative pressures, and in such conditions, $h_s$ measures the magnitude of the negative pressure head, and $h_c$ remains unchanged and represents the maximum height of the pipe.

3. Numerical Solution

In this work, the Godunov scheme was utilized to numerically solve the governing equations. The Godunov scheme is an explicit finite volume approach that is widely used in solving hyperbolic partial differential equations. The computational domain in this approach is discretized into some computation cells with a spatial length of Δx. The order of accuracy in the Godunov scheme depends on how the data are reconstructed in the computational cells. Piecewise constant data reconstruction provides the first order of accuracy, while a piecewise linear data reconstruction leads to a second-order accurate solution. Since first-order numerical schemes are generally more dissipative than second-order schemes, they better fit into the existing application because, as shown by Malekpour and Karney [14], lack of adequate numerical viscosity of the scheme during the pressurization of the conduit is the main cause of the production of the spurious numerical oscillations.

By discretizing Equation (1), unknowns at the current time level can explicitly be calculated based on the data retrieved from the previous timeline using the following equation:

$${\mathit{U}}_i^{n+1}={\mathit{U}}_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\left({\mathit{F}}_{i+\frac{1}{2}}^n-{\mathit{F}}_{i-\frac{1}{2}}^n\right)+\mathrm{\Delta }t\times {\mathit{S}}_i^n$$

(5)

where subscript $i$ is the computational cell number, $i+\frac{1}{2}$ and $i-\frac{1}{2}$ refer to the upstream and downstream boundaries of $i_{th}$ cell, respectively, $n$ and $n+1$ refer to the previous and current timelines, respectively, and $\mathrm{\Delta }t$ is the computational time step.

In the Godunov scheme, the fluxes at the cell boundaries are calculated by solving the Riemann problem, which is defined as a hyperbolic system of equations with a flow discontinuity. Although the exact Riemann solution is available for the current system of equations, considering the iterative nature of the solution, it compromises the efficiency of the numerical analysis. To resolve this problem, approximate Riemann solutions can also be utilized to speed up the calculations [23]. Several approximate Riemann solutions have been proposed, including Roe, HLL, and Harten–Lax–van Leer–contact (HLLC) [24], among which the HLL Riemann solution is one of the most efficient ones, and thus, we utilized it in this study.

The HLL Riemann solution assumes that the generated waves on either side of the discontinuity are both of shock wave type. Figure 1 describes the wave structure in the HLL solver.

The speed of the left and right waves can be easily approximated by the following equations [25]:

$$S_L=V_L-{\Omega }_L$$

$$S_R=V_R+{\Omega }_R$$

where indexes L, R refer to the left and right waves, respectively, and ${\Omega }_L$ and ${\Omega }_R$ can be calculating by the following equation:

$${\Omega }_{K\left(K=L, R\right)}=\{\begin{array}{c}\sqrt{\frac{g\left[Y_G{A}_G-{\left(h_s+h_c\right)}_K{A}_K\right]A_G}{A_K\left(A_G-A_K\right)}} if A_G>A_K \\ {\mathit{C}}_{\mathit{K}} if A_G>A_K\end{array}$$

(6)

${\mathit{C}}_{\mathit{K}}$ in Equation (6) is the gravity wave speed and calculated by the following equation and the index G, as discussed later, refers to a reference condition through which one can control the wave speeds:

$${\mathit{C}}_{\mathit{K}}=\sqrt{\mathit{g}{\mathit{D}}_{\mathit{K}}}$$

where $\mathit{D}$ is the hydraulic depth.

Having the $S_L$ and $S_R$ calculated, the fluxes can be approximated as follows:

$${\mathit{F}}_{\ast }=\{\begin{array}{cc}{\mathit{F}}_L& if S_L>0\\ \frac{S_R{\mathit{F}}_L-S_L{\mathit{F}}_R+S_L{S}_R\left({\mathit{U}}_R-{\mathit{U}}_L\right)}{S_R-S_L}& if S_L\le 0 and S_R\ge 0\\ {\mathit{F}}_R& if S_R<0\end{array}$$

(7)

As can be seen, when $S_L$ is positive, the flow regime is of supercritical type, and the flux at the star zone is equal to ${\mathit{F}}_L$. If the right wave moves to the left, the flow is supercritical, but it moves in the opposite direction. In such a condition, the upstream condition prevails, and the flux at the star zone becomes equal to ${\mathit{F}}_R$. In other conditions, the flow is subcritical, and the flux at the star zone is affected by both upstream and downstream flow conditions.

Although the fluxes presented in Equation (7) provide stable results over a wide range of flow conditions, it produces strong spurious numerical oscillations when the flow is being switched from an open channel to a pressurized flow. The induced oscillation is due to the significant change in the magnitude of the wave speed during the pressurization [14]. It is also found that during the pressurization, the numerical scheme fails to admit adequate numerical viscosity to suppress the numerical oscillation. Obviously, to control the numerical oscillation, artificial viscosity has to be added to the scheme when the pressurization occurs.

In the Godunov type finite volume scheme, the amount of artificial viscosity of the scheme can be controlled by changing the wave speeds based on which the fluxes are calculated. $Y_G$ in Equation (6) is to control the magnitude of the wave speed whenever required. Since the numerical oscillation is set up during the pressurization, it seems that the best place to increase the artificial numerical viscosity is the computational cell containing the pressurization front. However, increasing artificial viscosity only at the cell containing pressurization front gives rise to spurious numerical oscillation in the neighbor cells. To resolve this issue, it was found that the artificial viscosity should be distributed around the computational cell with a pressurization front.

If $Y_G$ is greater than the height of the conduit, the wave velocity calculated by Equation (6) does not differ from the gravity wave velocity except in the vicinity of the conduit roof. In other words, by using Equation (6), significant artificial velocity is admitted to the numerical scheme only when the water surface in the conduit becomes very close to the conduit roof, and the conduit is about to become pressurized. Extensive numerical experiments by Malekpour and Karney [14] resulted in the following formula for calculating $Y_G$:

$$Y_G=K_a\times MAX \left[d_{i-NS},d_{i-NS+1},\dots ,d_i,d_{i+1},\dots d_{NS}\right]$$

(8)

where $d=h_c+h_s$.

By using the above equation, the maximum $\mathrm{d}$ is calculated within a number of cells ($NS$) located on either side of the $i_{th}$ cell for which the wave velocity is being calculated. $Y_G$ is then calculated by multiplying the maximum $\mathrm{d}$ by the factor $K_a$. In this way, numerical viscosity is distributed within a number of cells rather than being injected just in one cell. If all$d$’s are greater than the conduit height, the system is pressurized, and a $K_a$ = 1.001 would provide reasonable results. When there is a pressurization front located within the cells, the numerical viscosity is further intensified by applying $K_a$ = 1.4. The number of cells ($NS$) considered depends on the resolution of the computational grid and should be selected such that the numerical viscosity is adequately distributed on either side of the computational cell for which the wave velocity is being calculated. Malekpour and Karney [14] suggested that the number of cells should cover a distance equal to at least three times as large as the conduit height, but in any case, it should not be less than three cells. If the $i_{th}$ computational cell is found near a boundary, we should also incorporate the flow depth at the boundary point into Equation (8).

3.1. Conduit Geometry

As mentioned earlier, the aim of this paper is to justify if the non-oscillatory TPA model proposed by Khani et al. [21] for calculating transient mix flow in circular and rectangular pipe cross sections performs equally well when other pipe shapes are used in the calculations. The most popular pipe shapes are assumed to be those supported by the program SWMM and presented in Figure 2. However, for the sake of generality, the custom shape whose dimensions can be defined by the users is also tested.

For the custom shape, the width of the conduit, W, at different heights needs to be defined; an example of such data is presented in

Table 1. Non-dimensional geometry information for the custom shapes considered in the study.

$\frac{{Y}}{{{Y}}_{\mathbf{f}\mathbf{u}\mathbf{l}\mathbf{l}}}$	$\frac{\mathbf{W}}{{\mathbf{Y}}_{\mathbf{f}\mathbf{u}\mathbf{l}\mathbf{l}}}$	$\frac{\mathbf{Y}}{{\mathbf{Y}}_{\mathbf{f}\mathbf{u}\mathbf{l}\mathbf{l}}}$	$\frac{\mathbf{W}}{{\mathbf{Y}}_{\mathbf{f}\mathbf{u}\mathbf{l}\mathbf{l}}}$
0.00	0.000	0.56	0.928
0.08	0.667	0.64	0.874
0.16	0.930	0.72	0.798
0.24	1.000	0.80	0.697
0.32	0.997	0.88	0.567
0.40	0.988	0.96	0.342
0.48	0.967	1.00	0.000

. Note that the values of the table are non-dimensional with respect to the height of the conduit, $Y_{full}$.

Shapes 1, 2, and 3 are all defined by the maximum height and maximum width of the pipe, while shapes 4 to 10 can be defined by just the maximum height of the pipe. Shape 12 is characterized by three parameters including maximum height, top width, and triangle height. Similarly, maximum height and top width define Shape 13 but instead of the triangle height, the bottom radius is used. Shape 14 is defined by maximum height, bottom width, and top radius.

The numerical solution uses the information related to the pipe geometry including flow area, top width, and hydraulic radius as a function of depth and flow depth as a function of flow area. For shapes 12, 13, and 14, all these parameters can be calculated analytically but for the other shapes, they should be interpolated from some precalculated tables. These tables are formed by calculating the geometry parameters at the incremental depths from the bottom to the top of the conduit. For the sake of generality, these data are stored in non-dimensional form with respect to the maximum height of the conduit. In this research, for the shapes 1 to 10, the data used by SWMM [26] were utilized. For the custom shape, the data defined by the user (an example of which is presented in Table 1) were employed to calculate the geometry-related parameters.

3.2. Model Verification

The validity of the proposed model, as well as its performance with regard to the suppression of the spurious numerical oscillations, were justified once before [21] through comparing the model results with the data obtained from a few experimental studies and analytical solutions but the investigation remained limited to the pipe systems with circular and rectangular shapes. Unfortunately, to the best of the authors’ knowledge, there is no experimental study on transient mixed flow in pipes with pipe shapes other than rectangular circular. Thus, in this research, an analytical solution was employed as a benchmark for verifying the model performance.

The analytical solution consists of a horizontal-frictionless conduit with the maximum height and the total length of 1 and 500 m, respectively. A reservoir at the upstream end of the pipe supplies the pipe and a reservoir at the downstream end collects the flow. The same water depth = 0.5 m at the upstream and downstream reservoirs made the pipe initially contains a stagnant water column with a depth of 0.5 m. By suddenly increasing the upstream reservoir level from 0.5 m to 6 m, a hydraulic bore is set up and moves along the pipe with a constant speed called wave speed. Over a specific period, the location of the wavefront can be easily calculated using the wave speed.

Different components of the moving hydraulic bore can be analytically calculated by applying the discrete form of continuity and momentum equations on a frame of reference moving with the same speed, as the hydraulic bore shown in Figure 3.

$$\overline{\mathit{Y}}\mathit{A}-\mathit{g}{\overline{\mathit{Y}}}_{\mathbf{0}}{\mathit{A}}_{\mathbf{0}}=\left(\mathit{V}-{\mathit{V}}_{\mathit{W}}\right)\mathit{A}\left(\mathit{V}-\stackrel{=0}{\overbrace{{\mathit{V}}_{\mathbf{0}}}}\right) \mathrm{Momentum} \mathrm{Equation}$$

(9)

$$\left(\stackrel{=0}{\overbrace{{\mathit{V}}_{\mathbf{0}}}}-{\mathit{V}}_{\mathit{W}}\right){\mathit{A}}_{\mathbf{0}}=\left(\mathit{V}-{\mathit{V}}_{\mathit{W}}\right)\mathit{A} \mathrm{Continuity} \mathrm{Equation}$$

(10)

where $\overline{\mathit{Y}}$ = distance between the HGL and the centroid of the pipe cross-sectional area, ${\overline{\mathit{Y}}}_{\mathbf{0}}$ = distance between the free surface and the centroid of the flow cross-sectional area, $\mathit{A}$ = pipe cross-sectional area, $\mathit{V}$ = flow velocity in the pressurized section, and ${\mathit{V}}_{\mathit{W}}$ = moving hydraulic bore speed.

The above equations contain three unknowns of $\mathit{V}$, ${\mathit{V}}_{\mathit{W}}$, and $\mathit{Y}$ and for being in the closed form, they need an extra equation, which is the energy balance on the upstream side of the system.

$$\mathit{H}=\mathit{Y}+\frac{{\mathit{V}}^{\mathbf{2}}}{\mathbf{2}\mathit{g}} \mathrm{The} \mathrm{energy} \mathrm{balance} \mathrm{at} \mathrm{the} \mathrm{upstream} \mathrm{end} \mathrm{of} \mathrm{the} \mathrm{conduit}$$

(11)

To calculate the three unknowns of the problem, Equations (9)–(11) can be solved simultaneously by any iterative method such as the Newton–Raphson method. Having the speed of the hydraulic bore calculated, the location of the pressurization front can be calculated at any given time after the formation of the bore. This method was applied to the hypothetical pipe system with all pipe shapes shown in Figure 1, and the results are summarized in

Table 2. Analytical solution for different pipe shapes.

Pipe Shape	Dimensions (in meter)	Y (m)	HGL (m)	V (m/s)	Wv (m/s)	Location of the Bore (m)
1	I = 1; II = 2	4.20	104.20	5.94	11.88	356.54
2	I = 1; II = 0.5	4.20	104.20	5.94	11.88	356.33
3	I = 1; II = 2	4.52	104.52	5.39	12.39	371.66
4	I = 1	3.91	103.91	6.40	11.54	346.10
5	I = 1	4.34	104.34	5.71	12.06	361.89
6	I = 1	4.32	104.32	5.73	12.06	361.69
7	I = 1	4.50	104.50	5.43	12.34	370.21
8	I = 1	4.59	104.59	5.26	12.53	375.97
9	I = 1	4.45	104.45	5.52	12.26	367.88
10	I = 1	4.50	104.50	5.43	12.34	370.16
11	I = 1	4.44	104.44	5.53	12.26	367.76
12	I = 1; III = 1; IV = 0.3	3.69	103.69	6.73	11.37	341.22
13	I = 1; III = 1; VII = 2	4.15	104.15	6.02	11.79	353.83
14	I = 1; V = 1; VI = 10	4.22	104.22	5.91	11.87	356.13

Note: I = maximum height; II = maximum width; III = top width; IV = triangle height; V = bottom width; VI = top radius; VII = bottom radius.

.

4. Numerical Results

For all test cases shown in Table 2, the model was run, and the hydraulic bore motion was numerically traced in 30 s. Since the spurious numerical oscillations are intensified at higher wave speeds, to make sure that the model provides oscillation-free solutions, the most possible acoustic speed that can occur in pipes (1400 m/s) was considered in the calculations. The numerical calculations were performed using 200 computational cells and the Courant number of 0.5, which resulted in a time step of 0.000911 s. Note that the fine computational grid was selected herein to better capture the hydraulic bore interface. The sharp interface can also give rise to small numerical wiggles in the solution regardless of the dissipative fluxes proposed in this research. Numerical exploration shows that the wiggles’ amplitude is reduced as the Courant number decreases. At Courant 0.5, it is found that the wiggles disappeared.

Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 compare the numerical results with the analytical solutions. The figures on the left compare the hydraulic grade line calculated by the model with that of the analytical solution at 30 s after the start of the simulation. The solid and dashed lines represent analytical and numerical solutions, respectively, for the figures showing the HGL along the pipe. In the figures demonstrating the HGL and flow velocity-time histories, the dashed and solid lines refer to HGL and velocity, respectively. As can be seen for all pipe shapes, the numerical results are in excellent agreement with the analytical solution, and there is no sign of numerical oscillation in the results. The figures on the right represent the time histories of the HGL and flow velocity at the middle of the conduit. It is obvious that no numerical oscillation is produced during the pressurization of the pipe for all pipe shapes considered. However, the time histories of the flow velocity contain tiny wiggles, which are significant. The calculated velocities and the pressure heads of the water column behind the hydraulic bore are compared with those obtained from the analytical solution, and the results are summarized in

Table 3. The comparison of numerical results with analytical solutions.

Pipe Index	Y (Analytical) (m)	Y (Numerical) (m)	Error (%)	V (Analytical) (m/s)	V (Numerical) (m/s)	Error (%)
1	4.2003	4.1974	0.0689	5.9423	5.9484	0.1032
2	4.2024	4.1784	0.5709	5.9388	5.9447	0.0998
3	4.5198	4.5004	0.4278	5.3891	5.3930	0.0734
4	3.9104	3.9417	0.7990	6.4029	6.4121	0.1433
5	4.3367	4.3455	0.2025	5.7125	5.7175	0.0861
6	4.3238	4.3253	0.0350	5.7348	5.7399	0.0882
7	4.4973	4.5160	0.4149	5.4298	5.4336	0.0702
8	4.5904	4.6066	0.3532	5.2590	5.2655	0.1242
9	4.4480	4.4647	0.3769	5.5183	5.5224	0.0746
10	4.4977	4.5161	0.4079	5.4291	5.4330	0.0724
11	4.4403	4.4571	0.3783	5.5319	5.5362	0.0772
12	3.6913	3.7362	1.2149	6.7303	6.7130	0.2567
13	4.1505	4.1523	0.0445	6.0239	6.0305	0.1092
14	4.2187	4.1926	0.6194	5.9118	5.9171	0.0898

. The table shows that the model captures the flow velocity for all pipe shapes quite accurately with an error ranging between 0.0702% and 0.1433%. Likewise, the calculated pressure heads are in excellent agreement with the analytical solution and the errors change from 0.0350% to 1.2149%.

Summary and Conclusions

A TPA-based numerical model was proposed for capturing transient mixed flows in a variety of pipe shapes. The model utilizes the first-order Godunov-type finite volume scheme to explicitly solve the governing equations. An HLL Riemann solver was proposed to calculate the numerical fluxes at the cells’ boundaries. To suppress the spurious numerical oscillations, the solver automatically increases the amount of the numerical viscosity of the computational cells located in the vicinity of the pressurization front.

Due to the lack of experimental results, the model performance was evaluated by using a hypothetical pipe system for which an analytical solution exists. Since the numerical oscillations intensify with increasing the pace of pressurization of the pipe, the hypothetical system was intentionally designed to give rise to a fast pace filling to make sure that the model is verified in the worst condition. Comparing the numerical results with the analytical solutions show that the model enables to successfully capture the hydraulics of the transient mixed flows quite accurately and can efficiently suppress the numerical oscillations even at as high acoustic pipe speed as 1400 m/s.

The model can perform well even at the worst condition, which is the combination of using very high pipe acoustic speeds and fast rates of pipe filling, confirming that it can be safely used for calculating transient mixed flows in real pipe systems, which may consist of a verity of pipe shapes.