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

#### **Zhonghao Mao, Guanghua Guan and Zhonghua Yang \***

State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China; 2011301580373@whu.edu.cn (Z.M.); ggh@whu.edu.cn (G.G.)

**\*** Correspondence: yzh@whu.edu.cn

Received: 21 February 2020; Accepted: 21 April 2020; Published: 27 April 2020

**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, *P*<sup>a</sup> and *P*b, and their values can be determined easily. It can sufficiently suppress numerical oscillations under an acoustic wave speed of 1000 ms<sup>−</sup>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.

**Keywords:** flow regime transition; finite volume methods; numerical oscillations; numerical viscosity; Preissmann slot model

#### **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–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–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–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–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–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.

**Figure 1.** A rectangular conduct with a slot on its top: (**a**) free surface flow; (**b**) pressurized flow.

Under the framework of PSM, the governing equations of one-dimensional flow regime transition with uniform cross-sections can be written as [26]:

$$\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}(\mathbf{U})}{\partial \mathbf{x}} = \mathbf{S}(\mathbf{U})$$

$$\mathbf{U} = \begin{pmatrix} A \\ Q \end{pmatrix} \mathbf{F}(\mathbf{U}) = \begin{pmatrix} Q \\ \frac{Q^2}{A} + gI(h) \end{pmatrix} \mathbf{S}(\mathbf{U}) = \begin{pmatrix} 0 \\ gA \{\mathcal{S}\_b - \mathcal{S}\_f\} \end{pmatrix} \tag{1}$$

$$I(h) = \int\_0^h (h - \xi) I(x, \xi) d\xi$$

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 = \frac{n^2 \mu |u|}{R^{4/3}}\tag{2}$$

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) := \begin{cases} B & h \le H \\\ B\_{sl} & h > H \end{cases} \tag{3}$$

$$A(h) = \begin{cases} Bh & h \le H \\ BH + B\_{sl}(h - H) & h > H \end{cases} \tag{4}$$

$$I(h) = \begin{cases} 0.5Bl^2 & h \le H\\ BH(h - 0.5H) + 0.5B\_{sl}(h - H)^2 & h > H \end{cases} \tag{5}$$

$$c(h) = \begin{cases} \sqrt{gh} & h \le H \\ \sqrt{g \frac{A}{B\_{sl}}} & h > H \end{cases} \tag{6}$$

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

$$\mathbf{U}\_{i}^{\*} = \mathbf{U}\_{i}^{n} - \frac{\Delta t\_{i}}{\Delta \mathbf{x}\_{i}} (\mathbf{F}\_{i+1/2} - \mathbf{F}\_{i-1/2}) \tag{7}$$

$$\mathbf{U}\_{i}^{n+1} = \mathbf{U}\_{i}^{\*} + \Delta t\_{i} \mathbf{S}\left(\mathbf{U}\_{i}^{\*}\right) \tag{8}$$

#### **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−<sup>1</sup> and the slot width is 9.8 <sup>×</sup> <sup>10</sup>−<sup>6</sup> 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.

**Figure 2.** Front view of the benchmark model.

Under initial condition, flow states **U***<sup>L</sup>* in the reservoir and **U***<sup>R</sup>* in the conduct are discontinuous:

$$
\begin{bmatrix} h\_L \\ u\_L \end{bmatrix} = \begin{bmatrix} 4 \text{ m} \\ 0 \text{ m} \text{s}^{-1} \end{bmatrix} \begin{bmatrix} h\_R \\ u\_R \end{bmatrix} = \begin{bmatrix} 0.6 \text{ m} \\ 0 \text{ m} \text{s}^{-1} \end{bmatrix} \tag{9}
$$

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 **U**<sup>1</sup> and **U**2, respectively. Then **U**<sup>2</sup> can be obtained by solving the following equations iteratively [24]:

$$h\_1 = h\_2 + \frac{u\_2^2}{2g} \tag{10}$$

$$
\mu\_2 = \mu\_\mathbb{R} + \sqrt{\frac{[gI(A\_\mathbb{R}) - gI(A\_2)](A\_\mathbb{R} - A\_2)}{A\_\mathbb{R}A\_2}}\tag{11}
$$

