Next Article in Journal
Bayesian Modeling of Flood Frequency Analysis in Bangladesh Using Hamiltonian Monte Carlo Techniques
Next Article in Special Issue
Scenario-Based Economic Impact Analysis for Bridge Closures Due to Flooding: A Case Study of North Gyeongsang Province, South Korea
Previous Article in Journal
Evaluating the Effects of Watershed Size on SWAT Calibration
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Hybrid Numerical Scheme of Preissmann Slot Model for Transient Mixed Flows

1
Department of Agricultural & Rural Engineering, Chungnam National University, Daejeon 34134, Korea
2
Climate Change Research Team, APEC Climate Center, Busan 48058, Korea
3
Department of Civil Engineering, University of Texas at Arlington, Arlington, TX 76019, USA
4
Software Development Center, K-water Convergence Institute, Deajeon 34045, Korea
*
Author to whom correspondence should be addressed.
Water 2018, 10(7), 899; https://doi.org/10.3390/w10070899
Submission received: 15 June 2018 / Revised: 3 July 2018 / Accepted: 4 July 2018 / Published: 6 July 2018
(This article belongs to the Special Issue Influence of the Urban Fabric on the Risks of Floods)

Abstract

:
The Preissmann slot model is one of the most widely used models to conceptualize both free-surface and pressurized flows in urban drainage systems. Despite its simplicity and wide range of applications, numerical solutions of the Preissmann slot model suffer from the spurious oscillations, especially when flow conditions switch from free-surface flow to pressurized flow. To overcome this problem, a new hybrid numerical flux solver of the Preissmann slot model is proposed herein, in which the upwind flux solver is combined with the centered flux solver. Numerical experiments are conducted for multiple flow conditions such as typical filling, pipe-filling, and transition-flow conditions. The numerical results indicate that the proposed scheme generally outperforms the conventional flux schemes for various hydraulic conditions and wave velocities. The proposed scheme should be useful to further enhance integrated urban water modeling in which transient mixed flow conditions significantly impact the simulation accuracy during extreme floods.

1. Introduction

Numerical modeling of sewer-pipe flow is of significant importance for a wide range of urban water problems such as urban inundation analysis, detecting the damage of a sewer pipe, or design of the sewer system. The type of flows in a sewer-pipe system is a free-surface flow but, during extreme flood events, it often becomes a partly free surface, partly pressurized flow (mixed flow) or fully pressurized flow. The numerical modeling of transient mixed flows is associated with infrastructure damage or operational problems. However, it is difficult to find a straightforward way because two different flow regimes should be considered in the same domain. When using two governing equations, it is difficult to track the mixed flow interfaces and switch smoothly between the two systems. Previous research [1,2] has shown that a compressible one-dimensional fluid equation with a linearized pressure law is an attractive alternative for modeling of transient mixed flows in pipe; however, a compressible fluid system is more complex to solve than a one-dimensional shallow water system usually used to represent open channel flows. The Preissmann slot model can avoid this difficulty and has been widely used to solve multiple flow conditions in sewer systems [3,4,5,6,7,8,9].
Unfortunately, the Preissmann slot model suffers from the spurious oscillations when the flow condition switches from open-channel flow to pressurized flow [10]. The conventional remedy for this problem is to reduce the impact of the acoustic-wave velocity by increasing the slot width [11,12]. However, this remedy can reduce the accuracy of the numerical solution because of the wider slot width and lower wave speed compared with the physical values. [13] presents an exact Riemann flux solver of the Preissmann slot model for highly transient mixed flows. However, the exact Riemann flux solver is not practical because it requires a large amount of computation for the iteration procedure. In addition, this method could not fully remove the spurious oscillations [6]. Ref. [3] proposed two methods based on the Roe’s flux scheme to prevent spurious oscillations. The first method is to introduce the artificial viscosity, which gradually increases with respect to the wave velocity. The second method uses the numerical filter to smooth variables such as water depth and momentum across the cells. Although both two methods can reduce the numerical oscillation within a certain level and stabilize the numerical results, oscillations can still be significant even with a comparatively moderate acoustic-wave speed of 100 m/s [14]. Leon [15] proposed a funnel-shaped slot to alleviate the abrupt change of wave velocity and the modified-shaped slot resulted in nonoscillatory solutions. However, this method is specifically tailored for a pipe with circular cross section. Malekpour and Karney [10] investigated the source of the spurious oscillation and proposed the modified Harten–Lax–van Leer (HLL) flux estimation based on the HLL flux proposed by Leon [15]. Their method removes the spurious oscillations even with acoustic-wave speed as high as 1000 m/s and produces accurate results for a wide range of cases. However, it requires tuning two additional parameters for each case, which may hinder the general application of this method.
The formulation of boundary conditions is also not a trivial issue in numerical modeling of sewer-pipe systems. The treatment of boundary conditions is strongly associated with the stability and accuracy of the numerical solution. Recently, ref. [16] presented a robust formulation that can represent a wide range of boundary conditions such as drop shafts, reservoirs, junctions, and dead ends. They considered common boundaries as junctions connected to an arbitrary number of pipes, and their method works well for a wide range of problems. However, equations for the implicit nonlinear system must be solved at each junction for each time step, which is not an easy task and may be an obstacle to applying this method.
In a general sense, the numerical flux solver for hyperbolic conservation laws can be categorized as either an upwind scheme or a centered (symmetric) scheme [17,18]. The upwind flux solvers, also called “Godunov-type schemes” or “approximated Riemann flux solvers,” are based on Riemann approaches or the equivalent and are suitable for capturing sharp water profiles such as shocks. Because of their superior performance, the upwind flux solver has been preferred over general hyperbolic conservation systems such as compressible flow, hydrodynamics systems, and shallow water systems, etc. [3,13,16,19,20,21]. However, for the Preissmann slot model, the monotonicity is not guaranteed by the widely used Riemann flux solvers because of the highly irregular shapes of pipes connected to the hypothetical slot, as presented in numerous previous works [5,10,22], whereas it is preserved for a regular-shaped flow domain. Conversely, the centered flux schemes are more diffusive than the upwind flux schemes, although they often over-smooth the water profile. The present study is based on the simple idea that a more diffusive flux solver may be required when flow is near the conduit roof [10,14].
The objective of this study is to develop a new hybrid numerical scheme of the Preissmann slot model and, through numerical experiments, demonstrate its applicability for a wide range of transient mixed-flow conditions. To achieve a more stable and accurate numerical scheme, the hybrid scheme combines both the upwind flux solver and the centered flux, although the proposed scheme is more than a combination of existing schemes. We also propose a modified formulation for boundary conditions, originally developed by [16,23], to explicitly compute the junction boundary. Numerical experiments are conducted to demonstrate the improved performance of the proposed method over conventional schemes for multiple flow conditions such as typical filling (Test 1), pipe-filling (Test 2), and transition-flow conditions (Test 3).

