Abstract

A well-balanced scheme with total variation diminishing Runge-Kutta discontinuous Galerkin (TVD-RK DG) method for solving shallow water equations is presented. Generally, the flux function at cell interface in the TVD-RK DG scheme is approximated by using the Harten-Lax-van Leer (HLL) method. Here, we apply the weighted average flux (WAF) which is higher order approximation instead of using the HLL in the TVD-RK DG method. The consistency property is shown. The modified well-balanced technique for flux gradient and source terms under the WAF approximations is developed. The accuracy of numerical solutions is demonstrated by simulating dam-break flows with the flat bottom. The steady solutions with shock can be captured correctly without spurious oscillations near the shock front. This presents the other flux approximations in the TVD-RK DG method for shallow water simulations.

1. Introduction

Hyperbolic balanced law for one-dimensional problem iswhere , , and represent solution vector, flux function, and source terms, respectively.

In this work, the hyperbolic equation is the shallow water equations which can be expressed bywhere is the water depth, is the discharge, is the flow velocity in the -direction, is the acceleration due to gravity, and is the bottom function. Equation (2) can be rewritten by setting

The total variation diminishing Runge-Kutta discontinuous Galerkin (TVD-RK DG) method can be applied to solve the shallow water equations; see [15]. This scheme has several advantages. For instance, it can be used to handle complex geometries, and in the same time, adaptive strategies are easily applied. The accuracy of numerical solutions can be improved by increasing the polynomial degree of approximating polynomial and by the method for estimating the flux function at cell interface.

By the concept of discontinuous Galerkin method, numerical solutions need not to be continuous at cell interface. So, we require an efficient flux approximation. Generally, there are several types of approximations. The most famous one is the Harten-Lax-van Leer flux (HLL), while the higher order approximation is the weighted average flux (WAF); see [69]. This approach approximates the flux functions by averaging flux in each direction along the wave solutions at the half-time step. It originates from random flux scheme, which was shown to be second-order accuracy in space and time in statistical sense by Toro [6]. From previous work, the weighted average flux has been successfully applied to solve various types of problems, especially in the finite volume method [7, 10, 11]. It improves the accuracy of the finite volume method to be second-order without reconstruction process. But the weighted average flux is rarely applied in the discontinuous Galerkin method. So the main objective of this work is the application of the WAF method in the TVD-RK DG method. We have also shown that the present scheme has consistent property with the WAF approximation.

The steady solutions of (3) can be obtained by settingFlux gradients are nonzero, and it must be balanced with the bottom gradient. Usually, numerical schemes without balancing these two quantities produce oscillate steady solution. So a balancing numerical scheme is needed, and it is called the well-balanced scheme. This method is designed to preserve energy at steady state. Generally, a numerical scheme is said to be well-balanced if it satisfies the exact C-property which was first introduced by Bermudez and Vazquez [12]. The exact C-property means that the numerical solution must satisfy still water condition at steady state given byThus, to obtain a well-balanced scheme, we must design the method that its steady solutions satisfy (5).

Recently, the second-order Runge-Kutta discontinuous Gelerkin method with well-balanced scheme is proposed by Kesserwani and Liang [13]. They presented the wetting and drying algorithms for solving one-dimensional problem. Xing and Shu [14] proposed a well-balanced finite volume method based on the weighted essentially nonoscillatory (WENO) scheme for solving one- and two-dimensional problems. They applied the decomposing source term technique into the sum of several terms to ensure the numerical balance between flux difference and source terms. Audusse et al. [15] presented a well-balanced finite volume with hydrostatic reconstruction to solve one-dimensional problem. Later, Noelle et al. [16] extended the well-balanced finite volume schemes proposed by [15] to any that desired the orders of accuracy.

In order to solve the shallow water equations with source terms, we modified the well-balanced discontinuous Galerkin scheme proposed by [15] and related work by Xing and Shu [14]. The weighted average flux (WAF) method [6, 7] has been applied in the modified DG scheme. We will show that the modified TVD-RK DG scheme with the WAF has consistent property. Various numerical experiments for both steady and unsteady flows are demonstrated to show the capability and accuracy of our numerical scheme.

Numerical scheme and its consistency property with WAF for the TVD-RK DG method have been shown in Section 2. The modified well-balanced TVD-RK DG scheme with source terms is presented in Section 3. Numerical examples are demonstrated in Section 4. Conclusions are finally made in Section 5.

2. Numerical Scheme for Shallow Water Equations without Source Term

