Journal of Applied MathematicsVolume 2012 (2012), Article ID 489269, 14 pagesdoi:10.1155/2012/489269

Research Article

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

Szu-Hsien Peng

Department of Spatial Design, Chienkuo Technology University, Number 1, Chieh Shou N. Road, Changhua City 500, Taiwan

Received 6 January 2012; Accepted 13 February 2012

Academic Editor: Shuyu Sun

Copyright © 2012 Szu-Hsien Peng. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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.

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

2. Numerical Model

2.1. 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 is the depth of flow above the rigid bed, is the unit width discharge, and is the mean velocity in the longitudinal flow direction. Letting denote the bed elevation above a reference datum, the slope can be written as The friction slope can also be also expressed with a relationship established for uniform flow, by using the Manning-Strickler formula [16] as follows: where is the Manning coefficient, recalling that, for a very wide channel, the hydraulic radius is equal to the flow depth.

2.2. 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 operator-splitting approach [17] is used to separately treat the hyperbolic and source components. The hyperbolic operator solves the homogeneous equations: The source operator, on the other hand, deals with the nonhomogeneous part in the absence of flux terms: 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.

2.2.1. Roe Scheme

Consider first the hyperbolic operator. The corresponding equations can be cast in the following matrix form: where

Fluxes at the interfaces between finite volumes are evaluated using the Roe scheme. Let and denote the cell states to the left and right of a given interface. A decomposition of the flux difference is sought and simultaneously satisfies where the , , and are the surge strengths, eigenvalues, and right eigenvectors of the Roe linearisation of Jacobian matrix . Following the Roe-Pike procedure, the eigenvalues and the eigenvectors are first written in terms of the averaged variables , , and : The surge strengths are then obtained by linearising (2.8):

Finally, the substitution of these expressions in (2.9) yields a set of algebraic equations which can be solved for the averages of , , and :

This decomposition is exploited as follows in a finite volume framework. Let be the time step and the size of each cell of a one-dimensional grid. Denoting by the state of cell at time , the state at the next time step is computed using the finite volume statement: where and 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 hence the hyperbolic operator is fully specified by the Roe scheme.

2.2.2. 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 is discretised into piecewise constant segments over finite intervals of constant length . The corresponding discharges are represented by fluxes sampled at the boundaries of the intervals. For the continuity equation, time step from to is achieved using the classical finite volume statement where is the standard HLL flux function [20, 21]. In the above formula, the left and the right surge speeds and are estimated from where and are the surge speed bounds.

Besides, the LHLL scheme [19] is used to discretize the momentum equation, which associates with the nonconservative product : where are lateralised corrections to the standard HLL flux: in which . In the above formulas, the surge speeds and 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].

2.3. Source Term

Consider now the source term operator. Using the Manning-Strickler formula to specify the friction slope for the computation of clear water, the equation for the momentum source term can be written as: Using an implicit backward Euler scheme, (2.22) is discretised as: Using the first component of the source operator, hence . 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 is the wave celerity and 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.

2.4. 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 -sweeps and then applies the result as the initial conditions to calculate -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.

3. Results and Discussions

3.1. 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 () and the right () of the dam are set equal to and , respectively; the total length of the channel is 1000 m; the discretisation interval is m; the Courant number is set to CFL = 0.9. The results for the Roe and HLL schemes at time  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  m and  m, respectively; the total length of the channel is 2000 m; m, and CFL = 0.9. Figure 2 indicates the reasonable results for  s, 60 s, and 90 s, respectively, in nondimensional form. The accuracy can be increased by decreasing the grid size from these two figures.

(a)

(b)

Figure 1: Idealized dam-break problem: (a) initial depth ratio ; (b) initial depth ratio .

Figure 2: Increasing accuracy by decreasing the grid size to  m.

