1D and 2D Numerical Modeling for Solving Dam-Break Flow Problems Using Finite Volume Method

The purpose of this study is to model the flow movement in an idealized dam-break configuration. One-dimensional and two-dimensional motion of a shallow flow over a rigid inclined bed is considered. The resulting shallow water equations are solved by finite volumes using the Roe and HLL schemes. At first, the one-dimensional model is considered in the development process. With conservative finite volume method, splitting is applied to manage the combination of hyperbolic term and source term of the shallow water equation and then to promote 1D to 2D. The simulations are validated by the comparison with flume experiments. Unsteady dam-break flow movement is found to be reasonably well captured by the model. The proposed concept could be further developed to the numerical calculation of non-Newtonian fluid or multilayers fluid flow.


Introduction
Earthquakes or heavy rainfall usually caused more than a dozen landslide dams to form across Taiwan streams, temporarily impounding large volumes of water after the Chichi Earthquake in 1999 and Typhoon Morakot in 2009 1 . Once formed, these natural dams were highly exposed to catastrophic failure. Partial or complete failure can lead to severe flooding downstream and possibly trigger further floods or debris flows such as Shiaolin landslide events in 2009 2 . For realizing the dam formation process and evaluate the potential consequences of subsequent failure, it is important to be able to model the dynamics of dam-break flows, such as computational hydraulics and laboratory experiments.
Computational Hydraulics is regarded as an important technology, which utilizes numerical methods for solving the governing equations and discusses the relationship between the flow field and the change of water depth. The commonly used numerical methods such as finite difference method, finite element method, method of characteristics, and finite volume method have been studied.
The first computer-based simulation model for shallow water flows was finite difference method FDM , which is still widely applied at present 3 . According to Taylor series, FDM is an approximate numerical solution directly turning shallow water equations into algebra questions. Based on the typical numerical theory, FDM was developed early and presented high processing efficiency that it was simple and easily accepted. In order to enhance the calculation accuracy, FDM requires more simulated time for calculation and is rather unstable. According to the past research applying various numerical methods to solving such problems, Liao et al. 4 also applied the commonly used finite difference method FDM . When catching shock wave, discontinuous numerical oscillation appeared in the numerical processing that the total variation diminishing scheme TVD was utilized as the research method 5, 6 . It remained two-order accuracy of time and space for the optimal solution for unsteady flows.
Finite element method FEM was first applied to structural mechanics. With the development of computers in 1970s, it was applied to computational hydraulics 3 . In finite element method, the computing zone is divided into several nonoverlapping and connected individuals; basis functions are selected from each element for linear combinations so as to approach the true solution of elements. A large system of linear equations is required for each time-point standard FEM is implicit so that the explicit FEM could be utilized for enhancing the efficiency. In spite of the fact that finite element method could solve irregular zones, it requires more time to solve matrix equations. In this case, parallel computing or specific solutions are required for saving the calculation time. FEM therefore has not been widely applied to hydraulics computing. Idelsohn et al. 7 suggested applying meshless finite element method MFEM and particle finite element method PFEM to the approximate partial differential equation of fluid, where MFEM covered irregular shapes to approach the real situation. It therefore continuously disperses the moving particles due to gravity and the surface energy owing to the interaction with the contiguous particles , as well as density, viscosity, conductivity, and so on. The changes of particle velocity and position are also defined. For this reason, PFEM is considered as an advantageous and effective model to solve the surface problem and simplify the interaction between fluid structures in the papers 7, 8 . Wu and Chen 9 applied method of characteristics MOC to the calculation of kinematic wave equation. Method of characteristics does not show the drawback of numerical diffusion as in finite difference method, and the accuracy and convenience are better than other numerical methods. But the mathematical deduction of MOC is more complicated and it merely shows more restrictions on the side flow term. However, it makes more definitely physical meaning of the algorithm for method of characteristics. Shi and Liu 10 regarded method of characteristics as the most accurate numerical solution for hyperbolic partial differential equations. The basic theory was the one-order simulation of linear hyperbolic partial differential equations with two-dimensional characteristics space. The curves of the two characteristics and the correspondent characteristic arithmetic expressions are deducted. In terms of characteristic arithmetic expressions, they are proceeded by the numerical solutions such as velocity and water depth to discretize characteristic equations. The physical concept of MOC is definite, and the calculation accuracy of numerical analyses is high. The discontinuity of dam-break flows is rather difficult to solve with general difference method. However, the characteristic line still exists Journal of Applied Mathematics 3 and can be solved with method of characteristics. Nevertheless, the ratio of time-space is restricted by stable conditions that the minimum time is obtained. When the analysis is abundant, it would require longer calculation time.
Finite volume method FVM used to be applied to aviation and aerodynamics. For hydraulics, finite difference method and finite element method were more widely applied 3, 11 . However, unstructured grids are utilized for the grid computing with both finite volume method and finite element method that they are acceptable in irregularly natural channels. Besides, finite volume method and finite difference method present similar calculating speed, but faster than finite element method. The application of finite volume method has therefore been emphasized in recent years. With distinct directions, characters, and coordinates or different grids being the numerical calculation, each method would show different advantages. Both FVM and FEM divide calculation zone into several regular or irregular shapes of elements or control volume, but the calculating speed of FVM is faster than it of FEM. Bellos et al. 12 utilized finite volume method for calculation simulation and verified by the dam-break test and observe the surges generated by the dam-break in an extreme-wide flume.
In this paper, a numerical model by using the finite volume method is presented for simulating the dam-break flows. Meanwhile, the model uses a splitting to deal with the source term. The advantage of this approach is that it can develop more other computations, for example, mudflows, debris flows, or the aggradation and degradation of sediment-laden flows 13 , in the future.