The TVD-RK DG with weighted average flux method for the case of without source terms is presented in this section.

2.1. TVD-RK Discontinuous Galerkin Method

Consider the one-dimensional hyperbolic conservation law,The computational domain is divided into cells. We denote the th cell by , for , with uniform cell size . The cell center is that , where and are the left and the right of cell boundaries, respectively.

Approximate solution is denoted by .

Multiplying (6) by a test function, , where is the polynomial space degree on the interval , and replacing by and then taking the integration by parts over , we obtain a weak form of numerical scheme aswhere the flux function at the cell interfaces is approximated by , which is the function of and at asHere and denote the approximate solutions at the left and the right of cell boundaries, respectively. If we apply the Legendre polynomials to be a local basis function, the approximate solution can be written bywhere is called the temporal coefficient. Now, the basis function is defined by the Legendre polynomial of degree over .

The test function is typically chosen to be the basis function; that is, . So after applying Legendre’s properties, (7) is simplified tofor , and .

The time derivative term in (10) can be approximated by applying the high-order TVD Runge Kutta (TVD-RK) method; see [1, 3, 17]. Note that, when the polynomial degree is applied, the TVD-RK method at least the order of must be applied to obtain the accuracy of order for smooth flows.

The TVD-RK DG method can be used to simulate shallow water flows with moving shocks. Unphysically oscillate behaviours are usually produced near the shock fronts. The slope limiter techniques can be applied to remove the oscillations. In this work, we apply the Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) limiter (see [1, 3, 6, 18]) in the TVD-RK DG method. This approach limits the present solution slope by comparing the slope with neighbor cells.

2.2. Weighted Average Flux (WAF)

The HLL numerical flux (Harten-Lax-van Leer) is usually used in the TVD-RK DG method. Another choice but higher order is the weighted average flux (WAF). It was first introduced by Toro; see [6, 7]. The WAF approximation is second-order accurate in both space and time in statistical sense [6]. This approximation has been extensively applied in the finite volume method, but it is rarely used in the TVD-RK DG method. So the main objective of this work is to try to modify the WAF in the TVD-RK DG method.

The weighted average flux, at the interface, , is defined by the integral average of a flux function at the half-time step,It can be written in the wave structure form aswhere is the number of wave solutions in the Riemann problem and is the flux of the Riemann problem. For one-dimensional shallow water flows, we have , where and . Flux component can be obtained from the HLL approach [6]. Weighted values, , are defined by , where is the Courant number of wave , , , and is the speed of wave .

To avoid spurious oscillations near a shock front, the WAF method will be modified by enforcing a total variation diminishing (TVD) scheme [6, 7, 9, 11, 19]. The TVD-WAF version is thatwhereHere is a WAF limiter function. There are various choices; see more details in [6, 7, 9, 11, 19]. In this work, we choose the basic one of minmod type.

2.3. Consistency of WAF with TVD-RK DG Method

In this section, we show the consistency of the weighted average flux in the TVD-RK discontinuous Galerkin method.

Lemma 1. The TVD Runge-Kutta discontinuous Galerkin method is consistent with the weighted average flux for smooth flows.

Proof. Considering the weak form of RKDG method for the conservation law,Solution is approximated by . The test function is estimated by . Flux function at cell interfaces is approximated in terms of numerical flux, , so we obtain a numerical scheme of the following form:Substituting by in (16) and subtracting with (15), we obtain a truncation error term asWhen the solution is smooth, it must be continuous at or .

The TVD weighted average flux is defined bywhere , , and the component which can be obtained from the HLL method [4, 6, 7]. Thus,Similarly, we have . Then, the weighted average flux is consistent.

For the TVD-RK DG method, the approximate solution in (16) is defined bywhere is the Legendre polynomial degree . Since the considering solution is smooth, we know from Theorem (3.1) in [1] thatThis norm is defined as a truncation error term due to approximating polynomial, denoted by .

After substituting the approximate solution into the numerical scheme (16), we obtain an ODE system,Next, we apply the TVD Runge-Kutta scheme for integrating in time. We approximate at each time step by . The order of accuracy for time integration depends on the TVD-RK order [1, 3, 17]. If the TVD-RK order is applied, a truncation error term is given byCombining all of the truncation error terms together, the total truncation error term is estimated byThe total truncation error term is approaching zero as and . Hence, the TVD-RK DG with the WAF method is consistent.

3. Well-Balanced TVD-RK DG with WAF Scheme

