A Modified Shock Tracking–Capturing Finite-Volume Approach for Transient Mixed Flows in Stormwater Systems | Journal of Hydraulic Engineering | Vol 148, No 11

AbstractThis paper proposes a modification to a previous shock tracking–capturing approach used for simulating transient flows in stormwater systems. Based on the Godunov-type finite-volume (FV) method, this modified approach uses two governing equations for free-surface and pressurized flows discretized. The free-surface flow is simulated using an explicit FV method. However, the pressurized flow is discretized on one control volume and simulated using an implicit FV method. Different examples will be solved to examine the abilities of the proposed modified approach. The examples are pressurized and partially pressurized flows, including the water hammer, the pipe-filling bore advancement, and air pocket entrapment problems. The partially pressurized flow following air pocket entrapment is analyzed under different conditions including various ranges of driving pressure, air and water lengths, as well as dry-bed and with and without air-release conditions. It will be shown that the proposed modification can result in improving the overestimated peak values and the underestimated attenuation of pressure distribution. In addition, it will be shown that the proposed modified approach exhibits less oscillatory behavior, particularly in simulating pipe-filling bore advancement.IntroductionThe 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.STC ApproachThe STC approach discretizes conservative forms of the governing equations of the pressurized and free-surface flows using an explicit FV method and treats the interface between these flows using a set of equations depending on either the interface moves toward free-surface or pressurized flow. These equations will be subsequently explained in detail. It is worth mentioning that the convergence and stability of the Godunov scheme was studied in the literature (e.g., Bressan and Jenssen 2000).However, it is known that the main operational problems in sewer systems may occur during rapid filling events caused by air entrapment. Therefore, it seems unnecessary to simulate the pressurized flow using a fine grid mesh, while a coarse grid mesh and even just one grid mesh may be sufficient. In addition, as subsequently shown, the STC approach, similar to the TPA approach, exhibits an oscillatory behavior in the bore advancement problem. In addition, as indicated by Malekpour and Karney (2011), if the water column is not internally disturbed, the RC model is the cheapest model, which has been shown that it can capture the fundamental features of the transient flow in sewer systems.ObjectiveThis study proposes a modification to the STC approach, hereafter called the MSTC approach, to take advantage of the STC approach and the rigid column theory. The MSTC approach includes the entire pressurized zone in one control volume with the size of pressurized flow length. Therefore, this modified approach allows the pressure and density of the pressurized flow to vary with respect to the time and particularly, along the pipe axis because the pressure and density at the boundaries and center of the control volume are different. In addition, the MSTC approach calculates the governing equations of the pressurized flow using an implicit FV method, while the discretized governing equations of the free-surface flow, similar to León et al. (2010), are calculated by an explicit FV method. Note that these FV methods have the first order of accuracy near shocks and discontinuities and have the second order of accuracy away from shocks. Note that Karney (1990) showed that the RC model due to the water incompressibility assumption is only valid when the potential energy of the pressurized flow is negligible compared with the kinetic energy. Thus, because the MSTC approach considers water compressibility, as subsequently discussed, it could be expected that this approach is applicable to different transient flows in sewer systems.Governing EquationsFollowing the literature (e.g., León et al. 2010), the conservative form of the governing equations are presented as (1) where x and t are the space and time variables. The vector of variables (U), the flux vector (F), and the vector of source term (S) for pressurized flow is (2) U=[ΩQm],F=[QmQm2Ω+Ap],andS=[0(S0−Sf)gΩ]where Ω is the mass per unit length, Qm is the mass discharge, p is the absolute pressure at the gravity center, A is the cross-sectional area of the pipe, S0 is the bed slope, Sf is the friction loss, and g is the gravitational acceleration equal to 9.81  m s−2. Note that the friction loss is calculated using the Darcy–Weisbach equation with steady-state friction factor as (3) where 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 (Ω,Qm,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 as (4) where a is the acoustic wave speed. Taking integral of Eq. (4) yields (5) P=ArefAPref+a2A(Ω−Ωref)where 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 yref=0.95D (León et al. 2010; Yuan 1984). Note that the parameter Aref is the hydraulic area below yref. It should also note that, during the computation, for stability reasons the pipe cross-sectional area (A) in Eq. (5) was replaced by Aref. In addition, Pref, and ρref are set to 101 kPa and 1,000  kg/m3, respectively, and Ωref is calculated as Ωref=ρrefAref.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): (6) U=[ρfAfρfQf],F=[ρfQf(ρfQf)2ρfAf+Afp ¯+Aρgha],andS=[0(S0−Sf)gρfAf]where ρf is the water density in the free-surface flow, assumed constant, Af is the cross-sectional flow area, and Qf is the flow discharge. The friction loss (Sf) is calculated using the Chézy formula as (7) where Rh is the hydraulic radius, n is the Manning roughness coefficient, and V is the velocity. The variable ha is the absolute air pressure head, and p ¯ is the average pressure of water column over the cross-sectional area. Note that, following León (2006), the term Afp ¯ is calculated as (8) Afp ¯=112ρfg[(3D2−4Dy+4y2)y(D−y)−3D2(D−2y)tan−1(yD−y)]where 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 as (9) where k is the polytropic exponent, and ∀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 (10) d∀adt=A[(QmΩ)r−1/2−(QmΩ)j+1/2]Numerical DiscretizationThe 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 (U) at a new time step is calculated as (11) Uin+1=Uin−ΔtΔx(Fi+1/2n−Fi−1/2n)+Δt[(S0−Sf)gρfAf]inwhere “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 Δt and Δx are the time and space increments, respectively. The fluxes at the interfaces (Fi±1/2n) 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 as (12) U1n+1=U1n−ΔtLn+1(F3/2n+1−F1/2n)+Δt[(S0−Sf)gΩ]1n+1where the index “1” indicates that there is only one control volume, in which the source term and the flux at the moving interface (F3/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 as (13) where 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 step (14) to 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).Positive Interface Moving DownstreamThe 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 (F3/2n+1), the variables at the interface (U3/2n+1) is needed, which is equal to the flow variables at the left of the moving interface at time level n+1 (UP). The flow variables at the left of the moving interface is related to the flow variables at the right of this interface (UF) as well as to wn+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 (uR+cR

Source link