A Modified Shock Tracking–Capturing Finite-Volume Approach for Transient Mixed Flows in Stormwater Systems

The gravity flow in stormwater systems can be altered to a transient, partially pressurized flow by rapid filling during intense rainfall events. The study of these transient flows is motivated by operational problems such as geysering through drop shafts, structural damages, as well as serious public health and environmental issues. The rapid filling process may not properly push the air above the free surface flow toward the boundaries so that an air pocket can be entrapped in the systems. As indicated in the literature (e.g., Malekpour et al. 2016; Wylie and Streeter 1993), in the presence of an air pocket, the maximum induced transient pressure head can be significantly larger than the driving pressure head. In addition, based on observations, Wright et al. (2017) showed that the pressurized air pocket is mainly responsible for the geyser formations, which is contrary to previous studies such as Guo and Song (1990), in which the geyser was only linked to the inertial surge.

Because experimental investigations of full-scale systems are impractical due to large sizes, numerical simulations are desirable to efficiently calculate the relevant variables of transient flows in terms of accuracy and cost. For this purpose, the focus of numerical investigations is on one-dimensional (1D) models, in which the air pocket is assumed intact and simulated using the thermodynamic relationship of an ideal gas. The rigid column (RC) model is probably the simplest 1D model, which, based on the mass and momentum balances, solves ordinary differential equations (ODEs) associated with the governing equations of the pressurized flow zone. This model either neglects the effects of the free-surface flow zone or treats this zone as a rigid column as well. Although the RC model is able to reproduce the fundamental features of this kind of transient flow and is easy to implement, it overestimates the peak values of the transient flow variables, which is due to incompressibility assumption and due to the presence of large air volume in cases with large velocity during pressurization and low acoustic wave speed (Malekpour and Karney 2014). In addition, as shown by Zhou and Liu (2013), the initial depth of the free-surface flow impacts the maximum air pressure. The studies of the RC model can be exemplified by the works of Hamam and McCorquodale (1982), Li and McCorquodale (1999), Zhou et al. (2002), and Vasconcelos and Wright (2003).

The shock-capturing approach only solves a single set of the Saint-Venant equations for both flow regimes, which is based on similarities between the governing equations of the pressurized and free-surface flows. Thus, there is no need to track the moving interface between pressurized and free-surface flows because this interface is captured as a shock wave. The Preissmann slot is the most common shock-capturing approach, in which the pressurized flow is simulated by a hypothetical narrow slot on top of the conduit. However, the Preissmann slot is unable to calculate subatmospheric pressures, which may occur in the pressurized flow zone. This approach, which was originally proposed by Preissmann (1961), can be exemplified by the works of Cunge (1980), Garcia-Navarro et al. (1994), Capart et al. (1997), Ji (1998), and Trajkovic et al. (1999). To remove the restriction of the Preissmann slot approach, Vasconcelos et al. (2006) proposed the Two-component Pressure Approach (TPA), in which the elastic behavior of the pipe allows simulating overpressurizations. The main problem with the TPA approach is the oscillatory behavior of the solutions near the pressurized and free-surface flows interface. Thus, in the calculation, small acoustic wave speeds need to use, which may compromise the accuracy of the solutions.

The shock-fitting approach is based on using different governing equations to best mimic physics. Thus, in this approach, the Saint-Venant equations are applied to the free-surface flow, and either the water hammer or RC model is applied to the pressurized flow. The combination of the water hammer and Saint-Venant equations can be exemplified by the works of Cardle (1984) and Fuamba (2002), in which the method of characteristics (MOC) was used for the discretization and the works of León (2006) and León et al. (2010), in which a finite-volume (FV) method was used for the discretization. The MOC method is a finite-difference method, which does not conform to the conservation law. In addition, when the free-surface flow becomes supercritical this method is unable to calculate the solutions. Therefore, an FV method is more desirable, which is applied to the conservative form of the governing equations so that it is able to properly conserve the mass and momentum. In addition, this method is able to automatically handle the supercritical conditions that may occur in the free-surface flow zone because interface fluxes are calculated by solving a Reimann problem (León et al. 2010).