In this section, we develop a well-balanced scheme for the TVD-RK DG with the WAF method. The main purpose is to present a modified scheme for solving the shallow water equations with source term. Also, this scheme must preserve exactly stationary solution when bottom slope exists. Let us start by considering the standard TVD-RK DG method,with initial conditionwhere are the weighted average flux in (13). Source term is given by

We will derive a well-balanced scheme based on [14], but the weighted average flux is applied in the TVD-RK DG instead of using the Lax-Friedrichs (LF) flux. The main modification is the source term treatments in the WAF method. Assume that the source terms can be written by the summation of some functions aswhere and are some functions which will be determined later.

If the solution at steady state is stationary, then in (28) can be decomposed by and , such thatSince is zero, we consider only ,From (28), we have that To balance the flux gradient and the source term approximation at steady state, it is required thatorwhere is arbitrary constant.

The integration of source term in (25) can be approximated byFunctions and on the RHS of (34) can be approximated by and , where and are the projection of and into the space of . Then, we obtain and . Thus,Then,To satisfy the weighted average flux in (12), must be modified. Here, we assume thator if the TVD version in (13) is applied, this term can be written byNote that are weighted values in the WAF approximations and is the WAF limiter function. So can be defined similar to .

For one-dimensional shallow water equations, we have , where and . The approximation of in the intermediate region can be obtained by the HLL approach,Here and are signal velocities in the Riemann problem.

At steady state, the solution is assumed to be stationary. It implies that , where is constant and . From (35), we have thatThis shows the suitable choices of and in the TVD-RK DG with the WAF approximation.

Next, we will show that the TVD-RK DG with the WAF method preserves the well-balanced property.

Proposition 2. The TVD-RK DG with the WAF scheme has preserved exactly stationary solutions at steady state.

Proof. Since and are equal to the same constant at each Gauss-Lobatto point over cell , thusHence, we obtain a truncation error term asUsing in (12), in (37), and , we have thatAfter applying condition (35) and rearranging terms, yields where is arbitrary constant. Hence, the TVD-RK DG with the WAF scheme preserves exactly the stationary solution at steady state.

Remark 3. One can show that the developed well-balanced TVD-RK DG with WAF scheme in the case of existing source term is consistent, by showing that the approximation of the source term is also consistent.

4. Numerical Results

In this section, various test cases have been investigated to demonstrate the accuracy of the present scheme, not only steady, but also unsteady flows.

4.1. Dam Break Flow

It has been shown in the previous section that the weighted average flux is consistent with the TVD-RK DG method. We apply this modified scheme to solve the shallow water equations without source terms in this subsection. The accuracy of numerical solutions is shown and compared with the standard TVD-RK DG when using the HLL method.

The computational domain is that . The initial water depth is given byThe initial velocity is assumed to be zero. The boundary conditions are transmissive boundaries. We perform 50, 100, and 200 cells in the numerical experiments. Polynomial degree zero, one, and two are applied as a local basis in the TVD-RK DG method. Simulation time is , with . The root mean squared errors (RMS) are shown in Table 1.

When and are fixed, the accuracy of numerical solution obtained from the WAF is higher than those obtained from the HLL method. The RMS errors decrease as increases.

If we fix and vary , the RMS error decreases as the polynomial degree increases.

The water depth profiles using the HLL and the WAF methods at when and are shown in Figures 1 and 2, respectively. The front of moving shock can be captured correctly by the HLL and the WAF methods. But the scheme using the WAF can capture shock and rarefaction wave more precisely than those using the HLL method. This investigation shows the accuracy of the TVD-RK DG with the WAF method. It has been generalized from the finite volume method.

4.2. Flow over Irregular Bed

The uniform channel is length of 1500 m. The bottom elevation is irregular that is shown in Figure 3. This problem is proposed by [5] for testing the accuracy of numerical scheme at stationary state. The boundary conditions are transmissive. The initial water depth is that , with zero initial velocity. We set and run simulation until .

It can be seen from Figure 3 that the well-balanced scheme (dot) gives exactly the stationary solution while the non-well-balanced scheme (dash line) gives solution error especially in the high gradient area of bottom elevation.

4.3. Steady Flow over a Bump

We consider the shallow water flows over a bump in a rectangular channel with length of 25 m. The bump elevation is given by

At steady state, classical flows can be characterised by subcritical flow, transcritical flow with shock, and transcritical flow without shock. The TVD-RK DG with the WAF method is performed to solve these problems. The accuracy of numerical solutions can be checked by comparing with the existing analytical solutions; see [20].

