Abstract

A numerical scheme based on discontinuous Galerkin method is proposed for the two-dimensional shallow water flows. The scheme is applied to model flows with shock waves. The form of shallow water equations that can eliminate numerical imbalance between flux term and source term and simplify computation is adopted here. The HLL approximate Riemann solver is employed to calculate the mass and momentum flux. A slope limiting procedure that is suitable for incompressible two-dimensional flows is presented. A simple method is adapted for flow over initially dry bed. A new formulation is introduced for modeling the net pressure force and gravity terms in discontinuous Galerkin method. To validate the scheme, numerical tests are performed to model steady and unsteady shock waves. Applications include circular dam break with shock, shock waves in channel contraction, and dam break in channel with bend. Numerical results show that the scheme is accurate and efficient to model two-dimensional shallow water flows with shock waves.

1. Introduction

Free surface flows that take place in rivers, oceans, and estuaries are of great importance to human activities. Numerical models for the open channel flows have been of great interests to hydraulic researchers and engineers. In practice, numerical models for flows with shock waves present difficulties in capturing shock wave and preserving conservative properties of the flow equations. Many shock-capturing methods have been developed in previous studies [1, 2]. The first-order upwind approximate Riemann solvers are commonly used to deal with discontinuous flows. However, the first-order schemes are known to generate oscillation around discontinuities. To deal with the unphysical oscillations, these methods are extended to high-resolution schemes based on the concept of Total Variation Diminishing (TVD) [3, 4].

In recent years, discontinuous Galerkin (DG) finite element method has gained popularity in modeling shallow water flows [5, 6]. The DG method is a combination of the finite volume method and the finite element method. In DG method, higher-order spatial accuracy can be achieved by employing higher-order interpolation functions. Since DG method allows discontinuous solutions across element boundaries, it provides better solution strategy for problems involving shocks and discontinuities. In addition, the explicit solution for one element at a time in the DG method is efficient in practical problems as calculation for a large global matrix is avoided, especially in case of nonlinear problems that require iterative solution [7].

Various upwind schemes can be applied to solve for the numerical flux across the boundaries of an element with discontinuities. These include Roe’s flux function [8], HLL (Harten-Lax-van Leer) flux function introduced by Harten et al. [9], HLLE (Harten-Lax-van Leer-Einfeldt) flux proposed by Einfeldt [10], and HLLC (Harten-Lax-van Leer-Contact) flux presented by Toro et al. [11]. In this study, the HLL flux function is used for two-dimensional shallow water flows.

In this paper, a numerical model based on DG method is presented for two-dimensional (2D) shallow water flows. In contrast to the previous studies based on the DG scheme, the net pressure force and the gravity terms are combined in this paper. A discretization scheme is adapted for the combined term that eliminates numerical imbalance between the pressure force term and gravity term and simplifies computation. The flux terms are approximated using HLL flux function. A modified HLL flux function is used to handle flow over dry bed. It is well known that when higher-order numerical schemes are used, nonphysical oscillations are generated around discontinuities and steep gradients. TVD slope limiters are widely used to minimize these oscillations and stabilize numerical schemes. To achieve oscillation-free solutions, a four-step slope-limiting procedure in DG method is adapted for the two-dimensional incompressible flows. The governing equations are first introduced, and then the DG formulation is briefly described. The implementation of HLL flux function and treatment of the source term are discussed. The Total Variation Diminishing (TVD) method for time integration scheme and slope limiter is presented. Finally, several numerical examples are presented to test the present numerical scheme for shallow water flows with shock waves and flow over dry bed.

2. Governing Equations

The vector form of the depth-averaged, two-dimensional, shallow water flow equations can be written as where the vectors of conserved variables, fluxes in the and directions, and the sources term can be written, respectively, as follows: where = water depth, = water surface level, , and are the unit width flow rates in the and directions, respectively, and are the velocity components in the and directions, respectively, = gravitational acceleration and = Manning’s roughness coefficient. Since the net hydrostatic pressure is included in the source term, the numerical treatment of this new source that accounts for an accurate estimate of numerical flux will be given later.

The two-dimensional shallow water equations are derived by integrating the Navier-Stokes equations along the depth of the fluid body. Several assumptions are made such as hydrostatic pressure distribution and uniform velocity profile in the vertical direction. The advantage is that free surface location is determined as part of the solution. The two-dimensional shallow water flow equations can be applied in situations where vertical acceleration may be neglected, and the horizontal extent is much greater than the depth of flow [12, 13].

3. Formulation of Discontinuous Galerkin Method

For as a function of and , (1) can be written in vector form as For the main element (0) and the three surrounding elements (1, 2, and 3), a typical mesh configuration of triangular elements in DG method is shown in Figure 1. A node is composed of multiple vertices. The number of vertices at a node depends on the number of elements joining the node.