Rokhzadi and Fuamba (2020), using the shock-fitting approach, which is a combination of the RC model and the Saint-Venant equations, calculated an air pocket entrapment problem in a partially pressurized transient flow. They showed that, compared to the fourth-order Runge–Kutta explicit time integration scheme, the dissipation property of the backward Euler implicit scheme can take the extra energy, which is left due to the RC model assumptions (Malekpour and Karney 2014) as well as due to neglecting the heat transfer process between the air and pipe wall (Lee 2005). Therefore, it was shown that this shock-fitting approach can improve overpredicted peak values and it can predict the attenuation behavior observed in the solutions. In addition, Rokhzadi and Fuamba (2021) analyzed the contribution of numerical dissipation on the RC model’s ability to solve the air pocket entrapment problem in transient flows. They showed that the time step size directly controls the numerical dissipation of implicit schemes.

However, the RC model is not able to calculate the pressure inside the pressurized zone, in which the pressure may exceed the system limit (Vasconcelos et al. 2006). In addition, the RC model cannot fully describe interactions between multiple waves.

For numerical simulation of transient flows in stormwater systems, a mathematical model and the discretization method used need to be able to calculate the flow variables during normal operations as well as during unusual events such as rapid filling, which may cause the free-surface flow to change to a partially pressurized flow. Among different models proposed for simulating this transient flow, the shock-fitting approach proposed by León et al. (2010) that was called a shock tracking–capturing method, and hereafter called the STC approach, was shown to properly handle different conditions including supercritical flows and a wide range of acoustic wave speeds.

Following the literature (e.g., León et al. 2010), the conservative form of the governing equations are presented aswhere $x$ and $t$ are the space and time variables. The vector of variables ($\mathit{U}$), the flux vector ($\mathit{F}$), and the vector of source term ($\mathit{S}$) for pressurized flow iswhere $\mathrm{\Omega }$ is the mass per unit length, $Q_m$ is the mass discharge, $p$ is the absolute pressure at the gravity center, $A$ is the cross-sectional area of the pipe, $S_0$ is the bed slope, $S_f$ is the friction loss, and $g$ is the gravitational acceleration equal to $9.81\mathrm{m}{\mathrm{s}}^{-2}$. Note that the friction loss is calculated using the Darcy–Weisbach equation with steady-state friction factor aswhere $f$ is the friction factor, and $D$ is the pipe diameter. Note that the absolute sign ensures that the head loss is applied in an appropriate direction.

As can be seen in Eq. (2), because the unknown variables ($\mathrm{\Omega },Q_m,p$) are more than the equations, an auxiliary equation is required to close the system of equations. Following León (2006), León et al. (2010), and Guinot (2003), the definition of the acoustic wave speed is used aswhere $a$ is the acoustic wave speed. Taking integral of Eq. (4) yieldswhere the subscript ($ref$) refers to the reference state, which is the location where pressurization of a control volume of the free-surface flow next to the moving interface occurs and all flow parameters in both free-surface and pressurized flow regimes are identical (León et al. 2010). The criterion for pressurization is defined when the piezometric depth of the free-surface flow ($y$) exceeds $y_{ref}=0.95D$ (León et al. 2010; Yuan 1984). Note that the parameter $A_{ref}$ is the hydraulic area below $y_{ref}$. It should also note that, during the computation, for stability reasons the pipe cross-sectional area ($A$) in Eq. (5) was replaced by $A_{ref}$. In addition, $P_{ref}$, and ${\rho }_{ref}$ are set to 101 kPa and $\mathrm{1,000}\mathrm{kg}/{\mathrm{m}}^3$, respectively, and ${\mathrm{\Omega }}_{ref}$ is calculated as ${\mathrm{\Omega }}_{ref}={\rho }_{ref}{A}_{ref}$.

The conservative form of the Saint-Venant equations applied to the free-surface flow can be formed by the following components, which are substituted into Eq. (1):where ${\rho }_f$ is the water density in the free-surface flow, assumed constant, $A_f$ is the cross-sectional flow area, and $Q_f$ is the flow discharge. The friction loss ($S_f$) is calculated using the Chézy formula aswhere $R_h$ is the hydraulic radius, $n$ is the Manning roughness coefficient, and $V$ is the velocity. The variable $h_a$ is the absolute air pressure head, and $\overline{p}$ is the average pressure of water column over the cross-sectional area. Note that, following León (2006), the term $A_f\overline{p}$ is calculated aswhere $y$ is the piezometric depth of the free-surface flow.

