(1)
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
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
Abstract: This paper describes a finite volume TVD-MacCormack scheme for the computation of two-dimensional open channel flows with abrupt changes. The algorithm modified the widely-used MacCormack scheme by implementing a conservative dissipation step to avoid the unphysical oscillation in the numerical solution. A series of simulation, such as oblique hydraulic jump, circular dam-break and two-dimensional dam-break are carried out to demonstrate its robustness and stability in capturing strong gradients and discontinuities in open channel flows.
Keywords: finite volume, MacCormack scheme, open channel flows
The prediction of shallow water flows with abrupt changes, such as hydraulic jumps, bores and surges are of great interest to hydraulic engineers. The variations of water depths and velocities in these extreme events are important parameters for the design of hydraulic systems and for flood control operations. However, due to the strong gradients inherent in the problems, most of the traditional simulation models display spurious oscillations at the shock fronts.
More recently, a number of shock-capturing schemes for solving the two-dimensional Saint Venant equation have been proposed (Fennema and Chaudhry, 1990; Alcrudo and Garcia-Navarro, 1993; Zhao et al., 1996; Mingham and Causon, 1998; Louaked and Hanich, 1998). But all these numerical models restrain the oscillatory behavior in the solution at the expense of increased algorithm complexity. In this study, the TVD-MacCormack scheme proposed by Garcia-Navarro et al. (1992) is modified and extended to a two-dimensional finite volume approach.
The governing equations for the unsteady two-dimensional shallow water flows can be written in conservative form as
(1)
where
(2)
The source term S includes the effects of bottom slopes and friction
(3)
where h is the water depth, u and v are the velocity components in the x- and y-directions, respectively, g is the gravitational acceleration, Sox and Soy are the slopes of channel bottom in the x- and y-directions, n is the Manning roughness coefficient.
By using the mean value theorem and divergence theorem, the governing equations represented in Eq. (1) can be integrated around the surface area s of each grid cell in the following way
(4)
where Qm and Sm are the mean quantity referred to the center of the volume V. In order to discretize Eq. (4), the surface integral is approximated as the sum of numerical fluxes over the four sides of the grid cell.
(5)
In order to avoid the numerical oscillation, a three-step procedure to include the total variation diminishing (TVD) dissipation term is employed
(6)
(7)
(8)
where the
superscript p and c denote the variables at predictor and
corrector steps, respectively;
,
,
are the flux vectors at the center
of the upstream, mid-stream and downstream volumes, respectively;
,
represent the surface area vectors
at the upstream and downstream sides of the volume.
The D term is defined as
(9)
The subscript
denotes the intermediate state
between grid points (i, j) and (i+1, j). The mean values
of velocity and water depth can be calculated as
(10)
(11)
(12)
The component
of
in Eq. (9) is defined as
(13)
where
is the mesh ratio. The
characteristic variable
is
(14)
The purpose of
the flux limiter function
in Eq. (13) is to supply
artificial dissipation when there is a discontinuity or a strong gradients,
while adding very little or no dissipation in regions of smooth variation.
In order to satisfied the numerical stability requirement, the time step is determined based on the Courant-Friedrich-Lewy criterion
(15)
where Cr is the Courant number. In this study, the value of Cr can be set to 1.0, usually 0.8.
An oblique hydraulic jump is a
phenomenon generated by the interaction between a supercritical flow and a
converging wall. The design of simulation condition in this study is similar to
that of Alcrudo and Garcia-Navarro (1993) and Zhao et al. (1996). The angle
between the converging wall and the stream-wise direction is set to be
= 8.95o. The upstream
condition is a uniform supercritical flow with a water depth of 1.0 m, flow
velocity u = 8.57 m/s and v = 0.0 m/s. This yields a upstream Froude no. Fr =
2.74. The flow domain is discretized into a 40
30 non-rectangular mesh. Based on the geometry and upstream condition, an
analytic solution can be found (Chaudhry, 1993). The solution are water depth hd
= 1.5 m, velocity
= 7.95 m/s at the downstream and
= 30o for the angle of the oblique jump.
Figure 1 shows the water depth contours predicted by the TVD-MacCormack scheme. In order to illustrate the non-oscillatory feature of proposed scheme, the result predicted by a finite volume MacCormack scheme is also plotted in Figure 1 for comparison. The contrast between them clearly indicates that the shock-capturing ability of the TVD-MacCormack scheme is better than that of traditional MacCormack scheme.
Another interesting test case for the proposed scheme is the breaking of a circular dam and the time evolution of the shock waves. This problem can be used to examine the capability of the numerical schemes to conserve symmetry. The initial condition is comprised of two regions of still water separated by a cylindrical wall of radius 11 m. The water depth inside the wall is 10 m and outside the dam 1.0 m. A circular mesh system consisting of 50 cells following the tangential direction and 25 cells of 1.0 m length along the radial direction is used. The simulation result obtained by the TVD-MacCormack scheme at time t = 0.69 sec is plotted in Figure 2. On the other hand, the traditional MacCormack scheme breaks down after time t = 0.22 sec for this case.
In this section, opening of a
non-symmetrical sluice gate (or partial breach of a dam) is simulated to examine
the robustness of the TVD-MacCormack scheme. The computation domain comprised of
a 200 m
200 m channel and is discretized into 40
40 rectangular mesh. At the instant of gate-opening, water is released into the
downstream side through a 75 m-wide gate. The water depths in the reservoir
= 10 m and in the tailwater
= 5.0 m are imposed as the initial
condition. The water depth contours at time t = 7.2 sec predicted by the MacCormack and TVD-MacCormack schemes
are compared in Figure 3. As can be seen in this comparison, the TVD-MacCormack
scheme does not induce numerical oscillation like that of the MacCormack scheme.
The three-dimensional view of the water surface elevation is shown in Figure 4.
Figure 5 shows the water depth
contours (at time t = 7.2 sec) simulated by the TVD-MacCormack scheme for three
cases of depth ratio
= 10, 100, 1000. This figure indicates that the shock front of depth ratio
= 1000 has the fastest propagation speed and the least amount of lateral
spreading. This phenomenon also can be seen in Figure 6, which compares the
water elevations at the centerline of the gate (at time t = 5.0 sec) for depth ratio
= 10, 100, 1000. This is primarily because large depth ratio has massive head
water momentum, in turn it can generate a fast shock wave and has less time for
lateral spreading as it propagates downstream.
A two-dimensional TVD-MacCormack scheme based on the finite volume approach is developed to simulate unsteady open channel flows. This algorithm modified the MacCormack scheme by implementing a conservative dissipation step to avoid the spurious oscillation. The extra step was devised according to the theory of total variation diminishing (TVD) schemes.
A series of simulation, such as oblique hydraulic jump, circular dam-break, and opening of a non-symmetrical gate are carried out to demonstrate its ability in capturing strong gradients and discontinuities in open channel flows. Good agreement between the simulated results with the analytic solution has been obtained, thus indicating the accuracy of the TVD-MacCormack scheme and the efficiency of a finite volume implementation.
Acknowledgment
Financial support for the research reported in this paper was provided by the National Science Council, (grant No. NSC 88-2611-E008-001).
[1] Alcrudo, F. and Garcia-Navarro, P. (1993) “A high resolution Godunov-type scheme in finite volumes for the 2D shallow water equations”, Inter. J. Numerical Methods Fluids, Vol. 16, pp.489-505.
[2] Chaudhry, M.H. (1993) Open-Channel Flow, Prentice-Hall Inc., London.
[3] Fennema, R.J. and Chaudhry, M.H. (1990) “Explicit methods for 2D transient free surface flows”, J. of Hydraulic Engineering, ASCE, Vol. 116, No. 11, pp.1013-1034.
[4] 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.
[5] Harten , A. and Hyman, P. (1983) “Self-adjusting grid methods for one-dimensional hyperbolic conservation laws”, J. Computation Physics, 50, pp.235-269.
[6] Louaked, M. and Hanich, L. (1998) “TVD scheme for the shallow water equations”, J. of Hydraulic Research, IAHR, Vol.36, No.3, pp.363-378.
[7] Mingham, C.G. and Causon, D.M. (1998) “High-resolution finite-volume method for shallow water flows”, J. Hydraulic Engineering, ASCE, Vol.124, 6, pp.605-614.
[8] Zhao, D.H., Shen, H.W., Lai, J.S. and Tabios, G.Q. (1996) “Approximate Riemann solvers in FVM for 2D hydraulic shock wave modeling”, J. Hydraulic Engineering, ASCE, Vol. 122, 12, pp. 692-702.
Fig.1 Water depth contours of oblique hydraulic jump. (a) MacCormack scheme; (b) TVD-MacCormack scheme
Fig.2 Water depth contours of the circular dam-break case
Fig.3 Water depth contours for the partial dam-break case at time t = 7.2 sec (a)MacCormack(b)TVD-acCormack
Fig.4 Water surface elevation of the partial dam-break case
Fig.5 Water
depth contours for the partial dam-break case at time t = 7.2 sec. (a)
= 10, (b)
= 100, (c)
= 1000
Fig.6 Water depth at time t = 5.0 sec for the partial dam-break case