The Discontinuous Galerkin formulation is written for each element. The variation of variables within an element is represented by the values of the variables at the vertices and shape functions , where is a diagonal matrix of basis function. Equation (3) is multiplied by the weight function, which is the same as basis function for the Galerkin method, and the resulting equation is integrated over an element. The integration by parts of the flux term results in the following equation where   is the numerical flux across the boundaries of an element and is the outward unit normal vector at the element boundary. Since discontinuous elements are connected through the numerical flux across a common boundary, the accuracy with which the numerical flux is calculated is crucial to the DG method.

4. Evaluation of Numerical Flux

Since the discontinuous elements are allowed in discontinuous Galerkin method, a generalized local Riemann problem can be solved for the numerical flux. The numerical flux in (4) can be evaluated using upwind numerical flux functions [14, 15]. In this study, the HLL flux function is adopted.

For and , as the components of the unit normal vector in the and directions, respectively, the rotational invariance of the flux yields where is the rotation matrix. Defining , the numerical flux can be obtained through the evaluation of numerical flux . The follows the same functional relationship between and as, is that, between and .

The numerical flux, , computed from HLL flux function, is given by where subscripts and stand for the left- and right-hand side of the element boundary. To determine the wave speed, the average velocities along the left and right boundaries under consideration are determined. The normal component, , of these velocities is determined. The wave speeds are given by [16]: For flow over an initially dry bed, the wave speed for right-hand dry bed and left-hand dry bed boundaries are given, respectively, as:

5. Source Term Treatment

The water level slopes in the source term can be determined with Green’s theorem as follows [17]: where is the water level at the boundary of element 0 and element and is the area of the element 0. The value of is determined by first calculating the average water surface elevation at the center of the element under consideration and its surrounding elements and then interpolating the water surface elevation at the boundary using distance weighting. Numerical results show that this treatment can account for the net pressure force terms efficiently and accurately.

6. Time Integration

The Total Variation Diminishing (TVD) scheme is chosen to eliminate numerically generated oscillations near shocks and steep gradients. The TVD Runge-Kutta time integration and the slope limiter are used here to achieve the TVD property. Former studies have shown that, to conserve the TVD property, the Runge-Kutta time integration scheme should be one order higher than the shape functions [1820]. For the linear shape functions adopted in this work, the second-order two-stage TVD Runge-Kutta scheme is employed. Thus, for the linear elements used in this study, the scheme is second order accurate in space and time [21]. Equation (3) can be written in the following form To advance the solution from time step to , the second-order TVD Runge-Kutta scheme, as given by Gottlieb and Shu [22], can be written as:

7. Slope-Limiting Procedure

A slope-limiting procedure developed by Tu and Aliabadi [23] for compressible flows is modified for application to incompressible flows. The limiting procedure includes the following four steps. Firstly, the average solution of the conserved variables at the element centroid is computed by determining the arithmetic mean of the solutions at each vertex of an element as follows: Secondly, the unlimited gradient in each element is computed using the interpolation functions and the vertex solutions as given below: Thirdly, limited gradient in the element under consideration is calculated by taking the weighted average of the unlimited gradient of the surrounding elements as follows: where the weighted factors are given by: and is a small number introduced to prevent indeterminacy. The parameters , , and are the square of the norm of the unlimited element gradients. Lastly, the limited conservative variables at vertices of each element are calculated from the limited gradient and average values of the variables. The requirements for the reconstructed solution are to satisfy each component of limited gradient and preserve the average at the element centroid.

8. Numerical Results

Tests are performed in this section to examine the accuracy of the proposed numerical scheme to model shallow water flows with shock waves. Numerical tests include circular dam break, shock wave in circular dam break, steady shock wave in channel contraction, and dam break in channel with 45° bend. Numerical results are compared with an exact solution or measured experimental data, if available.

8.1. Circular Dam Break

To test the symmetric shock-capturing capability of the scheme, the idealized circular dam break problem is used [2426]. The problem domain with horizontal bed is shown in Figure 2. The radius of the dam is 11 m. Initially, the water depth inside the dam is set to 10 m, and water depth outside is 1 m. The circular dam is removed instantaneously, and the flow in the domain is computed. Numerical results at 0.8 seconds after the removal of the dam are shown in Figures 3 and 4. The corresponding velocity field is shown in Figure 5. The symmetry of the forward moving wave is well preserved. The initial water volume is 5909 m3, and the water volume at 0.8 seconds is 5910 m3, showing the mass is well conserved.

8.2. Shock Wave in Circular Dam Break

The same domain as used in the previous test is adopted here with different initial conditions. The initial water depth is 1 m inside the dam and 10 m outside the dam. After removing the dam, the circular shock moves inwards, passes through the singularity, and then expands outwards. The shock at 2 seconds is shown below in Figures 6 and 7 as water surface in 3D view and water depth contour, respectively. The velocity field is shown in Figure 8. Water surface in Figure 6 shows the scheme is oscillation-free, and the diffusion effects are minimal. Water depth contours in Figure 7 show that the flow symmetry is well preserved with unstructured elements. The initial water volume is 21591 m3, and the water volume at 2 seconds is 21590 m3, showing the mass conservation is well preserved in this model.

8.3. Shock Wave in Channel Contraction

