# Consistent Weighted Average Flux of Well-Balanced TVD-RK Discontinuous Galerkin Method for Shallow Water Flows

Authors: Thida Pongsanguansin; Montri Maleewong; Khamron Mekchay
Journal: Modelling and Simulation in Engineering
 (2015)

---

## 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.


---

## Body

## 1. Introduction

Hyperbolic balanced law for one-dimensional problem is(1)Ut+FxU=GU,where U, F, and G represent solution vector, flux function, and source terms, respectively.In this work, the hyperbolic equation is the shallow water equations which can be expressed by(2)ht+qx=0,qt+q2h+gh22x=-ghbx,where h is the water depth, q=uh is the discharge, u is the flow velocity in the x-direction, g is the acceleration due to gravity, and b is the bottom function. Equation (2) can be rewritten by setting(3)U=hq,FU=qq2h+gh22,GU=0-ghbx.The total variation diminishing Runge-Kutta discontinuous Galerkin (TVD-RK DG) method can be applied to solve the shallow water equations; see [1–5]. 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 [6–9]. 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 setting(4)hu=constant,12u2+gh+b=constant.Flux 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 by(5)u=0,h+b=constant.Thus, 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 Section2. 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,(6)Ut+FxU=0.The computational domain 0,L is divided into K cells. We denote the jth cell by Ij=xj-1/2,xj+1/2, for j=1,…,K, with uniform cell size Δj=xj+1/2-xj-1/2. The cell center is that xj=xj+1/2+xj-1/2/2, where xj+1/2 and xj-1/2 are the left and the right of cell boundaries, respectively.Approximate solution is denoted byUh=hhqh.Multiplying (6) by a test function, vhx∈PNIj, where PNIj is the polynomial space degree N on the interval Ij, and replacing U by Uh and then taking the integration by parts over Ij, we obtain a weak form of numerical scheme as(7)∫Ij∂tUhvhdx-∫IjFUh∂xvhdx+F^j+1/2vhxj+1/2--F^j-1/2vhxj-1/2+=0,where the flux function F at the cell interfaces is approximated by F^, which is the function of Uh+ and Uh- at xj±1/2 as(8)F^j±1/2=F^Uhj±1/2-,Uhj±1/2+.Here Uh|j±1/2- and Uh|j±1/2+ 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 Uh can be written by(9)Uhx,t=∑m=0NUjmtφmx,where Ujmt is called the temporal coefficient. Now, the basis function φmx is defined by the Legendre polynomial Pm(x) of degree m over [-1,1].The test functionvhx is typically chosen to be the basis function; that is, vhx=φlxl=0N. So after applying Legendre’s properties, (7) is simplified to(10)dUjltdt=2l+1Δj∫IjFUh∂xφlxdx+2l+1Δj-1lF^j-1/2-F^j+1/2,for j=1,…,K, and l=0,…,N.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 N is applied, the TVD-RK method at least the order of N+1 must be applied to obtain the accuracy of order O(ΔxN+1) 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,F^WAF at the interface, xj+1/2, is defined by the integral average of a flux function F(U) at the half-time step,(11)F^j+1/2WAF=1Δx∫-Δx/2Δx/2FUj+1/2x,Δt2dx.It can be written in the wave structure form as(12)F^j+1/2WAF=∑k=1Nc+1ωkFj+1/2k,where Nc is the number of wave solutions in the Riemann problem and Fj+1/2(k) is the kth flux of the Riemann problem. For one-dimensional shallow water flows, we have Nc=2, where Fj+1/2(1)=FUj+1/2- and Fj+1/2(3)=FUj+1/2+. Flux component Fj+1/2(2) can be obtained from the HLL approach [6]. Weighted values, ωk, are defined by ωk=1/2ck-ck-1, where ck=SkΔt/Δx is the Courant number of wave k, c0=-1, cNc+1=1, and Sk is the speed of wave k.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 that(13)F^j+1/2WAF-TVD=12FUj+1/2-+FUj+1/2+-12∑k=1Ncsgn⁡ckϕj+1/2kΔFj+1/2k,where(14)ΔFj+1/2k=Fj+1/2k+1-Fj+1/2k.Here ϕj+1/2(k) 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,(15)∫Ij∂∂tUvxdx-∫IjFUv′xdx-FUj+1/2vxj+1/2-+FUj-1/2vxj-1/2+=0.Solution U is approximated by Uh∗. The test function v is estimated by vh. Flux function at cell interfaces is approximated in terms of numerical flux, F^j±1/2, so we obtain a numerical scheme of the following form:(16)∫Ij∂∂tUh∗vhxdx-∫IjFUh∗vh′xdx-F^j+1/2Uh∗j+1/2-,Uh∗j+1/2+vhxj+1/2-+F^j-1/2Uh∗j-1/2-,Uh∗j-1/2+vhxj-1/2+=0.Substituting Uh∗ by U in (16) and subtracting with (15), we obtain a truncation error term T1 as(17)
T
1
=
∫
I
j
∂
∂
t
U
-
U
v
h
x
d
x
-
∫
I
j
F
U
-
F
U
v
h
′
x
d
x
-
F
U
j
+
1
/
2
-
F
^
j
+
1
/
2
U
j
+
1
/
2
-
,
U
j
+
1
/
2
+
v
h
x
j
+
1
/
2
-
+
F
U
j
-
1
/
2
-
F
^
j
-
1
/
2
U
j
-
1
/
2
-
,
U
j
-
1
/
2
+
·
v
h
x
j
-
1
/
2
+
.When the solution is smooth, it must be continuous at xj+1/2 or Uj+1/2+=Uj+1/2-=Uj+1/2.The TVD weighted average flux is defined by(18)F^j+1/2Uj+1/2-,Uj+1/2+=12FUj+1/2-+FUj+1/2+-12∑k=12sgn⁡ckϕkΔFk,where F(1)=FUj+1/2-=FUj+1/2, F(3)=FUj+1/2+=FUj+1/2, and the component F(2)=FUj+1/2 which can be obtained from the HLL method [4, 6, 7]. Thus,(19)F^j+1/2Uj+1/2-,Uj+1/2+=12FUj+1/2+FUj+1/2-12∑k=12sgn⁡ckϕkΔFk=FUj+1/2.Similarly, we have F^j-1/2Uj-1/2-,Uj-1/2+=F(Uj-1/2). Then, the weighted average flux is consistent.For the TVD-RK DG method, the approximate solutionUh∗ in (16) is defined by(20)Uh∗x,t=∑k=0NUktφkx,where φk(x) is the Legendre polynomial degree k. Since the considering solution is smooth, we know from Theorem (3.1) in [1] that(21)U-Uh∗L20,1≤OΔxN+1/2,forN≥0.This norm is defined as a truncation error term due to approximating polynomial, denoted by T2.After substituting the approximate solution into the numerical scheme (16), we obtain an ODE system,(22)ddtUh∗=LUh∗.Next, we apply the TVD Runge-Kutta scheme for integrating in time. We approximate Uh∗ at each time step by U^h. The order of accuracy for time integration depends on the TVD-RK order [1, 3, 17]. If the TVD-RK order N+1 is applied, a truncation error term T3 is given by(23)T3=Uh∗-U^h≤OΔtN+1.Combining all of the truncation error terms together, the total truncation error term T is estimated by(24)T=T1+T2+T3≤OΔxN+1/2,ΔtN+1.The total truncation error term is approaching zero as Δx→0 and Δt→0. Hence, the TVD-RK DG with the WAF method is consistent.