2. Numerical Model

2.1. Governing Equation

The hyperbolic conservative form of the continuity and momentum equations of free surface flow in uniform open channels and pipes can be written as follows:
U t + F x = S ,
where U , F , and S are vectors representing flow variables, fluxes, and source terms, respectively. These terms are written as follows:
U = [ A Q ] , F = [ Q Q 2 A + A p ¯ ρ ] , S = [ 0 g A ( S 0 S f ) ] ,
where A is the flow cross-sectional area, Q is the flow discharge, p ¯ is the average pressure over the flow cross-sectional area, ρ is the density of water, g is the acceleration due to gravity, S 0 is the slope of the pipe, and S f is the energy slope, which can be expressed by Manning’s formula as follows:
S f = n m 2 Q | Q | A 2 R 4 / 3 ,
where R is hydraulic radius and n m is Manning’s coefficient. Note that the friction surface on a slot is ignored following a basic assumption of Preissmann slot model [15].
The pressure head of pressurized flow is conceptualized as the open-channel flow depth in the Preissmann slot model. Because the acoustic-wave velocity is hypothetically equivalent to the velocity of a gravity wave, the slot width is T s = g A f / a 2 , where A f is the conduit cross-sectional area and a is the acoustic-wave velocity of the pipe.
Equation (1) can be discretized by the finite volume and forward Euler method as follows:
U i n + 1 = U i n Δ t Δ x ( F i + 1 / 2 n F i 1 / 2 n ) + Δ t S n ,
where the superscript n is the time-step index, the subscript i is the index for the ith computational cell, Δ t is the time step, Δ x is the length of a computational cell, and F i + 1 / 2 is the numerical flux between cells i and i + 1.
The numerical flux function F i + 1 / 2 ( U L , U R ) is estimated by using U L = U i and U R = U i + 1 in the spatially first-order finite-volume numerical scheme with piecewise constant variable in the cell, where U L and U R represent the left and right variable vectors at cell-interface i + 1/2. The values of U L and U R can be reconstructed for the higher-order spatial estimate. Here, we consider the total variation diminishing scheme with piecewise linear distribution of variables. The U L and U R are reconstructed as
U L = U i + 1 2 φ ( r i ) Δ U i 1 / 2 , U R = U i + 1 1 2 φ ( r i + 1 ) Δ U i + 1 / 2 ,
where
Δ U i 1 / 2 = U i U i 1 ,    Δ U i + 1 / 2 = U i + 1 U i ,    r i = Δ U i + 1 / 2 Δ U i 1 / 2 ,    r i + 1 = Δ U i + 3 / 2 Δ U i + 1 / 2 ,
and φ is a slope-limiter function. Among the variety of slope-limiter functions, the minmod limiter is used herein for stability. The minmod limiter is
φ ( r ) = max [ 0 , min ( r , 1 ) ] .
Note that the minmod limiter is the most diffusive and stable limiter.