where **U**<sup>1</sup> = **U***<sup>L</sup>* by assuming that flow velocity inside the reservoir is negligible. In this case, *h*<sup>2</sup> and *u*<sup>2</sup> are 3.167 m and 4.0334 ms<sup>−</sup>1, a ghost cell is set at the upstream boundary adopting *h*<sup>2</sup> and *u*2. Since **U**<sup>2</sup> is connected to **U***<sup>R</sup>* through a right shock, the flow states inside the conduct ultimately will take place by **U**2. The propagation speed of the filling-bore is given by the Rankine–Hugoniot condition, which is 10.067 ms<sup>−</sup>1. Then we can construct the analytical result in the benchmark model at *t*0:

$$h(\mathbf{x}) = \begin{cases} 3.167 \,\mathrm{m} & \mathbf{x} < 10.067 t\_0 \\ 0.6 \,\mathrm{m} & \mathbf{x} > 10.067 t\_0 \end{cases}, \mu(\mathbf{x}) = \begin{cases} 4.0334 \,\mathrm{ms}^{-1} & \mathbf{x} < 10.067 t\_0 \\ 0 \,\mathrm{ms}^{-1} & \mathbf{x} > 10.067 t\_0 \end{cases} \tag{12}$$

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 **U***i*+1/<sup>2</sup> at *xi*<sup>+</sup>1/<sup>2</sup> satisfy the following equations:

$$u\_{i+1/2} = \begin{cases} u\_i - \sqrt{g \frac{\left[I(A\_i) - I\left(A\_{i+1/2}\right)\right]\left(A\_i - A\_{i+1/2}\right)}{A\_i A\_{i+1/2}}} & A\_{i+1/2} > A\_i\\ u\_i + \int\_0^{A\_L} \sqrt{\frac{\mathcal{S}}{ab}} d\alpha - \int\_0^{A\_{i+1/2}} \sqrt{\frac{\mathcal{S}}{ab}} d\alpha & A\_{i+1/2} < A\_i \end{cases} \tag{13}$$

$$u\_{i+1/2} = \begin{cases} u\_{i+1} + \sqrt{g \frac{\left[l(A\_{i+1}) - l\left(A\_{i+1/2}\right)\right]\left(A\_{i+1} - A\_{i+1/2}\right)}{A\_{i+1}A\_{i+1/2}}} & A\_{i+1/2} > A\_{i+1} \\\ u\_{i+1} - \int\_0^{A\_{i+1}} \sqrt{\frac{g}{ab}} da + \int\_0^{A\_{i+1/2}} \sqrt{\frac{g}{ab}} da & A\_{i+1/2} < A\_{i+1} \end{cases} \tag{14}$$

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:

$$\mathbf{U}\_{i}^{n+1} = (1 - 2\varepsilon)\mathbf{U}\_{i}^{n+1} + \varepsilon(\mathbf{U}\_{i-1}^{n+1} + \mathbf{U}\_{i+1}^{n+1})\tag{15}$$

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.

**Figure 3.** Comparison of results simulated by the numerical filtering method (filtered) and the analytical result (AR): (**a**) piezometric head; (**b**) flow velocity.

#### *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 λ*<sup>j</sup> <sup>i</sup>*+1/2 as eigenvalues of the linearized Jacobian matrix, **R***<sup>j</sup> <sup>i</sup>*+1/2 as the corresponding right eigenvectors and <sup>α</sup>*<sup>j</sup> <sup>i</sup>*+1/2 as the wave strengths across the cell boundary:

$$
\lambda\_{i+1/2}^1 = \frac{Q\_{i+1/2}}{\overline{A}\_{i+1/2}} - \overline{c}\_{i+1/2}, \\
\lambda\_{i+1/2}^2 = \frac{Q\_{i+1/2}}{\overline{A}\_{i+1/2}} + \overline{c}\_{i+1/2} \tag{16}
$$

$$\mathbf{R}\_{i+1/2}^{1} = \begin{bmatrix} 1, \lambda\_{i+1/2}^{1} \end{bmatrix}^{T}, \mathbf{R}\_{i+1/2}^{2} = \begin{bmatrix} 1, \lambda\_{i+1/2}^{2} \end{bmatrix}^{T} \tag{17}$$

