Next Article in Journal
Decision-Making Tools to Manage the Microbiology of Drinking Water Distribution Systems
Next Article in Special Issue
Hydraulic Simulation and Analysis of an Urban Center’s Aqueducts Using Emergency Scenarios for Network Operation: The Case of Thessaloniki City in Greece
Previous Article in Journal
Variability of Trends in Precipitation across the Amazon River Basin Determined from the CHIRPS Precipitation Product and from Station Records
Previous Article in Special Issue
A Stochastic Model to Predict Flow, Nutrient and Temperature Changes in a Sewer under Water Conservation Scenarios
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Suppress Numerical Oscillations in Transient Mixed Flow Simulations with a Modified HLL Solver

State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China
*
Author to whom correspondence should be addressed.
Water 2020, 12(5), 1245; https://doi.org/10.3390/w12051245
Submission received: 21 February 2020 / Revised: 20 April 2020 / Accepted: 21 April 2020 / Published: 27 April 2020
(This article belongs to the Special Issue Advances in Modeling and Management of Urban Water Networks)

Abstract

:
Transition between free-surface and pressurized flows is a crucial phenomenon in many hydraulic systems. During simulation of such phenomenon, severe numerical oscillations may appear behind filling-bores, causing unphysical pressure variations and computation failure. This paper reviews existing oscillation-suppressing methods, while only one of them can obtain a stable result under a realistic acoustic wave speed. We derive a new oscillation-suppressing method with first-order accuracy. This simple method contains two parameters, Pa and Pb, and their values can be determined easily. It can sufficiently suppress numerical oscillations under an acoustic wave speed of 1000 ms−1. Good agreement is found between simulation results and analytical results or experimental data. This paper can help readers to choose an appropriate oscillation-suppressing method for numerical simulations of flow regime transition under a realistic acoustic wave speed.

1. Introduction

In water conveyance systems, water flows under free-surface flow condition or pressurized flow condition. Under certain circumstances, transition between the two flow regimes may occur (i.e., the flow regime transition phenomenon). Following the transition, force exerting on structures changes violently and causes structural damage [1,2,3]. Numerical simulation of flow regime transition can provide substantial information for the design and management of river-crossing bridges, tunnels, conducts and culverts [4,5,6,7,8,9,10].
The complexity of flow regime transition lies in the presence of free-surface and pressurized flows, which are governed by different equations. This problem can be avoided by adopting one set of governing equations for two flow regimes. Based on this idea, the Preissmann slot model (PSM) is proposed [11]; it was adopted by many researchers and commercial software packages [12,13,14,15,16,17,18]. The strong gradient in piezometric head at the interface between two flow regimes forms a discontinuity in flow. Finite volume methods can capture the discontinuity in the flow implicitly, which makes them popular in simulating flow regime transition [19].
Despite of all the fine properties that finite volume methods have, the numerical oscillations in a flow regime transition simulation have troubled many hydraulic engineers [12]. These numerical oscillations have the same origin with “post-shock oscillations” in gas dynamics [20,21,22,23]. In analytical results, the thickness of the filling-bore is infinitely small, and the flow states at the two sides of a filling-bore satisfy the Rankine–Hugoniot condition. In numerical simulations, a filling-bore spreads over several computational cells, and the flow states at the two adjacent cells do not satisfy the Rankine–Hugoniot condition. This causes trivial discrepancies in the mass and momentum fluxes, which are amplified in simulation results due to the large acoustic wave speed. High-order finite volume methods cause more numerical oscillations because of low dissipation away from the shocks [21]. First-order upwind finite volume methods fail to prevent numerical oscillations without compromising the representation of the filling-bore [23,24,25].
A lot of effort has been spent to obtain a stable and accurate result of flow regime transition [12,23,24]. In this paper, some existing oscillation-suppressing methods are tested on a benchmark model. Considering the lack of an efficient and simple method, the authors derive a new method, which can suppress numerical oscillations and capture the filling-bore nicely under a high acoustic wave speed. The structure of this paper is as follows: Section 2 introduces the governing equations and the discretization method. Section 3 reviews the existing oscillation-suppressing methods, and their effects are evaluated on the benchmark model that was adopted by Malekpour and Karney [24]. Section 4 proposes a new and simple modified HLL solver to suppress numerical oscillations. Its accuracy and robustness are tested against the analytical results and experiment data in Section 5. Conclusions are drawn in the last section.

2. Governing Equations and Discretization Method