Following León et al. (2010) and Martin (1976), among others, the air pocket behavior is assumed to be simulated by the time derivative of the polytropic thermodynamic relationship aswhere $k$ is the polytropic exponent, and ${\forall }_a$ is the air volume. The air volumetric change caused by differences between discharge flowrates at the air pocket boundary (Fig. 1) is calculated as

The governing equations are discretized using an FV method, which for pure free-surface flow zone is the same as León et al. (2010) and León (2006). Thus, the variable ($\mathit{U}$) at a new time step is calculated aswhere “$i$” is the control volume index, for which the interfaces are located at “$i+1/2$” and “$i-1/2$”.

The superscript $n$ represents the current time level, and $\mathrm{\Delta }t$ and $\mathrm{\Delta }x$ are the time and space increments, respectively. The fluxes at the interfaces $(F_{i\pm 1/2}^n)$ are calculated using the Harten–Lax–van Leer (HLL) solver. The details of the HLL solver and the explicit FV method can be found in León (2006).

In the MSTC approach, the pressurized flow zone is placed in one control volume. Therefore, one of the cell interfaces is calculated as a boundary condition. To calculate the flux at this interface, the solutions at the time level $n+1$ are updated using one of the characteristics originating from the cell center of the pressurized flow and cuts through this boundary, as explained in León (2006). The other cell interface of the pressurized flow is the moving interface between pressurized and free-surface flows, in which the flux is calculated, as explained in the following sections. Therefore, the discretized governing equations are aswhere the index “1” indicates that there is only one control volume, in which the source term and the flux at the moving interface $({\mathit{F}}_{3/2})$ are calculated implicitly, while the flux at the other boundary is calculated explicitly. It is worth mentioning that Eqs. (2) and (12) show that the MSTC approach takes the time variation of the water compressibility of pressurized flow into account, which is the main difference from the RC model. Note that $L$ is the length of the pressurized flow, which is calculated aswhere $w$ is the speed of the moving interface between pressurized and free-surface flows. Note that the implicit scheme needs an iterative approach, in which the residuals must converge up to a specific tolerance, which was set to ${10}^{-16}$. It is worth mentioning that the condition presented in Eq. (14) was imposed on the time stepto ensure that the length variation of the pressurized flow at each time step does not exceed one length of the control volume of the free-surface flow. Therefore, the time step is calculated as the minimum time step associated with the Courant–Friedrichs–Lewy (CFL) conditions of the free-surface and pressurized flows, and the time step calculated in Eq. (14).

The positive interface was defined by Cardle (1984) as an interface moving toward the free-surface flow. The flow across a positive interface is equivalent to a flow across a hydraulic bore, which, in the relative coordinate system, the flow changes from supercritical to subcritical flow (Cardle 1984). Fig. 2 illustrates the positive interface moving downstream. As explained by León et al. (2010), to calculate the flux $({\mathit{F}}_{3/2}^{n+1})$, the variables at the interface $({\mathit{U}}_{3/2}^{n+1})$ is needed, which is equal to the flow variables at the left of the moving interface at time level $n+1$ (${\mathit{U}}_P$). The flow variables at the left of the moving interface is related to the flow variables at the right of this interface $({\mathit{U}}_F)$ as well as to $w^{n+1}$. As shown by Cardle (1984), in the relative coordinate system, a transition from a supercritical flow to a subcritical flow across the positive interface implies that the positive characteristic of the control volume of the free-surface flow next to the moving interface always cuts through the interface $(u_R+c_R<w)$, in which $c$ represents the gravity wave speed calculated as (León 2006)where $\theta =2{\cos }^{-1}(1-2y/D)$.

In addition, it is obvious that the negative characteristic of the same control volume of the free-surface flow always cuts through the interface. Thus, as also indicated by León et al. (2010), ${\mathit{U}}_F$ is equal to ${\mathit{U}}_R$. Therefore, to calculate ${\mathit{U}}_P$ and $w^{n+1}$, three equations are required, including two equations of the Rankine–Hugoniot condition applied to the mass and momentum conservations, andwhich holds along the positive characteristic of the control volume of the pressurized flow $(u_L+a)$ next to the moving interface (Fig. 2). More details of these auxiliary equations of a positive moving interface can be found in León et al. (2010).

