Simulation of Dam-Break Flows by an Improved TVD-MacCormack Scheme

Ming Hseng Tseng

Department of Information Management, Ling Tung College, Taiwan, China

Tel & Fax: 886-4-386-1709, E-mail: mhtseng@mail.ltc.edu.tw

 Chia R. Chu

Department of Civil Engineering, National Central University, Taiwan, China.

Tel & Fax: 886-3-427-9127, E-mail: crchu@cc.ncu.edu.tw

Abstract: A finite difference predictor-corrector TVD scheme is developed for the simulation of unsteady one-dimensional dam-break flows. The accuracy and robustness of the numerical scheme is verified with an analytic solution. Furthermore, a sensitivity study is carried out to investigate the efficiency and robustness of four different versions of the predictor-corrector schemes. It is found that the numerical scheme will have less computational error and higher efficiency when the direction of the predictor-corrector step is the same as the direction of the shock wave propagation.

Keywords:  finite difference, predictor-corrector TVD scheme, dam-break flows

1    Introduction

It is well known that numerical simulations for dam-break flows are difficult, because of the steep gradient inherent in the problem. Early works on dam-break simulations involved using the characteristics method and shock fitting techniques.  More recently, a number of shock-capturing schemes for solving the St. Venant equations have been proposed.  For instance, Fennema and Chaudhry (1987) used the MacCormack model along with the Beam-Warming scheme to simulate one-dimensional and two-dimensional dam-break flows.  Rahman and Chaudhry (1998) combined the MacCormack (1969) scheme and an adaptive grid technique to calculate such flows.  However, their models require an empirical artificial viscosity term to damp out the oscillation around the discontinuities.  The choice of this artificial viscosity requires both fine-tuning and adept judgement. 

In this study, the predictor-corrector TVD scheme proposed by Garcia-Navarro et al. (1992) is extended to include four versions of the predictor-corrector steps.  Also, the entropy correction function suggested by Harten and Hyman (1983) is used to eliminate the trial procedure for the entropy inequality condition.  Simulation results from the present model are verified with an analytical solution and experimental data.  Then, a sensitivity study is carried out to investigate the accuracy and efficiency of the predictor-corrector schemes. 

2    Governing Equation

Based on the assumptions of hydrostatic pressure distribution, unsteady open channel flows can be described by the St. Venant equations as

                                 (1)

                                                                (2)

where h is the water depth, and u is the depth-averaged velocity, respectively, and g is the gravitation acceleration.  The source term S includes the effect due to bottom slopes and friction, S_o is the slope of the channel bottom, and S_f is the slopes of the energy grade line,

                                                       (3)

n is the Manning roughness coefficient, and R is the hydraulic radius. 

         Eq. (1) can be expressed in quasi-linear form as

                                                       (4)

where A is the Jacobian matrix having two real eigenvalues

                                                                                                 (5)

where  is the wave celerity. 

3    Predictor-Corrector TVD Scheme

For the numerical solution of Eq. (1), the domain of integration is discretized as , , where is the size of a uniform mesh, and  is the time increment.  The present scheme is based on a predictor-corrector three-step procedure.  The first two steps of the MacCormack scheme are employed, and the third step equips the scheme with the total variation diminishing (TVD) dissipation property to keep the solution oscillation-free. 

The scheme can be written as

                                                                                (6)

                                                                                (7)

                                                        (8)

where  is the mesh ratio and the subscript p and c stand for the predictor and corrector step, respectively.  The  component in Eq. (10) is defined as

                                                      (9)

The subscript  denotes the intermediate state between grid points (i) and (i+1).  The entropy correction function  is defined as

                                               (10)

 is a small positive number whose value has to be determined for each individual problem.  For this study, a formula suggested by Harten and Hyman (1983) is used to calculate  in order to cut down on the trial and error process

                                                                   (11)

The characteristic variable  in Eq. (9) is defined as

                                               (12)

The purpose of the flux limiter function  in Eq. (9) is to supply artificial dissipation when there is a discontinuity or a strong gradient, while adding very little or no dissipation in regions of smooth variation. 

         Following the technique suggested by Roe (1981), the mean values for the velocity and the water depth can be calculated as follows:

                                                         (13)

