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

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]:𝜕ℎ𝜕𝑞𝜕𝑡+𝜕𝜕𝑥𝑞=0,𝜕𝑡+𝜕𝜕𝑥𝑞ℎ²+12𝑔ℎ²=𝑔ℎ𝑆⁰−𝑆_𝑓,(2.1) 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𝜕𝑧_𝑏𝜕𝑥=−𝑆⁰.(2.2) The friction slope 𝑆_𝑓 can also be also expressed with a relationship established for uniform flow, by using the Manning-Strickler formula [16] as follows:𝑆_𝑓=𝑞²𝑛²ℎ^10/3,(2.3) 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:𝜕ℎ𝜕𝑞𝜕𝑡+𝜕𝜕𝑥𝑞=0,𝜕𝑡+𝜕𝜕𝑥𝑞ℎ²+𝜕𝑧21_𝑏𝑔ℎ²+𝑔ℎ𝜕𝜕𝑧𝑥^𝑏=0,𝜕𝑡=0.(2.4) The source operator, on the other hand, deals with the nonhomogeneous part in the absence of flux terms:𝜕ℎ𝜕𝑞𝜕𝑡=0,𝜕𝑡=−𝑔ℎ𝑆^𝑓𝜕𝑧_𝑏,𝜕𝑡=0.(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.

2.2.1. Roe Scheme

Consider first the hyperbolic operator. The corresponding equations can be cast in the following matrix form:𝜕𝐔𝜕𝑡+𝜕𝜕𝑥𝐅(𝐔)+𝐇(𝐔)𝜕𝜕𝐔𝑥=0,(2.6) where 𝐔=⎡⎢⎢⎢⎢⎣𝑧ℎ𝑞_𝑏⎤⎥⎥⎥⎥⎦,𝐅(𝐔)=⎡⎢⎢⎢⎢⎢⎣𝑞₂𝑞ℎ+12𝑔ℎ²0⎤⎥⎥⎥⎥⎥⎦,𝐇(𝐔)=⎡⎢⎢⎢⎢⎣000000𝑔00ℎ⎤⎥⎥⎥⎥⎦.(2.7)

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Δ𝐔=³_𝑘=1Δ𝐅+𝛼_𝑘𝐊^(𝑘),(2.8)𝐇Δ𝐔=³_𝑘=1𝛼_𝑘𝜆̃_𝑘𝐊^(𝑘),(2.9) 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 𝑔ℎ:𝐊𝜆̃₁=𝑢̃−𝑔ℎ,𝜆̃₂=𝑢̃+𝑔ℎ,𝜆̃₃=0,₍₁₎=⎡⎢⎢⎢⎢⎣𝑢̃−10𝑔ℎ⎤⎥⎥⎥⎥⎦,𝐊⁽²⁾=⎡⎢⎢⎢⎢⎣𝑢̃+10𝑔ℎ⎤⎥⎥⎥⎥⎦,𝐊⁽³⁾=⎡⎢⎢⎢⎢⎢⎣𝑢₂̃10−𝑔ℎ𝑔ℎ⎤⎥⎥⎥⎥⎥⎦.(2.10) The surge strengths 𝛼_𝑘 are then obtained by linearising (2.8):𝛼₁=12⎡⎢⎢⎢⎣Δℎ−𝑔ℎ𝑢̃−𝑔ℎΔ𝑧^𝑏−ℎ𝑔ℎΔ𝑢⎤⎥⎥⎥⎦,𝛼²=12⎡⎢⎢⎢⎣Δℎ+𝑔ℎ𝑢̃+𝑔ℎΔ𝑧^𝑏+ℎ𝛼₃𝑔ℎΔ𝑢⎤⎥⎥⎥⎦,=𝑔ℎ𝑢²̃−𝑔ℎΔ𝑧^𝑏.(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 ℎ, 𝑢̃, and 𝑔ℎ:ℎ=√ℎ_𝐿ℎ_𝑅,𝑢̃=√ℎ_𝑅√ℎ_𝑅+√ℎ_𝐿𝑢_𝑅+√ℎ_𝐿√ℎ_𝑅+√ℎ_𝐿𝑢_𝐿,𝑔ℎ=12𝑔ℎ^𝑅+𝑔ℎ_𝐿.(2.12)

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 𝐔^𝑖𝑛+1 at the next time step is computed using the finite volume statement:𝐔^𝑖𝑛+1=𝐔^𝑛𝑖+Δ𝑡Δ𝑥𝐅^𝑖−1/2−𝐅_𝑖+1/2+Δ𝑡Δ𝑥14𝐇^𝑖−1+𝐇_𝑖+1𝐔_𝑖−1−𝐔_𝑖+1,(2.13) where 𝐅_𝑖−1/2 and 𝐅_𝑖+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𝐅_𝑖+1/2=12𝐅^𝑖+𝐅_𝑖+1−12³_𝑘=1𝛼_𝑘||𝜆̃_𝑘||𝐊^(𝑘),(2.14) 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 𝑥_𝑖−1/2<𝑥<𝑥_𝑖+1/2 of constant length Δ𝑥. The corresponding discharges 𝑞(𝑥) are represented by fluxes 𝑞_𝑖+1/2 sampled at the boundaries 𝑥_𝑖+1/2 of the intervals. For the continuity equation, time step from 𝑡 to 𝑡+Δ𝑡 is achieved using the classical finite volume statementℎ^{𝑖𝑡+Δ𝑡}=ℎ^𝑡𝑖+Δ𝑡Δ𝑥𝑞^{HLL𝑖−1/2}−𝑞^HLL_𝑖+1/2,(2.15) where 𝑞^HLL_𝑖+1/2=𝑆^𝑅𝑆_𝐿−𝑆_𝑅𝑞^𝑡𝑖−𝑆^𝐿𝑆_𝐿−𝑆_𝑅𝑞^𝑡𝑖+1+𝑆^𝐿𝑆_𝑅𝑆_𝐿−𝑆_𝑅ℎ^𝑡𝑖+1−ℎ^𝑡𝑖(2.16) is the standard HLL flux function [20, 21]. In the above formula, the left and the right surge speeds 𝑆_𝐿 and 𝑆_𝑅 are estimated from𝑆_𝐿=min𝑆_min,𝑖,𝑆_min,𝑖+1,0,𝑆_𝑅=max𝑆_max,𝑖,𝑆_max,𝑖+1,0,(2.17) where 𝑆_min and 𝑆_max are the surge speed bounds.

Besides, the LHLL scheme [19] is used to discretize the momentum equation, which associates with the nonconservative product 𝑔ℎ𝜕𝑧_𝑏/𝜕𝑥:𝑞^{𝑖𝑡+(1/2)Δ𝑡}=𝑞^𝑡𝑖+Δ𝑡Δ𝑥𝜎^{𝑅𝑖−1/2}−𝜎^𝐿𝑖+1/2,(2.18) where 𝜎^𝐿𝑖+1/2=𝜎^HLL_𝑖+1/2−𝑆^𝐿𝑆_𝑅−𝑆_𝐿𝑔ℎ^𝑡𝑖+ℎ^𝑡𝑖+12𝑧^𝑏𝜎+ℎ^𝑡𝑖+1−𝑧_𝑏+ℎ^𝑡𝑖,(2.19)_𝑅𝑖+1/2=𝜎^HLL_𝑖+1/2−𝑆^𝑅𝑆_𝑅−𝑆_𝐿𝑔ℎ^𝑡𝑖+ℎ^𝑡𝑖+12𝑧^𝑏+ℎ^𝑡𝑖+1−𝑧_𝑏+ℎ^𝑡𝑖(2.20) are lateralised corrections to the standard HLL flux:𝜎^HLL_𝑖+1/2=𝑆^𝑅𝑆_𝐿−𝑆_𝑅𝜎^𝑡𝑖−𝑆^𝐿𝑆_𝐿−𝑆_𝑅𝜎^𝑡𝑖+1+𝑆^𝐿𝑆_𝑅𝑆_𝐿−𝑆_𝑅𝑞^𝑡𝑖+1−𝑞^𝑡𝑖,(2.21) in which 𝜎=(𝑞²/ℎ)+(1/2)𝑔ℎ². 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:𝜕𝑞𝜕𝑡=−𝑔ℎ𝑞²𝑛²ℎ^10/3.(2.22) Using an implicit backward Euler scheme, (2.22) is discretised as:𝑞^{𝑖𝑡+Δ𝑡}=𝑞^{𝑖𝑡+(1/2)Δ𝑡}−Δ𝑡𝑔ℎ^−7/3𝑛²𝑞^{𝑖𝑡+Δ𝑡}².(2.23) Using the first component of the source operator, 𝜕ℎ/𝜕𝑡=0 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:Δ𝑡=CFLΔ𝑥max(|𝑢|+𝑐),(2.24) 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 𝐻_𝑟/𝐻_𝑡=2 and 𝐻_𝑟/𝐻_𝑡=1000, respectively; the total length of the channel is 1000 m; the discretisation interval is Δ𝑥=10m; the Courant number is set to CFL = 0.9. The results for the Roe and HLL schemes at time 𝑡=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 ℎ=10 m and ℎ₀=1 m, respectively; the total length of the channel is 2000 m; Δ𝑥=2m, and CFL = 0.9. Figure 2 indicates the reasonable results for 𝑡=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.

(a)

(b)

Figure 1 

Idealized dam-break problem: (a) initial depth ratio 𝐻𝑟/𝐻𝑡=2; (b) initial depth ratio 𝐻𝑟/𝐻𝑡=1000.

Figure 2 

Increasing accuracy by decreasing the grid size to Δ𝑥=2 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 𝐻₁=0.305 m and downstream of the dam 𝐻₂=0.0 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 𝑡=10 sec; (b) time evolution of water depth at 𝑥=70.1 m; (c) time evolution of water depth at 𝑥=85.4 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 𝑥=0.4 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 𝑥=1.2 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 𝑥=0.4 m, 0.8 m, and 1.2 m. The simulated grid number is 160 × 60 (Δ𝑥=Δ𝑦=0.01 m) that there are 9600 cells. The Manning coefficient is 𝑛=0.01. 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 𝑥=0.4 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 𝑥=1.2 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 (𝑡=0.725 sec) also corresponds to the simulated result, Figures 5(e) and 5(f). Finally, the water flow reaches the boundary of the tank at 𝑡=0.925 sec, which then forms surges crashing the columns and returning to the gate at 𝑡=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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 5 

Comparisons of experiments and numerical simulations—(a) and (b) when 𝑡=0.56 sec, dam-break flow arrives to the columns; (c) and (d) when 𝑡=0.625 sec, dam-break flow moves to the middle of the columns; (e) and (f) when 𝑡=0.725 sec, dam-break flow surrounds the columns and forms vortexes; (g) and (h) when 𝑡=0.925 sec, dam-break flow touches the boundary of tank; (i) and (j) when 𝑡=1.75 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.