Governing Equations
The Saint Venant equations are used to describe unsteady one-dimensional open-channel flow. Continuity and momentum balance are, respectively, written as 14, 15 : where h is the depth of flow above the rigid bed, q hu is the unit width discharge, and u is the mean velocity in the longitudinal flow direction. Letting z b denote the bed elevation above a reference datum, the slope S 0 can be written as The friction slope S f can also be also expressed with a relationship established for uniform flow, by using the Manning-Strickler formula 16 as follows: where n is the Manning coefficient, recalling that, for a very wide channel, the hydraulic radius is equal to the flow depth.

Hyperbolic Term
In this study, the above shallow water equations are solved numerically using a finite volume approach, well suited for transient problems such as dam-break flows. An operatorsplitting approach 17 is used to separately treat the hyperbolic and source components. The hyperbolic operator solves the homogeneous equations: ∂h ∂t The source operator, on the other hand, deals with the nonhomogeneous part in the absence of flux terms:

2.5
For both operators, the procedure outlined by documents 18-20 is adapted for the computation of geomorphic dam-break surges. In the hyperbolic operator, two schemes, including Roe and HLL, are used to deal with the partial differential equations, and an implicit backward Euler scheme to treat the source term which will be illustrated in the next section.

Roe Scheme
Consider first the hyperbolic operator. The corresponding equations can be cast in the following matrix form: Journal of Applied Mathematics 5 Fluxes at the interfaces between finite volumes are evaluated using the Roe scheme. Let U L and U R denote the cell states to the left and right of a given interface. A decomposition of the flux difference ΔF F R − F L is sought and simultaneously satisfies where the α k , λ k , and K k are the surge strengths, eigenvalues, and right eigenvectors of the Roe linearisation of Jacobian matrix ∂F/∂U. Following the Roe-Pike procedure, the eigenvalues and the eigenvectors are first written in terms of the averaged variables h, u, and gh:

2.10
The surge strengths α k are then obtained by linearising 2.8 :

2.11
Finally, the substitution of these expressions in 2.9 yields a set of algebraic equations which can be solved for the averages of h, u, and gh: This decomposition is exploited as follows in a finite volume framework. Let Δt be the time step and Δx the size of each cell of a one-dimensional grid. Denoting by U n i the state of cell i at time nΔt, the state U n 1 i at the next time step is computed using the finite volume statement: Journal of Applied Mathematics where F i−1/2 and F i 1/2 are the fluxes across the left and right interfaces of the cell. Based on the Roe wave decomposition derived previously, these fluxes are evaluated using the expression 14 hence the hyperbolic operator is fully specified by the Roe scheme.

HLL Scheme
Another scheme of hyperbolic operator adopted in the present work is an extension of the HLL scheme 20 widely used for shallow flows. Whereas the original HLL scheme applies to the equations in full conservation form, the momentum equation feature nonconservative product associated with pressure along the slope bed. This term is treated following the approach 19 . The source term associated with friction along the bed is further treated in the next section. Adopting a finite volume point of view, each depth h x is discretised into piecewise constant segments h i over finite intervals x i−1/2 < x < x i 1/2 of constant length Δx. The corresponding discharges q x are represented by fluxes q i 1/2 sampled at the boundaries x i 1/2 of the intervals. For the continuity equation, time step from t to t Δt is achieved using the classical finite volume statement is the standard HLL flux function 20, 21 . In the above formula, the left and the right surge speeds S L and S R are estimated from where S min and S max are the surge speed bounds. Besides, the LHLL scheme 19 is used to discretize the momentum equation, which associates with the nonconservative product gh∂z b /∂x: Journal of Applied Mathematics are lateralised corrections to the standard HLL flux: in which σ q 2 /h 1/2 gh 2 . In the above formulas, the surge speeds S L and S R are again estimated from 2.17 . The "lateralised flux correction" approach leading to the statements 2.18 -2.20 is presented by Fraccarollo et al. 19 .