3.2. 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. The WES’ experiments were conducted in a rectangular channel 122 m long and 1.22 m wide, with a bottom slope of 0.005 and Manning’s coefficient of 0.009. The dam was placed at the middle of the channel, giving the initial water depth upstream of the dam  m and downstream of the dam  m. In this simulation, the uniform grid spacing is 1.0 m and the Courant number CFL = 0.9. The agreement between the simulation and experimental results is satisfactory. The results also indicate that the finite volume method can capture the surges well and the numerical model has a good ability to deal with the “contact points” at locations where the depth reaches zero.

(a)

(b)

(c)

Figure 3: Comparisons with WES experiment: (a) water surface profile along the channel at  sec; (b) time evolution of water depth at  m; (c) time evolution of water depth at  m.

3.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.6 m (Length) × 0.6 m (Width) × 0.6 m (Height). A gate, located on the longitudinal  m, could be quickly removed for simulating the dam-break flow. The initial depth before gate is 0.15 m and the downstream depth is 0.01 m. Two square columns (of the length and width being 0.1 m and the height being 0.3 m) are placed on  m and are paralleled placed in the middle (with the distance 0.1 m) shown as Figure 4. Besides, in consideration of the natural channel not being as flat as the artificial construction, six small obstacles are placed on both sides for simulating the irregular banks. They are placed on the side walls at longitudinal  m, 0.8 m, and 1.2 m. The simulated grid number is 160 × 60  m) that there are 9600 cells. The Manning coefficient is . The boundary conditions are considered as closed boundary, that is, there is no flow-out at the sides. With the total simulation time 2 sec, it is simulated the water flow crashes the boundary after gate opening and returns to the initial location.

Figure 4: Initial conditions of experiment and 2D numerical simulation. The gate is located at  m, the dam-break modeling initial water depth 0.15 m, the downstream water depth 0.01 m, and 2 square columns (0.1 m × 0.1 m × 0.3 m) paralleled placed at  m with the distance 0.1 m.

The experimental results show that it takes about 0.56 sec that dam-break flow touches the square columns when the distance between the columns and the gate is 0.8 m. It presents the same time as the simulation time of the numerical model in Figures 5(a) and 5(b). With the flow movement of the small obstacles on both sides, it could also be simulated in the numerical model. When the time is 0.625 sec, the time of dam-break flow moving to the center of two columns is similar to it simulated in the model (see Figures 5(c) and 5(d)). The time of dam-break flow surrounding the two columns and forming vortexes  sec) also corresponds to the simulated result, Figures 5(e) and 5(f). Finally, the water flow reaches the boundary of the tank at  sec, which then forms surges crashing the columns and returning to the gate at  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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 5: Comparisons of experiments and numerical simulations—(a) and (b) when  sec, dam-break flow arrives to the columns; (c) and (d) when  sec, dam-break flow moves to the middle of the columns; (e) and (f) when  sec, dam-break flow surrounds the columns and forms vortexes; (g) and (h) when  sec, dam-break flow touches the boundary of tank; (i) and (j) when  sec, the flow returns to the location of gate.

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

Acknowledgments

Thanks for the parts of support from the Project (NSC 98-2622-B-270-001-CC3) of National Science Council. Special thanks are due to Dr. H. Capart, Professor of Civil Engineering at National Taiwan University, for reading the parts of the paper and making a number of helpful suggestions.

References

S. C. Chen, “Failure mechanism and disaster mitigation on landslide-dammed lakes,” Journal of Chinese Soil and Water Conservation, vol. 30, no. 4, pp. 299–311, 1999.

Z. Y. Feng, “The seismic signatures of the surge wave from the 2009 Xiaolin landslide-dam breach in Taiwan,” Hydrological Processes, 2011. In press

W. Y. Tan, Shallow Water Hydrodynamics, Elsevier, New York, NY, USA, 1992.

C. B. Liao, M. S. Wu, and S. J. Liang, “Numerical simulation of a dam break for an actual river terrain environment,” Hydrological Processes, vol. 21, no. 4, pp. 447–460, 2007.