The negative interface, illustrated in Fig. 3, is defined as the interface moving toward the pressurized flow. To calculate the flow variables at the interface $({\mathit{U}}_{3/2}^{n+1})$, the flow variables ${\mathit{U}}_F$, ${\mathit{U}}_{*}$, and ${\mathit{U}}_R$ are required to solve a free-surface flow Riemann problem. The details of this problem and how to solve it can be found in Toro (2001), and León (2006).

Like the positive interface moving, ${\mathit{U}}_F$ is related to ${\mathit{U}}_p$ and $w^{n+1}$, which are five unknown variables. Two equations of Rankine–Hugoniot conditions applied to the mass and momentum conservations across the interface can be used. The other equation iswhich holds along the positive characteristic of the pressurized flow $(u_L+a)$. Following León et al. (2010), the speed of the moving interface is calculated from the previous time step ($t_n$) as

One extra equation is required to close the system of equations. If the flow in the cell of the free-surface flow at the right of the moving interface is subcritical, then, following Cardle (1984), the negative characteristic originating from this cell cuts through the moving interface $(u_R-c_R<-w)$. Therefore, this characteristic can be used as the extra equationwhere $\phi =\int (c/A_f)d{A}_f$. Following León (2006), $\phi$ is calculated as

Fig. 3 shows that if the flow in the control volume of the free-surface flow at the right of the moving interface is supercritical, the negative characteristic has no interaction with the interface, thus, Eq. (19) is no longer valid. Following León et al. (2010), for supercritical flows, the energy equation across the negative interface is used, which is aswhere $h$ is the absolute pressure head relative to the pipe invert either at pressurized or free-surface flow zone. Note that for an air pocket problem, in which the air can get pressurized, the absolute air pressure head needs to be included in $h$ on the left-hand side of Eq. (21). As can be seen in Eq. (21), the pressure loss is neglected across the negative interface because the negative interface is similar to an expansion wave, which contains negligible head loss (León et al. 2010).

An additional remark is necessary to explain regarding the negative interface moving in test cases with the supercritical free-surface flow. As previously explained, the flow variables at the star zone (i.e., $U^{*}$ in Fig. 3) need to calculate. Using the MSTC approach, first, the flow variable $A_f$ in the star zone is calculated, by which $\theta$ is calculated iteratively using the Newton–Raphson method and the equation of $A_f=D^2/8(\theta -\sin \theta )$. Then, other parameters will be calculated. In some iteration, $A_f$ may become larger than $A$ so that the Newton–Raphson method diverges. In these situations, it is necessary to replace $A_f$ with $A_{ref}$.

Following León et al. (2010) and Yuan (1984), the threshold for pressurization of a control volume is when the piezometric depth ($y$) exceeds $y_{ref}=0.95D$. Note that except for the first cell of the free-surface flow, which is next to the moving interface, when pressurization occurs in a cell “$i$”, the following modifications need to be done:

When pressurization occurs at the first cell of the free-surface flow next to the moving interface, this cell will be included in the pressurized flow. In addition, the length of the pressurized flow needs to be updatedwhere, before computation begins, $L^{*}$ is equal to the initial water column length. However, as previously explained, the water column length in each time step is calculated by Eq. (13), implying that the water column length does not match a cell’s interface. Therefore, to be sure that, during each pressurization, only one control volume is included in the pressurized flow zone, immediately after Eq. (23), $L^{*}$ will be updated as

As explained by León et al. (2010), during pressurization, numerical instability may occur, for which the piezometric depths of the free-surface flow need to be modified. Further details of these modifications can be found in León et al. (2010), Section 4, “Threshold levels for piezometric depths.”

