Augmented versions of the HLL and HLLC Riemann solvers including source terms in one and two dimensions for shallow flow applications

https://doi.org/10.1016/j.jcp.2012.06.031Get rights and content

Abstract

Shallow water flows are found in a variety of engineering problems always dominated by the presence of bed friction and irregular bathymetry. These source terms determine completely the possible evolution of the flooded area in time. It is well known that appropriate numerical schemes for this type of flows must be well-balanced. Well-balanced numerical schemes are based on the preservation of cases of quiescent equilibrium over variable bed elevation. Commonly they are formulated as an adaptation of numerical solvers defined for cases without source terms. This procedure is insufficient when applied to real situations. Then, it is possible to argue that appropriate numerical schemes cannot arise directly from those derived from the simplest homogeneous case without source terms. New solutions are presented in this work by defining weak solutions that include the presence of source terms. To do that, the solvers presented in this work extend the number of waves in the well known HLL and HLLC solvers involving a stationary jump in the solution. This is done without modifying the original solution vector of conserved quantities. The resulting approximate Riemann solvers include variable bed level surface and friction. Solvers are systematically assessed via a series of test problems with exact solutions for one and two dimensions, including steady and unsteady flow configurations, variation of the flooded area in time and comparisons with experimental data. The obtained results point out that the new method is able to predict faithfully the overall behavior of the solution and of any type of waves.

Introduction

There is a wide range of physical situations where the finite volume method in combination with conservative schemes has proved to be a successful tool. When modeling flow in open channels and rivers, if using shallow water assumptions, the physics of the problem can be mathematically represented by a system of conservation laws of hyperbolic character. The development of conservative schemes for the solution of conservation laws by Godunov [1] was the starting point for the development of most well known numerical schemes. In general, Godunov type methods use solutions to the Riemann problem at cell interfaces.

Among Godunov type methods, the HLL developed by Harten et al. [2] is maybe the most disseminated approximate Riemann solver after Roe’s method [3] (which preceded HLL). The HLL method solves the flux at cell interfaces using an intermediate or star region between cells. The numerical flux arises directly from the governing equations. To do that, a control volume centered on a cell interface is integrated over space and time. By simply bounding the star region by two propagating waves, the flux across the interface is mathematically determined. Other approximate Riemann solvers, as Roe’s method [3], provide solutions to the Riemann problem at cell interfaces by decomposing the flux difference vector into positive and negative components. These components are defined according to the sign of the eigenvalues of the Jacobian matrix of the system. In contrast, the HLL method gives approximate solutions with less effort than flux-splitting methods.

The HLL method considers only two propagating waves, which is a valid assumption for two variable hyperbolic systems, associated to purely one-dimensional unsteady flow. When moving to two-dimensional unsteady flow, the set of governing equations increases, and this assumption is inappropriate [2] due to the necessity of two intermediate states, separated by a middle propagating wave. The solution to this problem was proposed by Toro et al. [4] by the introduction of a contact discontinuity, leading to the well known HLLC method where C stands for Contact wave.

The system of equations in realistic shallow water models includes source terms, that is, terms that are functions of the vector of unknowns. It is possible to find cases, where even naive discretizations of the source terms work reasonably well, but there are some well documented situations in which only adequate schemes can perform adequately. When solving real problems one is likely to encounter all sorts of situations, and naive schemes completely fail. In this framework many efforts have been reported on the search for the best methods able to handle source terms. They have resulted on a wide family of well-balanced numerical schemes. Their keystone is based on the preservation of motionless steady state over irregular geometries. Bermudez and Vázquez-Cendón [5] showed that the upwind treatment of source terms produces better results than the pointwise method in a rectangular channel with bed slope variation. Since then, the upwind approach for source term discretization has been widely used. The upwind approach was extended to to nonprismatic channels and high-resolution schemes García-Navarro and Vázquez-Cendón [6], Hubbard and García-Navarro [7]. Lee and Wright [8] integrated in the numerical flux the representation of the source terms by a straightforward modification of the governing equations. In contrast, other authors [9], [10] ensured quiescent equilibrium by using a centered discretization of the bed slope source term. Zhou et al. [9] ensured the well-balanced property by using the water surface gradient as the basis for data reconstruction in combination with the HLL solver. Liang and Borthwick [10] applied the HLLC solver [4] writing the set of equations in a deviatoric formulation and used the two-step unsplit MUSCL-Hancock method to formulate a second-order scheme. It is worth remarking that other alternatives are possible. For instance, Yulong and Chi-Wang [11] designed a finite difference high order WENO Lax–Friedrichs scheme, keeping well-balanced without reducing the high order accuracy of the scheme.