2.2. Upwind Flux Solver

The Riemann flux solver for the Godunov-type numerical scheme is constructed based on an analysis of the eigenstructure of the system. Since an exact Riemann flux solver is impractical because of its high computational cost, an approximate Riemann flux solver should be used in practical applications [24]. In this study, we consider Roe’s flux and HLL flux schemes.
Roe’s flux scheme determines the approximate solution by solving a constant-coefficient linear system instead of the original nonlinear system. The system of Equation (1) is linearized as
U t + J ^ ( U L , U R ) U x = S ,
where J ^ is a Jacobian matrix at inter-cells and is assumed constant between two cells. The Riemann problem can then be solved as a linear hyperbolic system at each cell interface as
F i + 1 / 2 Roe = 1 2 ( F L F R ) | J ^ ( U L , U R ) | ( U R U L ) ,
J ^ ( U L , U R ) = R ( U L , U R ) | Λ ( U L , U R ) | R 1 ( U L , U R ) ,
where R is the eigenvector matrix and Λ is the diagonal corresponding eigenvalue matrix. The final stage of the algorithm is to find a suitable average interface state to determine J ^ . As long as the space of the numerical solution is smooth, any reasonable average would probably give the right results. However, when discontinuities or shocks are present in the solution, the estimate of the average interface states becomes very important to find the right solution. Here, we use the average interface state as follows:
A ¯ = A L A R ,   Q ¯ = A R Q L + A L Q R A R + A L ,   c ¯ = g I R I L A R A L
where c is the gravity wave velocity and I = A y .
In the HLL approach, the solution of the Riemann problem is approximated by an intermediate region U , which is the constant state separated from the left and right states U L and U R by two infinitely thin shock waves. The numerical flux F i + 1 / 2 is computed by the HLL flux solver as follows:
F i + 1 / 2 HLL = { F L if   S L > 0 , S R F L S L F R + S L S R ( U R U L ) S R S L if   S L 0   and   S R 0 , F R if   S R < 0 ,
where F L and F R represent left and right flux vectors at cell-interface i + 1/2, respectively, and S L and S R left and right wave speeds, which are calculated as follows [15]:
S L = u L M L , S R = u R M R
where u L and u R are left and right flow velocities, respectively, and M K ( M = L , R ) is
M K = { ( A p ¯ ρ A K p ¯ K ρ ) A A K ( A A K ) if   A > K c K otherwise
where the subscript “*” represents the star (intermediate) region and both A and p ¯ are the functions of flow depth of the star region ( y ). The variables of the star region can be estimated from the characteristic of the system. Ref. [1] suggests several different approaches for approximating y based on two shock-wave approximations, two rarefaction-wave approximations, and the linearization of the governing equations. For example, the estimate based on the linearized governing equation is
A = A R + A L 2 + A ¯ 2 c ¯ ( u L u R ) ,
where A ¯ = ( A R + A L ) / 2 and c ¯ = ( c R + c L ) / 2 . However, in our experience all three approaches mentioned above result in spurious oscillations when the flow switches from open-channel flow to pressurized flow, which is consistent with the observation of [4]. They investigated the effect of Equation (15) and concluded that the intermediate state could be adjusted to admit optimal numerical viscosity. They proposed following the adjusted estimate of star region flow depth as follows:
y = K a max [ y i N S , y i N S + 1 , , y i 1 , y i , y i + 1 , , y i + N S 1 , y i + N S ] ,
where K a and N S are parameters that must be tuned empirically. The problem of Equation (16) is that both parameters significantly affect the performance of the numerical scheme and should be calibrated for each case. Unlike the experiment, there usually is no reference data in an actual sewer network, which can significantly hinder the practical application of Equation (16). Therefore, the method proposed by [4] is not considered herein.

2.3. Centered Flux Solver

The main feature of the centered scheme lies in its simplicity. This scheme does not require explicit knowledge of the eigenstructure of the system, nor the availability of a Riemann solver. The classical centered schemes include the Lax–Wendroff (LW) scheme and the Lax–Friedrichs (LF) scheme. The two-step version of the Lax–Wendroff scheme is
F i + 1 / 2 LW = F ( Q i + 1 / 2 LW ) , Q i + 1 / 2 LW = 1 2 ( Q i n + Q i + 1 n ) 1 2 Δ t Δ x ( F ( Q i + 1 n ) F ( Q i n ) ) .
The Lax–Friedrichs scheme is
F i + 1 / 2 LF = 1 2 ( F ( Q i n ) + F ( Q i + 1 n ) ) 1 2 Δ x Δ t ( Q i + 1 n Q i n ) .
Note that the two-step version of LW scheme was derived by estimating intermediate variables at n + 1/2 from LF scheme. Therefore, it is a second-order method in temporal and different with the first-order forward Euler method. Although these two schemes are pioneering works in the centered scheme, they are not suitable for a modern finite-volume numerical scheme. The LW scheme does not preserve monotonicity and often results in spurious oscillations. The LW scheme, however, usually results in excessively diffusive numerical solutions. Vasconcelos et al. [25] pointed out that the LF scheme is not appropriate for flow regime transition problems because a huge difference in the wave velocities before and after the shock front results in very small Courant numbers in the free-surface portion of the flow, which in turn exacerbates the diffusive nature of the LF scheme and compromises the representation of the pipe-filling bore.
The first-order centered (FORCE) flux, proposed by Toro and Billett [26], is another option for a centered scheme. The FORCE flux is derived as the deterministic version of the staggered-grid random choice method, and the FORCE flux is an arithmetic mean of the LF and LW fluxes:
F i + 1 / 2 FO = 1 2 ( F i + 1 / 2 LW + F i + 1 / 2 LF ) .
The FORCE flux is less diffusive than the conventional LF scheme and it preserves the monotonicity. It was applied to Euler and molecular hydrodynamics equation, and the high-order extension of the FORCE scheme was also presented on two-dimensional unstructured meshes by Toro et al. [17].
Because the centered flux schemes are more diffusive than the upwind flux schemes, we expect the spurious oscillations to be reduced by the centered flux schemes. However, the diffusive flux schemes may affect the numerical accuracy whereas the upwind schemes usually do not result in spurious oscillations when the flow depth is lower than the conduit roof. Therefore, the following simple hybrid flux is proposed:
F i + 1 / 2 Hybrid = { F i + 1 / 2 FO if   ( y L D ) ( y R D ) < 0 , F i + 1 / 2 HLL otherwise .
where D is the height of pipe. The strategy of the hybrid scheme is straightforward: FORCE flux is implemented only if the flow switches from the open-channel flow to the pressurized flow; otherwise, HLL flux is implemented.

2.4. Boundary Computation

Various types of boundaries should be considered for a typical storm-sewer system. These include drop shafts, reservoirs, junctions, dead ends, and control gates, etc. Leon et al. [27] generalized the various types of boundary conditions as an N-way junction problem. In general, 2N + 1 variables are unknown for an N-way junction; namely, the piezometric depth yb and the flow discharge Qb at each pipe boundary, and the water depth yd at the junction pond. Therefore, 2N + 1 equations must be solved for an N-way junction. The first equation is obtained from the Riemann invariants between the pipe boundary and the neighboring cell as d u ± ( c / A ) d A = 0 . The Riemann invariants are solved as follows Leon et al. [27]:
u b j n + 1 u b j n ± ( c b j n + 1 + c j n ) A b j n + 1 A j n A b j n + 1 A j n = 0 ,
where the positive (negative) sign is used for inflowing (outflowing) pipes, the subscript b refers to the boundary, and j is the pipe index.
The second equation is usually given by the following energy equation:
y b n + 1 y p n + 1 + ( u b j n + 1 ) 2 2 g k j u b n + 1 | u b n + 1 | 2 g = 0 ,
where k is the local head-loss coefficient and the subscript p is the depth of the junction pond. When the inflowing pipe is decoupled from the junction, the different boundary conditions can be used as presented in Leon et al. [27].
The last equation is given by the mass-balance equation of the junction pond:
A p y p t = j N A b j u b j ,
where Ap is the vertical cross-sectional area of the junction pond. Leon et al. [27] used the backward Euler method to discretize Equation (24) as follows:
A p y p n + 1 y p n Δ t = j N ( A b j u b j ) n + 1
The combination of Equations (21), (22) and (24) can represent wide range of boundary conditions such as drop shafts, reservoirs, junctions, and dead ends, as has been verified through several test problems. However, this approach requires 2N + 1 nonlinear equations to be solved at a junction in an implicit manner. The convergence of the nonlinear system could be difficult to achieve for a large number of connected junctions. To avoid this difficulty, the above approach is modified in this study. The explicit form of Equation (25) is
A p y p n + 1 y p n Δ t = j N ( A b j u b j ) n
Note that Equation (25) is derived from Equation (23) by using the forward Euler method. This simple change makes the problem much easier. Equation (25) can be solved separately, so the boundary conditions at each pipe can be solved in separate manner.

3. Numerical Results

3.1. Test 1: Typical Filling Bore Problem

For Test 1, a typical filling bore problem with an analytical solution, introduced by Malekpour and Karney [10], is used to understand the cause of the numerical oscillations. Figure 1 shows the description of typical filling bore problem and details regarding the experiments process can be found in Malekpour and Karney [10]. With the condition y R   =   0.6   m and u R = 0 , the analytical solution can be obtained by using the traveling-wave method as S W = 10.077 m / s , u L = 4.044 m / s , and y L = 3.167 m .
The Courant–Friedrichs–Lewy (CFL) number is set to 0.5 and two acoustic-wave velocities ( a = 100 ,   1000   m / s ) are considered to estimate the performance of numerical schemes with relatively low and high values of wave velocity. The number of numerical cells is 400, as used in the study of Malekpour and Karney [10]. The values at intermediate region for HLL scheme are estimated by Equation (15) for all test problems in this paper.
Figure 2 compares the analytical solution with the numerical solutions computed by HLL, Roe, FORCE, and hybrid flux schemes using the spatially first-order piecewise constant variables in cells. Both upwind schemes, HLL and Roe, produce spurious oscillations near the shock-wave front whereas FORCE and hybrid schemes do not. In particular, the Roe scheme produces the most significant oscillations for both cases with a = 100 ,   1000   m / s , presumably because it is the least diffusive of the four flux schemes while preserving the sharp water profile near the shock-wave front. Conversely, the FORCE scheme produces a smoothed water profile compared with the other flux schemes. Over-smoothing of the water profile becomes more significant when a higher acoustic-wave velocity, a   =   1000   m / s , is applied. Smoothed hydraulic-head profiles are also found in the other three schemes with high acoustic-wave velocity near the shock-wave front, although the magnitude of over-smoothing is smaller than those of FORCE. This result implies that the high acoustic-wave velocity in the Preissmann slot model may cause not only the computational instability but also the reduced accuracy in a practical aspect. Overall, of the four flux schemes, the proposed hybrid scheme performs the best.

3.2. Test 2. Pipe-Filling Bore Experiment

For Test 2, a laboratory experiment performed by Vasconcelos and Wright [25] is used to compare the performance of the numerical schemes. Figure 3 shows the temporal variations of the observed and simulated water levels at the downstream tank. All numerical flux schemes reproduce the observed phenomena of the experiment except for the Roe flux scheme. The Roe flux scheme fails to reproduce this experiment because of the instability caused by spurious oscillations. Figure 4 compares the hydraulic-head profiles obtained from different flux schemes at 6 s. Apparently, the HLL scheme causes the most spurious oscillations of the three flux schemes. The spurious oscillation is more significant with a   =   300   m / s than with a   =   100   m / s . The FORCE scheme produces the most stable results but the profile is over-smoothed near pipe-filling bores as compared with the other two flux schemes. The hybrid scheme has the best result with a   =   100   m / s whereas it suffers from spurious oscillations with a   =   300   m / s . Note that the FORCE scheme, which is the most diffusive flux scheme of the three schemes, also produces spurious oscillations with a   =   300   m / s . The results of this study are comparable to Figure 7(A) in Vasconcelos et al. [14]. Their results, obtained by using a modified Roe scheme, are more oscillatory than the results presented herein. Although the spurious oscillations are not perfectly removed in the case of high acoustic-wave velocity, the hybrid scheme performs satisfactorily for all cases in Test 2.

3.3. Test 3: Transition Flows in Pipe Experiment

The laboratory experiment, conducted at the Hydraulics Laboratory of the University of Calabria by Trajkovic et al. [12], is used for Test 3. This experiment evaluates how rapid changes affect the opening or closing of the sluice gates. In this study, we simulate the type-A set of experiments of Trajkovic et al. [12]. The initial condition for the experiment is an inflow rate of 0.0013 m3/s. The downstream sluice gate is totally opened, and the upstream sluice gate is opened by 0.014 m. A transient flow is generated after the downstream sluice gate is rapidly closed. Thirty seconds after the gate closes, the gate is partially reopened by 0.008 m. The pressure head appears 4.6 and 0.6 m upstream from the downstream sluice gate. The CFL number is set to 0.5 and two acoustic-wave velocities ( a   =   100 ,   300   m / s ) are considered to estimate the performance of the numerical schemes with relatively low and high wave velocity. The number of numerical cells is 100.
Figure 5 shows the observed and simulated pressure head 4.6 and 0.6 m upstream from the downstream sluice gate. The results of the HLL and hybrid schemes match well with the observed data overall whereas the FORCE scheme does not reproduce the observed pressure heads, even with a relatively low acoustic-wave speed a   =   100   m / s . We hypothesize that the FORCE scheme is too diffusive to capture a filling bore propagation to a satisfactory level for this problem. The effect of the diffusive numerical flux is more significant when high acoustic-wave speed is adopted, which is also shown in Test 1. Figure 6 shows the hydraulic-head profile at 30 s. The HLL and hybrid scheme produce very similar profiles with a   =   100   m / s . However, the HLL scheme causes significant spurious oscillations with a   =   300   m / s at x = 7.5 m, whereas the hybrid scheme does not. For both cases, the FORCE scheme produces diffusive filling bore profiles. Although the acoustic-wave speed of 1000 m/s is also tested, the results are not presented here because all numerical schemes fail to simulate due to spurious oscillations. Note that the hybrid scheme proposed herein can, to a certain degree, alleviate the instability from spurious oscillations, but the very high acoustic-wave speed of over 1000 m/s still can cause numerical instability, depending on the situation.

4. Conclusions

We propose a hybrid flux solver of the Preissmann slot model and evaluate it for the improved numerical analysis of transient mixed flows in sewer-pipe systems. The new hybrid scheme combines the upwind flux solver and the centered flux solver to exploit the benefits of the two different approaches. A modified treatment scheme for general boundary conditions, such as junctions, which was originally developed by Leon et al. [27], is also proposed and applied.
Three numerical test problems are investigated to verify the proposed hybrid flux scheme. The proposed scheme is usually more stable than the conventional HLL scheme and less diffusive than the central flux schemes such as the Force scheme. For Test 1, which is a typical filling bore problem with an analytical solution in a rectangular pipe, the hybrid scheme performs the best of the four different flux schemes with high and low wave celerity (a = 100, 1000 m/s). Tests 2 and 3 are experimental simulations of transient mixed flows in a circular pipe connected between two water tanks. The hybrid scheme performs generally better than the other conventional flux schemes. In particular, the Roe scheme, which is known as the least diffusive scheme, often results in spurious oscillations and fails the computations.
In addition to improved numerical accuracy and stability, another advantage of the hybrid scheme is its simplicity of application compared with the other methods presented in Vasconcelos et al. [26] and Malekpour and Karney [10]. Regardless of the shape of the cross section, the hybrid method can be applied to any type of pipe without adjusting additional parameters. The hybrid method is usually stable with a wave velocity ranging from ~100 to ~300 m/s through various test simulations. There remains, however, a limitation in terms of spurious oscillations, especially with a very high value of the wave velocity, such as a = 1000 m/s for two test problems. It remains difficult to present a threshold for the wave velocity that causes a numerical instability because it depends strongly on the problem. A more rigorous analysis of the relation between numerical stability and wave velocity is left for future work. The proposed scheme would be useful to further enhance integrated urban water modeling in which transient mixed flow conditions impact the simulation accuracy, especially for extreme flooding. We are currently pursuing a demonstration of the proposed method for real-world applications in a large urban area, and the results will be reported in the near future.

Author Contributions

Conceptualization, H.A. and S.L.; Methodology, H.A.; Software, H.A.; Validation, H.A., S.L. and Y.K.; Formal Analysis, H.A.; Investigation, S.L.; Resources, S.J.N.; Data Curation, S.L.; Writing-Original Draft Preparation, H.A.; Writing-Review & Editing, S.L., S.J.N., Y.K., J.N.; Visualization, H.A.; Supervision, J.N.; Project Administration, H.A.; Funding Acquisition, H.A.

Funding

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (2017R1C1B2005750).

Acknowledgments

S.L. acknowledges the support from the APEC Climate Center.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bourdarias, C.; Gerbi, S. A Finite Volume Scheme for a Model Coupling Unsteady Flows in Open Channels and in Pipelines. J. Comput. Appl. Math. 2007, 209, 109–131. [Google Scholar] [CrossRef]
  2. Bourdarias, C.; Gerbi, S. A Conservative Model for Unsteady Flows in Deformable Closed Pipes and Its Implicit Second Order Finite Volume Discretisation. Comput. Fluids 2008, 37, 1225–1237. [Google Scholar] [CrossRef]
  3. Bertsch, R.; Glenis, V.; Kilsby, C. Urban Flood Simulation Using Synthetic Storm Drain Networks. Water 2017, 9, 925. [Google Scholar] [CrossRef]
  4. Djordjević, S.; Prodanović, D.; Walters, G.A. Simulation of Transcritical Flow in Pipe/Channel Networks. J. Hydraul. Eng. 2004, 130, 1167–1178. [Google Scholar] [CrossRef]
  5. Kerger, F.; Archambeau, P.; Erpicum, S.; Dewals, B.J.; Pirotton, M. An exact Riemann solver and a Godunov scheme for simulating highly transient mixed flows. J. Comput. Appl. Math. 2011, 235, 2030–2040. [Google Scholar] [CrossRef]
  6. Kerger, F.; Archambeau, P.; Erpicum, S.; Dewals, B.J.; Pirotton, M. A fast universal solver for 1D continuous and discontinuous steady flows in rivers and pipes. Int. J. Numer. Meth. Fluids 2011, 66, 38–48. [Google Scholar] [CrossRef] [Green Version]
  7. Leandro, J.; Chen, A.S.; Djordjević, S.; Savić, D.A. Comparison of 1D/1D and 1D/2D Coupled (Sewer/Surface) Hydraulic Models for Urban Flood Simulation. J. Hydraul. Eng. 2009, 135, 495–504. [Google Scholar] [CrossRef]
  8. Lee, S.; Nakagawa, H.; Kawaike, K.; Zhang, H. Urban inundation simulation considering road network and building configurations. J. Flood Risk Manag. 2015. [Google Scholar] [CrossRef]
  9. Noh, S.J.; Lee, S.; An, H.; Kawaike, K.; Nakagawa, H. Ensemble urban flood simulation in comparison with laboratory-scale experiments: Impact of interaction models for manhole, sewer pipe, and surface flow. Adv. Water Resour. 2016, 97, 25–37. [Google Scholar] [CrossRef]
  10. Malekpour, A.; Karney, B.W. Spurious Numerical Oscillations in the Preissmann Slot Method: Origin and Suppression. J. Hydraul. Eng. 2015, 142. [Google Scholar] [CrossRef]
  11. Capart, H.; Sillen, X.; Zech, Y. Numerical and experimental water transients in sewer pipes. J. Hydraul. Res. 1997, 35, 659–672. [Google Scholar] [CrossRef]
  12. Trajkovic, B.; Ivetic, M.; Calomino, F.; D’Ippolito, A. Investigation of transition from free surface to pressurized flow in a circular pipe. Water Sci. Technol. 1999, 39, 105–112. [Google Scholar] [CrossRef]
  13. Beljadid, A.; Mohammadian, A.; Kurganov, A. Well-balanced positivity preserving cell-vertex central-upwind scheme for shallow water flows. Comput. Fluids 2016, 136, 193–206. [Google Scholar] [CrossRef] [Green Version]
  14. Vasconcelos, J.G.; Wright, S.J.; Roe, P.L. Numerical oscillations in pipe-filling bore predictions by shock-capturing models. J. Hydraul. Eng. 2009, 135, 296–305. [Google Scholar] [CrossRef]
  15. Leon, A.S. Improved Modeling of Unsteady Free Surface, Pressurized and Mixed Flows in Storm-Sewer Systems. Ph.D. Thesis, University of Illinois at Urbana-Champaign, Urbana-Champaign, IL, USA, 2007. [Google Scholar]
  16. Cai, L.; Feng, J.H.; Xie, W.X. A CWENO-type central-upwind scheme for ideal MHD equations. Appl. Math. Comput. 2005, 168, 600–612. [Google Scholar] [CrossRef]
  17. Toro, E.F.; Hidalgo, A.; Dumbser, M. FORCE Schemes on Unstructured Meshes I: Conservative Hyperbolic Systems. J. Comput. Phys. 2009, 228, 3368–3389. [Google Scholar] [CrossRef]
  18. Vázquez-Cendón, E.; Hidalgo, A.; Navarro, P.G.; Cea, L. Numerical Methods for Hyperbolic Equations; CRC Press: London, UK, 2012. [Google Scholar]
  19. Causon, D.M.; Ingram, D.M.; Mingham, C.G.; Yang, G.; Pearson, R.V. Calculation of shallow water flows using a Cartesian cut cell approach. Adv. Water Resour. 2000, 23, 545–562. [Google Scholar] [CrossRef]
  20. Ducros, F.; Laporte, F.; Soulères, T.; Guinot, V.; Moinat, P.; Caruelle, B. High-Order Fluxes for Conservative Skew-Symmetric-like Schemes in Structured Meshes: Application to Compressible Flows. J. Comput. Phys. 2000, 161, 114–139. [Google Scholar] [CrossRef]
  21. Koren, B. Upwind discretization of the steady Navier–Stokes equations. Int. J. Numer. Methods Fluids 1990, 11, 99–117. [Google Scholar] [CrossRef]
  22. Vasconcelos Jose, G.; Wright Steven, J.; Roe Philip, L. Improved Simulation of Flow Regime Transition in Sewers: Two-Component Pressure Approach. J. Hydraul. Eng. 2006, 132, 553–562. [Google Scholar] [CrossRef]
  23. León, A.S.; Liu, X.; Ghidaoui, M.S.; Schmidt, A.R.; García, M.H. Junction and drop-shaft boundary conditions for modeling free-surface, pressurized, and mixed free-surface pressurized transient flows. J. Hydraul. Eng. 2010, 136, 705–715. [Google Scholar] [CrossRef]
  24. Toro, E.F. Shock-Capturing Methods for Free-Surface Shallow Flows; John Wiley: New York, NY, USA, 2001. [Google Scholar]
  25. Vasconcelos, J.G.; Wright, S.J. Experimental Investigation of Surges in a Stormwater Storage Tunnel. J. Hydraul. Eng. 2005, 131, 853–861. [Google Scholar] [CrossRef]
  26. Toro, E.F.; Billett, S.J. Centred TVD schemes for hyperbolic conservation laws. IMA J. Numer. Anal. 2000, 20, 47–79. [Google Scholar] [CrossRef]
  27. Leon, A.S.; Ghidaoui, M.S.; Schmidt, A.R.; Garcia, M.H. A robust two-equation model for transient-mixed flows. J. Hydraul. Res. 2010, 48, 44–56. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Description of typical filling bore problem (Test 1).
Figure 1. Description of typical filling bore problem (Test 1).
Water 10 00899 g001
Figure 2. Comparison of results computed by different flux schemes for a typical filling bore problem (Test 1) with (a) a = 100 m/s; (b) a = 1000 m/s.
Figure 2. Comparison of results computed by different flux schemes for a typical filling bore problem (Test 1) with (a) a = 100 m/s; (b) a = 1000 m/s.
Water 10 00899 g002aWater 10 00899 g002b
Figure 3. Temporal variations of water level at downstream tank (Test 2) with (a) a = 100 m/s; (b) a = 300 m/s.
Figure 3. Temporal variations of water level at downstream tank (Test 2) with (a) a = 100 m/s; (b) a = 300 m/s.
Water 10 00899 g003aWater 10 00899 g003b
Figure 4. Hydraulic-head profiles simulated by different flux schemes at 6 s (Test 2) with (a) a = 100 m/s; (b) a = 300 m/s.
Figure 4. Hydraulic-head profiles simulated by different flux schemes at 6 s (Test 2) with (a) a = 100 m/s; (b) a = 300 m/s.
Water 10 00899 g004
Figure 5. Observed and simulated pressure heads at 4.6 m (P3) and 0.6 m (P7) from the downstream sluice gate (Test 3) with (a) a = 100 m/s; (b) a = 300 m/s.
Figure 5. Observed and simulated pressure heads at 4.6 m (P3) and 0.6 m (P7) from the downstream sluice gate (Test 3) with (a) a = 100 m/s; (b) a = 300 m/s.
Water 10 00899 g005
Figure 6. Hydraulic-head profiles simulated by different flux schemes at 30 s (Test 3). The two black lines represent the location of the sewer pipe with (a) a = 100 m/s; (b) a = 300 m/s.
Figure 6. Hydraulic-head profiles simulated by different flux schemes at 30 s (Test 3). The two black lines represent the location of the sewer pipe with (a) a = 100 m/s; (b) a = 300 m/s.
Water 10 00899 g006

Share and Cite

MDPI and ACS Style

An, H.; Lee, S.; Noh, S.J.; Kim, Y.; Noh, J. Hybrid Numerical Scheme of Preissmann Slot Model for Transient Mixed Flows. Water 2018, 10, 899. https://doi.org/10.3390/w10070899

AMA Style

An H, Lee S, Noh SJ, Kim Y, Noh J. Hybrid Numerical Scheme of Preissmann Slot Model for Transient Mixed Flows. Water. 2018; 10(7):899. https://doi.org/10.3390/w10070899

Chicago/Turabian Style

An, Hyunuk, Seungsoo Lee, Seong Jin Noh, Yeonsu Kim, and Jaekyoung Noh. 2018. "Hybrid Numerical Scheme of Preissmann Slot Model for Transient Mixed Flows" Water 10, no. 7: 899. https://doi.org/10.3390/w10070899

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop