Advanced Numerical Solver for Dam-Break Flow Application

In this paper, a HLL (Harten Lax van Leer) approximate Riemann solver with MUSCL scheme (Monotonic Upwind Schemes for Conservative Laws) is implemented in the presented FV (Finite Volume) model. The presented model is used to simulate different dam-break flow events to verify its capability. Four test cases are presented in this paper. In the first test case, a 1-Dimensional (1D) dambreak flow is simulated over a rectangular channel with different slope limiters of the FV model (namely Godunov, Superbee, Minmod, van Leer, and van Albada). The second test case consists of a simulation of shallow water discontinuous dam-break flow over a dry-downstream bed channel. The third test simulates the shallow water dam-break flow with the existence of bed slope and bed shear stress. Finally, in the last test, the HLL-MUSCL model used in this paper and some other solver models used in literature are compared against the referred exact solution in dam-break flow application. The presented HLL-MUSCL scheme is found to give the best agreement to the exact solution.


Introduction
In the rise of the awareness on hazards that are caused by dam-breaking events on different hydraulic flow systems, studies using computer modelling environment of such events have been given serious attention in the recent years. The virtual-prototyping prediction tools, that are modelled, give failure analyses to dams in different size and geometry, and a warning system can then be equipped by combination of those prediction tools with satellite data capturing. In statistical record, there were more than 200 dam failures reported in large dam-systems around the world in the last century, and they had caused a loss of more than 8,000 lives besides millions dollar of loss [1].
In the numerical simulation of dam-break flow, there are a lot of schemes and solvers that have been employed for its prediction. However, the most popular schemes are usually associated with the use of Riemann solvers, as they were found to be more consistent in handling the shock-capturing event on most rapidly varied free-surface flows [2].
Jha et al. [3] proposed the use of flux-differencesplitting (FDS) Lax-Wendroff scheme in solving the sudden closure of sluice gate and flood wave on the dam-break flow system. Seaid [4] adopted the total-variation-diminishing (TVD) Runge-Kutta scheme to solve the wet bed sluice gate dam-break and circular dam-break problems. Tseng and Chu [5], and Vincent et al. [6] had both used the TVD MacCormack scheme to investigate wet bed dambreak application; besides, Tseng [2] had also look into the same problem using Roe and essentially non-oscillatory (ENO) scheme. All of them were utilising approximate Riemann solver in their solution scheme formulation, and showed reasonable results on simulating dam-break flow.
In this paper, the shallow water equations are used as they show supremacy in representing the large water flow system; although their numerical simulation can be difficult due to their requirement of special treatment for achieving the conservative properties to avoid any spurious oscillation during their simulation [7,22]. The solution scheme employed by this paper is the HLL approximate Riemann solver with MUSCL scheme. Compared with the exact Riemann solver that considers multiwave constant stages in the smallest unit cell of finite volume method (FVM), the approximate _________________________ *corresponding author. Email: jhpu@nu.edu.kz Riemann solvers are more practical in terms of computational cost [8][9][10]. The use of the MUSCL scheme in the HLL approximate Riemann solver has two distinctive advantages, 1) it always keeps the HLL solver in its conservative form, and 2) it allows the HLL solver to achieve a second order of accuracy in its computation. A two-stage Hancock (predictor-corrector) scheme has been cooperated into the HLL-MUSCL solver in this paper. This two-stage Hancock scheme was proven to be efficient in numerical modelling of dam-break flow, as showing by literature [10][11][12].

Governing Equations
By assuming a water flow in the depth-averaged form, and by ignoring the effects of wind shear and Coriolis force, the governing equations of shallow water approach can be constructed as equations (1)(2)(3)(4)(5)(6)(7)(8). In this paper, both the St Venant and shallow water equations are used for the test cases simulation.

St Venant Equations
Before the existence of complicated computer schemes and applications for dam-break flow simulation, St Venant equations are the most popular approach used to understand and investigate those applications due to its simplicity and effectiveness to represent the flow in a single streamwise flow direction [6,13].
, and, In equation (2), φ = geopotential, where g h φ = ⋅ ; h = water height; g = gravitational acceleration; and u = velocity in x-direction. ox S and fx S are the bed slope and bed shear stress respectively in xdirection, and they are expressed as: and,

Shallow Water Equations
In the recent years, research of different hydraulic system, like oceans, rivers and dams, are often related to the use of shallow water approach as the vertical properties of those systems are usually very small and ignorable when compared with the properties of the other two directions of vector.
In shallow water formulation, y G % , oy S , and fy S represent the numerical flux vector, bed slope and bed shear stress in the y-direction. The oy S and fy S can respectively be expressed as:

Riemann Solver Scheme
In the early development of the representation of shallow water flow, studies were often carried out on the analytical basic. Their solutions are usually hard to comprehend or explain. In 1981, Roe [14] had proposed an idea for CFD (computational fluid dynamics), in which he proposed that for any conservative fluid flow, its solution can be found by assuming a numerical flux on its flow properties. The Riemann solver was used in his study by an approximate scheme, which was believed to be accurate in different flow conditions. In the more recent research, Toro [9] has proven that HLL-type approximate Riemann solver is one of the most accurate approximate Riemann solvers.
Godunov [16] has suggested that numerical fluxes of an inviscid fluid flow simulation can be obtained by solving the Riemann problem of the flow locally in each divided numerical cell or unit of consideration under a meshed flow system. He suggested that the difference properties at each left and right side of the divided numerical cell is the flux that changes the flow properties locally. Even through, the idea of Godunov to solve a flow by considering the flux difference on the left and righthand side of the FVM solution cell was theoretically founded, Godunov scheme is too idealised to constraint the fluxes into the actual value of right or left side solution of a numerical cell. Inspiring by the weakness of Godunov scheme, Harten et al. [17] had suggested a new Godunov-type scheme, called Godunov-type HLL-Riemann approximation scheme. The scheme is proven to be robust in most of the flow simulation. Besides, the HLL-approximate Riemann solver with Godunov-type scheme of numerical flux also possess with an upgradeable order of accuracy according to different applications.

HLL-Approximate Riemann Solver
In HLL approximate Riemann solver, wave speeds, L S and R S are used to separate the HLL approximate Riemann solution cells into three regions (left, right and * regions), hence 3 U % can be getting from this solution ( L U % , R U % , and * U % ), see Figure 1.

Fig. 1. HLL Approximate Riemann Solution Cell
The in viscid solution of the indicated cell can be written in conservation law as And by applying the in viscid solution of equation (9) into the control volume of AOBCEDA (in Figure 1), we can get where, subscripts L and R represent left and right sides of the solution cell, and F % = numerical flux vector of the solution cell.
By assuming the control volume of AOEDA and OBCEO are acting in the same behaviour as the control volume of AOBCEDA, we can reform equation (10) into two separate equations for the left cell (AOEDA) and the right cell (OBCEO) using the common numerical flux at * region, * F % .
( ) and, ( ) By eliminating * U % at equations (11)(12), we can get a useful expression for * F % , and it can be expressed in the conditions of: where, L S and R S can be expressed respectively as: and, and, In order for HLL approximate Riemann solver to compute a dry bed dam-break flow, extra conditions is needed for R S and L S , and the conditions are For the left dry bed ( For the right dry bed (

MUSCL Scheme
The robustness of MUSCL-type (monotonic upwind schemes for conservative laws) scheme is demonstrated by Marques and Pereira [15] ( ) and, ( ) where, and, i represents the space step.
There are 5 types of different slope limiter, Ψ , which is presented in this paper, Godunov, Minmod, Superbee, van Leer and van Albada slope limiters (see Table 1). Those slope limiters are tested in Section 5.1 for their accuracy to represent a shock-capturing dam-break event.

Time Integration
In this paper, a Hancock-two stage predictorcorrector scheme is utilised to update i U % in time step. It has the ability to maintain stability and to achieve second order of accuracy in time domain for the flux-limiting scheme.

Predictor
Step Step ( ) where, Ω = cell volume, and N represents the time step.
For t Δ , a stability criterion is followed to ensure the utilised t Δ does not exceed its maximum allowable limit.

t CFL c
where, ⋅ q s = where, m = last space step in computation.

Source Terms
Hu et al. [10,12], Mingham, and Causon [11] have consistently suggested that the source terms of any flow simulation cause no much of numerical instability in its effect towards the solutions of shallow water equations (for new source terms treatment study refer to Pu et al. [22]). Therefore, the complex solution of source terms is not necessary here. In this paper, a more direct approach compare to the inviscid terms solver, is used for solving source terms.
In the source terms simulation of a flow, HLL approximate Riemann solver is more powerful than most of the other approaches, like Roe approximate Riemann solver, which Jacobian matrix is involved in their solvers. It is because Jacobian matrix is having difficulty in handling distorted mesh of a flow system, and hence the more complicated solver of the source terms is required (see Delis [20]).

Numerical Simulations and Results
In this paper, four test problems are investigated to verify the numerical scheme developed in the previous sections. In the first test problem, a 1D dam-break flow is tested with different slope limiters of the FV model (namely Godunov, Superbee, Minmod, van Leer, and van Albada) to determine their accuracy. The second test problem is used to investigate the shallow water discontinuous dam-break flow over a drydownstream bed channel. As a continuation, in the third test the shallow water dam-break flow is simulated with the existence of source terms (bed slope and bed shear stress). Finally in the last test, the proposed HLL-MUSCL solver and two other solvers proposed by literature , and Petrov-Galerkin [21] solvers) are compared against the exact solution in a dam-break flow application.
The CFL number used in all 1D test problems is 0.95. It demonstrates the numerical stability of the proposed solver. The CFL number for the shallow water test problems is set to be about 0.3-0.4, in order to avoid numerical difficulties during simulation.

Test 1: Slope Limiters Comparison
In this test, a 1.0m long channel is simulated with upstream reservoir headwater, h US = 1.0 m, and downstream tail water, h DS = 0.01 m, separated by a thickness-free dam located at the middle of the channel. A schematic of the dam-break problem is shown in Figure 2. The channel is flat, straight and frictionless. The water in the channel, which initially at rest, is allowed to flow through the dam when it is instantly and completely removed at 0 t = . In this simulation, 100 computational cells with simulation time of 0.25s are used.

. Water Height Simulation with Different Slope Limiters
The simulation is conducted by using different slope limiters on the FV model. Besides, an exact solution derived from the analytical approach is also used to compare with the different slope limiter results, since it is commonly used as the benchmark solution in numerical research. By referring to Figures 3-4, one can see that, in general, all slope limiters (Godunov, Superbee, Minmod, van Leer, and van Albada) are simulated results of good agreement with the exact solution. However, van Leer and van Albada slope limiters are showing the best agreement with the exact solution (refers to the enlarged portions at Figures  3-4). Further comparing those 2 slope limiters, one can identify that van Leer slope limiter is showing a better result than the van Albada slope limiter (refers to the enlarged portion of Figure 4). Due to the supremacy of the van Leer slope limiter over other slope limiters, it is used in the rest of the test problems that are discussed in the coming sections.

Test 2: Shallow Water Discontinuous Dry Dam-Break Flow
Shallow water shock wave is a natural turbulent water flow event to test the stability and capability of a proposed numerical model, as it requires highresolution shock-capturing numerical schemes to avoid any numerical instability. In this test problem, a 200 m × 200 m rectangular channel is considered with a 15m thick dam at the middle. The dam is having a partial dam-break or rapidly opening of sluice-gap at the starting of the simulation (refers to Figure 5). Initially, the upstream headwater is set to 10.0m, whereas downstream tailwater is set to dry. In the numerical simulation, 42 × 42 computational cells are used, where Δ x = Δ y = 5 m. The simulation time is 6.0 s.

Fig. 5. Dam Break Layout for Inviscid Shock Wave
The results of water height, u-momentum and vmomentum are presented at Figures 6-8 respectively. At the instant of dam breaking, water is released through the 75m-wide sluice-gate. A forward bore wave is formed at downstream of the channel; while negative depression wave is generated at the upstream. From Figure 6, one can observe that both bore and depression waves are well-formed in the simulation of the proposed FV model. The results of this simulation are having good agreements with those presented in various literatures like, Valiani et al. [18], Caleffi et al. [19], and Mingham, and Causon [11], who used this same test case in their numerical studies. From Figures 7-8, one can observe that the u-and vmomentum are at their maximum absolute magnitudes at the position of the sluice-gate. It is caused by the changing from bore to depression wave of water flow in the sluice opening. During the process of changing depression to bore wave, the water flow is associated with a great deal of velocities change in x and y directions, hence causing the u-momentum and v-momentum to be in their maximum absolute magnitudes (the results are also agreeable to the numerical tests done by Valiani et al. [18], and Seaïd [4]).

Test 3: Discontinuous Dam-Break Flow with Source Terms
As a continuation of Test 2 in Section 5. conditions are similar to those presented at Figure  5, except the 75m sluice-gate opening is situated 30m from the left of the dam. The change is able to test the ability of meshing boundary conditions to assign the dam walls. In Figure 9, the simulation results of sloping frictional dam-break channel over a dry downstream bed are presented in different times. Compared to the inviscid dam-break flow in Test 2, where the water flow spreads laterally and radially when it reaches downstream of the channel; in this test problem, water flows rapidly in an decline manner at downstream of the channel. Caused by the bed slope and friction inputted into the simulation, water flow is not spreading wide at the downstream, but being pulled by gravitational force to flow directly to the end of the slope (refers to t =10 s diagram in Figure 9).  Figures 10-12, we can observe that in general the proposed HLL-MUSCL model is presenting a closer result to the exact solution compared to Roe-HLLE [20] and Petrov-Galerkin [21] models. The water height diagram presented in Figure 10 shows the Roe-HLLE solution is depleted away from the exact solution at the start of depression wave and the end of bore wave formation, and the Petrov-Galerkin solution shows a huge numerical oscillation at the depression wave generated. In contrast, the proposed HLL-MUSCL model shows the best agreement to the exact solution. In Figure 11, the water discharge distribution through the dam-break channel is presented, and the result from the proposed HLL-MUSCL model is compared with the results from the Roe-HLLE model and the exact solution. Again, the bore and depression wave formation of the Roe-HLLE model simulation is depleted away from the exact solution if compared to the proposed HLL-MUSCL model simulation. If one refers to Figure 12, the single peaking jump results are detected in the Froude number distribution across flow domain. The single peaking jump is caused by the supercritical condition at the end of the bore wave formation. The location of the peaking jump at the Froude number for the proposed HLL-MUSCL solution (which happens at the location about x = 1750 m) is also presented more accurately than the Roe-HLLE solution (which happens at the location about x =1600 m), if compare with the exact solution (which happens at the location about x = 1800 m). In short, the proposed HLL-MUSCL model is well agreed with the exact solution, and it is a highly accurate and robust approximate Riemann solver even a high CFL number is used.

Conclusions
In this study, a HLL-MUSCL approximate Riemann solver is successfully implemented in a FV model. The proposed approximate Riemann solver is robust, efficient and able to avoid inaccurate numerical oscillations or depletions when tested against the exact solution. The proposed model is accurate in simulating the shallow water discontinuous shock wave across any dam-break channel with dry sections. It is also capable in handling the unsteady supercritical condition of a flow. Besides, the proposed model is also having the capability in simulating any frictional channel with steep slope. The capabilities of the proposed model are demonstrated by the simulations of different flow test problems in dambreak application.
"Environmental Assessment of Sediment Pollution Impact on Hydropower Plants". The authors acknowledge the research support from Prof Stefaan Simons and Prof Sergey Mikhalovsky at Nazarbayev University, Kazakhstan during the preparation of the manuscript. The first author would also like to thank Dr Khalid Hussain and Prof Simon Tait at the University of Bradford, UK for their guidance on the SWE modelling during his PhD study.