Numerical treatment in resonant regime for shallow water equations with discontinuous topography

Abstract

This paper deals with numerical treatments for the shallow water equations with discontinuous topography when the initial data belong to both supersonic region and subsonic region. This kind of data are present in both engineering and rivers, but they are not always well-treated in existing schemes. Our goal is to improve the well-balanced scheme constructed earlier in our work by introducing a computing corrector into the construction of the scheme. First, a further study in the construction of the well-balanced scheme reveals that the errors could make the approximate states near the critical surface that ought to be in one side of the critical surface fall into the other side. This qualitative change, though small, may cause much larger errors following stationary hydraulic jumps formed from these approximate states due to the jump of the bottom. Then, we introduce a corrector in the computing algorithm that selects the equilibrium states in the construction of the well-balanced scheme such that the approximate stationary hydraulic jumps always remain in the right region. Numerical tests show that the well-balanced method using an underlying numerical flux such as Lax–Friedrichs flux, FORCE, GFORCE, or Roe fluxes can approximate very well the exact solution even when the initial data are on both supercritical region and subcritical region.

Introduction

We are interested in numerical approximations in the resonant regime of solutions of the one-dimensional shallow water equations with discontinuous topography

$$\begin{array}{cc}\hfill & {\partial }_{t}h+{\partial }_x(\mathit{hu})=0\text{,}\hfill \\ \hfill & {\partial }_t(\mathit{hu})+{\partial }_x\left(h\left(u^2+\frac{\mathit{gh}}{2}\right)\right)+\mathit{gh}{\partial }_{x}a=0\text{,}\hfill \end{array}$$

where h is the height of the water from the bottom to the surface, u is the water velocity, g is the gravity constant, and $a=a(x)$ (with $x\in \mathbb{R}$) is the height of the bottom from a given level.

Supplementing the system (1.1) by the trivial equation

$${\partial }_{t}a=0\text{,}$$

we can transform the system of balance laws with source term (1.1) into a system of balance laws in nonconservative form. Then, the theory of shock waves, in particular, the Riemann problem, can be considered, see [28], [29]. The technique of supplementing an equation like (1.2) was proposed in [25], [26] for the model of fluid flows in a nozzle with variable cross-sections.

When the Riemann data belong to a strictly hyperbolic domain, that is, the Riemann data are either in the subsonic region or in the supersonic region, the well-balanced scheme [36], and the Godunov-type scheme [29] both can give good approximations to the solutions of (1.1). However, if the Riemann data belong to both supersonic region and subsonic region, then the Godunov-type scheme fails to approximate the exact solution. Moreover, as seen later, tests with the well-balanced scheme in [36] for the same initial data give unsatisfactory results as well. Mathematically, on one hand, it is worth to improve a numerical method such that it could be “robust” and could work properly for any data. On the other hand, in engineering structures and rivers, water at high speed may discharge into a zone of the river or engineered structures which can only sustain a lower speed. This causes stationary hydraulic jumps to occur, like outfalls of dams and irrigation works. In these hydraulic jumps the water is supercritical (the Froude number larger than 1) and the water height is lower before the jump, and the water is subcritical (the Froude number less than 1) and the water height is higher after the jump. Stationary hydraulic jumps can eventually occur when a rapid flow encounters a submerged object which pushes the water upwards. See [9], for example. Thus, it is important to have a numerical method that can work well for data in both supercritical and subcritical regions and the bottom suffers a discontinuity. The simplest model of discontinuous bottoms is the one with piece-wise constant levels:

$$a(x)=\left\{\begin{array}{cc}{a}_L\text{,}\hfill & x<0\text{,}\hfill \\ a_R\text{,}\hfill & x>0\text{,}\hfill \end{array}\right.$$

where $a_L\ne a_R$ are constants.

Our purpose in this work is to present a numerical treatment for shallow water equations with discontinuous topography (1.1), (1.2) when the initial data belong to both supercritical region and subcritical region by improving the well-balanced scheme [36]. We first carry out an investigation into the computing algorithm in the construction of this scheme which selects the admissible equilibrium state resulted by an admissible stationary jump. Then we realize that a computer selection procedure may be different from the theoretical procedure, probably due to the propagation of errors for states near the critical surface. The errors could force the approximate states near the critical surface that ought to be in one side of the critical surface fall into the other side. This initial qualitative change, possibly small, may cause much larger errors by admissible stationary discontinuities jumping from these states that always remain in the same region. In this case, the well-balanced scheme cannot give a satisfactory result and could eventually provide an approximate solution having a big difference with the exact solution. To deal with this, we build a computing corrector into the computing algorithm selecting admissible equilibrium states in the construction of the well-balanced scheme. Then, tests are presented with the Lax–Friedrichs scheme, FORCE, GFORCE, and Roe schemes. Let us recall that the FORCE and GFORCE schemes are constructed by using a convex linear combination of the Lax–Friedrichs and the Lax–Wendroff numerical fluxes, see [33], [34], [10]. These tests show that the well-balanced method using an underlying scheme among the one mentioned above can give good approximations of the solutions of (1.1).

There are many works for the shallow water Eqs. (1.1). Theoretically, solutions of (1.1) can be understood in the sense of nonconservative product by [11]. See [26], [30], [27], [35], [19], [18], [13], [2], [3], [28] for the works concerning the study of shock waves in various systems of balance las in nonconservative form. Numerical treatment for shallow water equations was considered by an early work [16], and then developed in [8], [36], [20], [21], [12], [31]. A well-balanced high-order scheme for the shallow water equations with topography and dry areas was presented in [14]. Well-balanced numerical schemes for a single conservation law with source term were studied in [17], [6], [7], [15], [4]. Well-balanced schemes for fluid flows in a nozzle with variable cross-section were constructed in [23], [22], [20], [21]. Well-balanced schemes for multi-phase flows and other models were studied in [5], [24], [32], [1], [37], [38], [39], etc.

The organization of this paper is as follows. Section 2 provides us with preliminaries of the system (1.1), (1.2), together with the selecting procedure of the admissible stationary waves. In Section 3 we recall the construction of the well-balanced scheme. We provide an analysis into the construction of the scheme and predict a cause of the failure for the well-balanced scheme when data are in the resonance regime. Then we build a computing corrector. Section 4 deals with numerical test cases, which address all possibilities of the Riemann data on both sides of the critical surface. The well-balanced scheme equipped by the computing corrector using an underlying numerical flux among the Lax–Friedrichs, FORCE, GFORCE, and Roe numerical fluxes is tested with different mesh sizes. All tests show that the well-balanced scheme equipped by the computing corrector can work very well for data in the resonance regime. Finally, in Section 5 we draw some conclusions.