All these numerical schemes were constructed to ensure good balance in cases of quiescent equilibrium retaining as valid the approximate solutions described in earlier works [3], [2], even though the presence of source terms modifies the characteristics of the exact solution [12], [13]. Their presence changes the solution of the Riemann problem, and the characteristic constant U solution defined for the homogeneous case in [3], [2] is no longer valid. In fact, this situation can be clearly noticed in situations of quiescent equilibrium over a step in the bottom, where the velocity must be nil on both sides of the RP, but the solution for the water depth is discontinuous. The formulation of this basic RP situation is impossible with a single U value for both water depth and discharge. Furthermore, in dam break problems over a bottom step, although the unit discharge must be constant across the discontinuity, neither the water depth nor the surface elevation are.

Therefore, even ensuring the discrete equilibrium formulated in well-balanced schemes, the direct application of the method derived for the homogeneous case to cases with source terms leads to important difficulties in real applications. One of them is the appearance of negative values of water depth in wet/dry fronts. Simple methods can be used to modify the local bed slope so as not to induce spurious flow in dry areas [14], [10], [15] in particular cases. But unphysical results can also appear in initially wet/wet Riemann problems [16]. If these cases are ignored by setting a nil water depth, mass conservation is violated. In general, numerical solvers in [5], [6], [7], [8], [9], [10] are unable to explain the effect of the source terms in the stability region.

As a consequence numerical stability seems controlled by a heavier restriction than the classical Courant–Friedrichs–Lewy (CFL) condition [17] on the time step size that may lead to inefficient computations. It is possible to find cases where unphysical results appear independently of the time step selected. Even when trying to overcome partially these issues by selecting a tuning parameter over the variables, in most cases defining tolerances for the water depth [10], the results shown differ depending on the value selected. This is a major problem for practitioners. Difficulties may also arise from the presence of relatively large friction terms that reduce considerably the stability region of the explicit schemes [18]. In [15] the HLLC solver was extended to avoid this problem using a splitting implicit technique over the friction source term. Although implicit formulations strongly help to relax the friction related stability restrictions, they are only feasible in the context of separated discretization [18], and do not allow a perfect balance in non-zero velocity steady states.

It is worth mentioning that, even in the absence of source terms, transcritical flows require entropy fixes in order to provide accurate results [19]. While an adequate selection of wave speed estimates is sufficient for the HLLC solver, the Roe solver requires further modifications of the wave speed estimates given by the definition of the approximate Jacobian matrix [3]. In any case, this option is derived once the internal structure of the associate Riemann problem is completely defined. It can be argued that the presence of source terms justifies the construction of new weak solutions appropriate to the nature of the equations, rather than the use of those constructed for the simple, homogeneous case. Recently augmented approximate Riemann solvers have been presented [20], [21], [22] in order to include the presence of the source terms in the weak solution. In [22] problems with m equations were presented using an augmented Roe’s Riemann solver defined by m+1 states, by means of an extra stationary wave linked to the source term. The augmented Roe’s Riemann solver in [22] completely explained the previously reported problems and has proved an accurate and robust tool. In [24] it was shown how the introduction of bed slope and friction source terms in the approximate solution allows a correct evaluation of the stability region even in cases of complex rheology and complex geometry.

Following the ideas in [22], the work presented in this paper tries to recover the simplicity and clarity of the HLL and HLLC solvers by means of introducing the presence of the source term in the weak solution. The resulting numerical schemes, designed to be first order accurate, will be referred to as HLLS and HLLCS.

Section snippets

One-dimensional flow

In this section, the mathematical model for one-dimensional flow is briefly presented. Fig. 1 considers a generic control volume for a horizontal flow over a bottom step where the velocity is depth averaged. The relevant integral formulation of the model derives from the depth-averaged equations of mass and momentum. In a Cartesian coordinate system (x,z) with a horizontal x-axis and a vertical z-axis, the integral formulation per unit width is written astxLxRUdx+FxR-FxL-xLxRSτdx-Sz=0withU=h

One-dimensional flow with shear waves and species equations

The differential formulation can be expanded allowing for shear waves and extra species-like equations by enlarging the set of equations as followsU=hhuhvhϕF=huhu2+12gh2huvhuϕS=0S200with v representing the normal velocity to plane (x,z) and ϕ representing the scalar depth-averaged volumetric concentration.

In presence of contact discontinuities, shear waves and material interfaces, the HLLS solver smears the numerical solutions, as it does not consider the multiple eigenvalue u, that appears in