The PSM places an infinitely high narrow slot on the top of the conduct, so that it becomes an open-channel with a composite cross-section. The water depth in the slot represents the piezometric head of the pressurized flow inside the original conduct, as shown in Figure 1. The slot width needs to be very small so that the gravity wave speed inside it is identical to the acoustic wave speed.
Under the framework of PSM, the governing equations of one-dimensional flow regime transition with uniform cross-sections can be written as [26]:
U t + F ( U ) x = S ( U )
U = ( A Q ) , F ( U ) = ( Q Q 2 A + g I ( h ) ) , S ( U ) = ( 0 g A ( S b S f ) )
I ( h ) = 0 h ( h ξ ) l ( x , ξ ) d ξ
where Q is the volume flow rate, A is the wetted area, g is the acceleration of gravity, h is the water depth, I is the moment of inertia, l is the cross-sectional width and l(x, h) = b(x) the water surface width. The external force acting on the flow is accounted in the source term, where Sb is the bed slope and Sf is the friction slope, which can be computed using the Manning relation:
S f = n 2 u | u | R 4 / 3
where u is the flow velocity, n is the Manning coefficient and R is the hydraulic radius. For a rectangular cross-section with a slot on its top, the parameters b, A, I and wave speed c can be expressed as functions of h:
b ( h ) = { B h H B s l h > H
A ( h ) = { B h h H B H + B s l ( h H ) h > H
I ( h ) = { 0.5 B h 2 h H B H ( h 0.5 H ) + 0.5 B s l ( h H ) 2 h > H
c ( h ) = { g h h H g A B s l h > H
where H and B are the cross-sectional height and width, and Bsl is slot width. In order to make the gravity wave speed inside the slot equal to the acoustic wave speed a, the slot width Bsl = gAf a−2, where Af is the full cross-sectional area of the conduct [27]. Using the Godunov-type finite volume methods with first-order accuracy and assuming a piecewise constant data construction, the governing equations are discretized as
U i * = U i n Δ t i Δ x i ( F i + 1 / 2 F i 1 / 2 )
U i n + 1 = U i * + Δ t i S ( U i * )

3. Review of Current Oscillation-Suppressing Methods

The benchmark model was proposed by Malekpour and Karney [24], it consists of a conduct with square-unit cross-sections that are connected to a reservoir at the upstream end. The acoustic wave speed is 1000 ms−1 and the slot width is 9.8 × 10−6 m. Under the initial condition, 0.6 m-deep stagnant water is in the conduct while the water level inside the reservoir is constantly 4 m, as shown in Figure 2.
Under initial condition, flow states UL in the reservoir and UR in the conduct are discontinuous:
[ h L u L ] = [ 4   m 0   ms 1 ] , [ h R u R ] = [ 0.6   m 0   ms 1 ]
Since hL is larger than hR, a shock wave (filling-bore) belonging to the second characteristic field is formed at the conduct inlet and propagates downstream. Consider 1–1 as a cross-section in the reservoir and 2–2 as a cross-section at conduct inlet: flow sates at 1–1 and 2–2 are U1 and U2, respectively. Then U2 can be obtained by solving the following equations iteratively [24]:
h 1 = h 2 + u 2 2 2 g
u 2 = u R + [ g I ( A R ) g I ( A 2 ) ] ( A R A 2 ) A R A 2
where U1 = UL by assuming that flow velocity inside the reservoir is negligible. In this case, h2 and u2 are 3.167 m and 4.0334 ms−1, a ghost cell is set at the upstream boundary adopting h2 and u2. Since U2 is connected to UR through a right shock, the flow states inside the conduct ultimately will take place by U2. The propagation speed of the filling-bore is given by the Rankine–Hugoniot condition, which is 10.067 ms−1. Then we can construct the analytical result in the benchmark model at t0:
h ( x ) = { 3.167   m x < 10.067 t 0 0.6   m x > 10.067 t 0 , u ( x ) = { 4.0334   ms 1 x < 10.067 t 0 0   ms 1 x > 10.067 t 0
As customary, x denotes the distance to the conduct inlet. In a numerical simulation, the size of each computational cell is 1 m, the time step is 0.008 s and the Courant number is 0.8. When the acoustic wave speed adopted in the simulation exceeds 100 ms−1, the magnitude of the numerical oscillations become so large that the simulated piezometric head become negative, and the simulation will not proceed [24]. In the remaining part of this section, the readers will see that only one method can get a satisfactory result under a high acoustic wave speed, while its performance rely on two parameters which must be well tuned. This shows the importance of devising an alternative method, which is stable and convenient.

3.1. Numerical Filtering Method

In this method, the exact Riemann solver is adopted to solve the Riemann problem at each cell boundary. Flow states Ui+1/2 at xi+1/2 satisfy the following equations:
u i + 1 / 2 = { u i g [ I ( A i ) I ( A i + 1 / 2 ) ] ( A i A i + 1 / 2 ) A i A i + 1 / 2 A i + 1 / 2 > A i u i + 0 A L g α b d α 0 A i + 1 / 2 g α b d α A i + 1 / 2 < A i
u i + 1 / 2 = { u i + 1 + g [ I ( A i + 1 ) I ( A i + 1 / 2 ) ] ( A i + 1 A i + 1 / 2 ) A i + 1 A i + 1 / 2 A i + 1 / 2 > A i + 1 u i + 1 0 A i + 1 g α b d α + 0 A i + 1 / 2 g α b d α A i + 1 / 2 < A i + 1
Equations (13) and (14) can be solved iteratively, then the wave structure in the Riemann problem can be determined and utilized to compute the flux; see Kerger et al. [26] for detail. Although the exact solver can obtain the complete wave structure, serious numerical oscillations appear in the simulation result. Vasconcelos et al. [23] proposed to suppress the numerical oscillations by averaging the flow states among the three conjunct cells at each time step:
U i n + 1 = ( 1 2 ε ) U i n + 1 + ε ( U i 1 n + 1 + U i + 1 n + 1 )
The authors suggest ε to be between 0.025 and 0.050. This method will increase the spreading length of the filling-bore front and remove any physical oscillations that appear in the solution. The simulation result of this method using ε = 0.04 is drawn in Figure 3.

3.2. Hybrid Flux Method

The hybrid flux method uses two types of numerical fluxes alternatively. The first one is based on the Roe solver [28] and the LxF scheme [23]. Define λ i + 1 / 2 j as eigenvalues of the linearized Jacobian matrix, R i + 1 / 2 j as the corresponding right eigenvectors and α i + 1 / 2 j as the wave strengths across the cell boundary:
λ i + 1 / 2 1 = Q ˜ i + 1 / 2 A ˜ i + 1 / 2 c ˜ i + 1 / 2 , λ i + 1 / 2 2 = Q ˜ i + 1 / 2 A ˜ i + 1 / 2 + c ˜ i + 1 / 2
R i + 1 / 2 1 = [ 1 , λ i + 1 / 2 1 ] T , R i + 1 / 2 2 = [ 1 , λ i + 1 / 2 2 ] T
α i + 1 / 2 1 = λ i + 1 / 2 2 ( A i + 1 A i ) ( Q i + 1 Q i ) λ i + 1 / 2 2 λ i + 1 / 2 1 , α i + 1 / 2 2 = ( Q i + 1 Q i ) λ i + 1 / 2 1 ( A i + 1 A i ) λ i + 1 / 2 2 λ i + 1 / 2 1
where Q ˜ i + 1 / 2 , A ˜ i + 1 / 2 and c ˜ i + 1 / 2 are Roe averages:
A ˜ i + 1 / 2 = ( A i A i + 1 ) 1 / 2
Q ˜ i + 1 / 2 = Q i + 1 ( A i ) 1 / 2 + Q i ( A i + 1 ) 1 / 2 ( A i ) 1 / 2 + ( A i + 1 ) 1 / 2
c ˜ i + 1 / 2 = { [ g I ( A i + 1 ) I ( A i ) A i + 1 A i ] 1 / 2 A i + 1 A i [ g A i + A i + 1 b i + b i + 1 ] 1 / 2 A i + 1 = A i
Then the fluxes obtained by the Roe solver are written as
F i + 1 / 2 Roe = 1 2 [ F ( U i + 1 n ) + F ( U i n ) ] 1 2 j = 1 2 | λ i + 1 / 2 j | α i + 1 / 2 j R i + 1 / 2 j
The Roe solver is known to be vulnerable to numerical oscillations [21,29], while the LxF scheme is robust against numerical oscillations but causes too much diffusion of the filling-bore:
F i + 1 / 2 LxF = 1 2 [ F ( U i + 1 n ) + F ( U i n ) ] Δ x 2 Δ t ( U i + 1 U i )
Compare Equations (22) with (23): The difference between the Roe solver and LxF scheme is the choice of eigenvalues. Vasconcelos et al. [23] proposed a new method to determine the eigenvalues:
| λ i + 1 / 2 1 , 2 | = min [ Δ x Δ t , | Q ˜ i + 1 / 2 A ˜ i + 1 / 2 c ˜ i + 1 / 2 | + L ( Δ c i + 1 / 2 ) L Δ x Δ t ]
Δ c i + 1 / 2 = | c ˜ i + 1 c ˜ i | max ( | c ˜ i + 1 c ˜ i | ) i = 1 N 1
This method is referred to as the hybrid one from here on further. At the cell boundary between the free-surface and pressurized flows, Δci+1/2 = 1, as L changes from 0 to 1, the eigenvalues switches from those obtained by the Roe solver to those adopted in the LxF scheme. At the other cell boundaries, Δci+1/2 is approximately 0, thus the eigenvalues in Equation (24) remain close to those obtained by the Roe solver. In this way, numerical viscosity is added at the cell boundary where flow condition transition happens, and its amount increases with L. The simulation results of the hybrid one with L = 0.6 are drawn in Figure 4. This method overestimates the spreading length of the filling-bore, and it fails to suppress the numerical oscillations under a high acoustic wave speed.
Hyunuk et al. [12] chose different fluxes; they used the FORCE scheme [30] at the cell boundary between the free-surface and pressurized flows, and the HLL solver [30] elsewhere:
F i + 1 / 2 hybrid = { F i + 1 / 2 FORCE ( h L H L ) ( h R H R ) < 0 F i + 1 / 2 HLL otherwise
This method is referred to as hybrid two from here on further. The HLL solver starts by estimating the largest wave speed SwR and smallest wave speed SwL in the Riemann problem. The fluxes are computed based on the signs of SwL and SwR:
F i + 1 / 2 HLL = { F ( U i ) S wL > 0 S wR F ( U i ) S wL F ( U i + 1 ) + S wR S wL ( U i + 1 U i ) S wR S wL S wL 0 and S wR 0 F ( U i + 1 ) S wR < 0
The choices of SwL and SwR follow Toro [31]:
S w L = u i Ω i , S w R = u i + 1 + Ω i + 1
Ω K ( K = i , i + 1 ) = { g [ I ( A i + 1 / 2 ) I ( A K ) ] A i + 1 / 2 A K ( A i + 1 / 2 A K ) A i + 1 / 2 > A K c K A i + 1 / 2 A K
where Ai+1/2 is an estimate of the wetted area at xi+1/2; we adopt the one proposed by Leno et al. [32], which admits the minimum amount of numerical viscosity:
A i + 1 / 2 = A i + A i + 1 2 ( 1 + u i u i + 1 c i + c i + 1 )
The FORCE flux can be written as the algebraic average of the LxF scheme and Lax–Wendorff scheme:
F i + 1 / 2 LW = F ( U i + 1 / 2 LW ) , U i + 1 / 2 LW = 1 2 ( U i n + U i + 1 n ) Δ x 2 Δ t [ F ( U i + 1 n ) F ( U i n ) ]
F i + 1 / 2 FORCE = 1 2 ( F i + 1 / 2 LW + F i + 1 / 2 LxF )
The FORCE scheme is a centred scheme; thus, it is robust against numerical oscillations. It is less diffusive than the LxF scheme, which can reduce the over-smearing at strong gradients [33]. The simulation results of hybrid two are depicted in Figure 5.
The hybrid flux methods have added enough numerical viscosity at the cell boundary between the two flow regimes but, still, serious numerical oscillations appear in the simulation results. The reason is that the numerical oscillations appear as soon as the flow regime transition happens, and while the two methods start to add numerical viscosity one time-step falls behind of it. If one method can foresee the happening of the flow regime transition, and admits numerical viscosity ahead of it, it will achieve a more stable result.

3.3. Modified HLL Solver

Malekpour and Karney [24] pointed out that the amount of numerical viscosity in the HLL fluxes increases with the magnitude of SwL and SwR. In fact, the HLL fluxes equal the LxF fluxes when |SwL| and |SwR| equal Δxt. To increase the amount of numerical viscosity, they proposed a modified HLL solver (referred to as the M_HLL solver). In the M_HLL solver, Ai+1/2 in Equation (29) is computed according to a reference depth hG:
A i + 1 / 2 = A ( h G ) , h G = K a max ( h i N S , h i N S + 1 , , h i 1 , h i , h i + 1 , , h i + N S 1 , h i + N S )
The depth hG is defined as Ka, multiplying the highest piezometric head in the 2NS+1 cells surrounding cell i, while Ka > 1 and NS ≥ 3. The solution of Equation (33) produces a larger magnitude of SwL and SwR; thus, increasing the numerical viscosity before the flow regime transition happens. The simulation results of M_HLL with Ka = 1.4 and NS = 5 are drawn in Figure 6.
The M_HLL solver can suppress numerical oscillations in the benchmark model, and it only slightly increases the spreading length at the filling-bore front. However, the values of Ka and NS can affect the diffusion and accuracy of the M_HLL solver to a great extent; see Figure 7.
Thus, the values of Ka and NS must be well-tuned. Meanwhile, the way to determine hG makes the HLL solver hard to use in parallelized computation. In the next section, the authors present an alternative method, which is equally efficient as the M_HLL solver.

4. A New Modified HLL Solver

In this section, we present a new modified HLL solver (referred to as the P_HLL solver) to suppress numerical oscillations. In the P_HLL solver, the solution to Ai+1/2 depends on the water depths at cell i and i + 1. When hi and hi+1 are below PbH, Ai+1/2 is computed using Equation (30) to admit the minimum amount of numerical viscosity, otherwise Ai+1/2 is computed according to a constant depth PaH:
A i + 1 / 2 = A ( P a H ) , h i > P b H or h i + 1 > P b H
Pb must be smaller than one, and a preferable value is between 0.6 and 0.8. PaH must be larger than the piezometric head peak during the transition to admit enough numerical viscosity.
To illustrate the effect of Pa and Pb, we study the Riemann problem at xi+1/2 in the benchmark model. Suppose hi and hi+1 is 3.167 m and 0.6 m, respectively, and ui and ui+1 is 4.0334 ms−1 and 0 ms−1, respectively. The solution of Equation (30) lies between Ai and Ai+1, and after substituting it into Equation (28), one will get SwL (noted as SwL1) as the speed of the left rarefaction wave head, and SwR (noted as SwR1) as the speed of right shock wave:
S wL 1 = u i c i , S wR 1 = u i + 1 + g [ I ( A i + 1 / 2 ) I ( A i + 1 ) ] A i + 1 / 2 A i + 1 ( A i + 1 / 2 A i + 1 )
A sketch of two waves in the Riemann problem is shown in Figure 8.
Because cell i is under a pressurized flow condition, it is easy to see that |SwL1| nearly equals the acoustic wave speed. The entropy condition of right shock wave requires
u i + 1 / 2 + c i + 1 / 2 > S wR 1 > u i + 1 + c i + 1
The magnitude of SwR1 equals the propagation speed of the filling-bore, which is 10.067 ms−1, and it is much smaller than the acoustic wave speed.
The solution of Equation (34) is unconditionally larger than Ai and Ai+1, which produces two shock waves in the Riemann problem. Consequently, SwR2 is the speed of the right shock wave, and SwL2 is the speed of the left shock wave; see Figure 9.
The entropy condition of the left shock wave requires
u i + 1 / 2 c i + 1 / 2 < S wL 2 < u i c i
Since cell i is under a pressurized flow condition, ui << ci; thus, |SwL2| > |SwL1| and they are both close to the acoustic wave speed. The speed of the right shock wave increases with Ai+1/2; thus, SwR2 > SwR1. This larger magnitude of SwR admits more mass and momentum into cell i + 1 before it becomes pressurized. The loci of Ui+1 simulated by HLL and P_HLL (Pa = 5, Pb = 0.7) are drawn in Figure 10.
Vertexes appear in the locus simulated by HLL after cell i + 1 becomes pressurized, which represent numerical oscillations in the simulation result. In contrast, P_HLL achieves a smooth and stable transition between the free-surface flow and pressurized flow. The locus simulated by P_HLL converges to a 3.159 m depth and 4.033 ms−1 velocity, which is very close to the flow states at the entrance of the conduct. The discrepancy in water depth is more pronounced due to the small value of the slot width. Therefore, P_HLL preserves the conservation in mass and momentum.
A larger magnitude of SwR2 admits more mass and momentum into cell i + 1 before it is pressurized; thus, it increases the diffusion of the filling-bore.The magnitude of SwL2 and SwR2 are related to the value of Ai+1/2, which consequently depends on the value of Pa; see Figure 11 and Figure 12:
Figure 11 shows that SwR2 increases with Pa, while|SwL2| barely changes with Pa and stays close to the acoustic wave speed. However, SwR2 does not increase significantly when Pa changes from 1 to 10, which denotes that the diffusion of the filling-bore does not increase significantly when Pa changes from 1 to 10. The simulation results using the P_HLL solver with Pb = 0.7 and different values of Pa are drawn in Figure 13.
Although Pa = 10 produces a more diffusive filling-bore, the spreading length of the filling-bore does not increase significantly compared to that using Pa = 5. During a realistic transition phenomenon, the piezometric head peak seldom exceeds 10 times the cross-sectional depth. Therefore, one can always start by adopting a large Pa (for example 10) in the P_HLL solver and do not worry about significantly compromising the representation of the filling-bore.
Compared to Pa, the value of Pb has a more significant effect on the numerical oscillations, for it determines the threshold depth where numerical viscosity starts to increase. Pb must be smaller than one so that the numerical viscosity is added before the flow regime transition happens. A smaller Pb leads to more stable result, but it may cause more diffusion. The simulation results using Pa = 5 and Pb = 0.7 or 0.9 are drawn in Figure 14.
When Pb = 0.9, a smooth transition between the free-surface and pressurized flows cannot be guaranteed. Therefore, we suggest a Pb between 0.6 and 0.8 to avoid causing two much diffusion of the filling-bore. This is also supported by the numerical tests in the next section.
In conclusion, at a free-surface cell, P_HLL admits numerical viscosity once the water depth is higher than a threshold value PbH. Thus, a smooth transition from the free-surface flow to pressurized flow can be obtained. Meanwhile, P_HLL causes diffusion of the filling-bore. In realistic applications, a Pa of 10 is large enough to suppress numerical oscillations without significantly increasing the spreading-length of the filling-bore. The value of Pb is suggested to be between 0.6 and 0.8.

5. Numerical Tests

5.1. Two Filling-Bores

This test evaluates the accuracy of P_HLL under the presence of two filling-bores. It consists of a 200 m-long conduct with square-unit cross-sections and two reservoirs connected to it at the upstream end and downstream end; the acoustic wave speed is 1000 ms−1. At initial conditions, water in the conduct is stationary with a depth of 0.6 m, while water depth at the upstream and downstream reservoir is 4 m and 3 m, respectively. The model set up is sketched in Figure 15.
At t = 0 s, the conduct inlet and outlet are opened immediately, forming two filling-bores that propagate in the opposite direction. The upstream filling-bore is identical with that in the benchmark model. The downstream filling-bore can be derived analogously; its propagation speed is 8.429 ms−1, and the flow states behind it are 2.42 m and −3.3717 ms−1. Boundary conditions adopt ghost cells at the conduct inlet and outlet. Before the two filling-bores collide with each other, the analytical result at t0 is
h ( x ) = { 3.167   m x < 10.067 t 0 0.6   m 10.067 t 0 < x < 200 8.429 t 0 2.42   m 200 8.429 t 0 < x u ( x ) = { 4.0334   ms 1 x < 10.067 t 0 0   ms 1 10.067 t 0 < x < 200 8.429 t 0 3.3717   ms 1 200 8.429 t 0 < x
In P_HLL, we choose Pa = 5 and Pb = 0.8, optionally. In M_HLL, we choose Ka = 1.4 and NS = 5 as suggested by Malekpour and Karney. The computational cell is 1 m, time step is 0.0008 s and the Courant number is 0.8. The simulation results in the two tests at t = 6 s are drawn in Figure 16. In this paper, an error indicator based on the L2-norm [34] is used to evaluate the accuracy of P_HLL and M_HLL. In the following equation, φi stands for the simulation result at cell i, while φref stands for the analytical result.
L 2 = ( 1 N i = 1 N ( φ i φ ref ) 2 ) 0.5
The L2 in the piezometric head of P_HLL and M_HLL are 0.2963 and 0.2913, respectively; the L2 in the velocity of P_HLL and M_HLL are 0.2879 and 0.2873, respectively. In P_HLL, the spreading length of the right filling-bore is slightly longer than that in M_HLL. This denotes that P_HLL is more diffusive than M_HLL in there. At the same time, P_HLL has eliminated some minor numerical oscillations while M_HLL does not. Both solvers are very robust and stable.

5.2. Negative Pressure Flow

In this test, P_HLL and M_HLL are adopted to simulate a water hammer phenomenon in a 600 m-long circular pipe with a 0.5 m diameter and acoustic wave speed is 1200 ms−1. The pipe is horizontal and frictionless; a steady flow rate of 0.477 m3s−1 is initially in it. It connects to a reservoir at the downstream end, and water depth inside it is 45 m; see Figure 17.
At t = 0 s, the inflow rate is decreased to 0.4 m3s−1, which triggers a water hammer phenomenon; the water hammer pressure is 48.05 m according to Kerger et al. [26]. In P_HLL, the values of Pa and Pb are 100 and 0.8. In M_HLL, the values of Ka and NS are 1.2 and 12, as suggested by Malekpour and Karney. The computational cell is 1.2 m, the time-step is 0.0008 s and Courant number is 0.8. A ghost cell is set at the upstream boundary, and flow rate in it is constantly 0.4 m3s−1, while the piezometric head adopts the transmissive condition. Another ghost cell is set at the downstream boundary in which the piezometric head is constantly 45 m and the flow rate adopts the transmissive condition.
The L2 indicator is adopted to evaluate the accuracy of the two solvers; it is defined as
L 2 = ( 1 N i = 1 N ( φ i φ ref ) 2 ) 0.5
where Nt is the number of time step, φi is the simulation result at the midpoint of the pipe, and φref is the analytical result, which is given as
h i = { 45   m n < i Δ t < n + 0.25 3.05   m n + 0.25 < i Δ t < n + 0.75 45   m n + 0.75 < i Δ t < n + 1.25 93.05   m n + 1.25 < i Δ t < n + 1.75 45   m n + 1.75 < i Δ t < n + 2 , n = 0 , 1 , 2 , 3 u i = { 2.4293   ms 1 n < i Δ t < n + 0.25 2.0377   ms 1 n + 0.25 < i Δ t < n + 0.75 1.6461   ms 1 n + 0.75 < i Δ t < n + 1.25 2.0377   ms 1 n + 1.25 < i Δ t < n + 1.75 2.4293   ms 1 n + 1.75 < i Δ t < n + 2 , n = 0 , 1 , 2 , 3
The history of the flow states at the pipe midpoint in the simulation and analytical results are drawn in Figure 18. Both solvers nicely capture the reflection of the water-hammer wave in the pipe. The L2 in the piezometric head of P_HLL and M_HLL are 6.3965 and 6.3970, respectively, while L2 in the velocity of P_HLL and M_HLL are 0.1332 and 0.1333, respectively.
Interestingly, in this test, the P_HLL solver can obtain almost the same result with an arbitrary value of Pa; for example, the simulation result using Pa = 1 is drawn in Figure 19.
In Section 4, we have proven that under the framework of PSM, any non-negative value of Pa will produce a wave speed that is close to the acoustic wave speed, provided that the cell is under pressurized flow condition; see Figure 12 for detail. In this test, all the computational cells are under a pressurized flow condition, which makes the simulation result of Pa = 100 and Pa = 1 almost the same. In the flow regime transition simulation, PaH must be larger than the highest piezometric head.

5.3. Vasconcelos’s Experiment

This experiment was designed by Vasconcelos et al. [35] to study the filling process in a realistic stormwater storage tunnel. It is a 14.6 m-long horizontal tunnel with circular cross-sections of 9.4 cm in diameter; initially, stagnant water of 7.3 cm depth is established in the tunnel. A fill box with a 25 cm × 25 cm bottom and 31 cm height connects to the tunnel inlet. A surge tank with a constant diameter of 19 cm connects to the tunnel outlet. The experiment starts by constantly supplying 3.1 Ls−1 water into the fill box, and when water level inside the fill box reaches its top, water is simply spilled away. A gate is installed at the tunnel outlet; its opening is smaller than the initial water depth. When the filling bore collides with the gate, it triggers a water-hammer phenomenon. A ventilation tower is fixed upstream of the gate so that no air is trapped in the tunnel. The experiment setup is drawn in Figure 20.
The manning coefficient is 0.016 m1/6, acoustic wave speed is 300 ms−1 and head loss coefficient associated with the partially opened gate is 12, as suggested by Malekpour and Karney. In P_HLL, the values of Pa and Pb are 5 and 0.8. In M_HLL, the values of Ka and NS are 1.2 and 12. The computational cell is 0.1 m, and the time-step is set for a Courant number of 0.8. At the upstream end, the three unknowns are the discharge, the wetted area at the tunnel inlet and the water level in the fill box. At the downstream end, the three unknowns are discharge, the wetted area at tunnel outlet and the water level in the surge tank. The continuity, energy and characteristic equations are applied to obtain the three unknowns at each boundary [36]. The loci of flow states at x = 9.9 m in the simulation results of P_HLL and M_HLL are shown in Figure 21 and Figure 22.
The simulation results are in good agreement with the experimental data, and the two solvers have correctly computed the arrival time of the filling-bore front at x = 9.9 m. P_HLL has computed a piezometric head peak slightly larger than M_HLL. Most importantly, P_HLL and M_HLL do not smear the physical pressure oscillations.

5.4. Aureli’s Experiment

This experiment is conducted on a 12.12 m-long pipe with a 19.2 cm-diameter and 4 mm-wall thickness [25]. We used x to present the longitudinal distance to the pipe inlet. There is a sharp turn at about x = 7 m; the downward part has a slope of about 8.4% and the upward part has a slope of −27.7%. Six piezometric transducers were installed in the bottom at x = 1 m, 3 m, 4.5 m, 6.8 m, 7.06 m and 8.52 m. A sluice gate was installed at approximately x = 5 m; it is closed at the initial conditions. Stagnant water with a 22.5 cm piezometric head at transducer 1 was established behind the sluice gate; the rest of the pipe was empty. As the experiment began, the gate was lifted within 0.2 s, setting flush into the pipe. The pipe inlet was blocked so that no water flows through it, while the outlet was totally open. The experimental setup is shown in Figure 23.
In P_HLL, the values of Pa and Pb are 4 and 0.7; in M_HLL, the values of Ka and NS are 1.4 and 5. The computational cell is 0.04 m, acoustic wave speed is 200 ms−1 and time-step is set for a Courant number of 0.8. A reflective boundary condition was set at the upstream end, while a transmissive boundary condition was set at the downstream end. At the wet/dry interface, the estimates of wave speed followed Leon et al. [27]. The loci of the flow states at x = 6.8 m simulated by P_HLL and M_HLL are drawn in Figure 24 and Figure 25.
The simulation results of P_HLL and M_HLL show certain discrepancies: M_HLL overestimates the wave speed. It is because P_HLL adds numerical viscosity only if the depth is higher than the threshold value of PbH, while M_HLL adds numerical viscosity at any free-surface cells. Nonetheless, the simulation results of the two solvers are in good agreement with the experimental data.

6. Conclusions

Numerical oscillation is a critical problem in transient mixed flow simulations. These numerical oscillations arise from the ambiguity about the location of the filling-bore within one computational cell, which cannot be captured even with high-order finite volume methods. First order finite volume methods have failed to suppress them while capturing the filling-bore front.
Four oscillation-suppressing methods were tested on a benchmark model, with three of them failing to get a satisfactory result under a high acoustic wave speed. The key is to admit numerical viscosity before the flow regime transition occurs. Numerical viscosity can be added by artificially increasing the magnitude of the wave speed in the HLL solver. Following this idea, this paper presents a P_HLL solver that has two parameters: Pa and Pb. Pa needs to be larger than the highest piezometric head while Pb needs to be between 0.7 and 0.9. This solver adds numerical viscosity when the water depth is above PbH so that a smooth transition between the free-surface and pressurized flows can be achieved. The amount of numerical viscosity increases with Pa, while a large Pa does not smear the simulation result significantly. Therefore, one can always start by adopting a Pa that is large enough under realistic applications, for example 10.
The P_HLL solver is tested in several numerical tests, in which a good agreement between the simulation results and analytical results or experimental data is found. In the test where the wet/dry interface is presented, the P_HLL solver achieves a more accurate result than M_HLL. Meanwhile, P_HLL solver has sufficiently suppressed the numerical oscillations, and accurately captured the propagation of the filling-bore. Compared to the M_HLL solver proposed by Malekpour and Karney, the P_HLL is more convenient to use as the parameters in this solver are easier to determine.
The result in this paper can provide useful information for readers to design an oscillation-suppressing method. It provides an alternative method to the M_HLL solver, which may be used in the parallelization of computation.
This method is proposed under the framework of PSM; it is limited to the simulation of flow regime transition. Besides, since the air pressure is not counted, this method cannot be applied to simulate flow regime transition where air pockets are present, which is out of the scope of this paper.

Author Contributions

Conceptualization, Z.M., G.G., Z.Y.; resources, Z.M., G.G., Z.Y.; writing—original draft, Z.M.; writing—review and editing, Z.M. All authors have read and agreed to the published version of the manuscript.

Funding

The authors acknowledge the support of the NSFC grant 51979202 and NSFC grant 51879199.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Cong, J.; Chan, S.N.; Lee, J.H.W. Geyser formation by release of entrapped air from horizontal pipe into vertical shaft. J. Hydraul. Eng. 2017, 143, 04017039. [Google Scholar] [CrossRef]
  2. Liu, L.J.; Shao, W.Y.; Zhu, D.Z. Experimental study on stormwater geyser in vertical shaft above junction chamber. J. Hydraul. Eng. 2020, 146. [Google Scholar] [CrossRef]
  3. Chegini, T.; Leon, A.S. Numerical investigation of field-scale geysers in a vertical shaft. J. Hydraul. Res. 2019. [Google Scholar] [CrossRef]
  4. Cataño-Lopera, Y.A.; Tokyay, T.E.; Martin, J.E.; Schmidt, A.R.; Lanyon, R.; Fitzpatrick, K.; Scalise, C.F.; García, M.H. Modeling of a transient event in the tunnel and reservoir plan system in Chicago, Illinois. J. Hydraul. Eng. 2014, 140, 05014005. [Google Scholar] [CrossRef]
  5. Fuamba, M.; Bouaanani, N.; Marche, C. Modeling of dam break wave propagation in a partially ice-covered channel. Adv. Water Resour. 2007, 30, 2499–2510. [Google Scholar] [CrossRef]
  6. Holly, F.M.; Merkley, G.P. Unique problems in modeling irrigation canals. J. Irrig. Drain. Eng. 1993, 119, 656–662. [Google Scholar] [CrossRef]
  7. Guo, Q.; Song, C.C.S. Surging in urban storm drainage systems. J. Hydraul. Eng. 1990, 116, 1523–1537. [Google Scholar] [CrossRef]
  8. Liu, G.Q.; Guan, G.H.; Wang, C.D. Transition mode of long distance water delivery project before freezing in winter. J. Hydroinformatics 2013, 15, 306–320. [Google Scholar] [CrossRef]
  9. Mao, Z.H.; Guan, G.H.; Yang, Z.H.; Zhong, K. Linear model of water movements for large-scale inverted siphon in water distribution system. J. Hydroinformatics 2019, 21, 1048–1063. [Google Scholar] [CrossRef] [Green Version]
  10. Guan, G.; Clemmens, A.J.; Kacerek, T.F.; Wahlin, B.T. Applying water-level difference control to central arizona project. J. Irrig. Drain. Eng. 2011, 137, 747–753. [Google Scholar] [CrossRef]
  11. Cunge, J.A.; Wegner, M. Intégration numérique des équations d’écoulement de barré de Saint-Venant par un schéma implicite de différences finies. Houille Blanche Rev. Int. De L Eau 1964, 1, 33–39. [Google Scholar] [CrossRef] [Green Version]
  12. 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. [Google Scholar] [CrossRef] [Green Version]
  13. Dazzi, S.; Maranzoni, A.; Mignosa, P. Local time stepping applied to mixed flow modelling. J. Hydraul. Res. 2016, 54, 145–157. [Google Scholar] [CrossRef]
  14. Fernandez-Pato, J.; Garcia-Navarro, P. A pipe network simulation model with dynamic transition between free surface and pressurized flow. In 12th International Conference on Computing and Control for the Water Industry, Ccwi2013, Perugia, Italy, 2–4 September 2013; Brunone, B., Giustolisi, O., Ferrante, M., Laucelli, D., Meniconi, S., Berardi, L., Campisano, A., Eds.; Elsevier Science Bv: Amsterdam, The Netherlands, 2014; Volume 70, pp. 641–650. [Google Scholar]
  15. Kerger, F.; Erpicum, S.; Dewals, B.J.; Archambeau, P.; Pirotton, M. 1D unified mathematical model for environmental flow applied to steady aerated mixed flows. Adv. Eng. Softw. 2011, 42, 660–670. [Google Scholar] [CrossRef] [Green Version]
  16. Ferreri, G.B.; Freni, G.; Tomaselli, P. Ability of preissmann slot scheme to simulate smooth pressurisation transient in sewers. Water Sci. Technol. 2010, 62, 1848–1858. [Google Scholar] [CrossRef] [PubMed]
  17. Maranzoni, A.; Mignosa, P. Numerical treatment of a discontinuous top surface in 2D shallow water mixed flow modeling. Int. J. Numer. Methods Fluids 2018, 86, 290–311. [Google Scholar] [CrossRef]
  18. Maranzoni, A.; Dazzi, S.; Aureli, F.; Mignosa, P. Extension and application of the Preissmann slot model to 2D transient mixed flows. Adv. Water Resour. 2015, 82, 70–82. [Google Scholar] [CrossRef]
  19. Toro, E.F.; Garcia-Navarro, P. Godunov-type methods for free-surface shallow flows: A review. J. Hydraul. Res. 2007, 45, 736–751. [Google Scholar] [CrossRef]
  20. Kitamura, K.; Shima, E. Numerical experiments on anomalies from stationary, slowly moving, and fast-moving shocks. AIAA J. 2019, 57, 1763–1772. [Google Scholar] [CrossRef]
  21. Johnsen, E. Analysis of numerical errors generated by slowly moving shock waves. AIAA J. 2013, 51, 1269–1274. [Google Scholar] [CrossRef]
  22. Zaide, D.W.; Roe, P.L. Flux functions for reducing numerical shockwave anomalies. In Proceedings of the Seventh International Conference on Computational Fluid Dynamics, Big Island, Hawaii, 9–13 July 2012. [Google Scholar]
  23. 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] [Green Version]
  24. Malekpour, A.; Karney, B.W. Spurious numerical oscillations in the preissmann slot method: Origin and suppression. J. Hydraul. Eng. 2016, 142, 04015060. [Google Scholar] [CrossRef] [Green Version]
  25. Aureli, F.; Dazzi, S.; Maranzoni, A.; Mignosa, P. Validation of single- and two-equation models for transient mixed flows: A laboratory test case. J. Hydraul. Res. 2015, 53, 440–451. [Google Scholar] [CrossRef]
  26. 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] [Green Version]
  27. Leon, A.S.; Ghidaoui, M.S.; Schmidt, A.R.; Garcia, M.H. Application of Godunov-type schemes to transient mixed flows. J. Hydraul. Res. 2009, 47, 147–156. [Google Scholar] [CrossRef]
  28. Vasconcelos, J.G.; Wright, S.J.; Roe, P.L. Improved simulation of flow regime transition in sewers: Two-component pressure approach. J. Hydraul. Eng. 2006, 132, 553–562. [Google Scholar] [CrossRef]
  29. Kim, S.S.; Kim, C.; Rho, O.H.; Hong, S.K. Cures for the shock instability: Development of a shock-stable Roe scheme. J. Comput. Phys. 2003, 185, 342–374. [Google Scholar] [CrossRef]
  30. Toro, E.F.; Billett, S.J. Centred TVD schemes for hyperbolic conservation laws. IMA J. Numer. Anal. 2000, 20, 47–79. [Google Scholar] [CrossRef]
  31. Toro, E. Shock Capturing Methods for Free Surface Shallow Water Flows; John Wiley: New York, NY, USA, 2001. [Google Scholar]
  32. Leon, A.S.; Ghidaoui, M.S.; Schmidt, A.R.; Garcia, M.H. Godunov-type solutions for transient flows in sewers. J. Hydraul. Eng.-Asce 2006, 132, 800–813. [Google Scholar] [CrossRef] [Green Version]
  33. 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]
  34. Bertaglia, G.; Ioriatti, M.; Valiani, A.; Dumbser, M.; Caleffi, V. Numerical methods for hydraulic transients in visco-elastic pipes. J. Fluids Struct. 2018, 81, 230–254. [Google Scholar] [CrossRef] [Green Version]
  35. Vasconcelos, J.; Wright, S. Experimental Investigation of Surges in a Stormwater Storage Tunnel. J. Hydraul. Eng. 2005, 131, 853–861. [Google Scholar] [CrossRef] [Green Version]
  36. Vasconcelos, J.; Wright, S. Numerical Modeling of the Transition between Free Surface and Pressurized Flow in Storm Sewers. J. Water Manag. Model. 2004. [Google Scholar] [CrossRef] [Green Version]
Figure 1. A rectangular conduct with a slot on its top: (a) free surface flow; (b) pressurized flow.
Figure 1. A rectangular conduct with a slot on its top: (a) free surface flow; (b) pressurized flow.
Water 12 01245 g001
Figure 2. Front view of the benchmark model.
Figure 2. Front view of the benchmark model.
Water 12 01245 g002
Figure 3. Comparison of results simulated by the numerical filtering method (filtered) and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Figure 3. Comparison of results simulated by the numerical filtering method (filtered) and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Water 12 01245 g003
Figure 4. Comparison of the results simulated by hybrid one and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Figure 4. Comparison of the results simulated by hybrid one and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Water 12 01245 g004
Figure 5. Comparison of the results simulated by hybrid two and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Figure 5. Comparison of the results simulated by hybrid two and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Water 12 01245 g005
Figure 6. Comparison of the results simulated by M_HLL, with Ka = 1.4 and NS = 5, and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Figure 6. Comparison of the results simulated by M_HLL, with Ka = 1.4 and NS = 5, and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Water 12 01245 g006
Figure 7. Comparison of the results simulated by M_HLL, with Ka = 1.4 and NS = 3, and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Figure 7. Comparison of the results simulated by M_HLL, with Ka = 1.4 and NS = 3, and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Water 12 01245 g007
Figure 8. A left rarefaction wave with a right shock wave at xi+1/2.
Figure 8. A left rarefaction wave with a right shock wave at xi+1/2.
Water 12 01245 g008
Figure 9. A left shock wave with a right shock wave at xi+1/2.
Figure 9. A left shock wave with a right shock wave at xi+1/2.
Water 12 01245 g009
Figure 10. (a) Locus of the flow states in cell i + 1 simulated by HLL and P_HLL (Pa = 5, Pb = 0.7); (b) history of the flow states in cell i + 1 simulated by P_HLL (Pa = 5, Pb = 0.7).
Figure 10. (a) Locus of the flow states in cell i + 1 simulated by HLL and P_HLL (Pa = 5, Pb = 0.7); (b) history of the flow states in cell i + 1 simulated by P_HLL (Pa = 5, Pb = 0.7).
Water 12 01245 g010
Figure 11. SwR2 under different values of Pa.
Figure 11. SwR2 under different values of Pa.
Water 12 01245 g011
Figure 12. |SwL2| under different values of Pa.
Figure 12. |SwL2| under different values of Pa.
Water 12 01245 g012
Figure 13. Comparison of the results simulated by P_HLL, with Pb = 0.7 and Pa = 5 or Pa = 10, and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Figure 13. Comparison of the results simulated by P_HLL, with Pb = 0.7 and Pa = 5 or Pa = 10, and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Water 12 01245 g013
Figure 14. Comparison of the results simulated by P_HLL, with Pa = 5 and Pb = 0.7 or Pb = 0.9, and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Figure 14. Comparison of the results simulated by P_HLL, with Pa = 5 and Pb = 0.7 or Pb = 0.9, and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Water 12 01245 g014
Figure 15. A sketch of the reservoir–conduct–reservoir model.
Figure 15. A sketch of the reservoir–conduct–reservoir model.
Water 12 01245 g015
Figure 16. Comparison of the results simulated by P_HLL (Pa = 5 and Pb = 0.8), M_HLL (Ka = 1.4 and NS = 5) and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Figure 16. Comparison of the results simulated by P_HLL (Pa = 5 and Pb = 0.8), M_HLL (Ka = 1.4 and NS = 5) and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Water 12 01245 g016
Figure 17. A sketch of the conduct–reservoir model.
Figure 17. A sketch of the conduct–reservoir model.
Water 12 01245 g017
Figure 18. Locus of the flow states at the midpoint of the pipe simulated by P_HLL (Pa = 100 and Pb = 0.8), M_HLL (Ka = 1.2 and NS = 12) and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Figure 18. Locus of the flow states at the midpoint of the pipe simulated by P_HLL (Pa = 100 and Pb = 0.8), M_HLL (Ka = 1.2 and NS = 12) and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Water 12 01245 g018
Figure 19. Locus of the flow states at the midpoint of pipe simulated by P_HLL (Pa = 1 and Pb = 0.8) and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Figure 19. Locus of the flow states at the midpoint of pipe simulated by P_HLL (Pa = 1 and Pb = 0.8) and the analytical result (AR): (a) piezometric head; (b) flow velocity.
Water 12 01245 g019
Figure 20. A sketch of the experimental setup.
Figure 20. A sketch of the experimental setup.
Water 12 01245 g020
Figure 21. Locus of flow states at x = 9.9 m simulated by P_HLL (Pa = 5 and Pb = 0.8) and experimental data: (a) piezometric head; (b) flow velocity.
Figure 21. Locus of flow states at x = 9.9 m simulated by P_HLL (Pa = 5 and Pb = 0.8) and experimental data: (a) piezometric head; (b) flow velocity.
Water 12 01245 g021
Figure 22. Locus of flow states at x = 9.9 m simulated by M_HLL (Ka = 1.2 and NS = 12) and experimental data: (a) piezometric head; (b) flow velocity.
Figure 22. Locus of flow states at x = 9.9 m simulated by M_HLL (Ka = 1.2 and NS = 12) and experimental data: (a) piezometric head; (b) flow velocity.
Water 12 01245 g022
Figure 23. A sketch of the experimental setup, drawn by Aureli et al. [25].
Figure 23. A sketch of the experimental setup, drawn by Aureli et al. [25].
Water 12 01245 g023
Figure 24. Locus of flow states at x = 6.8 m simulated by P_HLL (Pa = 4 and Pb = 0.7) and experimental data: (a) piezometric head; (b) flow velocity.
Figure 24. Locus of flow states at x = 6.8 m simulated by P_HLL (Pa = 4 and Pb = 0.7) and experimental data: (a) piezometric head; (b) flow velocity.
Water 12 01245 g024
Figure 25. Locus of flow states at x = 6.8 m simulated by M_HLL (Ka = 1.4 and NS = 5) and experimental data: (a) piezometric head; (b) flow velocity.
Figure 25. Locus of flow states at x = 6.8 m simulated by M_HLL (Ka = 1.4 and NS = 5) and experimental data: (a) piezometric head; (b) flow velocity.
Water 12 01245 g025

Share and Cite

MDPI and ACS Style

Mao, Z.; Guan, G.; Yang, Z. Suppress Numerical Oscillations in Transient Mixed Flow Simulations with a Modified HLL Solver. Water 2020, 12, 1245. https://doi.org/10.3390/w12051245

AMA Style

Mao Z, Guan G, Yang Z. Suppress Numerical Oscillations in Transient Mixed Flow Simulations with a Modified HLL Solver. Water. 2020; 12(5):1245. https://doi.org/10.3390/w12051245

Chicago/Turabian Style

Mao, Zhonghao, Guanghua Guan, and Zhonghua Yang. 2020. "Suppress Numerical Oscillations in Transient Mixed Flow Simulations with a Modified HLL Solver" Water 12, no. 5: 1245. https://doi.org/10.3390/w12051245

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