The steady supercritical shock wave due to channel contraction is simulated to test the numerical scheme. The plan view of shock wave in a symmetric channel contraction is illustrated in Figure 9. In the figure, is the length of the channel contraction, is the angle of wall deflection, and and are the shock front angles. The flow velocities in regions 1, 2, and 3 are , , and , respectively. Ippen and Dawson [27] showed that a proper width ratio could minimize the disturbance in region 3 and limit the standing shock waves within the contraction part. In the following, the channel contraction suggested by Lin et al. [28] is adopted, so that results can be compared with exact solutions and other numerical schemes. The geometry and computational mesh of the channel are shown in Figure 10. The channel width at the upstream end is 20 m and the width at the downstream end is 10.548 m, the angle of wall deflection () is 12 degrees, and the length of contraction () is 22.234 m. For the inflow boundary conditions, the Froude number is 2.7, the water depth is 1 m, the longitudinal velocity is 8.4566 m/s, and the lateral velocity is zero. The exact solution of shock wave angles is found to be as and , while the water depth in region 2 and region 3 are 1.868 m and 2.562 m, respectively.

Steady flow solutions for water depth are shown in Figures 11 and 12 with a coarse mesh (3632 elements, see Figure 10). Figures 13 and 14 show a comparison of the exact solution and the simulated water depth along the dash line and the solid line (see Figure 10). Results from the coarse mesh and a refined mesh (14528 elements) are compared to show that the numerical scheme can perform adequately even for a coarse mesh. Results using the refined mesh show better resolution at the shock and weaker oscillation after the shock. Moriasi et al. [29] suggested several methods to evaluate/quantify the accuracy of the simulated results with respect to measured data. The current simulated results are evaluated using the Nash-Sutcliffe efficiency (NSE), which indicates how well the plot of simulated versus observed data fits the line of unit slope, percent bias (PBIAS), and ratio of the root mean square error to the standard deviation of measured data (RSR). NSE ranges between –∞ and 1, with 1 being the optimal value. PBIAS and RSR are error index with zero being the optimal value. Table 1 reports the accuracy of the simulated results with respect to the analytical solution along dashed and solid lines. The statistics confirm that the model can accurately simulate the shock waves due to supercritical flow in a channel contraction.

8.4. Dam Break in a Channel with 45° Bend

Physical models were built in the Civil Engineering Department Laboratory, Université Catholique de Louvain (UCL, Belgium) to model dam break and strong transient flows in sharp bends. Experimental data were collected and used to validate numerical models developed by the CADAM group [30].

The plan view of the channel with horizontal bed and 45° bend is shown in Figure 15. The gauge points are also shown in the figure, and their positions are listed in Table 2. The dam is represented by a gate at the outlet of the reservoir. The gate is pulled up rapidly to simulate instantaneous failure of the dam. The initial water level in the upstream reservoir is 0.25 m above the horizontal channel bed, and the channel downstream is dry. The Manning’s roughness coefficients of 0.0095 s/m1/3 for bottom and 0.0195 s/m1/3 for wall, as suggested by Frazão et al. [30], are adopted. These values for Manning’s roughness are based on the steady uniform flow.

After the removal of the gate, water flows rapidly into the channel and reaches the bend. The water reflects against the wall, and a shock forms and moves upstream. The velocity field at 3 seconds is shown in Figure 16, and the water surface at 10 seconds is shown in Figure 17. The velocity field shows that the flow is two-dimensional at the inlet and in the bend region. The reflected shock wave can be clearly seen in the water surface profile. The simulated hydrographs at 9 gauging points are compared with measured data in Figures 18, 19, 20, 21, 22, 23, 24, 25, and 26. Accuracy evaluation of the numerical results is shown in Table 2. Numerical results are in good agreement with the measured data, except at G2 (low NSE, high PBIAS and RSR) which is located at the exit of the reservoir. At G2, the magnitude of the reflected wave and its arrival time are predicted accurately. The difference in the water level drop immediately after the gate opening may be related to difference in the manner in which the gate is actually opened and simulated. It should be mentioned that the simulated results are similar to or better than that reported by previous studies.

9. Conclusions

A new numerical model is developed for two-dimensional shallow water flows with shock waves. In this model, the two-dimensional shallow water equations are solved by discontinuous Galerkin (DG) finite element method. The form of shallow water equations that simplifies computation and eliminate numerical imbalance between flux term and source term is used here. In this formulation, the net pressure force term is combined with the source term due to gravity. A proper treatment of the new source term is given. The HLL approximate Riemann solver is employed to calculate the mass and momentum flux. A slope-limiting procedure that is suitable for incompressible two-dimensional flows in DG solver is presented. The performance of the numerical scheme is tested by simulating circular dam break, shock in channel contraction, and dam break in channel with 45° bend. These tests include steady and unsteady shock waves, flow over initially dry bed, and subcritical and supercritical flow. The circular dam break tests show the shock capturing and symmetry preservation capabilities of the proposed scheme. The comparisons of the numerical results with analytical solutions and measured data demonstrate that the scheme is capable of simulating two-dimensional shallow water flows with shock waves as well as flow over dry bed.