The threshold for depressurization is when the piezometric head of the pressurized flow drops below $0.84D$ (Yuan 1984). Following León et al. (2010), the piezometric head is calculated aswhere ${\mathrm{\Omega }}_d$ in the STC approach is the mass per unit length at the last cell of the pressurized flow next to the moving interface, while, in the MSTC approach, it is the mass per unit length at the cell center. Note that, besides this threshold, an extra condition needs to be satisfied to ensure that depressurization occurs. The extra condition is that the absolute pressure head relative to the pipe invert at the first cell of free-surface flow must be larger than the same pressure head at the adjacent cell of the pressurized flow. When depressurization occurs, a new control volume with the length of $\mathrm{\Delta }x$ is excluded from the pressurized zone, for which the variable $({\rho }_f{A}_f)$ is calculated aswhere in the STC approach, ${\mathrm{\Omega }}_{bR}$ is the mass per unit length of the last cell of the pressurized flow, while in the MSTC approach, it is equal to ${\mathrm{\Omega }}_{3/2}$. Then, the other variables associated with this new cell will be calculated. It should be noted that in the test case with small driving pressure, subsequently presented, first, the water depth of the new control volume is set to the piezometric head presented in Eq. (25), then, other variables will be calculated accordingly. The reason is that in this test case, the water density does not change noticeably. Therefore, Eq. (26) causes the whole pressurized zone to depressurize. In addition, the length of the pressurized flow, after a depressurization, will be updated as

Note that, after each depressurization, the parameter $L^{*}$ needs to be updated as in Eq. (24). Similar to León et al. (2010), in the present study, no numerical instabilities were observed during the depressurization process.

The abilities of the MSTC approach in simulating transient flows under different conditions are examined and compared with the STC approach. In the first two examples that were experimentally studied by Zhou (2000), the air pocket entrapment occurs in a partially pressurized flow in a reservoir–pipe system with a large driving pressure and with different acoustic wave speeds and different initial air lengths. Then, both STC and MSTC approaches are used to solve another test case, also experimentally studied by Zhou (2000), with the dry-bed condition and with air release and with a larger driving pressure than the previous two examples. Note that in all these examples, a supercritical flow regime develops in the free-surface zone. Afterward, another partially pressurized flow following air pocket entrapment, which was experimentally studied by Hatcher et al. (2015), will be solved in which the driving pressure is slightly greater than the atmospheric pressure so that a subcritical flow develops in the free-surface zone. Then, both STC and MSTC approaches are used to predict the bore advancement following a rapid filling process, studied by Vasconcelos et al. (2006), which pushes the air toward a downstream manhole. Finally, a pressurized flow in a reservoir–pipe system, in which sudden obstructing downstream end causes the water hammer phenomenon, will be solved by both MSTC and STC approaches.

Zhou (2000) experimentally studied transient flows in a horizontal galvanized steel pipe with a diameter of 35 mm connected to a pressure tank at the upstream end, with a constant pressure head.

Three different test cases with different initial air and water column lengths were taken from this experiment, in which the transient flow with and without air entrapment was analyzed. A valve separates two columns, and it takes 0.06–0.08 s to be opened. The upstream column was pressurized, and the downstream column was either fully filled with air or was a free-surface water flow contacting with the atmospheric air. The downstream end is either closed to study the air entrapment problem or partially open with a cap containing a centered, sharp-edged orifice to study the air release and dry-bed condition.

For the examples with air entrapment, the total pipe length is $L_t=8.96\mathrm{m}$, and the nondimensional initial pressurized flow lengths $({\lambda }_0=L/L_t)$ are 0.56 and 0.89, and the initial water depth in the free-surface flow is $y/D=0.4$. The absolute pressure head in the reservoir is $H_0^{*}/H_b^{*}=2.43$, in which $H_b^{*}$ is the atmospheric pressure head. Note that the Manning roughness and Darcy–Weisbach coefficients were taken as $0.012{\mathrm{m}}^{1/6}$ and 0.035, respectively. As Zhou (2000) indicated, the acoustic wave speed, in test cases with ${\lambda }_0=0.56$ and 0.89 were taken as $a=200$ and $700\mathrm{m}/\mathrm{s}$, respectively. Also, the polytropic coefficient was set to $k=1.4$ for both approaches.

In order to show that the results are mesh independent, three different mesh sizes are used to solve the test case with ${\lambda }_0=0.89$.