D. Liang, R. A. Falconer, and B. Lin, “Comparison between TVD-MacCormack and ADI-type solvers of the shallow water equations,” Advances in Water Resources, vol. 29, no. 12, pp. 1833–1845, 2006.

M. H. Tseng, “Improved treatment of source terms in TVD scheme for shallow water equations,” Advances in Water Resources, vol. 27, no. 6, pp. 617–629, 2004.

S. R. Idelsohn, E. Oñate, F. Del Pin, and N. Calvo, “Fluid-structure interaction using the particle finite element method,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 17-18, pp. 2100–2123, 2006.

T. W. H. Sheu, C. Chiao, and C. Huang, “Development of a particle interaction kernel function in MPS method for simulating incompressible free surface flow,” Journal of Applied Mathematics, Article ID 793653, 16 pages, 2011.

C. M. Wu and S. C. Chen, “Simulating runoff using the method of characteristics with unsteady rainfall events,” Journal of Mechanics, vol. 21, no. 3, pp. 171–178, 2005.

H. D. Shi and Z. Liu, “Review and progress of research in numerical simulation of dam-break water flow,” Advances in Water Science, vol. 17, no. 1, pp. 129–135, 2006.

C. Hirsch, Numerical Computation of Internal and External Flows, Fundamentals of Numerical Discretization, vol. 2, John Wiley & Sons, New York , NY, USA, 1988.

C. V. Bellos, J. V. Soulis, and J. G. Sakkas, “Experimental investigation of two-dimensional dam-break induced flows,” Journal of Hydraulic Research, vol. 30, no. 1, pp. 47–75, 1992.

G. Simpson and S. Castelltort, “Coupled model of surface water flow, sediment transport and morphological evolution,” Computers and Geosciences, vol. 32, no. 10, pp. 1600–1614, 2006.

M. B. Abbott, Computational Hydraulics Elements of the Theory of Free Surface Flows, vol. 1 of Monographs and Surveys in Water Resource Engineering, Pitman Publishing, London, UK, 1979.

M. H. Chaudhry, Open-Channel Flow, Prentice-Hall, Englewood Cliffs, NJ, USA, 1993.

W. H. Graf, Fluvial Hydraulics, John Wiley & Sons, Chichester, UK, 1998.

E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, Berlin, Germany, 2nd edition, 1999.

H. Capart, Dam-break induced geomorphic flows and the transition from solid- to fluid-like behaviour across evolving interfaces, Ph.D. thesis, University of Louvain, Belgium, 2000.

L. Fraccarollo, H. Capart, and Y. Zech, “A Godunov method for the computation of erosional shallow water transients,” International Journal for Numerical Methods in Fluids, vol. 41, no. 9, pp. 951–976, 2003.

A. Harten, P. D. Lax, and B. van Leer, “On upstream differencing and Godunov-type schemes for hyperbolic conservation laws,” SIAM Review, vol. 25, no. 1, pp. 35–61, 1983.

L. Fraccarollo and E. F. Toro, “Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems,” Journal of Hydraulic Research, vol. 33, no. 6, pp. 843–864, 1995.

R. J. LeVeque, Finite volume methods for hyperbolic problems, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 2002.

J. J. Stoker, Water Waves, Interscience, New York, NY, USA, 1957.

M. H. Tseng, C. A. Hsu, and C. R. Chu, “Channel routing in open-channel flows with surges,” Journal of Hydraulic Engineering, vol. 127, no. 2, pp. 115–122, 2001.

Z. Wang and H. T. Shen, “Lagrangian simulation of one-dimensional dam-break flow,” Journal of Hydraulic Engineering, vol. 125, no. 11, pp. 1217–1221, 1999.

USACE, Flood Resulting from Suddenly Breached Dams, US Army Engineer Waterways Experiment Station, Vicksburg, Miss, USA, 1960.