(1)
Chen Zuhua1
Wang Guangqian2
Wang Zhishi3
1 PhD. Candidate in Tsinghua University, Key Laboratory for Water and Sediment sciences of Ministry of Education, China
2 Prof., in Tsinghua University, Key Laboratory for Water and Sediment sciences of Ministry of Education, China
3
Prof., in University of Macau
Abstract:
A numerical model for the solution of the two-dimensional unsteady
depth-averaged flow is described. The model, which is based on the Roe’s
approximate Riemann solver, is used to solve the 2D depth-averaged equations and
solute advection-diffusion equation on an irregular structured grid. A modified
MUSCL technique is combined to preclude the restrictions on the Roe’s scheme
duo to bed slope as well as to obtain a second-order accuracy. The wetting and
drying process is considered by using an artificial critical depth. The model is
verified by comparing the model output with experimental data and the field data
of Macau.
Keywords: finite-volume method, Roe’s scheme, drying-wetting process, two-dimensional unsteady flows
The shallow-water flow model is popular in modeling free-surface flow in many engineering applications. They may be tidal movement in estuaries and coastal areas, dams breaking, as well as wind-driven circulations in lakes. The shallow water wave equations admit discontinuous as well as smooth or classical solutions due to the hyperbolic character of them. Even for the case in which the initial data is smooth, the nonlinear character combined with the hyperbolic type of the equations can lead to discontinuous solutions in finite time, such as hydraulic jumps, tidal surges, that make it difficult to obtain the numerical solutions to these equations. It is well known that common numerical schemes to deal with the nonlinear advection items suffer from serious numerical diffusion or unphysical spurious oscillations when they are used to solve convection-dominated flows, or problems with discontinuity solutions.
Recently, some Riemann solver-based techniques, namely, the Godunov-type schemes, which have been successfully applied to a wide class of problems involving the Euler equation of gas dynamics, have also been used for the numerical simulation of the shallow water flow, for example[1],[3],[10]. When applying those Godunov-type schemes to the simulations of the shallow water flow, one of important aspects differed from the applications in the gas dynamics is that the bed slopes and the frictions on the bottom, called source terms, are important influencing factors in shallow water flow. A careless treatment of the source terms can seriously damage the conservativity of the scheme.
The numerical method we are presenting in the paper is a high-resolution, oscillation-free, based on the Roe’s approximate Riemann solver, finite-volume method which solved the bidimensional depth-averaged equations on irregular structured grids. The modified Monotone Upstream Scheme for Conservation Laws (MUSCL) is used to deal with the bed level changes as well as to provide a second-order accurate solution. The wetting-drying process, which is an important phenomenon in many practical cases, is incorporated by using a critical water depth (Zou and Chen, 1989). The mass advection-diffusion equation is also included to predict solute concentration distributions in coastal and estuarine waters to meet the needs of environmental impact assessments in those regions.
The equations that form the basis of the model are depth-averaged, fluid, momentum and solute conservation equations. The vectorized form of the governing equations is given by
(1)
(2)
(3)
where h = water depth measured vertically, u and v = depth-averaged velocity components in x- and y-direction, respectively; c = solute concentration; r = water density; f = Coriolis coefficient; tax and tay = wind shear stress components in x- and y-direction, respectively; S0x and S0y = bed slope components in x- and y-direction, respectively; tbx and tby = bottom shear stress components in x- and y-direction, respectively; Sc = source or sink term of solute. The bottom stresses are modeled according to the following formulas:
(4)
where Cf is
the drag coefficient, Cf =
, in which n is the Manning’s coefficient.
The turbulent stresses and the solute turbulent diffusion flux terms are closured as below:
(5)
(6)
where nt = turbulent kinematic viscosity, and sc = the turbulent Schmidt number for solute. In this paper, sc is given a constant 1.0 and the following simple closure for nt (van Rijn 1987) is adopted:
(7)
in which k is von Karman constant; k = 0.4
The finite-volume method is adopted to solve the governing equations. A structured grids are used and every mesh can be irregular quadrilateral so that the realistic geometries can be represented as accurate as possible. If U is treated as constant in a cell and is equal to the value at the center of form figure of the cell, the finite-volume discrete form of (1) is expressed as:
(8)
where Dt denotes the time increment and n+1 denotes the time level of the solution, WC = the area of cell. The subscripts E, S, W, N indicate the midpoints of the borders of the cell while l is the cell face length, see Fig.1.
Fig. 1 Sketch of Typical Cell
The numerical flux
across the cell interface,
, is normal to the cell border, positive in the direction of outward normal
line, and composed of two parts, one is the nonlinear advection flux F^
，the
other is the diffusion flux
, that is,
(9)
As stated above, the evaluation of numerical advection flux is of most importance for the numerical solution of shallow water equations. In the present model, the Roe’s first-order upstream scheme based on the Roe’s approximate Riemann solver (Roe, 1981) is adopted. Firstly a local coordinates transformation is introduced to simplify the calculation of the numerical advection flux. Let us write
(10)
where (nx, ny ) are components of the outward unit normal to the cell border. The transformation V=TU will transform the vector of conservative variables U in global coordinates to one in the local normal coordinates of the border. Denote V = [h hU^ hU|| hc]T in which U^ = velocity normal to the cell border, and U|| = velocity parallel to the cell border. Then the advection flux F^ can be obtained by the inverse transformation
(11)
where n denotes the vector of the outward unit normal, n = [nx ny]. F(VL,VR,n) is the numerical flux in the local normal coordinates. According to the Roe’s scheme, it can be written as
(12)
in which the subscripts L and R denote the
left (the inverse direction of the outward normal) and right states with respect
to the interface. The matrix
is diagonal and contains approximations of the absolute value of the eigenvalues
of the Jacobian of F at the interface
of cell. It is of the following form
(13)
where
.
The matrix R in which the columns are the corresponding right eigenvectors is
written as
(14)
The characteristic intensities of the
differences of the conservative variables across the cell interface, defined as DV=R-1(VR-VL),
are as below:
(15)
The variables denoted with a ‘~’ are known
as Roe averages and are determined by
(16)
The averages
and
are computed in a manner analogous
to
.
One problem with the Roe scheme is that it cannot distinguish a rarefaction wave with a zero eigenvalue but treats it as a discontinuity. In the present model, a entropy modification proposed by Monthe L.A. et al (1999) is adopted so that the positivity of the scheme is preserved.
The MUSCL approach (Van Leer ,1979) is used to reconstruct VL and VR in order to obtain second-order accurate. Take an example of the east cell border: denotes the vector of the primitive variables W=[h u v c]T, then
(17a)
(17b)
where
and
are the components of limited gradient of W
in the cell (i, j) in x- and y-direction, respectively; x and y is coordinates
of a point; superscript E mean point E (Fig.1) and superscript G indicates the
center of form figure of cell. Similar expressions may be written for the other
cell interface. The limited gradient upon the cell is evaluated by the following
procedure. First,
and
are determined as the minimum
points of the following quadratic function:
(18)
where K(i, j) = {(i+1,j), (i,j+1), (i-1,j), (i,j-1)}, is the indices set of neighborhood cells of the cell (i,j). The sign f can be replaced by each component of W. The equation (18) can be solved by the least square method. Then the limited gradient is evaluated by limitation techniques to guarantee the monotonicity of the scheme. A general two-dimensional limiter is obtained by
(19)
where
.
is evaluated in the same way.
For shallow water flow , the first component in W, that is, water depth, is replaced by water surface level in the MUSCL procedure. After the left and the right surface levels with respect to the cell interface are evaluated, the corresponding water depths are calculated by subtract the bed level at the midpoint of cell border from the surface levels. This method is also used by Hu & Tan(1994).
The diffusion flux
can be computed by central
differencing scheme as below:
(20)
where avg(*) means average value at the cell interface and is evaluated, i.e., for txx,
(21)
in which the quantities
with an overbar means the arithmetic average of the values at the left and the
right cells of the interface. Similar expressions may be written for the other
components of
and
.
The pointwise approach to deal with the source terms in the shallow water equations is adopted for the advantage of simpleness. That is
(22)
Because the drying-wetting process is an important phenomenon for flows in estuarine and coastal regions, it is considered in the model. One simple way to handle this problem is to put the null value for depth and velocity in cells when the water depth goes under an artificially critical value (Zou and Chen, 1989). This method is adopted in this paper. However, the critical depth can be chosen very small so that the conservative error is reduced as small as possible. When the water depth is very small, the numerical stabilities are guaranteed by evaluated the friction velocity as
(23)
One physical explanation for (23) is that the velocity can diminish at most to zero on the action of the bottom friction and no invert velocity can be produced.
The implementation of boundary conditions is achieved by the use of a virtual cell at the boundary, but one outside the domain of interest. The flux at the boundary can be evaluated in the manner previously described by specifying the appropriate quantities in the virtual cell.
There are two types of boundaries in shallow water flow. These are solid-wall and open boundary. At a solid wall, the primitive variables in the virtual cell are specified in this manner: h and c are equal to those in the inner boundary cell, the velocity vector in the virtual cell and that in the inner boundary cell are symmetrical about the boundary wall so that the perpendicular velocity at the wall is zero while the parallel velocity remains unchanged, i.e., non-slip condition. These conditions must be kept when the MUSCL approach is applied.
For open boundary, usually one or more quantities must be given while the unknown quantities are determined by additional numerical boundary conditions on the assumption of local unidimensional flow. For example, for inflow boundary, the discharge of per unit width and the solute concentration are specified and the remaining quantities are extrapolated from the flow domain. For subcritical outflow boundary, the water level is specified while the velocities are zero-order extrapolated from the flow domain.
To demonstrate the performance of the method presented in this paper, several typical tests
were computed. In the first example we consider the test case developed as a laboratory experiment by J.M. Hiver from the Free University of Brussels, Belgium (Garcia-Navarro & Vazquez-Cendon, 2000). Additionally, the solute transport is included in our computation. The concentration of solute is assumed to be 10 g/l in the upstream of the gate. The diffusion terms and the source term of solute transport are omitted in this example. The gate is suddenly removed and the volume of water is released. Figs.2 display water and concentration profiles along the channel at time t=10s after the gate removal. The Manning roughness coefficient n = 0.0125 is chosen as in the original reference. Fig.3 shows a comparison of the results obtained by this paper against measured data at x = 19.5m.
The second example is the simulation of tidal currents in the surrounding waters of Macau. The computational domain is about 32 km long and 22.5 km wide and covers from Sanzao Island to the west part of Lingdingyang Bay as shown in Fig.4. From the calibration, the Manning roughness coefficient is chosen as 0.028.
Fig. 2 Dam-break wave evolution.water level. t =10 s
Fig. 3 Water level at x = 19.5 m. numerical results versus experimental data
Fig. 4 Bathymetric contours of computational domain
The south boundaries are specified by the observed simultaneous tidal profiles [4] and at the West River estuary the discharge is given as 1149 m3/s [6]. However, the data for the other boundaries are unavailable, then, numerical boundary conditions are supplemented. The depth-averaged velocity profile at P3 was also available for model calibration and shown in Fig.5. Good agreement with the measurements can be seen.
(a) Velocity Profile (b) Flow Direction
Fig. 5 Comparison between computed and observed data profiles of velocity at P3
A positivity preserving finite-volume numerical model for the numerical solution of the 2D depth-averaged flow and solute transport has been described. The Roe’s upstream scheme with the entropy modification is employed to treat the nonlinear advection parts of the governing equations in order to deal with discontinuities that may be produced in shallow water flows. The modified MUSCL technique in which water level is interpolated instead of water depth is adopted to make it possible to apply the Roe’s scheme to the situations that bed level changes greatly. In order to evaluate the gradients of variables more accurately in irregular grids, the least square method is used. And a general two-dimensional limiter is introduced to accomplish the monotonicity of the scheme. The wetting-drying process is straightforward incorporated in the model with recourse to a threshold limitation and the use of (23) combined the positivity of the scheme guaranteed the numerical stability so that the threshold can be much more smaller than that in previous use.
Verification of the model is made using two examples. The
model is shown to be conservative, accurate and robust. It is capable of
handling discontinuities, adapting the topography with steep bed slopes and
simulating the wetting-drying process. Though attractive results were shown in
the example of Macau, it is not proven in theory to interpolate water level
instead of depth in the MUSCL method.
Acknowledgements
This project was funded by the National Natural Science Foundation of P.R. China, No.49976025 and 59890200.
References
[1] Bradford S. F. and Katopodes N. D. (1999) Hydrodynamics of turbid underflows. I: Formulation and numerical analysis. Journal of Hydr.Eng. Vol.125, No.10, pp1006—1015.
[2] Garcia-Navarro P. & Vazquez-Cendon M.E. (2000) On numerical treatment of the source terms in the shallow water equations, Computers & Fluids,Vol.29, pp951-979.
[3] Hu Siyi, Tan Weiyan (1995) Numerical modeling of two-dimensional shallow water flows on unstructured grids, Advances in water science, Vol.6, n1, in Chinese, pp1—9.
[4] Instituto Hidrográfico (1992) The numerical hydromorphological model of the surrounding waters of Macao, Final report.
[5] Monthe L.A., Benkhaldoun F., Elmahi I. (1999) Positivity preserving finite volume Roe schemes for transport-diffusion equations. Comput.Methods Appl.Mech.Engrg, Vol.178, pp215-232.
[6] PRWRC (1990) Report on the physical tidal modeling of the surrounding waters of Macao.
[7] Roe, P..L. (1981) Approximate Riemann solvers, parameter vectors, and difference schemes, J. Computat. Phys., vol.43, pp357—372.
[8] Van Leer, B. (1979) Towards the ultimate consevative difference scheme. V. A second order sequel to Godunov’s method. J. Computat.Phys. Vol.32, pp101—136.
[9] Van Rijn, L.C.(1987) Mathematical modeling of morphological processes in the case of suspended sediment transport. Communication No.382, Delft Hydr., Delft, the Netherlands.
[10] C. Zoppou and S. Roberts (2000) Numerical solution of the two-dimensional unsteady dam break, Applied Mathematical Modeling, Vol.24, pp457—475.
[11] Zou GY, and Chen YS (1989) A numerical study of the motion of a wave running up a beach. Acta Mechanica sinica, Vol.21, n1, pp1—8, in Chinese.