Fig. 4 shows the air pressure distributions solved by both approaches using different control volume sizes. Note that the lengths of the control volumes $(\mathrm{\Delta }x)$ are 0.11, 0.071, and 0.052 m, which correspond to 10, 15, 20 control volumes set in the free-surface zone, respectively. Note that $\mathrm{\Delta }x=0.052\mathrm{m}$ is the finest grid size, which both approaches tolerate. Also, in the STC approach, these numbers of control volumes result in the total number of 84, 129, and 175 control volumes, respectively, applied to the entire flow. Also, $\mathrm{\Delta }x=0.052\mathrm{m}$ is the finest grid size, which the MSTC approach can tolerate. In addition, as previously explained, the MSTC approach uses one control volume for the pressurized zone. It is worth mentioning that for the mesh independence test, the CFL number was set to 0.7 for both STC and MSTC approaches. As can be seen in Fig. 4, the results are almost the same. It should also be noted that for all test cases subsequently presented, the number of control volumes set in the free-surface zone was selected such that the length of the control volume becomes $\mathrm{\Delta }x=0.071\mathrm{m}$. Thus, it can be expected that the results of the other test cases are mesh-independent as well. Note that, with this control volume length, both approaches in some test cases cannot tolerate $\mathrm{CFL}=0.7$. Therefore, in order to provide similar conditions, all test cases are solved by both STC, and MSTC approaches with $\mathrm{CFL}=0.5$.

Fig. 5 shows the air pressure distributions of both test cases with ${\lambda }_0=0.89$ on the left and ${\lambda }_0=0.56$ on the right. As can be seen, the peak values and the damping predicted by the MSTC approach are less deviated from the experimental results than those predicted by the STC approach, while the phase shift caused by both approaches is almost similar. Note that these improvements are more obvious in the test case with ${\lambda }_0=0.56$, in which the initial air length is larger. Also, in this test case, the phase shift between numerical and experimental results is less significant. As Zhou (2000) indicated, the phase shift may cause by assuming that the air pocket remains intact, while it could roll up and split into several small pockets. It seems plausible because in the test case with smaller air length (i.e., ${\lambda }_0=0.89$), a larger portion of the air pocket can split into smaller pockets and the implication of assuming an intact air pocket becomes more obvious.

In order to show that the MSTC approach can capture the depressurization of the pressurized zone, the variation of the pressure head relative to the threshold value (i.e., $H_p=h-0.84D$) calculated by both approaches is shown in Fig. 6. Note that, using the MSTC approach, the pressure head at the cell center of the pressurized flow and, by using the STC approach, the pressure head at the last cell of the pressurized flow are shown in Fig. 6. Also, the depressurization can be realized by sharp drops of the pressure to zero. As can be seen, the MSTC approach can capture the depressurization competitively to the STC approach.

Another test case of Zhou (2000) has been solved by both approaches to compare their abilities in a dry-bed condition and with air release. In this test case, the pipe length is 10 m, the reservoir absolute pressure head is $H_0^{*}/H_b^{*}=4.57$, the upstream pressurized column has the initial length ratio of ${\lambda }_0=0.8$, and the downstream column is full of air with the atmospheric pressure. The downstream end is partially open with a cap containing a centered, sharp-edged orifice plate. Following Zhou (2000), the discharge coefficient of the plate is taken as $C_d=0.65$ and the dimeter of the orifice plate is taken as $d_0=0.028D$. Also, the acoustic wave speed is taken as $a=\mathrm{1,000}\mathrm{m}/\mathrm{s}$. Note that a common practice to solve a dry-bed condition is to assume a very thin layer of water because the transient flow in sewer systems does not begin with a dry-bed condition (León et al. 2006; Abbott and Minns 1997). Thus, in the present paper, for both STC and MSTC approaches, the initial depth of free-surface flow is assumed as $y=0.01D$.

It should be noted smaller initial depths up to $y=0.005D$ were also examined. However, the results (not shown here) are not significantly different. Following León et al. (2010) and Zhou (2000), the air pressure variation is calculated aswhere ${\mathrm{m}}_a$ is the air mass and the air mass rate is calculated aswhere ${\rho }_a$ is the air density and $A_0$ is the orifice area. Note that the air volume change is calculated using Eq. (10). After calculation at each time step, the air density is updated as ${\rho }_a=m_a/{\forall }_a$.