## 2.1. TVD-RK Discontinuous Galerkin Method

Consider the one-dimensional hyperbolic conservation law,(6)Ut+FxU=0.The computational domain 0,L is divided into K cells. We denote the jth cell by Ij=xj-1/2,xj+1/2, for j=1,…,K, with uniform cell size Δj=xj+1/2-xj-1/2. The cell center is that xj=xj+1/2+xj-1/2/2, where xj+1/2 and xj-1/2 are the left and the right of cell boundaries, respectively.Approximate solution is denoted byUh=hhqh.Multiplying (6) by a test function, vhx∈PNIj, where PNIj is the polynomial space degree N on the interval Ij, and replacing U by Uh and then taking the integration by parts over Ij, we obtain a weak form of numerical scheme as(7)∫Ij∂tUhvhdx-∫IjFUh∂xvhdx+F^j+1/2vhxj+1/2--F^j-1/2vhxj-1/2+=0,where the flux function F at the cell interfaces is approximated by F^, which is the function of Uh+ and Uh- at xj±1/2 as(8)F^j±1/2=F^Uhj±1/2-,Uhj±1/2+.Here Uh|j±1/2- and Uh|j±1/2+ 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 Uh can be written by(9)Uhx,t=∑m=0NUjmtφmx,where Ujmt is called the temporal coefficient. Now, the basis function φmx is defined by the Legendre polynomial Pm(x) of degree m over [-1,1].The test functionvhx is typically chosen to be the basis function; that is, vhx=φlxl=0N. So after applying Legendre’s properties, (7) is simplified to(10)dUjltdt=2l+1Δj∫IjFUh∂xφlxdx+2l+1Δj-1lF^j-1/2-F^j+1/2,for j=1,…,K, and l=0,…,N.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 N is applied, the TVD-RK method at least the order of N+1 must be applied to obtain the accuracy of order O(ΔxN+1) 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,F^WAF at the interface, xj+1/2, is defined by the integral average of a flux function F(U) at the half-time step,(11)F^j+1/2WAF=1Δx∫-Δx/2Δx/2FUj+1/2x,Δt2dx.It can be written in the wave structure form as(12)F^j+1/2WAF=∑k=1Nc+1ωkFj+1/2k,where Nc is the number of wave solutions in the Riemann problem and Fj+1/2(k) is the kth flux of the Riemann problem. For one-dimensional shallow water flows, we have Nc=2, where Fj+1/2(1)=FUj+1/2- and Fj+1/2(3)=FUj+1/2+. Flux component Fj+1/2(2) can be obtained from the HLL approach [6]. Weighted values, ωk, are defined by ωk=1/2ck-ck-1, where ck=SkΔt/Δx is the Courant number of wave k, c0=-1, cNc+1=1, and Sk is the speed of wave k.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 that(13)F^j+1/2WAF-TVD=12FUj+1/2-+FUj+1/2+-12∑k=1Ncsgn⁡ckϕj+1/2kΔFj+1/2k,where(14)ΔFj+1/2k=Fj+1/2k+1-Fj+1/2k.Here ϕj+1/2(k) 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,(15)∫Ij∂∂tUvxdx-∫IjFUv′xdx-FUj+1/2vxj+1/2-+FUj-1/2vxj-1/2+=0.Solution U is approximated by Uh∗. The test function v is estimated by vh. Flux function at cell interfaces is approximated in terms of numerical flux, F^j±1/2, so we obtain a numerical scheme of the following form:(16)∫Ij∂∂tUh∗vhxdx-∫IjFUh∗vh′xdx-F^j+1/2Uh∗j+1/2-,Uh∗j+1/2+vhxj+1/2-+F^j-1/2Uh∗j-1/2-,Uh∗j-1/2+vhxj-1/2+=0.Substituting Uh∗ by U in (16) and subtracting with (15), we obtain a truncation error term T1 as(17)
T
1
=
∫
I
j
∂
∂
t
U
-
U
v
h
x
d
x
-
∫
I
j
F
U
-
F
U
v
h
′
x
d
x
-
F
U
j
+
1
/
2
-
F
^
j
+
1
/
2
U
j
+
1
/
2
-
,
U
j
+
1
/
2
+
v
h
x
j
+
1
/
2
-
+
F
U
j
-
1
/
2
-
F
^
j
-
1
/
2
U
j
-
1
/
2
-
,
U
j
-
1
/
2
+
·
v
h
x
j
-
1
/
2
+
.When the solution is smooth, it must be continuous at xj+1/2 or Uj+1/2+=Uj+1/2-=Uj+1/2.The TVD weighted average flux is defined by(18)F^j+1/2Uj+1/2-,Uj+1/2+=12FUj+1/2-+FUj+1/2+-12∑k=12sgn⁡ckϕkΔFk,where F(1)=FUj+1/2-=FUj+1/2, F(3)=FUj+1/2+=FUj+1/2, and the component F(2)=FUj+1/2 which can be obtained from the HLL method [4, 6, 7]. Thus,(19)F^j+1/2Uj+1/2-,Uj+1/2+=12FUj+1/2+FUj+1/2-12∑k=12sgn⁡ckϕkΔFk=FUj+1/2.Similarly, we have F^j-1/2Uj-1/2-,Uj-1/2+=F(Uj-1/2). Then, the weighted average flux is consistent.For the TVD-RK DG method, the approximate solutionUh∗ in (16) is defined by(20)Uh∗x,t=∑k=0NUktφkx,where φk(x) is the Legendre polynomial degree k. Since the considering solution is smooth, we know from Theorem (3.1) in [1] that(21)U-Uh∗L20,1≤OΔxN+1/2,forN≥0.This norm is defined as a truncation error term due to approximating polynomial, denoted by T2.After substituting the approximate solution into the numerical scheme (16), we obtain an ODE system,(22)ddtUh∗=LUh∗.Next, we apply the TVD Runge-Kutta scheme for integrating in time. We approximate Uh∗ at each time step by U^h. The order of accuracy for time integration depends on the TVD-RK order [1, 3, 17]. If the TVD-RK order N+1 is applied, a truncation error term T3 is given by(23)T3=Uh∗-U^h≤OΔtN+1.Combining all of the truncation error terms together, the total truncation error term T is estimated by(24)T=T1+T2+T3≤OΔxN+1/2,ΔtN+1.The total truncation error term is approaching zero as Δx→0 and Δt→0. 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,(25)∫Ij∂tUhvhdx-∫IjFUh∂xvhdx+F^j+1/2vhxj+1/2--F^j-1/2vhxj-1/2+=∫IjGUhvhdx,with initial condition(26)∫IjUhx,0vhxdx=∫IjU0xvhxdx,where F^j±1/2 are the weighted average flux in (13). Source term is given by(27)GU=G1U,xG2U,x=0-ghbx.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 as(28)GpU,x=∑isiaU,xti′x,p=1,2,where si and ti are some functions which will be determined later.If the solution at steady state is stationary, thenaU,x in (28) can be decomposed by a1U,x and a2U,x, such that(29)a1≡h+b=constant,a2≡u=0.Since G1 is zero, we consider only G2,(30)G2U,x=-ghbx=-gh+bbx+12gb2x.From (28), we have that (31)s1=s1a1=-gh+b,s2=12g,t1x=b,t2x=b2.To balance the flux gradient and the source term approximation at steady state, it is required that(32)ddxFUx-∑isiaU,xtix=0,or(33)FUx-∑isiaU,xtix=α,where α is arbitrary constant.The integration of source term in (25) can be approximated by(34)∫IjG2Uh,xvhdx=∫Ij∑isiaU,xti′xvhxdx=∑i12siaUh,xj-1/2++siaUh,xj+1/2-tixi+1/2vhxi+1/2--tixi-1/2vhxi-1/2+-∫Ijtixvh′xdx+∑i∫IjsiaUh,x-12siaUh,xj-1/2++siaUh,xj+1/2-ti′xvhxdx.Functions auh,x and ti on the RHS of (34) can be approximated by ah(Uh,x) and tih, where ah(Uh,x) and tih are the L2 projection of aUh,x and ti into the space of Vh. Then, we obtain t1h=bh and t2h=bh2. Thus,(35)FUhx-∑isiaUh,xtixh=α.Then,(36)∫IjG2Uh,xvhdx=∑i12siahUh,xj-1/2++siahUh,xj+1/2-·t^ih,j+1/2vhxi+1/2--t^ih,j-1/2xi-1/2vhxi-1/2+-∫Ijtihxvh′xdx+∑i∫IjsiahUh,x-12siahUh,xj-1/2++siahUh,xj+1/2-tih′xvhdx.To satisfy the weighted average flux in (12), t^ih,j+1/2 must be modified. Here, we assume that(37)t^ih,j+1/2=∑k=1Nc+1ωkth,j+1/2k,or if the TVD version in (13) is applied, this term can be written by(38)t^ih,j+1/2=∑k=1Nc+1sgn⁡ckϕkth,j+1/2k.Note that ωk are weighted values in the WAF approximations and ϕ(k) is the WAF limiter function. So th,j+1/2(k) can be defined similar to Fj+1/2(k).For one-dimensional shallow water equations, we haveNc=2, where th,j+1/2(1)=thxj+1/2- and th,j+1/2(3)=thxj+1/2+. The approximation of th,j+1/2(2) in the intermediate region can be obtained by the HLL approach,(39)th,j+1/22=thxj+1/2-,ifSL≥0,SRthxj+1/2--SLthxj+1/2++SRSLthxj+1/2+-thxj+1/2-SR-SL,ifSL<0<SR,thxj+1/2+,ifSR≤0.Here SL and SR are signal velocities in the Riemann problem.At steady state, the solution is assumed to be stationary. It implies thath+b=c, where c is constant and u=0. From (35), we have that(40)FUh-∑isiahU,xtix=12ghh2-12gbh2-ghh+bhbh=12ghh+bh2=12gc2=α.This shows the suitable choices of a(Uh,x) and ti 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.
Sinceah(Uh,x) and ah(Uh,x)i∓1/2± are equal to the same constant c at each Gauss-Lobatto point over cell Ij, thus(41)∑i∫IjsiahUh,x-12siahUh,xj-1/2++siahUh,xj+1/2-tih′xvhdx=0.Hence, we obtain a truncation error term T as(42)
T
=
∫
I
j
F
U
h
x
,
t
v
h
′
x
d
x
-
F
^
j
+
1
/
2
v
h
x
j
+
1
/
2
-
+
F
^
j
-
1
/
2
v
h
x
j
-
1
/
2
+
+
∑
i
1
2
s
i
a
h
U
h
,
x
j
-
1
/
2
+
+
s
i
a
h
U
h
,
x
j
+
1
/
2
-
t
^
i
h
,
j
+
1
/
2
v
h
x
i
+
1
/
2
-
-
t
^
i
h
,
j
-
1
/
2
x
i
-
1
/
2
v
h
x
i
-
1
/
2
+
-
∫
I
j
t
i
h
x
v
h
′
x
d
x
.Using F^ in (12), t^i in (37), and ahUh,xj∓1/2±=c, we have that(43)T=∫IjFUhx,t-∑isictihvh′xdx-∑k=13ωkFj+1/2k-∑isic∑k=13ωkth,j+1/2k·vhxj+1/2-+∑k=13ωkFj-1/2k-∑isic∑k=13ωkth,j-1/2k·vhxj-1/2+=∫IjFUhx,t-∑isictihvh′xdx-∑k=13ωkFj+1/2k-∑isicth,j+1/2kvhxj+1/2-+∑k=13ωkFj-1/2k-∑isicth,j-1/2kvhxj-1/2+.
After applying condition (35) and rearranging terms, yields (44)T=αvhxj+1/2--vhxj-1/2+-αvhxj+1/2-+αvhxj-1/2+=0,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-5≤x≤5. The initial water depth is given by(45)hx,t=1,ifx≤0,0.6,ifx>0.The 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 t=2, with Δt=0.005. The root mean squared errors (RMS) are shown in Table 1.Table 1
RMS errors whenN=0, N=1, and N=2.
N
Numerical flux
K = 50
K = 100
K = 200
N = 0
HLL
0.0280
0.0215
0.0165
WAF
0.0182
0.0127
0.0094
N = 1
HLL
0.0130
0.0102
0.0072
WAF
0.0107
0.0093
0.0067
N = 2
HLL
0.0103
0.0079
0.0053
WAF
0.0100
0.0076
0.0052WhenN and K 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 K increases.If we fixK and vary N, the RMS error decreases as the polynomial degree increases.The water depth profiles using the HLL and the WAF methods att=2 when N=1 and K=100 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.Figure 1
Exact solution and water depth profile obtained by the TVD-RK DG HLL method.Figure 2
Exact solution and water depth profile obtained by the TVD-RK DG WAF method.


