Chiang-An Hsu
Civil & Hydraulic Engineering Research Center, Sinotech Engineering Consultants Inc.
171 Nanking E. Rd. Sec. 5, Taipei 105, Taiwan, China
Fax: +886-2-27655010; E-mail: cahsu@sinotech.org.tw
Abstract: A general and
powerful numerical model, named SEC-HY21, to compute the open-channel flows by
solving the depth-averaged, two-dimensional unsteady shallow water equations is
presented. The spatial solution algorithm used is a cell-centered multi-block
finite-volume method in which the inviscid numerical fluxes are calculated by a
class of high-resolution shock-capturing TVD (Total Variation Diminishing) and
ENO (Essential Non-Oscillatory) schemes, and the viscous terms are evaluated
with usual second-order accurate central difference. Several time-marching
schemes including the second-order accurate Runge-Kutta explicit scheme and
LUSSOR (Lower-Upper Symmetric Successive Over-Relaxation) implicit scheme that
is distinguished for its computational efficiency as it needs only scalar
diagonal inversions and is completely vectorizable on oblique planes of sweep
are incorporated into the model. Two zero-equation eddy viscosity models, an
empirical formula introduced by Fischer and a depth-averaged subgrid-scale model
of Smagorinsky, and one depth-averaged standard
model are available for resolving
the turbulent flow. General and flexible boundary-condition setups and the
treatment of drying/wetting and hydraulic structures are implemented. The model
is capable of simulating smooth or discontinuous steady or unsteady flow
conditions on complex shaped domains with complex topography. To assess the
performance of the proposed model, both steady and unsteady smooth and
discontinuous flows are simulated to verify its accuracy and robustness.
Keywords:
two-dimensional shallow water equations, discontinuous flow, multi-block
finite-volume method, TVD/ENO shock-capturing scheme, LUSSOR, Runge-Kutta
Most open-channel flows can be treated approximately as shallow water problems, because the effect of vertical motions is usually insignificant. In this paper a general and powerful numerical model, named SEC-HY21, to compute such flow by solving the depth-averaged, two-dimensional unsteady shallow water equations is presented. Many popular and well-designed schemes originally developed in the CFD field for solving the Euler and Navier-Stokes equations have been transplanted and integrated into the present model. A semi-discrete approach is adopted to uncouple the solution algorithm for time and space. The basic flow solver in space is a cell-centered finite-volume Godunov-type upwind scheme in which the inviscid numerical fluxes are obtained on the edges of each quadrangular cell with a flux-difference-splitting algorithm in conjunction with a second-order accurate variation-bounded, TVD or ENO, reconstruction ensuring the oscillation-free shock capturing capacity. As usually, the viscous fluxes are evaluated with second-order accurate central-difference. In the temporal discretization, both explicit and implicit schemes are implemented. In summary, the present model, SEC-HY21, has the following features:
(1) It is capable of handling problem domains of arbitrary complexity by use of multi-block boundary-fitted structured mesh.
(2) It is capable of simulating subcritical or supercritical steady or unsteady flow.
(3) It is capable of computing discontinuous flows such as those associated with strong bores or oblique hydraulic jumps.
(4) Second-order accurate solution in both space and time can be achieved.
(5) Turbulent flows can be modeled using
Boussinesq’s eddy-viscosity concept. Two zero-equation models and one
depth-averaged standard
model are included.
(6) It has the flexibility to specify general boundary conditions, either steady or unsteady.
In the remainder of the paper, the governing equations are given, and the basic solution algorithms are described, and then two numerical studies applying the proposed model are presented and compared with corresponding experimental data.
Neglecting the wind shear effect and Coriolis force, the depth-averaged two-dimensional unsteady shallow water equations can be given in conservation form as
(1)
where
(2)
(3)
here
water depth;
depth-averaged velocity components in x- and y-directions
respectively;
dimensionless surface-stress coefficient;
gravitational acceleration;
density of water;
dimensionless bed-friction coefficient in which
is the Chezy constant and
is the Manning coefficient;
friction velocity;
depth-averaged Reynolds stresses; and
bed slopes in x- and y-directions respectively. Use of
Boussinesq’s eddy-viscosity concept, the Reynolds-stress terms can be given as
(4)
where
represent
and
respectively, and
is the eddy viscosity. Two
zero-equation eddy viscosity models, an empirical formula introduced by Fischer
[1] and a depth-averaged subgrid-scale model of Smagorinsky [2], and one
depth-averaged standard
model [3] are available in this
version of SEC-HY21 model.
The whole
flow solver is built on a semi-discrete, cell-centered, multi-block,
finite-volume framework. Taking a semi-discrete approach is advantageous since
both spatial and temporal discretization schemes can be more easily and flexibly
designed. Briefly, the difference equation approximating Eq. (1), applied to a
quadrangular computational cell labeled (
), can be given in semi-discrete form as
(5)
where
and
are the average
and
of cell (
) stored at the cell center,
is the area of cell (
) ,
is the length of the
-th edge of the cell, and
are the numerical fluxes
approximating the physical inviscid flux
and viscous flux
through the edge respectively in
which
are the components of the
outward-pointing unit vector
normal to the edge. As
is calculated in a conventional
central-difference way, only the schemes for
are addressed below.
Use of flux-difference-splitting technique [4]
as an approximate Riemann solver, the inviscid numerical flux
, at the cell face between cell (
) and cell (
), can be given by
(6)
where
and
are the reconstructed left and
right states respectively of the cell face, and
is the normal flux Jacobian
evaluated by Roe’s average of
and
.
Similar to our previous works [5-6], the left and right states generated by a class of second-order accurate, total variation bounded reconstruction are
(7)
where
is the slope limiter defined as
(8)
here function
is the so-called limiter function
such as the popular “minmod” limiter. In equation (8), if
one has a second-order TVD scheme,
and if
one has a uniformly second-order
ENO scheme.
Both explicit and implicit schemes solving Eq. (5) are available in the current model. For explicit schemes we use Runge-Kutta methods. For implicit scheme, the LUSSOR scheme developed by Yoon and Jameson [7] for Navier-Stokes equations are applied. To save space, only the implicit scheme is outlined below.
Use of the LUSSOR scheme, Eq. (5) can be approximately discretized as
(9)
where
(10)
here the inviscid and
viscous Jacobians in the direction normal to the cell faces
and
have been denoted as
and
respectively.
are the split inviscid flux
Jacobians defined as
where
are the spectral radius of
respectively. And
and
are the spectral radius of
respectively. The above LUSSOR
scheme is noted for its computation efficiency as only scalar diagonal inversion
is involved and is completely vectorizable on the i+j = constant oblique planes of sweep.
Two major kinds of boundaries are needed in a numerical simulation, i.e., open and closed boundaries. At closed (or wall) boundaries, the no-slip, full-slip and partial-slip velocity condition can be prescribed, and zero gradient of water-surface elevation is usually used to obtain the water depth at the boundaries. For open boundaries, a locally one-dimensional, normal to the boundaries, characteristics theory is used in conjunction with appropriate hydrodynamic conditions, such as specifying a hydrograph of water-surface elevation or discharge, or even a depth-discharge rating curve for outlet boundaries. The velocity gradient normal to the open boundaries can be set to be zero (for outlet boundaries) or non-zero (for inlet boundaries). The discharge hydrograph used as an upstream boundary condition can be specified in two manners, i.e., specifying discharge per unit width which is uniform along the boundary, or specifying total discharge which is then distributed over the cross section according to local conveyance along the boundary.
Two numerical examples with corresponding experiments are presented. The first one is an unsteady dam-break flow to demonstrate the supreme shock-capturing capability. The second one is a steady flow near the spur-dike to examine the solution accuracy for vortical flows.
The first numerical example presented is a
two-dimensional flood wave resulting from the instantaneous break of a dam
located at the throat of a converging-diverging flume, as shown in Fig. 1,
experimentally investigated by Bellos et al. [8]. The wet-bed case with
reservoir depth
, tail water depth
, bed slope
, and Manning
as suggested by Bellos et al. is
chosen since the resulting flow is featured with several bore waves propagation
and reflection. A weir-type depth-discharge rating-curve is specified at the
downstream end of the flume to simulate the effect of the vertical wall with
height equal to the tail water depth and the free overflow if the water depth is
higher than its height. At the upstream end and banks of the flume, the
full-slip wall boundary condition is applied. The computational mesh consists of
42 cells and 10 cells longitudinally and transversely respectively. Computed
results and experimental data measured along the midstream line of the flume are
compared in Fig. 2. The measured depth hydrographs show several bore waves
propagating and reflecting on the channel boundaries with consequent
superimposition. The computation reproduces the experimental data excellently at
all the measured points where
and
are inside the reservoir;
is just upstream of the dam,
is inside the enlargement, and
is inside the prismatic sections.
Test run A1 of the experimental data obtained
by Nawachukwu [9] is used to further validate the numerical model. For this
case, tests were conducted in a straight flume 37-m long and 0.915-m wide. The
spur-dike consisted of an aluminum plate, 3-mm thick and 152-mm long, projecting
well above the water surface. The computational mesh used is
, covering the streamwise range from the inlet boundary, located 1.8 m upstream
of the spur, to the outlet boundary located 3.6 m downstream of the spur. At the
inlet boundary, the velocities are fixed to
. At the outlet boundary the depth is fixed to
. The turbulent eddy viscosity
and Manning
were used as suggested by Molls and
Chaudhry [10]. The computed recirculation zone ends at
, as seen from the streamline plot in Fig. 3, which coincides with the
experimental data. Fig. 4 shows a comparison of the experimental and computed
velocities. The resultant velocity
was non-dimensionalized with reference to the inflow velocity
. The values are plotted at
in the y-direction and from
to
. Here b is the length of spur-dike.
The agreement between the experimental and numerical results is satisfactory
except at
where similar discrepancy was also
reported by previous numerical investigators [10-11].
A versatile and powerful numerical model, named SEC-HY21, for simulating the open channel flows by solving the depth-averaged, two-dimensional unsteady shallow water equations is presented. The basic flow solver is based on a cell-center, second-order accurate, multiblock finite-volume method combining the high-resolution TVD and ENO discretization and an LUSSOR implicit or a Runge-Kutta explicit time-marching algorithm. The model is capable of handling the full range of flow regimes including gradually varied and discontinuous flows, steady and unsteady flows, subcritical and supercritical flows on irregular domain and complex topology. Verifications of the current model are made by comparison of the computed results with the available experimental data for both the steady flow near spur-dike and the unsteady dam-break flow, and both results show satisfactory or excellent agreement.
References
[1] Fisher, H. B., “Longitudinal dispersion and turbulent mixing in open channel flow”, Annu. Rev. of Fluid Mech., 5, 59-78, 1973.
[2] Smagorinsky, J., “General circulation experiments with the primitive equations. I. the basic experiment”, Mon. Weather Rev., 91, 99-164, 1963.
[3]
Younus, M. and Chaudhry, M. H., “A depth-averaged
turbulence model for the
computation of free-surface flow”, J. Hydr. Res., 32(3), 415-444, 1994.
[4] Roe, P. L., “Approximate Riemann solvers, parameter vectors and difference schemes,” J. Comput. Phys., 43, 357-372, 1981.
[5] Hsu, C. A., “High resolution non-oscillatory schemes for hyperbolic conservation laws with applications to aerodynamics”, Ph.D. Dissertation, I.A.M., N.T.U., 1993.
[6] Yang, J. Y. and Hsu, C. A., “Computation of free surface flows. part II. two-dimensional unsteady bore diffraction”, J. Hydr. Res., 31(3), 403-414, 1993.
[7] Yoon, S. and Jameson, A., “An LUSSOR scheme for the Euler and Navier-Stokes equations”, AIAA paper 87-0600, 1987.
[8] Bellos, C. V., Soulis. J. V., and Sakkas, J. G., "Experimental investigation of two-dimensional dam-break induced flow", J. Hydr. Res., 40(1), 47-63, 1992.
[9] Nawachukwu, B. A., “Flow and erosion near groyne-like structures”, Ph.D. thesis, Univ. of Alberta, Edmonton, Canada, 1979.
[10] Molls, T. and Chaudhry, M. H., “Depth-averaged open-channel flow model”, J. Hydr. Engrg., ASCE, 121(6), 453-465, 1995.
[11] Tingsanchali, T. and Maheswaran, S., “2-D depth-average flow computation near groyne”, J. Hydr. Engrg., ASCE, 116(1), 71-86, 1990.
Fig. 1 Experimental configuration of Bellos et al.
Fig. 2 Experimental and numerical hydrographs for Bellos et al. Dam-break flow
Fig. 3 Streamline plot for flow near spur-dike
Fig. 4 Experimental and computed Resultant velocities