Fig. 7 shows the air pressure distribution calculated by both approaches, which are compared to the experimental data of Zhou (2000). The first noticeable feature is a large phase shift, for almost one period, between numerical results and the experimental data. However, as previously explained, this phase shift is due to the intact air pocket assumption in the calculation. Compared with Fig. 5, a larger phase shift is expected due to more severe conditions, including larger driving pressure for almost twice and the dry-bed condition. These conditions may cause the air to more strongly roll up and split to smaller air pockets. However, as can be seen in Fig. 7, the phase of the STC and MSTC approaches are almost the same, even though the MSTC approach seems to less slightly deviate from the experimental data. Also, the MSTC approach results in closer peak values and damping to the experimental results. Another possible reason for the phase shift is the discharge coefficient of the orifice plate. As Zhou (2000) indicated, the discharge coefficient depends on many factors and more importantly it oscillates with time. It should also be noted that León et al. (2010) did not find the same phase shift. However, their paper did not report the discharge coefficient’s value.

In order to show that the MSTC approach does not compromise the conservation law, for this test case, the log scale of the ratio of the difference between the water mass calculated at each time step and the initial water mass to the initial water mass is calculated and shown in Fig. 8. As can be seen, the MSTC approach can appropriately preserve mass conservation compared with the STC approach.

The experiment of Hatcher et al. (2015) includes a pipe section with adverse slopes, in which a reservoir supplies the water at the upstream in a pressurized flow regime. The water somewhere at the downstream detaches from the pipe wall, and a free-surface flow develops, which freely releases to the atmosphere at the downstream end. A valve was used to block the downstream end. One test case of this experiment has been solved in which the pipe diameter and length are $D=53\mathrm{mm}$, and $L_t=10.7\mathrm{m}$, respectively. The pipe is adversely sloped with $S_0=2.0\%$. The initial flow rate is $Q^{*}=0.15$, in which $Q^{*}=Q/\sqrt{g{D}^5}$, and the initial air pocket volume is $V_a^{*}=2.63$, in which $V_a^{*}=V_a/D^3$. The initial pressure head at the reservoir is assumed as $H_0^{*}=H_b^{*}+0.204\mathrm{m}$. To calculate the initial water column length, the free-surface flow with a uniform cross-sectional area and a constant depth, which is assumed as $y/D=0.5628$, is considered. Note that this flow depth corresponds to the water depth in a horizontal frictionless pipe (Alves et al. 1993; Benjamin 1968). Thuswhere $V_a^0$ is the initial air pocket volume. Note that for this test case ${\lambda }_0=L/L_t=0.96$. The Manning roughness coefficient and the Darcy–Weisbach friction factor were set to $0.012{\mathrm{m}}^{1/6}$ and 0.025, respectively. In addition, a local head coefficient $(K_{loss})$ was considered, which was set to 2.9. Note that the acoustic wave speed in this example was set to $a=\mathrm{1,200}\mathrm{m}/\mathrm{s}$. In addition, the polytropic coefficient was set to $k=1.2$.

Fig. 9 illustrates the distributions of the nondimensional air pocket pressure head $(H_a^{*}=H_a-H_b^{*}/D)$ in the left and the nondimensional flow rate in the right. As can be seen, the MSTC approach effectively calculates the peak values and damping compared with the STC approach.

Note that, in this example, the flow in the free-surface flow zone is subcritical, thus, the negative characteristic originating from the first cell of the free-surface flow cuts through the moving interface. Therefore, to treat a negative interface moving, this characteristic can be used for calculation. For further clarification, Fig. 10 illustrates the speed of the positive and negative characteristics of the first cell of the free-surface flow next to the moving interface and the speed of the moving interface. As shown in Fig. 10, the speed of the negative characteristic is lower than the moving interface speed, which implies interactions of this characteristic with the moving interface. Note that the main difference between this test case and Zhou (2000) is the reservoir pressure, which in this test case is much less. Comparing the left graph of Fig. 9 to Fig. 5 shows that in Fig. 9, the phase shift between numerical solutions and the experimental data is insignificant. Thus, the statement of Zhou et al. (2002), that linked the phase shift to the air pocket splitting to smaller air pockets, seems to be appropriate because in Fig. 9 the driving pressure is slightly different from the atmospheric pressure and the air pocket is more likely to remain intact.