$$a\_{i+1/2}^1 = \frac{\lambda\_{i+1/2}^2 (A\_{i+1} - A\_i) - (Q\_{i+1} - Q\_i)}{\lambda\_{i+1/2}^2 - \lambda\_{i+1/2}^1}, \\ a\_{i+1/2}^2 = \frac{(Q\_{i+1} - Q\_i) - \lambda\_{i+1/2}^1 (A\_{i+1} - A\_i)}{\lambda\_{i+1/2}^2 - \lambda\_{i+1/2}^1} \tag{18}$$

where *<sup>Q</sup>*2*i*<sup>+</sup>1/2, *<sup>A</sup>*2*i*+1/2 and <sup>2</sup>*ci*<sup>+</sup>1/2 are Roe averages:

$$
\widetilde{A}\_{i+1/2} = (A\_i A\_{i+1})^{1/2} \tag{19}
$$

$$\widetilde{Q}\_{i+1/2} = \frac{Q\_{i+1}(A\_i)^{1/2} + Q\_i(A\_{i+1})^{1/2}}{(A\_i)^{1/2} + (A\_{i+1})^{1/2}} \tag{20}$$

$$\widetilde{c\_{i+1/2}} = \begin{cases} \left[ \mathcal{S} \frac{l(A\_{i+1}) - l(A\_i)}{A\_{i+1} - A\_i} \right]^{1/2} & A\_{i+1} \neq A\_i\\ \left[ \mathcal{S} \frac{A\_i + A\_{i+1}}{b\_i + b\_{i+1}} \right]^{1/2} & A\_{i+1} = A\_i \end{cases} \tag{21}$$

Then the fluxes obtained by the Roe solver are written as

$$\mathbf{F}\_{i+1/2}^{\text{Roe}} = \frac{1}{2} \Big[ \mathbf{F} \Big( \mathbf{U}\_{i+1}^{\text{n}} \big) + \mathbf{F} \Big( \mathbf{U}\_{i}^{\text{n}} \big) \Big] - \frac{1}{2} \sum\_{j=1}^{2} \Big| \boldsymbol{\lambda}\_{i+1/2}^{j} \Big| \boldsymbol{\alpha}\_{i+1/2}^{j} \mathbf{R}\_{i+1/2}^{j} \tag{22}$$

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:

$$\mathbf{F}\_{i+1/2}^{\mathrm{LxF}} = \frac{1}{2} \Big[ \mathbf{F} \Big( \mathbf{U}\_{i+1}^{n} \Big) + \mathbf{F} \Big( \mathbf{U}\_{i}^{n} \Big) \Big] - \frac{\Delta x}{2 \Delta t} (\mathbf{U}\_{i+1} - \mathbf{U}\_{i}) \tag{23}$$

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:

