Augmented versions of the HLL and HLLC Riemann solvers including source terms in one and two dimensions for shallow flow applications
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 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 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 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 with a horizontal x-axis and a vertical z-axis, the integral formulation per unit width is written aswith
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 followswith v representing the normal velocity to plane () 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 is split into two components, and standing for friction and bed level variation respectively.
Under the hypothesis 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 byandaccording 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 conservationwherewithwith
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 byfor a channel of length 25 m. The initial conditions are taken asDepending 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 bythe velocity isand 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)
Approximate Riemann solvers, parameter vectors, and difference schemes
J. Comput. Phys.
(1981)- et al.
Upwind methods for hyperbolic conservation laws with source terms
Comput. Fluids.
(1994) - et al.
On the numerical treatment of the source terms in the shallow water equations
Comput. Fluids.
(2000) - et al.
Flux difference splitting and the balancing of source terms and flux gradients
J. Comput. Phys.
(2000) - et al.
The surface gradient method for the treatment of source terms in the shallow-water equations
J. Comput. Phys.
(2001) - et al.
Adaptive quadtree simulation of shallow flows with wet–dry fronts over complex topography
Comput. Fluids.
(2009) - et al.
Exact solutions to the Riemann problem of the shallow water equations with a bottom step
Comput. Fluids.
(2001) - et al.
Exact solution of the Riemann problem for the shallow water equations with discontinuous bottom geometry
J. Comput. Phys.
(2008) - et al.
Numerical resolution of well-balanced shallow water equations with complex source terms
Adv. Water Res.
(2009) - et al.
Conservative numerical simulation of multicomponent transport in two-dimensional unsteady shallow water flow
J. Comput. Phys.
(2009)
Self adjusting grid methods for one-dimensional hyperbolic conservation laws
J. Comput. Phys.
Augmented Riemann solvers for the shallow water equations over variable topography with steady states and inundation
J. Comput. Phys.
Weak solutions for partial differential equations with source terms: application to the shallow water equations
J. Comput. Phys.
On Godunov-type methods near low densities
J. Comput. Phys.
A finite difference method for the computation of discontinuous solutions of the equations of fluid dynamics
Mater. Sb.
On upstream differencing and Godunov type methods for hyperbolic conservation laws
SIAM Rev.
Restoration of the contact surface in the HLL Riemann solver
Shock Waves
Cited by (71)
Effects of bed sediment conditions on debris flow propagation from the two-phase flow modelling perspective
2024, Advances in Water ResourcesThe entropy fix in augmented Riemann solvers in presence of source terms: Application to the Shallow Water Equations
2023, Computer Methods in Applied Mechanics and EngineeringA family of well-balanced WENO and TENO schemes for atmospheric flows
2023, Journal of Computational PhysicsDevelopment of POD-based Reduced Order Models applied to shallow water equations using augmented Riemann solvers
2023, Computer Methods in Applied Mechanics and EngineeringAn efficient semi-implicit friction source term treatment for modeling overland flow
2023, Advances in Water Resources