Due to the fact that the MacCormack scheme incorporates forward and backward differences in separated predictor and corrector steps, an alternative would be to reverse the order of the predictor and corrector steps:

                                                                                (14)

                                                                                (15)

In this study, four possible combinations are adopted to investigate the effectiveness of the schemes for dam-break flow simulations.  The four different schemes are denoted as follows: the Forward-Backward (F-B) scheme combines Eqs. (6), (7) and (8); the Backward-Forward (B-F) scheme combines Eqs. (14), (15) and (8); the Forward-Backward-cycled (F-B-C) scheme combines Eqs. (6), (7), (8), (14), (15) and (8) in a cyclic manner; and the Backward-Forward-cycled (B-F-C) scheme combines Eqs. (6), (7), (8), (14), (15) and (8) in a cyclic manner.

In order to satisfy the numerical stability requirement, the time step is determined based on the Courant-Friedrich-Lewy criterion

                                                          (16)

where Cr is the Courant number.  A zero-gradient boundary condition is applied to the water depth at the inlet and outlet boundaries.   

4    Model Verification

In order to demonstrate that the proposed model is capable of describing dam-break flows, model predictions are compared with the laboratory dam-break experiments of Waterway Experiment Station (WES), U.S. Corps of Engineers (1960).  The experiments are conducted in a rectangular channel with a channel length of 122 m, a width of 1.22 m, a bottom slope of 0.005, and the Manning coefficient n = 0.0085. The initial water depth in the reservoir is  = 0.305 m and in the tail water  = 0 (dry bed). The flow domain is discretized into 122 grids with a uniform grid spacing = 1.0 m. Figure 1(a) shows a comparison of the computed and measured water depth variations along the centerline of the flume at time t = 10 sec.  Figures 1(b) and 1(c) compare the simulated and experimental data at the downstream distances of x = 70.1 m and x = 85.4 m, respectively.  The good agreement between the computed and measured water depth verifies that the proposed TVD scheme is capable of dam-break flow simulations.   

5    Results and Discussion

In this section, the efficiency and robustness of four versions of the predictor-corrector schemes are investigated. The design of simulation conditions is similar to the configuration used by Rahman and Chaudhry (1998).  The computation domain is comprised of a 2000 m long channel with a horizontal channel bottom.  The dam is located at the downstream distance x = 1025 m. The initial water depth in the reservoir is  = 10 m.  The flow domain is discretized into 80 uniform grids.  The time evolution of the water depth can be used to examine the shock-capturing capability of the numerical scheme. 

The water depths predicted by the four different schemes, namely the F-B, B-F, F-B-C, and B-F-C schemes, for a subcritical flow ( = 0.5 at time t = 60 sec) are shown in Figure 2.  The water depth, simulated by a finite-difference MacCormack scheme, is also plotted in Figure 2 for comparison.  This comparison clearly shows that all four TVD schemes can accurately predict the shock wave without spurious oscillations, while the traditional MacCormack scheme has an unphysical oscillation at the shock front.  Table 1 summarizes the relative error L₂ norm between the simulated and analytic solutions.  The L₂ norm is defined as

                                                                   (17)

where ,  are the simulated and exact solutions at grid point (i).  In general, the errors in all four TVD schemes are smaller than for the MacCormack scheme, but the differences among the former are not significant enough.  The Courant numbers Cr are set to be 0.95 in all the schemes for this simulation. The value of Cr is an indicator of the efficiency of the numerical scheme.  As the value of Cr increases, the efficiency of the numerical scheme increases.

A comparison of water depth for a supercritical dam-break flow (depth ratio  = 0.004) at time t = 60 sec is shown in Figure 3. Again, the agreement between the simulation results and the exact solution are satisfactory. The traditional MacCormack model however fails to calculate this same case. As can be seen in Table 1, the errors for the F-B and F-B-C schemes are smaller than the errors for the B-F and B-F-C schemes, in this case.  The highest values of Cr can be used for the B-F and B-F-C schemes are 0.25.  This implies that the F-B and F-B-C schemes have a better computational efficiency.