Integration of source terms

Once the structure of the wave solution is defined, we are interested in the analysis of the consequences of different discretizations of the source terms when inserted in the approximate Riemann solver before the updating step. In order to distinguish between frictional and gravity forces in the bed, the source term S¯i+1/2 is split into two components, S¯τ,i+1/2 and S¯z,i+1/2 standing for friction and bed level variation respectively.

Under the hypothesis (a) of smooth variation of the

Wave-speed estimates

When computing the numerical fluxes, the selection of estimates for the wave speeds is of utmost importance, as they can compromise the quality of the numerical results. The direct wave speed estimates are the simplest methods providing minimum and maximum signal velocities. The simplest of these is provided byλL=uL-cLλR=uR+cRandλL=min(uL-cL,uR-cR),λR=max(uL+cL,uR+cR)according to Davis [25]. In the context of gas dynamics, Einfeldt [26] combined the direct estimates in (87) with the Roe

Stability region

In the homogeneous case, if appropriate estimates for the celerity waves are not selected, unphysical values of water depth may appear. The estimates suggested by Einfeldt [26] lead to positively conservative methods [27]. In case of quiescent equilibrium over variable bed level no problem may be detected when applying the approximate solvers defined in this work, as integration of the source term in (19) can be done exactly.

When dealing with stationary cases with non-zero velocity, the

Extension to multi-dimensions

In this section the 2D extension of the model is presented. The numerical scheme in written in a general form, and definitions, written in terms of an unstructured mesh, can be easily written for a Cartesian mesh if desired. The 2D extension of the model is formulated by means of depth averaged equations expressing volume conservation and momentum conservationUt+F(U)x1+G(U)x2=SwhereU=hhu1hu2hϕF=hu1hu12+12gh2hu1vhu1ϕG=hu2hu2u1hu22+12gh2hu2ϕwithS=0,-ghzx1-cfu1u,-ghzx2-cfu2u,0with (u1

Steady flow over a hump

The purpose of this test case is to compare the performance of the HLLCS solver presented in this work with other solvers presented by other authors in previous works, as this case has been widely used to test numerical schemes for shallow water equations [9], [10], [11], [31]. The bed level is given byz(x)=0.2-0.05(x-10)2if8x120otherwisefor a channel of length 25 m. The initial conditions are taken ash(x,0)=0.5-z(x)u(x,0)=0Depending on the boundary conditions, the flow can be subcritical or

Long wave resonance in a circular parabolic frictionless basin

The first test case presented is motivated by the necessity of checking the performance of the HLLCS solver in a problem of transient boundaries dominated by the nonconservative terms and with analytical solution. The analytical solution of a long wave resonating in a circular, frictionless parabolic basin was presented in [35] for the shallow water equations. The free surface displacement is given byζ(r,t)=ζ0(1-A2)1/21-Acosωt-1-r2a21-A2(1-Acosωt)2-1the velocity isu(r,t)=sinωtAωr2-2A2and the

Conclusions

Based on the HLLC method, approximate Riemann solvers for the shallow water equations with source terms have been presented. They are defined by incorporating an extra wave that accounts for the presence of the source term. The definition of the new numerical fluxes that include explicitly this new wave, leads as a consequence to the desired well-balance property. Therefore, the method is not constructed departing, as usual, from numerical relationships that force equilibrium among fluxes and

Acknowledgment

This work has been funded by the Spanish Ministry of Science and Technology under research projects CGL2011-28590.

References (35)

  • A. Harten et al.

    Self adjusting grid methods for one-dimensional hyperbolic conservation laws

    J. Comput. Phys.

    (1983)
  • D.L. George

    Augmented Riemann solvers for the shallow water equations over variable topography with steady states and inundation

    J. Comput. Phys.

    (2008)
  • J. Murillo et al.

    Weak solutions for partial differential equations with source terms: application to the shallow water equations

    J. Comput. Phys.

    (2010)
  • B. Einfeldt et al.

    On Godunov-type methods near low densities

    J. Comput. Phys.

    (1991)
  • S.K. Godunov

    A finite difference method for the computation of discontinuous solutions of the equations of fluid dynamics

    Mater. Sb.

    (1959)
  • A. Harten et al.

    On upstream differencing and Godunov type methods for hyperbolic conservation laws

    SIAM Rev.

    (1983)
  • E.F. Toro et al.

    Restoration of the contact surface in the HLL Riemann solver

    Shock Waves

    (1994)
  • Cited by (71)

    View all citing articles on Scopus
    View full text