Vasconcelos et al. (2006) studied different experiments, one of which is the simulation of rapidly pipe filling events. The experiment was done by means of a horizontal pipe with a 9.4-cm diameter and a 14.33-m long transparent acrylic pipe. The initial water depth throughout the pipe is $y/D=0.78$ with zero velocity. An inflow $Q=3.1\mathrm{L}/\mathrm{s}$ is suddenly supplied by a reservoir connected to the upstream end with sides of $0.25\times 0.25\mathrm{m}$ and a height of 0.31 m relative to the pipe invert. A hydraulic bore is generated in the pipe and travels toward a surge tank, which is connected to the downstream pipe end with a diameter of 0.19 m. The spill level in the surge tank is much higher than that of the reservoir to ensure that no flow is lost through the surge tank. In addition, to ensure that the air will not escape through the surge tank, a gate was installed at the downstream end. Also, to avoid air pressurization during the bore advancement, a ventilation tower with a 3.2 cm diameter was installed upstream of the gate. When the bore arrives at the surge tank, the water level in the upstream reservoir increases and eventually starts spilling while inertial oscillations continue in the surge tank. Afterward, a fully pressurized flow develops, in which the oscillations have been damped by system losses. Note that the Manning roughness coefficient and the Darcy–Weisbach friction factor were set to $0.012{\mathrm{m}}^{1/6}$ and 0.025, respectively. In addition, the local head coefficient $(K_{loss})$ was considered, which was set to 2.9. The velocity and pressure at a distance of 9.9 m downstream of the reservoir were measured by Vasconcelos et al. (2006).

Following Vasconcelos et al. (2006), the acoustic wave speed was set to $25\mathrm{m}/\mathrm{s}$. Fig. 11 illustrates the velocity and pressure distributions calculated by both STC and MSTC approaches, which are compared to the TPA approach as well as the experimental data of Vasconcelos et al. (2006). Note that, in this test case, the minimum mesh size, which the MSTC approach tolerates, was $\mathrm{\Delta }x=0.7\mathrm{m}$. From the right graph, it can be seen that the MSTC approach predicted the experimental pressure distribution with less deviation than the other two approaches. In addition, less oscillatory behavior exhibited by the MSTC approach can be seen in both graphs. Also, the MSTC approach calculates the velocity peak values more accurately than the TPA approach, and it is competitive to the STC approach.

The last test case solved in the present study is the simulation of a pressurized flow in a reservoir–pipe system, in which the downstream end is suddenly closed by a valve. Following León (2006), the flow is assumed frictionless in a pipe with the diameter and length of 1 and 10,000.0 m, respectively. The reservoir pressure head is assumed constant and set to 200.0 m, in which the water is constantly supplied with an inflow of $2.0{\mathrm{m}}^3/\mathrm{s}$. Note that the acoustic wave speed is set to $a=\mathrm{1,000}\mathrm{m}/\mathrm{s}$. Also, CFL is set to 0.9, and by using the STC approach, the number of control volumes is set to 100. Fig. 12 shows the gauge pressure head calculated by the STC, and MSTC approaches. Note that, in the MSTC approach, the pressure head at the cell center is calculated and, in the STC approach, the pressure head at the last cell is calculated. As can be seen, the MSTC approach can capture the peak values and the frequency of the pressure variation as accurately as the STC approach.

In this paper, following the shock-fitting approach proposed by León et al. (2010), called STC, a similar approach has been proposed to further improve simulation of transient mixed flows in stormwater systems. The proposed approach, called MSTC, considers the whole pressurized flow zone in one control volume with the length of the pressurized flow. Following León et al. (2010), the governing equations have been discretized using an FV method and HLL solver to be able to handle supercritical and subcritical flows automatically. Note that the governing equations of the pressurized flow have been solved using an implicit FV method, while, similar to León et al. (2010), the governing equations of the free-surface flow are discretized using an explicit FV method. These FV methods have the first order of accuracy near shocks and discontinuities and have the second order of accuracy away from shocks The ability of the proposed approach has been examined by solving different transient flow problems, including pressurized flow and water hammer phenomenon, partially pressurized flow following air entrapment with and without air release, and rapidly pipe filling event following bore advancement. Also, different conditions have been considered, including dry-bed condition, super and subcritical flows, different air and water column lengths, different driving pressure heads, and a wide range of acoustic wave speeds. It was shown that in all test cases, the MSTC approach calculates the peak values and damping with less deviation than the STC approach. Especially, in simulating the bore advancement, the MSTC approach exhibits less oscillatory behavior.

All data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

The writers would like to express their gratitude to the Natural Sciences and Engineering Research Council of Canada (NSERC) for the financial support.