In order to investigate the directional dependence of the predictor-corrector schemes, the water depths in the reservoir and the tail water are reversed, i.e. the depth ratio  = 250 and the shock wave propagates from right to left. The computational results are shown in Figure 4 and Table 1. Table 1 indicates that the errors in the F-B and F-B-C schemes are larger than that in the B-F and B-F-C schemes for the reversed case.  Furthermore, the value of Cr suggests that the B-F and B-F-C schemes have a better efficiency. 

                         Table 1    Simulation parameters for dam-break flows.

This sensitivity study indicates that when the direction of the predictor-corrector step is the same as the shock wave propagation, the simulation will have less computational error and a higher efficiency.  Therefore, the choice of numerical scheme should depend on the direction of the shock wave propagation.  Table 1 also indicates that the error for the cyclic scheme is smaller that for the non-cyclic scheme.  However, contrary to the belief of most people, the F-B-C scheme performs better than the B-F-C scheme for a depth ratio of  = 0.004, and the B-F-C scheme has a better performance than the F-B-C scheme for a depth ratio of  = 250. This is because the accumulation of initial error caused by the directional dependence of the shock wave can not be attenuated by the cyclic procedure.   

6    CONCLUSIONS

A finite-difference predictor-corrector TVD scheme is developed to simulate one-dimensional dam-break flows. This model modifies the widely used MacCormack scheme by implementing a conservative dissipation step to avoid any unphysical oscillation in the vicinity of strong gradients.  The accuracy of the proposed model is verified with an analytic solution and experimental data.  Furthermore, the efficiency and robustness of four versions of the predictor-corrector schemes are investigated.  Based on the results of this study, several conclusions can be made: 

(1)         The accuracy and efficiency of the proposed TVD scheme compares favorably than that of the traditional MacCormack scheme for dam-break flow simulations. 

(2)         The numerical scheme has less computational error and a higher efficiency when the direction of the predictor-corrector step is in the same direction as the shock wave propagation.

(3)         The simulation results also demonstrate that the errors of the cyclic schemes are smaller that the non-cyclic schemes and that the accumulation of initial errors caused by the directional dependence can not be eliminated by the cycled step.     

Acknowledgment

Financial support for the research reported in this paper was provided by the National Science Council, Taiwan, China. (grant No. NSC 88-2611-E008-001)

References

[1]    Fennema, R.J., and Chaudhry, M.H. (1987) “Simulation of one-dimensional dam-break flows”, J. of Hydraulic Research, Vol. 25, pp.41-51

[2]    Garcia-Navarro, P., Alcrudo, F. and Saviron, J.M., (1992) “1-D open channel flow simulation using TVD-MacCormack scheme”, J. Hydraulic Engineering, ASCE, Vol.118, 10, pp.1359-1372

[3]    Harten, A. and Hyman, P. (1983) “Self-adjusting grid methods for one-dimensional hyperbolic conservation laws”, J. Comp. Physics, 50, pp.235-269

[4]    MacCormack, R.W. (1969) “The effect of viscosity in hyper-velocity cratering”, AIAA, Paper No.69-354, Washington, D.C.

[5]    Rahman, M. and Chaudhry, M.H. (1998) “Simulation of dam-break flow with grid adaptation”, Advances in Water Resources, Vol. 21, No.1, pp.1-9

[6]    Roe, P.L. (1981) “Approximate Riemann solvers, parameter vectors and different schemes.”, J. Computation Physics, 43, pp.357-372

[7]   US Corps of Engineers (1960) Flood resulting from suddenly breached dams, Miscellaneous paper, 2 (374), Report 1, US Army Engineer Waterways Experiment Station, Corps of Engineers, Vicksburg, Mississippi

Fig. 1    Comparison of computed and experimental results for a dam-break flow.
(a) at time t = 10 sec; (b) x = 70.1 m; (c) x = 85.4 m

Fig.2    Water depth variation for a depth ratio of    = 0.5

Fig. 3    Water depth variation for a depth ratio of    = 0.004

Fig. 4    Water depth variation for a depth ratio of    = 250