Subcritical Flow over a Bump. The upstream boundary is imposed by  m2/s, while the downstream boundary is set by m. The initial water depth is with zero initial velocity. Time step is that . The RMS errors obtained from the HLL and the WAF methods are shown in Table 2. This shows that the RMS errors obtained by the TVD-RK DG WAF are less than those obtained by the HLL for all .

The water depth and the bump profiles are shown in Figure 4. The numerical solution is very close to the analytical solution. We also found in this test case that the well-balanced scheme converges to the steady solution faster than the scheme that is not well-balanced.

Transcritical Flow with Shock over a Bump. The upstream boundary is given by /s, while the downstream boundary is set by m. The initial condition is that m. The comparison of water surfaces is shown in Figure 5. The numerical result is in good agreement with the analytical result. This shows the accuracy of the well-balanced scheme that can capture the shock front without any oscillations. It is also found that the well-balanced scheme converges to the steady solution faster than the non-well-balanced scheme.

Transcritical Flow without Shock over a Bump. The upstream boundary is prescribed by /s, while the downstream boundary is not specified. The initial conditions is that m with zero initial velocity.

The water depth profiles are shown in Figure 6. The numerical result agrees well with the analytical solution. These results show the accuracy of the well-balanced scheme for solving transcritical flow problem.

4.4. Small Perturbation of Steady State Water

This problem is first proposed by [18, 21, 22]. It can be used to study the capability of numerical scheme for solving small perturbation in shallow water flow. The bottom topography is defined byThe initial conditions are specified bywhere is a nonzero perturbation constant. The boundary conditions are transmissive boundaries. In this work, we consider the cases of and 0.01. The disturbance of initial water depth from small should split the initial wave into two waves. They propagate to the left and the right with characteristic speed at the early stage. A standard scheme which is not well-balanced usually faces with some difficulties to capture correctly the wave speed.

In our simulation, we use 400 uniform grid cells and polynomial degree one in the TVD-RK DG with the WAF method. The simulation time is that .

The comparison of water depths between our results and the Leveque’s solutions is shown in Figures 7 and 8 for and 0.2, respectively. They are in good agreement for both amplitude and wave speed. These test cases show the ability of our numerical scheme for solving the quasi-stationary flow with initial disturbance.

4.5. Flow over Nonhorizontal Bed

This test case is presented by [8]. The main propose is to study the ability of numerical scheme for solving unsteady flow over topography. The uniform channel is length of 30 m. The bed elevation is defined byThe initial conditions are given by Here, we set and use 200 uniform grid cells with polynomial degree one to simulate the unsteady flow. The boundary conditions are transmissive boundaries.

Our numerical results when comparing with Toro’s solutions [8] at  s and 4 s are shown in Figure 9. Wave speed and shock profiles are very close. These results show the accuracy of the present scheme for solving unsteady flow with source term.

5. Conclusions

In this work, we present the TVD-RK discontinuous Galerkin method (TVD-RK DG) for solving nonlinear shallow water equations. Most of the TVD-RK DG methods in the literatures usually approximate intercell flux by applying the HLL method. But here we apply another approach called the weighted average flux (WAF) in the TVD-RK DG. We have also shown the consistent property of the TVD-RK DG with the WAF approximation. Then the well-balanced TVD-RK DG scheme with the WAF approximation is developed. The present method can be used to simulate not only steady flows, but also unsteady flows. The accuracy of modified numerical scheme is demonstrated by various test cases, flow over irregular bed, steady flow over a bump, quasi-stationary, and flow over nonhorizontal bed. The well-balanced TVD-RK DG with the WAF method can be used to solve all the kinds of these problems. Moreover, if we restrict at steady state, the scheme using the WAF method converges to the steady solution faster than the scheme using the HLL method. In addition, the well-balanced scheme converges to the steady solution faster than the scheme that is not well-balanced. Due to its advantages of numerical accuracy, simplicity, and well-balanced property, the present scheme can be modified and extended to simulate two-dimensional problems. However, depending on the types of elements, for example, triangles or rectangles, it is not trivial to extend for two-dimensional problems due to the polynomial basis functions and the WAF fluxes at element interfaces and is not considered in this paper.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research is partially supported by the Science Achivement Scholarship Thailand (SAST) to the first author and financially supported by the Thailand Research Fund (TRF) under the Grant no. RSA5680038.