$$\left| \lambda\_{i+1/2}^{1,2} \right| = \min \left[ \frac{\Delta x}{\Delta t} \Big| \frac{\overline{Q}\_{i+1/2}}{\overline{A}\_{i+1/2}} \mp \overline{c}\_{i+1/2} \right| + L(\Delta c\_{i+1/2})^L \frac{\Delta x}{\Delta t} \Big| \tag{24}$$

$$\Delta c\_{i+1/2} = \frac{|\overleftarrow{c}\_{i+1} - \overleftarrow{c}\_i|}{\max\big( |\overleftarrow{c}\_{i+1} - \overleftarrow{c}\_i| \big)\_{i=1..N-1}} \tag{25}$$

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*<sup>+</sup>1/<sup>2</sup> = 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*<sup>+</sup>1/<sup>2</sup> 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.

**Figure 4.** Comparison of the results simulated by hybrid one and the analytical result (AR): (**a**) piezometric head; (**b**) flow velocity.

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:

$$\mathbf{F}\_{i+1/2}^{\text{hybrid}} = \begin{cases} \mathbf{F}\_{i+1/2}^{\text{FORCE}} & (h\_L - H\_L)(h\_R - H\_R) < 0\\ \mathbf{F}\_{i+1/2}^{\text{HLL}} & \text{otherwise} \end{cases} \tag{26}$$

*Water* **2020**, *12*, 1245

This method is referred to as hybrid two from here on further. The HLL solver starts by estimating the largest wave speed *S*wR and smallest wave speed *S*wL in the Riemann problem. The fluxes are computed based on the signs of *S*wL and *S*wR:

$$\mathbf{F}\_{i+1/2}^{\text{HLL}} = \begin{cases} \text{F}(\mathbf{U}\_{i}) & S\_{\text{wL}} > 0\\ \frac{S\_{\text{wR}}\text{F}(\mathbf{U}\_{i}) - S\_{\text{wL}}\text{F}(\mathbf{U}\_{i+1}) + S\_{\text{wR}}\text{S}\_{\text{wL}}(\mathbf{U}\_{i+1} - \mathbf{U}\_{i})}{S\_{\text{wR}} - S\_{\text{wL}}} & S\_{\text{wL}} \le 0\\ \text{F}(\mathbf{U}\_{i+1}) & S\_{\text{wR}} < 0 \end{cases} \quad \text{s}\_{\text{wR}} \le 0 \quad \text{and} \quad S\_{\text{wR}} \ge 0 \tag{27}$$

The choices of *S*wL and *S*wR follow Toro [31]:

$$S\_{wL} = \mu\_i - \Omega\_{i\prime} S\_{wR} = \mu\_{i+1} + \Omega\_{i+1} \tag{28}$$

$$\Omega\_{K(K=i;i+1)} = \begin{cases} \sqrt{8 \frac{[I\{A\_{i+1/2}\} - I(A\_{\mathcal{K}})]A\_{i+1/2}}{A\_{\mathcal{K}}(A\_{i+1/2} - A\_{\mathcal{K}})}} & A\_{i+1/2} > A\_{K} \\ \mathcal{c}\_{K} & A\_{i+1/2} \le A\_{K} \end{cases} \tag{29}$$

where *Ai*<sup>+</sup>1/<sup>2</sup> is an estimate of the wetted area at *xi*<sup>+</sup>1/2; we adopt the one proposed by Leno et al. [32], which admits the minimum amount of numerical viscosity:

$$A\_{i+1/2} = \frac{A\_i + A\_{i+1}}{2} \left( 1 + \frac{u\_i - u\_{i+1}}{c\_i + c\_{i+1}} \right) \tag{30}$$

The FORCE flux can be written as the algebraic average of the LxF scheme and Lax–Wendorff scheme:

$$\mathbf{F}\_{i+1/2}^{\text{LW}} = \mathbf{F} \{ \mathbf{U}\_{i+1/2}^{\text{LW}} \} \mathbf{U}\_{i+1/2}^{\text{LW}} = \frac{1}{2} \{ \mathbf{U}\_{i}^{\text{n}} + \mathbf{U}\_{i+1}^{\text{n}} \} - \frac{\Delta x}{2\Delta t} \Big[ \mathbf{F} \{ \mathbf{U}\_{i+1}^{\text{n}} \} - \mathbf{F} \{ \mathbf{U}\_{i}^{\text{n}} \} \Big] \tag{31}$$

$$\mathbf{F}\_{i+1/2}^{\text{FORCE}} = \frac{1}{2} (\mathbf{F}\_{i+1/2}^{\text{LW}} + \mathbf{F}\_{i+1/2}^{\text{L\times F}}) \tag{32}$$

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.

**Figure 5.** Comparison of the results simulated by hybrid two and the analytical result (AR): (**a**) piezometric head; (**b**) flow velocity.

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 *S*wL and *S*wR. In fact, the HLL fluxes equal the LxF fluxes when |*S*wL| and |*S*wR| equal Δ*x*/Δ*t*. 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*<sup>+</sup>1/<sup>2</sup> in Equation (29) is computed according to a reference depth *h*G:

$$A\_{i+1/2} = A(h\_{\mathbb{G}}), h\_{\mathbb{G}} = \mathsf{K}\_{\mathsf{a}} \max(h\_{i-\mathbb{NS}}, h\_{i-\mathbb{NS}+1}, \dots, h\_{i-1}, h\_i, h\_{i+1}, \dots, h\_{i+\mathbb{NS}-1}, h\_{i+\mathbb{NS}}) \tag{33}$$

The depth *h*<sup>G</sup> is defined as *Ka*, multiplying the highest piezometric head in the 2*NS*+1 cells surrounding cell *i*, while *Ka* > 1 and *NS* ≥ 3. The solution of Equation (33) produces a larger magnitude of *S*wL and *S*wR; 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.

**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.

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.

**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.

Thus, the values of *Ka* and *NS* must be well-tuned. Meanwhile, the way to determine *h*<sup>G</sup> 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*<sup>+</sup>1/<sup>2</sup> depends on the water depths at cell *i* and *i* + 1. When *hi* and *hi*<sup>+</sup><sup>1</sup> are below *PbH*, *Ai*<sup>+</sup>1/<sup>2</sup> is computed using Equation (30) to admit the minimum amount of numerical viscosity, otherwise *Ai*<sup>+</sup>1/<sup>2</sup> is computed according to a constant depth *PaH*:

$$A\_{i+1/2} = A(P\_{\text{d}}H), \ h\_{i} > P\_{\text{b}}H \quad \text{or} \quad h\_{i+1} > P\_{\text{b}}H \tag{34}$$

*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*<sup>+</sup>1/<sup>2</sup> in the benchmark model. Suppose *hi* and *hi*<sup>+</sup><sup>1</sup> is 3.167 m and 0.6 m, respectively, and *ui* and *ui*<sup>+</sup><sup>1</sup> is 4.0334 ms−<sup>1</sup> 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 *S*wL (noted as *S*wL1) as the speed of the left rarefaction wave head, and *S*wR (noted as *S*wR1) as the speed of right shock wave:

$$S\_{\rm w11} = u\_i - \mathbf{c}\_{i\prime} S\_{\rm wR1} = u\_{i+1} + \sqrt{g \frac{[I(A\_{i+1/2}) - I(A\_{i+1})]A\_{i+1/2}}{A\_{i+1}(A\_{i+1/2} - A\_{i+1})}} \tag{35}$$

A sketch of two waves in the Riemann problem is shown in Figure 8.

**Figure 8.** A left rarefaction wave with a right shock wave at *xi*<sup>+</sup>1/2.

Because cell *i* is under a pressurized flow condition, it is easy to see that |*S*wL1| nearly equals the acoustic wave speed. The entropy condition of right shock wave requires

$$u\_{i+1/2} + c\_{i+1/2} > S\_{\text{WR1}} > u\_{i+1} + c\_{i+1} \tag{36}$$

The magnitude of *S*wR1 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, *S*wR2 is the speed of the right shock wave, and *S*wL2 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\_{\text{w1.2}} < u\_i - c\_i \tag{37}$$

Since cell *i* is under a pressurized flow condition, *ui* << *ci*; thus, |*S*wL2|>|*S*wL1| and they are both close to the acoustic wave speed. The speed of the right shock wave increases with *Ai*<sup>+</sup>1/2; thus, *S*wR2 > *S*wR1. This larger magnitude of *S*wR admits more mass and momentum into cell *i* + 1 before it becomes pressurized. The loci of **U***i*+<sup>1</sup> simulated by HLL and P\_HLL (*Pa* = 5, *Pb* = 0.7) are drawn in Figure 10.

**Figure 9.** A left shock wave with a right shock wave at *xi*<sup>+</sup>1/2.

**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).

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−<sup>1</sup> 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 *S*wR2 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 *S*wL2 and *S*wR2 are related to the value of *Ai*<sup>+</sup>1/2, which consequently depends on the value of *Pa*; see Figures 11 and 12.