### 4.2. Flow over Irregular Bed

The uniform channel is length of 1500 m. The bottom elevation is irregular that is shown in Figure3. 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 h+b=16, with zero initial velocity. We set Δt=0.01 and run simulation until t=100.Figure 3
Stationary flow over irregular bed by well-balanced and non well-balanced TVD-RK DG with WAF method.It can be seen from Figure3 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(46)bx=0.2-0.05x-102,if8<x<12,0,otherwise.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 q=4.42 m2/s, while the downstream boundary is set by h=2m. The initial water depth is h+b=2 with zero initial velocity. Time step is that Δt=0.01. 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 K.Table 2
RMS errors whenN=1 for subcritical flow over a bump.
K = 50
K = 100
K = 200
K = 400
HLL
1.6954e − 03
4.5563e − 04
1.2546e − 04
3.5900e − 05
WAF
1.5816e − 03
4.1648e − 04
1.1474e − 04
3.3019e − 05The water depth and the bump profiles are shown in Figure4. 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.Figure 4
Water depth and bump profiles for subcritical flow.Transcritical Flow with Shock over a Bump. The upstream boundary is given by q=0.18m2/s, while the downstream boundary is set by h=0.33m. The initial condition is that h+b=0.33m. 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.Figure 5
Transcritical flow with shock over a bump.Transcritical Flow without Shock over a Bump. The upstream boundary is prescribed by q=1.53m2/s, while the downstream boundary is not specified. The initial conditions is that h+b=0.4m with zero initial velocity.The water depth profiles are shown in Figure6. The numerical result agrees well with the analytical solution. These results show the accuracy of the well-balanced scheme for solving transcritical flow problem.Figure 6
Transcritical flow without shock over a bump.


