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

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 correctlywithoutspuriousoscillationsneartheshockfront.ThispresentstheotherfluxapproximationsintheTVD-RKDGmethod forshallowwatersimulations.


Introduction
Hyperbolic balanced law for one-dimensional problem is where , , 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 by ℎ  +   = 0, where ℎ 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 ) ,  () = ( 0 −ℎ  ) . ( The total variation diminishing Runge-Kutta discontinuous Galerkin (TVD-RK DG) method can be applied to solve the shallow water equations; see [1][2][3][4][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 2 Modelling and Simulation in Engineering 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][7][8][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 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 Cproperty means that the numerical solution must satisfy still water condition at steady state given by 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 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.

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.

TVD-RK Discontinuous Galerkin
Method.Consider the one-dimensional hyperbolic conservation law, The computational domain (0, ) is divided into  cells.We denote the th cell by  6) by a test function, V ℎ () ∈   (  ), 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 as where the flux function  at the cell interfaces is approximated by F, which is the function of  + ℎ and  − ℎ at  ±1/2 as Here  ℎ | − ±1/2 and  ℎ | + ±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  ℎ can be written by where    () is called the temporal coefficient.Now, the basis function   () is defined by the Legendre polynomial   () of degree  over [−1, 1].
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  + 1 must be applied to obtain the accuracy of order (Δ +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.

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.
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 where Here  () +1/2 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.

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 V is estimated by V ℎ .Flux function at cell interfaces is approximated in terms of numerical flux, F±1/2 , so we obtain a numerical scheme of the following form: Substituting  * ℎ by  in ( 16) and subtracting with (15), we obtain a truncation error term  1 as When the solution is smooth, it must be continuous at The TVD weighted average flux is defined by where  (1) = ( − +1/2 ) = ( +1/2 ),  (3) = ( + +1/2 ) = ( +1/2 ), and the component  (2) = ( +1/2 ) 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 by where   () is the Legendre polynomial degree .Since the considering solution is smooth, we know from Theorem (3.1) in [1] that This norm is defined as a truncation error term due to approximating polynomial, denoted by  2 .
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  + 1 is applied, a truncation error term  3 is given by Combining all of the truncation error terms together, the total truncation error term  is estimated by The total truncation error term is approaching zero as Δ → 0 and Δ → 0. Hence, the TVD-RK DG with the WAF method is consistent.

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 condition where F±1/2 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 as where   and   are some functions which will be determined later.
If the solution at steady state is stationary, then (, ) in (28) can be decomposed by  1 (, ) and  2 (, ), such that From (28), we have that (31) To balance the flux gradient and the source term approximation at steady state, it is required that or where  is arbitrary constant.The integration of source term in (25) can be approximated by Functions ( ℎ , ) and   on the RHS of (34) can be approximated by  ℎ ( ℎ , ) and (  ) ℎ , where  ℎ ( ℎ , ) and (  ) ℎ are the  2 projection of ( ℎ , ) and   into the space of  ℎ .Then, we obtain ( 1 ) ℎ =  ℎ and ( 2 ) ℎ =  2 ℎ .Thus, Then, To satisfy the weighted average flux in ( 12), ( t ) ℎ,+1/2 must be modified.Here, we assume that or if the TVD version in ( 13) is applied, this term can be written by Note that   are weighted values in the WAF approximations and  () is the WAF limiter function.So  () ℎ,+1/2 can be defined similar to  () +1/2 .
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  ℎ ( ℎ , ) ± ∓1/2 are equal to the same constant  at each Gauss-Lobatto point over cell   , thus Hence, we obtain a truncation error term  as Using F in ( 12), ( t ) in (37), and  ℎ ( ℎ , ) ± ∓1/2 = , we have that After 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.

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.

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 ≤  ≤ 5.The initial water depth is given by 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  = 2, with Δ = 0.005.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  = 2 when  = 1 and  = 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.

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 ℎ +  = 16, with zero initial velocity.We set Δ = 0.01 and run simulation until  = 100.
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.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  = 4.42 m 2 /s, while the downstream boundary is set by ℎ = 2 m.The initial water depth is ℎ +  = 2 with zero initial velocity.Time step is that Δ = 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 .
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  = 0.18 m 2 /s, while the downstream boundary is set by ℎ = 0.33 m.The initial condition is that ℎ +  = 0.33 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 wellbalanced 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-wellbalanced scheme.Transcritical Flow without Shock over a Bump.The upstream boundary is prescribed by  1.53 m 2 /s, while the downstream boundary is not specified.The initial conditions is that ℎ +  = 0.4 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.

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  The initial conditions are specified by 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 ±√ℎ at the early stage.A standard scheme which is not well-balanced usually faces with some difficulties to capture correctly the wave 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  = 0.7.
The comparison of water depths between our results and the Leveque's solutions is shown in Figures 7 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 quasistationary flow with initial disturbance.

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 (50) Here, we set Δ = 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  = 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.

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

= 2 Figure 1 :Figure 2 :
Figure 1: Exact solution and water depth profile obtained by the TVD-RK DG HLL method.

Figure
Figure Water depth and bump profiles for subcritical flow.

5 :
Transcritical flow with shock over a bump.

Figure 6 :
Figure 6: Transcritical flow without shock over a bump.

Figure 9 :
Figure 9: Flow over non-horizontal bed at time 1 s (a) and 4 s (b).

Table 1 :
RMS errors when  = 0,  = 1, and  = 2.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.