Figure 11 shows that *S*wR2 increases with *Pa*, while|*S*wL2| barely changes with *Pa* and stays close to the acoustic wave speed. However, *S*wR2 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 *P*<sup>b</sup> = 0.7 and different values of *P*<sup>a</sup> 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 *P*<sup>b</sup> = 0.7 or 0.9 are drawn in Figure 14.

**Figure 11.** *S*wR2 under different values of *Pa*.

**Figure 12.** |*S*wL2| under different values of *Pa*.

**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.

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.

**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.

#### **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<sup>−</sup>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.

**Figure 15.** A sketch of the reservoir–conduct–reservoir model.

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<sup>−</sup>1, and the flow states behind it are 2.42 m and <sup>−</sup>3.3717 ms<sup>−</sup>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 *t*<sup>0</sup> is

$$\begin{aligned} h(\mathbf{x}) &= \begin{cases} 3.167 \text{ m} & \mathbf{x} < 10.067 t\_0 \\ 0.6 \text{ m} & 10.067 t\_0 < \mathbf{x} < 200 - 8.429 t\_0 \\ 2.42 \text{ m} & 200 - 8.429 t\_0 < \mathbf{x} \\ 4.0334 \text{ ms}^{-1} & \mathbf{x} < 10.067 t\_0 \\ 0 \text{ ms}^{-1} & 10.067 t\_0 < \mathbf{x} < 200 - 8.429 t\_0 \\ -3.3717 \text{ ms}^{-1} & 200 - 8.429 t\_0 < \mathbf{x} \end{cases} \end{aligned} \tag{38}$$

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 *L*2-norm [34] is used to evaluate the accuracy of P\_HLL and M\_HLL. In the following equation, ϕ*<sup>i</sup>* stands for the simulation result at cell *i*, while ϕref stands for the analytical result.