### 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 by(47)bx=0.25cos⁡10πx-1.5+1,if1.4<x<1.6,0.otherwise.The initial conditions are specified by(48)qx,0=0,hx,0=1-bx+ϵ,if1.1<x<1.2,1-bx,otherwise,where ϵ is a nonzero perturbation constant. The boundary conditions are transmissive boundaries. In this work, we consider the cases of ϵ=0.2 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 ±gh 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 thatt=0.7.The comparison of water depths between our results and the Leveque’s solutions is shown in Figures7 and 8 for ϵ=0.01 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.Figure 7
Quasi-stationary flow whenϵ=0.01.Figure 8
Quasi-stationary flow whenϵ=0.2.


### 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 by(49)bx=0,if0≤x≤10,0.1x-10,if10≤x≤20,1,if20≤x≤30.The initial conditions are given by (50)hx+bx=4,if0≤x<5,2,if5≤x≤30.Here, we set Δt=0.01 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 t=1 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.Figure 9
Flow over non-horizontal bed at time 1 s (a) and 4 s (b).
(a)
(b)


## 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-5≤x≤5. The initial water depth is given by(45)hx,t=1,ifx≤0,0.6,ifx>0.The 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 t=2, with Δt=0.005. The root mean squared errors (RMS) are shown in Table 1.Table 1
RMS errors whenN=0, N=1, and N=2.
N
Numerical flux
K = 50
K = 100
K = 200
N = 0
HLL
0.0280
0.0215
0.0165
WAF
0.0182
0.0127
0.0094
N = 1
HLL
0.0130
0.0102
0.0072
WAF
0.0107
0.0093
0.0067
N = 2
HLL
0.0103
0.0079
0.0053
WAF
0.0100
0.0076
0.0052WhenN and K 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 K increases.If we fixK and vary N, the RMS error decreases as the polynomial degree increases.The water depth profiles using the HLL and the WAF methods att=2 when N=1 and K=100 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.Figure 1
Exact solution and water depth profile obtained by the TVD-RK DG HLL method.Figure 2
Exact solution and water depth profile obtained by the TVD-RK DG WAF method.