Source Term
Consider now the source term operator. Using the Manning-Strickler formula to specify the friction slope S f for the computation of clear water, the equation for the momentum source term can be written as: ∂q ∂t −gh q 2 n 2 h 10/3 .

2.22
Using an implicit backward Euler scheme, 2.22 is discretised as:

2.23
Using the first component of the source operator, ∂h/∂t 0 hence h t Δt i h t i . Thus, one can solve 2.23 for the unit width discharge of clear water at the next time step.
To advance the solution at each time step, the hyperbolic operator is first applied to obtain a partial update. These results are then used as the initial conditions for the source operator to yield the complete update. Since the hyperbolic update is explicit, stability of the scheme is subject to the CFL condition on the time step: where c is the wave celerity and u the velocity at any given grid point, and CFL is the Courant number. The value CFL should be smaller than 1 and it is found to be satisfactory in our case.

Extend to Two-Dimensional Model
The numerical method of two-dimensional model in this study applies the explicit square grid to discrete governing equations in finite volume method. It also applies the approach of one-dimensional model which deals with the Riemann solutions and with the Godunov-type scheme 20 . The calculations of numerical fluxes apply to HLL scheme, and LHLL scheme is utilized for calculating the nonconservative term including the slope bed 19 . Finally, Strang splitting is applied to calculate the source term with frictions 17 . Since this study applies conservative finite volume method, the hyperbolic term and the source term in differential equations are separately processed so that the numerical model is largely improved. For example, the promotion of one-dimensional to two-dimensional grids utilizes dimensional splitting that it first calculates x-sweeps and then applies the result as the initial conditions to calculate y-sweeps 22 . Similarly, Strang splitting is also applied to the calculation of bed shear stress. In each time step, the hyperbolic term of shallow water equations is first calculated; the result is applied as the initial conditions for friction calculations. In this case, it is easy to expand, such as expanding from one-dimensional model to two-dimensional model, or expanding from clear water flows to mudflows or debris flows; merely the geometric mesh or source term is regulated.

Idealized Dam-Break Problem
To test the computation of the hyperbolic term, calculations for the clear water are first compared with the classical analytical solution of Stoker 23 for the sudden breach of a dam over a horizontal frictionless bed. The initial data used by Tseng et al. 24 are adopted: the ratios of water depths at the left H r and the right H t of the dam are set equal to H r /H t 2 and H r /H t 1000, respectively; the total length of the channel is 1000 m; the discretisation  interval is Δx 10 m; the Courant number is set to CFL 0.9. The results for the Roe and HLL schemes at time t 30 s are plotted in Figure 1. It is clearly shown in Figure 1 that the Roe and HLL schemes can capture the shock wave but not accurate enough. Hence, the second test refers to Wang and Shen 25 and uses the conditions: water depths to the left and the right of the dam are set equal to h 10 m and h 0 1 m, respectively; the total length of the channel is 2000 m; Δx 2 m, and CFL 0.9. Figure 2 indicates the reasonable results for t 30 s, 60 s, and 90 s, respectively, in nondimensional form. The accuracy can be increased by decreasing the grid size from these two figures.

Dam-Break Experiment in Rigid Channel
To test the numerical model for clear water, the simulation is compared with the laboratory experiments of the US Army Engineer Waterways Experiment Station 26 in Figure 3.

Comparison of 2D Numerical Simulation and Experiment
In order to prove the two-dimensional numerical model, this study designs a rectangular flume, which is a closed tank of 1. 0.925 sec, which then forms surges crashing the columns and returning to the gate at t 1.75 sec. In the entire simulation, the water crashing obstacles and forming surges as well as the flow movement due to small obstacles could be captured well by the present model.

Conclusions
This paper proposes the one-dimensional and two-dimensional numerical model with the finite volume method based on the shallow water equations. A key feature of the model is the use of an operator-splitting method to divide the governing equations into hyperbolic and source terms. This approach provides an easy method for using different numerical schemes such as the Roe and HLL schemes. With the conservative finite volume method, the model can be applied to various numerical methods or spatial dimensions. As described in the study, one-dimensional model could be easily developed into two-dimensional model so as to save the time for programming the codes. The comparisons between the experimental results and the model simulation are matched. In addition, it is very easy to develop another computation, for example, non-Newtonian fluid or multilayers fluid flows, from the present model in the future.