(ϕ*<sup>i</sup>* − ϕref)

2 ⎞

0.5

(39)

*L*<sup>2</sup> = ⎛

1

*N*

⎜⎜⎜⎜⎜⎝ *N i* = 1 ⎟⎟⎟⎟⎟⎠ (**a**) (**b**)

**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.

The *L*<sup>2</sup> in the piezometric head of P\_HLL and M\_HLL are 0.2963 and 0.2913, respectively; the *L*<sup>2</sup> 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−<sup>1</sup> 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.

**Figure 17.** A sketch of the conduct–reservoir model.

The *L*<sup>2</sup> indicator is adopted to evaluate the accuracy of the two solvers; it is defined as

$$L\_2 = \left(\frac{1}{N} \sum\_{i=1}^{N} \left(q\_i - q\_{\text{ref}}\right)^2\right)^{0.5} \tag{40}$$

where *Nt* is the number of time step, ϕ*<sup>i</sup>* is the simulation result at the midpoint of the pipe, and ϕref is the analytical result, which is given as

$$\begin{aligned} h\_{i} &= \begin{cases} 45 \text{ m} & n < i \Delta t < n + 0.25 \\ -3.05 \text{ m} & n + 0.25 < i \Delta t < n + 0.75 \\ 45 \text{ m} & n + 0.75 < i \Delta t < n + 1.25 \\ 93.05 \text{ m} & n + 1.25 < i \Delta t < n + 1.75 \\ 45 \text{ m} & n + 1.75 < i \Delta t < n + 2 \\ 2.4293 \text{ ms}^{-1} & n < i \Delta t < n + 0.25 \\ 2.0377 \text{ ms}^{-1} & n + 0.25 < i \Delta t < n + 0.75 \\ 1.6461 \text{ ms}^{-1} & n + 0.75 < i \Delta t < n + 1.25 \\ 2.0377 \text{ ms}^{-1} & n + 1.25 < i \Delta t < n + 1.75 \\ 2.4293 \text{ ms}^{-1} & n + 1.75 < i \Delta t < n + 2 \end{cases} \end{aligned} \tag{41}$$

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 *L*<sup>2</sup> in the piezometric head of P\_HLL and M\_HLL are 6.3965 and 6.3970, respectively, while *L*<sup>2</sup> in the velocity of P\_HLL and M\_HLL are 0.1332 and 0.1333, respectively.

**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.

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.

**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.

#### *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−<sup>1</sup> 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.

**Figure 20.** A sketch of the experimental setup.

The manning coefficient is 0.016 m1/6, acoustic wave speed is 300 ms−<sup>1</sup> 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 Figures 21 and 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.

**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 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.

#### *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.

**Figure 23.** A sketch of the experimental setup, drawn by Aureli et al. [25].

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−<sup>1</sup> 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 Figures 24 and 25.

**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 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.

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**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
