Non-dispersive conservative regularisation of nonlinear shallow water (and isothermal Euler) equations

TL;DR

This paper introduces a novel non-dispersive, conservative regularisation for shallow water and Euler equations that preserves key physical quantities and produces smooth, non-oscillatory shocks, demonstrated through dam-break simulations.

Contribution

It presents a new variational, non-dissipative regularisation method that maintains shock speeds and conserves mass, momentum, and energy in hyperbolic PDEs.

Findings

Regularised shocks propagate at the same speed as original shocks.

The method produces smooth, non-oscillatory hydraulic jumps.

Numerical tests confirm the effectiveness on dam-break problems.

Abstract

A new regularisation of the shallow water (and isentropic Euler) equations is proposed. The regularised equations are non-dissipative, non-dispersive and possess a variational structure. Thus, the mass, the momentum and the energy are conserved. Hence, for instance, regularised hydraulic jumps are smooth and non-oscillatory. Another particularly interesting feature of this regularisation is that smoothed `shocks' propagates at exactly the same speed as the original discontinuous ones. The performance of the new model is illustrated numerically on some dam-break test cases, which are classical in the hyperbolic realm.

Tables

Table 1. Table 1. Numerical and physical parameters used in the dam-break and shock wave simulations.

Parameter	Dam-break value	Shock wave value
Regularisation parameter, $ϵ$	${10}^{-3}$	${10}^{-2}$
Domain half-length, $\ell /d$	$25$	$25$
Number of Fourier modes, $N$	$1024$	$1024$
Initial transition length, $1/\delta d$	$1$	$1$
Final simulation time, $T\sqrt{g/d}$	$15$	$5$
Water depth on the left, $h_l/d$	$1.5$	$1.5$
Water depth on the right, $h_r/d$	$1$	$1$
Horizontal velocity on the left, $u_l/\sqrt{gd}$	$0$	$1$
Horizontal velocity on the right, $u_r/\sqrt{gd}$	$0$	$0.5435645\dots$

Equations

$$\displaystyle h_{t}\ +\ \partial_{x}\!\left[\,h\,u\,\right]\$$

$${\partial }_t\left[hu\right]+{\partial }_x\left[h{u}^2+\frac{1}{2}g{h}^2\right]$$

$${\mathcal{L}}_0\stackrel{\mathrm{def}}{=}\frac{1}{2}h{u}^2-\frac{1}{2}g{h}^2+\left\{h_t+{\left[hu\right]}_x\right\}\varphi ,$$

$${\mathcal{L}}_1\stackrel{\mathrm{def}}{=}\frac{1}{2}h{u}^2+\frac{1}{6}{h}^3{u}_x^2-\frac{1}{2}g{h}^2+\left\{h_t+{\left[hu\right]}_x\right\}\varphi .$$

$${\mathcal{L}}_2\stackrel{\mathrm{def}}{=}\frac{1}{2}h{u}^2+\left(\frac{1}{6}\text{/}+\text{/}\frac{1}{4}\text{/}\text{textbeta}\right)h^3{u}_x^2-\frac{1}{2}g{h}^2\left(1\text{/}+\text{/}\frac{1}{2}\text{/}\text{textbeta}\text{/}{h}_x^2\right)+\left\{h_t+{\left[hu\right]}_x\right\}\varphi ,$$

$${\mathcal{L}}_3\stackrel{\mathrm{def}}{=}\frac{1}{2}h{u}^2+\left(\frac{1}{6}\text{/}+\text{/}\frac{1}{4}\text{/}{\text{textbeta}}_1\right)h^3{u}_x^2-\frac{1}{2}g{h}^2\left(1\text{/}+\text{/}\frac{1}{2}\text{/}{\text{textbeta}}_2\text{/}{h}_x^2\right)+\left\{h_t+{\left[hu\right]}_x\right\}\varphi .$$

$$\delta \varphi :$$

$$\delta u:$$

$$\delta h:$$

$$+\left(\frac{1}{2}+\frac{3}{4}{\text{textbeta}}_1\right)h^2{u}_x^2-\frac{1}{2}{\text{textbeta}}_{2}gh{h}_x^2+\frac{1}{2}{\text{textbeta}}_{2}g[h^2{h}_x{]}_x,$$

$$\displaystyle h_{t}\ +\ \partial_{x}\!\left[\,h\,u\,\right]\$$

$${\varphi }_{xt}+{\partial }_x\left[u{\varphi }_x-\frac{1}{2}{u}^2+gh-\left(\frac{1}{2}+\frac{3}{4}{\text{textbeta}}_1\right)h^2{u}_x^2-\frac{1}{2}{\text{textbeta}}_{2}g(\text{/}{h}^2\text{/}{h}_{xx}\text{/}+\text{/}h\text{/}{h}_x^2\text{/})\right]$$

$$\displaystyle u\ -\ \left({\textstyle{1\over 3}}\ +\ {\textstyle{1\over 2}}\,\text{\textbeta}_{1}\right)\!\left(\,3\,h\,h_{x}\,u_{x}\,+\,h^{2}\,u_{xx}\,\right)\$$

$$u_t+u{u}_x+g{h}_x+\frac{1}{3}{h}^{-1}{\partial }_x\left[h^2\Gamma \right]=0,$$

$$\Gamma \stackrel{\mathrm{def}}{=}\left(1+\frac{3}{2}{\text{textbeta}}_1\right)h\left[u_x^2-u_{xt}-u{u}_{xx}\right]-\frac{3}{2}{\text{textbeta}}_{2}g\left[h{h}_{xx}+\frac{1}{2}{h}_x^2\right].$$

$$\displaystyle 0\$$

$$(\frac{1}{6}+\frac{1}{4}{\text{textbeta}}_1)h^2{u}_x^2+gh\left(1+\frac{1}{4}{\text{textbeta}}_2{h}_x^2\right)+\frac{1}{3}h\Gamma \}hu+\frac{1}{2}{\text{textbeta}}_{2}g{h}^3{h}_x{u}_x].$$

$$\displaystyle\text{\texteta}_{t}\ +\ d\,u_{x}\$$

$$\frac{c_0^2}{gd}=\frac{2+{\text{textbeta}}_2(kd{)}^2}{2+(\frac{2}{3}+{\text{textbeta}}_1)(kd{)}^2}.$$

$${\mathcal{L}}_{ϵ}\stackrel{\mathrm{def}}{=}\frac{1}{2}h{u}^2-\frac{1}{2}g{h}^2+\left\{h_t+{\left[hu\right]}_x\right\}\varphi +\frac{1}{2}ϵh^2\left(h{u}_x^2-g{h}_x^2\right),$$

$$\displaystyle h_{t}\ +\ \partial_{x}\!\left[\,h\,u\,\right]\$$

$$\displaystyle\partial_{t}\!\left[\,h\,u\,\right]\,+\ \partial_{x}\!\left[\,h\,u^{2}\,+\,{\textstyle{1\over 2}}\,g\,h^{2}\,+\,\upvarepsilon\,\mathscr{R}\,h^{2}\,\right]\$$

$$h\left(u_x^2-u_{xt}-u{u}_{xx}\right)-g\left(h{h}_{xx}+\frac{1}{2}{h}_x^2\right)$$

$${\partial }_t\left[\frac{1}{2}h{u}^2+\frac{1}{2}g{h}^2+\frac{1}{2}ϵh^3{u}_x^2+\frac{1}{2}ϵg{h}^2{h}_x^2\right]+$$

$${\partial }_x\left[\left\{\frac{1}{2}{u}^2+gh+\frac{1}{2}ϵh^2{u}_x^2+\frac{1}{2}ϵgh{h}_x^2+ϵh\mathcal{R}\right\}hu+\text{upvarepsilon}g{h}^3{h}_x{u}_x\right]$$

$${\partial }_t\left[h\mathcal{U}\right]+{\partial }_x\left[hu\mathcal{U}+\frac{1}{2}g{h}^2-\text{upvarepsilon}{h}^2\left(2h{u}_x^2+gh{h}_{xx}+\frac{1}{2}g{h}_x^2\right)\right]$$

$${\mathcal{U}}_t+{\partial }_x\left[u\mathcal{U}-\frac{1}{2}{u}^2+gh-\frac{3}{2}\text{upvarepsilon}{h}^2{u}_x^2-\text{upvarepsilon}gh\left(h{h}_{xx}+h_x^2\right)\right]$$

$$uh=\text{constant}\stackrel{\mathrm{def}}{=}-cd⟹u=-cd/h,$$

$$\displaystyle\frac{2\,\epsilon\,(\,h\,h_{\,xx}\,-\,h_{\,x}^{\,2}\,)}{g\,h\,/\,c^{\,2}}\ -\ \frac{\epsilon\,(\,2\,h\,h_{\,xx}\,+\,h_{\,x}^{\,2}\,)}{d^{2}\,/\,h^{2}}\ +\ \frac{2\,c^{\,2}}{g\,h}\ +\ \frac{h^{2}}{d^{\,2}}\$$

$$\displaystyle\frac{\epsilon\,(\,2\,h\,h_{\,xx}\,-\,h_{\,x}^{\,2}\,)}{2\,g\,h^{2}\,/\,c^{\,2}\,d}\ -\ \frac{\epsilon\,(\,h\,h_{\,xx}\,+\,h_{\,x}^{\,2}\,)}{d\,/\,h}\ +\ \frac{c^{\,2}\,d}{2\,g\,h^{2}}\ +\ \frac{h}{d}\$$

$$\text{upvarepsilon}{\left(\frac{\mathrm{d}h}{\mathrm{d}\text{/}x}\right)}^2=\frac{\mathcal{F}-(1+{\mathcal{C}}_1+2\mathcal{F})(h/d)+(2+2{\mathcal{C}}_2+\mathcal{F})(h/d{)}^2-(h/d{)}^3}{\mathcal{F}-(h/d{)}^3},$$

$$x=\text{/}{\int }_0^{\xi }{\left[\frac{\mathcal{F}\text{/}{d}^3-h^3({\xi }^{\prime })}{(\mathcal{F}-1)d^3}\right]}^{\frac{1}{2}}\text{/}\mathrm{d}\text{/}{\xi }^{\prime },\frac{\text{texteta}(\xi )}{d}=(\mathcal{F}-1){\mathrm{sech}}^2\left(\kappa \xi \right),(\kappa \text{/}d{)}^2=\frac{1}{\text{upvarepsilon}}.$$

$$\left(u-\dot{s}\right)[[u_{xx}]]+g[[h_{xx}]]=0,$$

$$\left(u-\dot{s}\right)[[h_{xx}]]+h[[u_{xx}]]=0,$$

$$\dot{s}(t)=u(x,t)\pm \sqrt{gh(x,t)}\text{ at }x=s(t).$$

$$[[u_{xt}]]+u[[u_{xx}]]+g[[h_{xx}]]=0,$$

Full text

Key words and phrases: Shallow water flows; Dispersionless; Conservative; Regularisation.

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Model equations

3. 3 Regularised shallow water equations

4. 3.1 Linear dispersion relation

5. 3.2 Regularised Saint-Venant (rSV) equations

6. 3.3 Steady motion

7. 3.4 Rankine–Hugoniot type conditions

8. 3.5 Remarks on the total energy

9. 4 Numerical illustrations

10. 4.1 Dam-break problem

11. 4.2 Shock wave

12. 5 Discussion

13. Acknowledgments

Introduction

In fluid mechanics, many phenomena can be described by hyperbolic partial differential equations, such as the inviscid Burgers equation [6], the isentropic Euler equations [26] and the shallow water (Airy or Saint-Venant [11]) equations. The latter, for flat seabeds in one horizontal dimension, are most often written as mass and momentum flux conservations

[TABLE]

where $u\ =\ u(x,\,t)$ is the depth-averaged horizontal velocity ($x$ the horizontal coordinate, $t$ the time), $h\ =\ d\ +\ \text{\texteta}\,(x,\,t)$ is the total water depth ($d$ the mean depth, \texteta the surface elevation from rest) and $g$ is the (downward) acceleration due to gravity (see the sketch in Figure 1). These equations describe nonlinear non-dispersive long surface gravity waves propagating in shallow water. They are equations of choice when one is interested in modelling large scale phenomena without resolving the details at the small scales, for instance, in tsunamis and tides modelling. It should be noted that equations (1.1) – (1.2) are mathematically identical to the isothermal Euler equations describing some compressible fluids [8]. Here, we focus on shallow water equations, but it is clear that our claims apply as well to the isothermal Euler and mathematically similar equations.

Hyperbolicity is a nice feature of equations (1.1) – (1.2) because they can be tackled analytically and numerically with powerful methods (e.g., characteristics, finite volumes, discontinuous Galerkin). A major inconvenient is that these equations admit non-unique weak solutions and entropy considerations have been proposed to ensure unicity [23]. In gas dynamics, these weak solutions correspond to shock waves and the loss of regularity can be problematic, in particular for computations using spectral methods (even if some spectral approaches have been developed for hyperbolic equations as well [17]). Several methods have then been introduced to regularise the equations and, in particular, to avoid the formation of sharp discontinuous shocks (replacing them by smooth $\tanh$-like profiles). Perhaps, the first regularisation was proposed by J. Leray [25] in the context of incompressible Navier–Stokes equations. His theoretical programme consisted in showing the existence of solutions in regularised equations, subsequently taking the limit $\upvarepsilon\ \to\ 0$ ($\upvarepsilon$ a small regularising parameter) to obtain weak solutions of Navier–Stokes.

A method of regularisation consists in first adding an artificial viscosity into the equations, and in taking the limit of vanishing viscosity in a second time. This method was introduced by von Neumann and Richtmyer [33]. It allows to generalise the classical concept of a solution and to prove eventually the uniqueness, existence and stability results for viscous regularised solutions [10, 4]. Due to the added viscosity, the energy is no longer conserved, that can be a serious drawback for some applications, for instance for long time simulations when the shocks represent (unresolved) small scale phenomena that are not dissipative.

Another regularisation consists in adding weak dispersive effects to the equations [22]. As shown by Lax [24], the dispersive regularisation is not always sufficient to obtain a reasonable limit to weak entropy solutions as the dispersion vanishes. Consequently, the most successful approach to study non-classical shock waves is to consider the combined dispersive-diffusive approximations [19]. Also, the added dispersion can generate high-frequency oscillations that must be resolved by the numerical scheme, resulting in a significant increase in the computational time. Nonlinear diffusive–dispersive regularisations for the scalar case were considered in [1]. The main goal was to obtain a regularised model which admits the existence of classical solutions globally in time.

Yet another, less known, regularisation inspired by Leray’s method [25], consists filtered the velocity field such that the resulting equations are non-dissipative and non-dispersive. Such regularisations have been proposed for the Burgers [2], for isentropic Euler [3, 28, 29] and other [7] equations. In the literature, this regularisation method appears with various denominations, such as Leray*-type regularisation*, $\upalpha-$regularisation and Helmholtz* regularisation*. A drawback of this method is that the regularised (then smooth) shocks propagate at a speed different than that of the original equations. This drawback, among other things, is addressed in the present paper.

In this paper, we propose a new type of regularisation which is both non-dissipative and non-dispersive. This regularisation preserves the conservations of mass, momentum and energy, that is an important feature for some physical applications, in particular for long time simulations. The derivation proposed below is based on a recent work [9] where a Lagrangian, suitable for long waves propagating in shallow water, was modified to incorporate one free parameter that can be used to improve the dispersion properties. Here, we make another step introducing two independent parameters. But, instead of improving the dispersion, these parameters are chosen to cancel the dispersion and thus to provide a regularisation of the classical shallow water (Saint-Venant) equations. In addition to being conservative and non-dispersive, this model yields regularised (smooth) ‘shocks’ that propagate exactly at the same speed as in the original model. The properties of the obtained model are discussed below, mostly via numerical evidences.

The present manuscript is organised as follows. In Section 2, a new two-parameter generalised shallow water model is introduced. In Section 3, the two parameters are related in a way to provide a non-dispersive non-diffusive and conservative regularised shallow water equations. The jump conditions of Rankine–Hugoniot type on both sides of a shock wave are discussed in Section 3.4. These equations admit regular travelling wave solutions, as shown in Section 3.3. In Section 4, several numerical examples are provided, demonstrating the efficiency of the method. Finally, some conclusions and perspectives are outlined in Section 5.

Model equations

For two-dimensional surface gravity waves propagating in shallow water over a horizontal seabed, the shallow water equations (1.1) – (1.2) can be derived from the Lagrangian density

[TABLE]

where $\phi\,(x,\,t)$ is, physically, a velocity “potential” introduced here as a Lagrange multiplier. Taking into account the next order of approximation, one can derive the very-well known fully-nonlinear weakly dispersive classical Serre–Green–Naghdi (cSGN) [18, 30, 32, 34, 35]. These equations can be derived in many ways, but the simplest derivation is from the Euler–Lagrange equations of the Lagrangian density

[TABLE]

Recently, a modified version of these equations was proposed in [9] in order to improve the dispersion properties, if needed. These improved Serre–Green–Naghdi (iSGN) can be derived from the Lagrangian density

[TABLE]

where \textbeta is a free parameter at our disposal. The Lagrangian densities $\mathscr{L}_{1}$ and $\mathscr{L}_{2}$ have the same order of approximation (see [9] for details) and, obviously, $\mathscr{L}_{1}$ is a special case of $\mathscr{L}_{2}$ when $\text{\textbeta}\ =\ 0\,$. The case $\text{\textbeta}\ =\ 2/15$ is the best choice for improving the dispersion properties of infinitesimal waves. The reader is referred to [9] for further details about the cSGN and iSGN models.

In the present paper, we consider a natural two-parameter extension of $\mathscr{L}_{2}$ in the form

[TABLE]

This Lagrangian density obviously generalises the ones above. The corresponding Euler–Lagrange equations are

[TABLE]

thence one obtains the generalised Serre–Green–Naghdi (gSGN) equations

[TABLE]

From these equations, a non-conservative momentum equation can be obtained as

[TABLE]

where

[TABLE]

The quantity $\Gamma$ plays the role of a relaxed version of fluid vertical acceleration at the free surface. Conservative equations for the momentum flux and the energy can also be obtained as

[TABLE]

Thus, for any choice of the parameters $\text{\textbeta}_{j}\,$, the gSGN equations conserve mass, momentum and energy.

It should be noted that $\mathscr{L}_{3}$ is asymptotically consistent with $\mathscr{L}_{1}$ and $\mathscr{L}_{2}$ only if $\text{\textbeta}_{1}\ =\ \text{\textbeta}_{2}\,$. When $\text{\textbeta}_{1}\ \neq\ \text{\textbeta}_{2}\,$, $\mathscr{L}_{3}$ remains asymptotically consistent with $\mathscr{L}_{0}$ for all choices of the parameters $\text{\textbeta}_{j}\,$, however. Thus, in the next section, we exploit this feature to derive a regularised version of the Saint-Venant equations.

Regularised shallow water equations

Here, we consider non-dispersive version of the gSGN equations above, obtaining thus a conservative regularised modification of the classical (i.e. dispersionless) shallow water (Saint-Venant) equations.

Linear dispersion relation

For infinitesimal waves, \texteta and $u$ being small, the gSGN equations can be linearised as

[TABLE]

Seeking for travelling wave solutions of the form $\text{\texteta}\ =\ a\,\cos k\,(x\ -\ c_{0}\,t)\,$, one obtains the linear dispersion relation

[TABLE]

A non dispersive model (i.e. $c_{0}$ independent of $k$) is obtained taking $\text{\textbeta}_{1}\ =\ \text{\textbeta}_{2}\ -\ 2/3\,$. It is this special choice that is investigated in this paper and that we call the regularised Saint-Venant (rSV) equations.

Regularised Saint-Venant (rSV) equations

Choosing $\text{\textbeta}_{1}\ \mathop{\stackrel{{\scriptstyle\mathrm{def}}}{{=}}}\ 2\/\upvarepsilon\ -\ 2/3$ and $\text{\textbeta}_{2}\ \mathop{\stackrel{{\scriptstyle\mathrm{def}}}{{=}}}\ 2\/\upvarepsilon$ (in order to obtain a non-dispersive model) and introducing $\Gamma\ \mathop{\stackrel{{\scriptstyle\mathrm{def}}}{{=}}}\ 3\/\upvarepsilon\/\mathscr{R}$ for convenience ($\upvarepsilon$ being a free parameter), the Lagrangian density $\mathscr{L}_{3}$ becomes

[TABLE]

and the resulting equations are

[TABLE]

If $\upvarepsilon\ =\ 0\,$, the classical Saint-Venant (cSV), or nonlinear shallow water equations (NSWE), are recovered. For $\upvarepsilon\ \neq\ 0\,$, $\mathscr{R}$ is a regularising term that prevents the formation of shocks, as shown below in Section 4. Of course, the rSV equations yield a conservative equation for the energy

[TABLE]

It should be noted that, though $\mathscr{R}$ involves high-order derivatives, the resulting rSV equations are not dissipative and not dispersive. For numerical resolutions, (3.4) can be advantageously replaced by one of the two equations

[TABLE]

where the new variable $\mathcal{U}\/\mathop{\stackrel{{\scriptstyle\mathrm{def}}}{{=}}}\/u\/-\/\upvarepsilon\/h\/(\/3\/h_{x}\/u_{x}\/+\/h\/u_{xx}\/)$ can be physically interpreted as an approximation of the tangential velocity at the free surface (see [9, Appendix B]). Note that $h\mathcal{U}\ =\ hu\ -\ \upvarepsilon\,\partial_{x}[h^{3}\,u_{x}]$ and that rSV is different from the regularised isentropic Euler equations proposed in [3] where both the mass and momentum equations are modified.

Steady motion

Consider now the special case of travelling waves with permanent form, studied in the frame of reference moving with the wave where the flow is steady. The functions $h$ and $u$ being then independent of the time $t\,$, the mass conservation can be integrated as

[TABLE]

so $c$ is the wave phase velocity observed in the frame of reference without mean flow. In the latter frame of reference, the wave travels toward the increasing $x-$direction if $c\ >\ 0\,$. With (3.9) and after some algebra, equations (3.7) and (3.8) give, respectively,

[TABLE]

where $\mathcal{C}_{j}$ are dimensionless integration constants ($\mathcal{C}_{1}\ =\ \mathcal{C}_{2}\ =\ 0$ if $h\ \to\ d$ as $x\ \to\ \infty$). Eliminating $h_{xx}$ between these two relations, one obtains the first-order ordinary differential equation

[TABLE]

where $\mathcal{F}\ \mathop{\stackrel{{\scriptstyle\mathrm{def}}}{{=}}}\ c^{2}/gd$ is a Froude number squared. If $\upvarepsilon\ =\ 0\,$, the equation (3.12) does not admit physically admissible regular smooth solutions.

If $\upvarepsilon\ \neq\ 0\,$, exact solutions can be easily obtained in parametric form $x\ =\ x(\xi)\,$, $h\ =\ h(\xi)\ =\ d\ +\ \text{\texteta}\,(\xi)$ and in terms of Jacobian elliptic functions. For brevity, we give here only the solitary wave (i.e. $\mathcal{C}_{1}\ =\ \mathcal{C}_{2}\ =\ 0$) solution

[TABLE]

This solution is admissible only if $\upvarepsilon$ is positive. It corresponds to a dispersionless solitary wave, as can be seen in the independence of the trend parameter $\kappa$ with respect of the amplitude. This is therefore not a suitable model for solitary surface gravity waves. However, is some media, there exist dispersionless solitary waves [27] and the rSV equations could be used as model (with, of course, different physical interpretations of the variables $h$ ans $u$ and of the parameters $g$ and $d$).

Note that the solitary wave solution shows that the phase velocity varies with the wave amplitude. The rSV waves (like the cSV ones) have thus amplitude dispersion though they have no frequency dispersion. In this paper dispersion always refers to frequency dispersion.

Rankine–Hugoniot type conditions

In the theory and practice of hyperbolic equations, it is well known that equations in the velocity $u$ or $\mathcal{U}$ such as (3.8) with $\upvarepsilon\ \to\ 0$ are not suitable for discontinuous solutions since they yield physically incorrect Rankine–Hugoniot jump conditions [15, 16]. However, this is not a problem if $\upvarepsilon\ \neq\ 0$ because no (discontinuous) shocks are formed, as illustrated numerically below (see Section 4).

Assuming that $h_{x}$ and $u_{x}$ are both continuous if $\upvarepsilon\ >\ 0$ and that discontinuities (if any) occur only in $h_{xx}$ and $u_{xx}$ (and thus in $\mathcal{U}$ too). Equations (3.7) and (3.8) (together with $\llbracket\/\mathcal{U}\/\rrbracket\ =\ -\upvarepsilon\/h^{2}\llbracket\/u_{xx}\/\rrbracket$) yield the same jump condition

[TABLE]

where $\dot{s}\ \mathop{\stackrel{{\scriptstyle\mathrm{def}}}{{=}}}\ \mathrm{d}s/\mathrm{d}t$ is the speed of the smoothed shock located at $x\ =\ s\,(t)$ and $\llbracket\/f\/\rrbracket\mathop{\stackrel{{\scriptstyle\mathrm{def}}}{{=}}}f(x\!=\!s^{+})-f(x\!=\!s^{-})$ denotes the jump across the shock for any function $f\,$. For brevity, we call ‘shock’ both the discontinuous (classical) and smoothed (regularised) shocks.

Differentiating twice with respect of $x$ the mass conservation (3.3), the jump condition of the resulting equation is

[TABLE]

and the elimination of $\llbracket u_{xx}\rrbracket$ (or $\llbracket h_{xx}\rrbracket$) between (3.14) and (3.15) yields at once

[TABLE]

The shock speed is thus independent of $\upvarepsilon$ and the shock propagates along the characteristic lines of the classical Saint-Venant equations.

Under the same regularity conditions, the jump condition for the momentum and energy equations, respectively (3.4) and (3.6), both yield $\llbracket\/\mathscr{R}\/\rrbracket\ =\ 0$ (i.e., $\mathscr{R}$ is continuous), thus from (3.5)

[TABLE]

thence, using (3.14), one gets

[TABLE]

The continuity of $\mathscr{R}$ is compatible with its definition (3.5) and with (3.14) and (3.17), i.e.

[TABLE]

$\mathscr{R}$ being continuous across the shock, applying one spatial derivative $\partial_{x}$ to the momentum equation (3.4), the jump condition for the resulting equation is $\llbracket\/\mathscr{R}_{x}\/\rrbracket=0\,$, so $\mathscr{R}_{x}$ too is continuous across the shock. Applying the same procedure to the energy equation (3.6), one gets the jump condition

[TABLE]

which, with (3.14) and (3.15), is identically fulfilled. Thus, the derivative of the energy equation does not provide any additional information. It should be noted that all the jump relations above are independent of the regularising parameter $\upvarepsilon\,$. Moreover, the regularising term vanishes identically on constant states; this property is necessary to preserve shock conditions of the original hyperbolic system.

Remarks on the total energy

For a domain $\Omega\,$, the energy density $\mathscr{H}_{\epsilon}$ and total energy $\mathcal{H}_{\epsilon}$ of the rSV equation are

[TABLE]

If $\Omega$ is periodic or if the flux of energy is constant at the boundaries of $\Omega\,$, then $\mathcal{H}_{\epsilon}$ is constant (i.e., $\mathrm{d}\mathcal{H}_{\epsilon}/\mathrm{d}\/t\ =\ 0$) because the energy equation (3.6) is conservative. However, the quantity $\mathcal{H}_{0}\,(t)$ (corresponding to the energy of the cSV equations) is not constant, in general.

If $\mathscr{H}_{0}$ represents the density of physical energy (kinetic plus potential) then $\mathscr{H}_{\epsilon}\ -\ \mathscr{H}_{0}$ can be interpreted as a density of ‘internal’ energy, $\mathcal{H}_{\epsilon}\ -\ \mathcal{H}_{0}$ being the ‘internal’ energy of the domain $\Omega\,$. This ‘internal’ energy can also be interpreted as an ‘entropy’. Indeed, when the temporal evolution of a smooth initial condition leads to stiffening of the free surface (see numerical simulations below), the quantities $\lvert\,h_{x}\,\rvert$ and $\lvert\,u_{x}\,\rvert$ increase in time in the vicinity of the forming shock, and so does $\mathscr{H}_{\epsilon}\ -\ \mathscr{H}_{0}\,$. We have then a transfert of energy between the physical one to the ‘internal’ one. This behaviour is consistent with the cSV equations where the energy decreases across shocks.

Numerical illustrations

In order to study the properties of the proposed regularisation method, we solve numerically the equations (3.3) – (3.5) using a Fourier-type pseudo-spectral method. It is totally fine since for any $\upvarepsilon\ >\ 0$ we are dealing with smooth solutions. The periodicity is enforced by symmetrising the solution. The time stepping is done with automatic time step selection using embedded Runge–Kutta methods [31]. The relative and absolute errors were set to $10^{-10}\,$. Our pseudo-spectral approach to non-hydrostatic equations is described in [13]. As anti-aliasing we use an eight-order Erfc–Log filter [5]. In order to compare the regularised solutions with entropic solutions to the original hyperbolic system, we employ the finite volume scheme described in some detail in [12]. Thus, we do not reproduce the numerical details in the present study. All numerical results are performed in dimensionless variables where $g\ =\ d\ =\ 1\,$.

Dam-break problem

The dam-break problem has become the standard test-case for shallow water equations [14]. Here, we consider the so-called wet dam-break problem where the water depth is positive everywhere and a wave system is generated due to the initial difference of water level on the left and on the right from a certain point, which is chosen to be the origin without loss of generality. The initial condition is chosen as a regularised wave front:

[TABLE]

The initial condition parameters are given in Table 1 (central column) and it is depicted in Figure 2a. The temporal evolution of this initial condition is shown in Figure 2b–d. In particular, one can see that during the propagation, the initially smooth front becomes steeper to produce a real shock wave propagating rightwards. On the left, we observe a rarefaction wave according to the classical Riemann problem solution. The most important point to notice is that we have no dissipation and no dispersion in the regularised solution according to our model construction. For the value of the regularisation parameter $\upvarepsilon\ =\ 10^{-3}$ the curves are non-distinguishable to graphical resolution.

An important feature of the rSV equations is that the smoothed shock speed is independent of the regularisation parameter $\upvarepsilon\,$. This can be clearly seen doing the same simulations as in Figure 2, but with $\upvarepsilon\ =\ 1$ (Fig. 3) and with $\upvarepsilon\ =\ 5$ (Fig. 4). In these two simulations, the discrepancies between the regularised and original shallow water equations are of course more important than with small $\upvarepsilon\ =\ 10^{-3}\,$, but it is clear that the shock speeds are not affected by the regularisation.

Shock wave

The dam-break problem considered above can be slightly refined by choosing thoroughly water levels on both sides of the dam. If it is done according to the classical Rankine–Hugoniot conditions, only one isolated shock wave will emerge to move rightwards with a constant speed, which is also given by the same conditions. Such initial condition parameters are given in Table 1 (rightmost column). The shape of the initial condition is the same as above (for the free surface elevation, see Figure 2a). The snapshot of the free surface profile at $t\/\sqrt{g/d}\ =\ 5$ is depicted in Figure 5. For the value $\upvarepsilon\ =\ 10^{-2}$ of the regularisation parameter, we can barely distinguish the two profiles. The most important is that the shock wave position is the same in both models, according to the theoretical predictions (see Section 3.4). We underline also the fact that there are no oscillations around the jump, showing again that the proposed regularisation is non-dispersive even in fully nonlinear simulations.

As the wave is stiffening, the ‘physical’ energy $\mathcal{H}_{0}$ decreases while the ‘internal’ energy $\mathcal{H}_{\epsilon}\ -\ \mathcal{H}_{0}$ increases (Fig. 6). The slope being limited in the rSV equations, these variations are plateauing as the steady regularised shock state is reached. This behaviour is expected as it is consistent with the energy jump in the cSV equations. This indicates that the solutions of the rSV equations should tend to the ones of the cSV equations as $\epsilon\ \to\ 0^{+}\,$. This claim is illustrated numerically but it remains to be proven rigorously.

Discussion

This paper presents several new developments for the non-dispersive shallow water equations. First, a new Lagrangian for long wave propagating in shallow water is proposed. This Lagrangian contains two free parameters that can be chosen independently. In the present study, our goal was to remove the frequency dispersion effects in order to obtain a generalised version of the celebrated hyperbolic Saint-Venant equations [11]. Tough the dispersion relation analysis is only linear, the analysis of steady flows shows that the non-dispersive feature remains in nonlinear solutions. In order to investigate the non-dispersive features for unsteady flows, we performed numerical studies of the fully nonlinear equations. The numerical results confirm the absence of dispersive effects in the unsteady nonlinear equations. As a result, we obtained a non-hydrostatic non-dispersive model for nonlinear shallow water waves. Second, for any value of the regularisation parameter $\upvarepsilon\ >\ 0\,$, we obtain smooth and monotonic solutions to the classical dam-break problem [14]. The regularising effect is not limited to this peculiar problem. For instance, the hyperbolic shock waves in Saint-Venant equations are replaced by smoothed kink-like fronts [20] without introducing any dissipative effects into the model. We remind that the rSV model is fully conservative and, by construction, it inherits a variational structure as well. In addition, the main advantage of the method proposed in the present study is that: (i) this regularisation can be supplied with a clear physical meaning; (ii) the whole variational structure of original equations is conserved; (iii) the regularised shock speed is independent of the regularising parameter and is identical to the original speed. To our knowledge, it is the first regularisation which seems to achieve these goals.

The numerical simulations of the regularised equations where performed with a pseudo-spectral scheme. We chose this method for its speed and accuracy but, more importantly, because it requires great regularity of the solution. The fact that it worked fine here demonstrates that other numerical methods should a fortiori work for the regularised equations.

An efficient regularisation is obtained only for a well-chosen regularising function $\mathscr{R}\,$. Leray-like regularisations are obtained introducing an ad hoc linear correction. In the present context, this means that we should have used $\mathscr{R}\ =\ -d\,u_{\,xt}\ -\ g\,d\,h_{\,xx}$ into the momentum equation, leading to the regularised velocity $\mathcal{U}\/=\/u\/-\/\epsilon\/d^{2}\/u_{xx}\,$. However, doing so, the speed of the regularised shocks is modified, so this regularisation is not optimal in that respect. Thus, the nonlinear terms in the definition of $\mathscr{R}$ play an important role and it is not at all trivial a priori to find out a regularising term $\mathscr{R}$ with all the desired properties. The complexity of the $\mathscr{R}$ used here shows that such term would have been very unlikely discovered by trail and error. Thanks to the relaxed variational formulation described in Section 2, a suitable choice for $\mathscr{R}$ appeared naturally. This is an illustration of the power of this variational approach.

The proposed rSV equations have to be studied deeper from the mathematical point of view. In particular, the limit of solutions when $\upvarepsilon\ \to\ 0^{+}$ has to be rigorously established along the lines of [24, 21, 22], to give a few references. Moreover, at the limit, one has to recover not only a weak solution to the original hyperbolic system, but an entropy weak solution. More importantly, a question is whether the proposed regularisation allows to obtain existence and uniqueness (or maybe even stability) results for the limiting hyperbolic system. We hope that this work will attract the mathematical community’s attention to these opportunities.

Moreover, we hope that this approach will be generalised to other important examples of conservation laws which arise in applications [15]. It is clear that the success of this operation will greatly depend on the existence of underlying variational structure of equations.

Acknowledgments

The authors would like to acknowledge the support of CNRS under the PEPS InPhyNiTi 2015 project FARA. D. Dutykh would like to acknowledge the hospitality of the Laboratory J. A. Dieudonné during his visits in 2016.