## 4.2. Flow over Irregular Bed

The uniform channel is length of 1500 m. The bottom elevation is irregular that is shown in Figure3. 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 h+b=16, with zero initial velocity. We set Δt=0.01 and run simulation until t=100.Figure 3
Stationary flow over irregular bed by well-balanced and non well-balanced TVD-RK DG with WAF method.It can be seen from Figure3 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(46)bx=0.2-0.05x-102,if8<x<12,0,otherwise.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 q=4.42 m2/s, while the downstream boundary is set by h=2m. The initial water depth is h+b=2 with zero initial velocity. Time step is that Δt=0.01. 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 K.Table 2
RMS errors whenN=1 for subcritical flow over a bump.
K = 50
K = 100
K = 200
K = 400
HLL
1.6954e − 03
4.5563e − 04
1.2546e − 04
3.5900e − 05
WAF
1.5816e − 03
4.1648e − 04
1.1474e − 04
3.3019e − 05The water depth and the bump profiles are shown in Figure4. 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.Figure 4
Water depth and bump profiles for subcritical flow.Transcritical Flow with Shock over a Bump. The upstream boundary is given by q=0.18m2/s, while the downstream boundary is set by h=0.33m. The initial condition is that h+b=0.33m. 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.Figure 5
Transcritical flow with shock over a bump.Transcritical Flow without Shock over a Bump. The upstream boundary is prescribed by q=1.53m2/s, while the downstream boundary is not specified. The initial conditions is that h+b=0.4m with zero initial velocity.The water depth profiles are shown in Figure6. The numerical result agrees well with the analytical solution. These results show the accuracy of the well-balanced scheme for solving transcritical flow problem.Figure 6
Transcritical flow without shock over a bump.


## 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 by(47)bx=0.25cos⁡10πx-1.5+1,if1.4<x<1.6,0.otherwise.The initial conditions are specified by(48)qx,0=0,hx,0=1-bx+ϵ,if1.1<x<1.2,1-bx,otherwise,where ϵ is a nonzero perturbation constant. The boundary conditions are transmissive boundaries. In this work, we consider the cases of ϵ=0.2 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 ±gh 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 thatt=0.7.The comparison of water depths between our results and the Leveque’s solutions is shown in Figures7 and 8 for ϵ=0.01 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.Figure 7
Quasi-stationary flow whenϵ=0.01.Figure 8
Quasi-stationary flow whenϵ=0.2.


## 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 by(49)bx=0,if0≤x≤10,0.1x-10,if10≤x≤20,1,if20≤x≤30.The initial conditions are given by (50)hx+bx=4,if0≤x<5,2,if5≤x≤30.Here, we set Δt=0.01 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 t=1 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.Figure 9
Flow over non-horizontal bed at time 1 s (a) and 4 s (b).
(a)
(b)


## 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.



---

*Source: 591282-2015-07-07.xml*