Hyperbolic Model for Free Surface Shallow Water Flows with Effects of Dispersion, Vorticity and Topography

TL;DR

This paper introduces a hyperbolic system to model dispersive shallow water flows, simplifying numerical computations and accurately capturing wave phenomena like undular bores and solitary waves.

Contribution

A novel hyperbolic model for dispersive shallow water flows that avoids elliptic problems and reduces computation time, validated against the original dispersive model.

Findings

Model closely matches dispersive solutions over time

Effectively simulates undular bores and wave breaking

Reduces numerical complexity and computation time

Abstract

We derive a hyperbolic system of equations approximating the two-layer dispersive shallow water model for shear flows recently proposed by Gavrilyuk, Liapidevskii \& Chesnokov (J. Fluid Mech., vol. 808, 2016, pp. 441--468). The use of this system for modelling the evolution of surface waves makes it possible to avoid the major numerical challenges in solving dispersive shallow water equations, which are connected with the resolution of an elliptic problem at each time instant and realization of non-reflecting conditions at the boundary of the calculation domain. It also allows one to reduce the computation time. The velocities of the characteristics of the obtained model are determined and the linear analysis is performed. Stationary solutions of the model are constructed and studied. Numerical solutions of the hyperbolic system are compared with solutions of the original dispersive model. It is shown that they almost coincide for large time intervals. The system obtained is applied to study non-stationary undular bores produced after interaction of a uniform flow with an immobile wall, non-hydrostatic shear flows over a local obstacle and the evolution of breaking solitary wave on a sloping beach.

Tables

Table 1. Table 1: Computational time T α subscript 𝑇 𝛼 T_{\alpha} and number of iterations (time steps) M α subscript 𝑀 𝛼 M_{\alpha} of model ( 12 ) for different values of the relaxation parameter α 𝛼 \alpha .

	$\alpha =0$	$\alpha =3$	$\alpha =6$	$\alpha =9$	$\alpha =12$
$T_{\alpha }$	5.8660	8.0554	10.5936	12.8388	14.2673
$M_{\alpha }$	2560	3483	4566	5446	6183

Equations

$$\begin{array}[]{l}\displaystyle u_{t}+uu_{x}+wu_{z}+\rho^{-1}p_{x}=0,\\[5.69054pt] \varepsilon^{2}\big{(}w_{t}+uw_{x}+ww_{z}\big{)}+\rho^{-1}p_{z}=-g,\quad u_{x}+w_{z}=0.\end{array}$$

$$E_{t}+\big{(}(E+\rho^{-1}p)u\big{)}_{x}+\big{(}(E+\rho^{-1}p)w\big{)}_{z}=0,$$

$$b_{t}+ub_{x}-w\big{|}_{z=b}=0,\quad Z_{t}+uZ_{x}-w\big{|}_{z=Z}=0,\quad p\big{|}_{z=Z}=0.$$

$$(b+h)_{t}+u(b+h)_{x}-w\big{|}_{z=b+h}=-M,$$

$$U=\frac{1}{h}{\int }_b^{b+h}udz,\bar{u}=\frac{1}{\eta }{\int }_{b+h}^{Z}udz,q^2=\frac{1}{\eta }{\int }_{b+h}^Z(u-\bar{u}{)}^{2}dz.$$

$$\begin{array}[]{l}\displaystyle h_{t}+(Uh)_{x}=-M,\quad(Uh)_{t}+(U^{2}h+P)_{x}=-MU+g\eta h_{x}-ghb_{x},\\[11.38109pt] \displaystyle\eta_{t}+(\bar{u}\eta)_{x}=M,\quad(\bar{u}\eta)_{t}+\bigg{(}(\bar{u}^{2}+q^{2})\eta+\frac{g\eta^{2}}{2}\bigg{)}_{x}=MU-g\eta h_{x}-g\eta b_{x},\\[11.38109pt] \displaystyle\bigg{(}\bigg{(}\frac{\bar{u}^{2}+q^{2}+g\eta}{2}+g(b+h)\bigg{)}\eta\bigg{)}_{t}+\bigg{(}\bigg{(}\frac{\bar{u}^{2}+3q^{2}}{2}+g(b+h+\eta)\bigg{)}\bar{u}\eta\bigg{)}_{x}=\\[11.38109pt] \displaystyle\quad\quad\quad=\bigg{(}\frac{U^{2}}{2}+g(b+h+\eta)\bigg{)}M+g\eta(b+h)_{t}-\frac{\sigma\kappa q^{3}}{2}.\end{array}$$

$$P=gh\eta +\frac{g{h}^2}{2}-\frac{{\epsilon }^2{h}^3}{3}({\dot{U}}_x-U_x^2),M=\sigma q.$$

$$\begin{array}[]{l}\displaystyle\bigg{(}\bigg{(}\frac{U^{2}+gh}{2}+gb+\frac{\varepsilon^{2}}{6}h^{2}U_{x}^{2}\bigg{)}h\bigg{)}_{t}+\bigg{(}\bigg{(}\frac{U^{2}+gh}{2}+gb+\frac{\varepsilon^{2}}{6}h^{2}U_{x}^{2}\bigg{)}Uh+UP\bigg{)}_{x}=\\[11.38109pt] \displaystyle\quad\quad\quad=-\bigg{(}\frac{U^{2}}{2}+g(b+h+\eta)+\frac{\varepsilon^{2}}{2}h^{2}U_{x}^{2}\bigg{)}M-g\eta h_{t}+ghb_{t}.\end{array}$$

$$q_{t}+(\bar{u}q)_{x}=\frac{\sigma}{2\eta}\big{(}(U-\bar{u})^{2}-(1+\kappa)q^{2}\big{)}$$

$$K=U-\frac{\varepsilon^{2}}{3h}\big{(}h^{3}U_{x}\big{)}_{x}$$

$$K_{t}+\bigg{(}KU+g(h+\eta)-\frac{U^{2}}{2}-\varepsilon^{2}\bigg{(}\frac{U_{x}^{2}h^{2}}{2}+\sigma qhU_{x}\bigg{)}\bigg{)}_{x}=-\frac{\varepsilon^{2}}{3}\sigma qhU_{xx}-gb_{x}.$$

$$\begin{array}[]{l}\displaystyle\rho_{t}+(u\rho)_{x}=0,\quad(u\rho)_{t}+(u^{2}\rho+P)_{x}=0,\\[8.53581pt] \displaystyle P=P\bigg{(}\rho,\frac{d\rho}{dt},\frac{d^{2}\rho}{dt^{2}}\bigg{)},\quad\frac{d}{dt}=\frac{\partial}{\partial t}+u\frac{\partial}{\partial x}\end{array}$$

$$P=P\bigg{(}\rho,\frac{d\tilde{\rho}}{dt},\frac{d^{2}\tilde{\rho}}{dt^{2}}\bigg{)}.$$

$$\tilde{\rho }(s)=\tilde{\rho }(t)+{\tilde{\rho }}^{\prime }(t)(s-t)+{\tilde{\rho }}^{\prime \prime }(t)(s-t{)}^2/2+o({\tau }^2),s\in (t-\tau ,t+\tau ).$$

$$\rho (t)=\frac{1}{2\tau }\underset{t-\tau }{\overset{t+\tau }{\int }}\tilde{\rho }(s)ds=\tilde{\rho }(t)+\frac{1}{6}{\tilde{\rho }}^{\prime \prime }(t){\tau }^2+o({\tau }^2)$$

$${\tilde{\rho }}^{\prime \prime }(t)=\alpha (\rho (t)-\tilde{\rho }(t))+O(\tau ),\alpha =6/{\tau }^2.$$

$$\begin{array}[]{l}\displaystyle\rho_{t}+(u\rho)_{x}=0,\quad(u\rho)_{t}+(u^{2}\rho+\bar{P})_{x}=0,\\[5.69054pt] \displaystyle\tilde{\rho}_{t}+u\tilde{\rho}_{x}=v,\quad v_{t}+uv_{x}=\alpha(\rho-\tilde{\rho}),\end{array}$$

$$U_x=-\frac{\dot{h}+\sigma q}{h},{\dot{U}}_x=-\frac{\ddot{h}+\sigma \dot{q}}{h}+\frac{(\dot{h}+\sigma q)\dot{h}}{h^2}$$

$$-\frac{\varepsilon^{2}h^{3}}{3}\big{(}\dot{U}_{x}-U_{x}^{2}\big{)}=\frac{\varepsilon^{2}h^{2}}{3}\ddot{h}+O(\sigma\varepsilon^{2}).$$

$$(Uh)_{t}+\bigg{(}U^{2}h+\frac{gh^{2}}{2}+gh\eta+\frac{\varepsilon^{2}h^{2}}{3}\ddot{h}\bigg{)}_{x}=g\eta h_{x}-\sigma qU-ghb_{x}.$$

$$\begin{array}[]{l}\displaystyle(h+\eta)_{t}+(Uh+\bar{u}\eta)_{x}=0,\quad\eta_{t}+(\bar{u}\eta)_{x}=\sigma q,\\[11.38109pt] \displaystyle(Uh)_{t}+\bigg{(}U^{2}h+\frac{gh^{2}}{2}+gh\eta+\frac{\alpha\varepsilon^{2}}{3}(h-\zeta)h^{2}\bigg{)}_{x}=g\eta h_{x}-\sigma qU-ghb_{x},\\[11.38109pt] \displaystyle(\bar{u}\eta)_{t}+\bigg{(}(\bar{u}^{2}+q^{2})\eta+\frac{g\eta^{2}}{2}\bigg{)}_{x}=-g\eta h_{x}+\sigma qU-g\eta b_{x},\\[11.38109pt] \displaystyle q_{t}+(\bar{u}q)_{x}=\frac{\sigma}{2\eta}\big{(}(U-\bar{u})^{2}-(1+\kappa)q^{2}\big{)},\\[11.38109pt] \displaystyle(h\zeta)_{t}+(Uh\zeta)_{x}=Vh-\sigma q\zeta,\quad(Vh)_{t}+(UVh)_{x}=\alpha(h-\zeta)h-\sigma qV.\end{array}$$

$$U_{t}+\bigg{(}\frac{U^{2}}{2}+gH\bigg{)}_{x}=-gb_{x},\quad Q_{t}+\bigg{(}U^{2}h+(\bar{u}^{2}+q^{2})\eta+\frac{gH^{2}}{2}\bigg{)}_{x}=-gHb_{x},$$

$$\begin{array}[]{l}\displaystyle\zeta_{t}+U\zeta_{x}=V,\quad V_{t}+UV_{x}=\alpha(h-\zeta),\\[11.38109pt] \displaystyle s_{t}+\bar{u}s_{x}=\frac{\sigma}{2\eta^{2}}\big{(}(U-\bar{u})^{2}-(3+\kappa)q^{2}\big{)}.\end{array}$$

$$\begin{array}[]{l}\displaystyle U_{t}+UU_{x}-\frac{\alpha h}{3}\zeta_{x}+\bigg{(}g+\alpha\bigg{(}h-\frac{2\zeta}{3}\bigg{)}\bigg{)}h_{x}+g\eta_{x}=-gb_{x},\\[11.38109pt] \displaystyle\bar{u}_{t}+\bar{u}\bar{u}_{x}+2\eta^{2}ss_{x}+gh_{x}+(g+3s^{2}\eta)\eta_{x}=\sigma s(U-\bar{u})-gb_{x}.\end{array}$$

$${\mathbf{U}}_t+\mathbf{A}{\mathbf{U}}_x=\mathbf{F}$$

$$\chi (\lambda )=(U-\lambda {)}^2(\bar{u}-\lambda )\hat{\chi }(\lambda )=0,$$

$$\begin{array}[]{l}\displaystyle\hat{\chi}(\lambda)=\big{(}(U-\lambda)^{2}-a_{1}\big{)}\big{(}(\bar{u}-\lambda)^{2}-a_{2}\big{)}-g^{2}h\eta,\\[8.53581pt] \displaystyle a_{1}=gh+\alpha h\bigg{(}h-\frac{2\zeta}{3}\bigg{)},\quad a_{2}=g\eta+3q^{2}.\end{array}$$

$${\lambda }_{1,2}=U\pm \sqrt{a_1},{\lambda }_{3,4}=\bar{u}\pm \sqrt{a_2}.$$

$$\Phi =(\bar{u}-\lambda )/\sqrt{a_2},\Psi =(U-\lambda )/\sqrt{a_1}.$$

$$({\Phi }^2-1)({\Psi }^2-1)=r,$$

Full text

††journal: Computers & Fluids

1 Introduction

The second-order approximation of the shallow water theory is commonly used in modelling of the interaction of nonlinear waves. There are different versions of the governing equations taking into account the influence of non-hydrostatic pressure distribution on the structure of nonlinear surface and internal waves in the long-wave approximation [14, 39, 12]. Alternative formulations of these dispersive models within the framework of hyperbolic systems of equations have been proposed in [24, 22, 8]. The effect of non-hydrostatic distributions of pressure is used in these models by using additional ‘internal’ variables in the equations. In [15, 26] hyperbolic techniques are applied to dispersive systems in a similar way. The main advantage of such hyperbolic approximation of the dispersive equations is essential simplification of the algorithms of numerical calculation and formulation of the boundary conditions. One example is the Serre–Green–Naghdi equations describing dispersive water waves [14, 32, 34]. In particular, the inversion of an elliptic operator is needed at each time step when the model is numerically solved [25, 20]. As a consequence, this drastically increases the calculation time. Another important numerical problem is how to impose artificial non-reflecting (transparent) conditions at the boundary of the calculation interval for dispersive equations [13, 1]. These conditions are crucial when one looks for waves passing through a bounded numerical domain. This is always an open problem for solving general dispersive equations. On the other hand, the description of full propagation of nonlinear waves from the shoaling zone to the surf zone has considerable practical importance for coastal water waves. However, dispersive models such as the Green–Naghdi equations or their hyperbolic approximation predict that no wave breaks, which is obviously incorrect. For this reason, it is needed the construction of more complex models of the theory of long waves taking into account both dispersion effects and the possibility of wave breaking.

Shallow water flows often exhibit turbulent structures such as coastal water waves in the surf zone and rollers of hydraulic jumps. The appearance of these structures is usually related to the mechanism of wave breaking. Several experimental and theoretical works have considered the flow of hydraulic jumps or spilling breakers and studied the mechanism of turbulence generation [33, 35, 27]. In particular, work [27] highlighted the importance of an intense shear layer which spreads downstream from the toe of the breaker below the spilling breaker and the recirculating region of the roller where the turbulent kinetic energy is the most intense. The development of the method for constructing depth-averaged shallow water equations for shear flows in hydrostatic approximation [35] made it possible to describe flows exhibiting vortex structures such as hydraulic jump rollers [31] and turbulent roll waves [30, 17]. Numerical modelling of multi-dimensional turbulent hydraulic jumps formed in convergent radial flow was recently performed and qualitatively compared to the experimental observations in [18]. However, these shallow water models for shear flows did not take into account dispersion effects.

The approach proposed in [2, 36] is based on coupling between the dispersive Green–Naghdi equations describing long waves far from the coast and the hyperbolic Saint-Venant equations describing wave breaking near the coast. The difficulty is to understand when we replace one model by the other. The search for a ‘switching criterion’ is not a well-defined problem even if several empirical criteria have been proposed in the literature. For instance, wave phase velocity becomes larger than the flow velocity or the wave slope attains the critical value. A unified model that is capable of describing both dispersion and breaking waves, as well as the transition between different flow regimes, was derived in [10]. This model is based on a two-layer long-wave approximation of the homogeneous Euler equations with a free surface evolving over a mild slope. The upper layer is turbulent and hydrostatic, while the lower one is almost potential and can be described by Green–Naghdi equations. The interface separating these two layers is considered as a discontinuity surface where turbulent mixing occurs. The model was validated, in particular, on the propagation of Favre waves [7] and shoaling of solitary waves [16]. Generalization of this model for flows of stratified fluid was proposed in [11]. Recently a one-layer dispersive model of shoaling and breaking waves with shear effects complemented by turbulent viscous dissipative terms was obtained and verified in [19]. A hierarchy of non-hydrostatic shallow water-type models approximating the incompressible Euler and Navier-Stokes systems with free surface was presented in [4, 9]. Suggested in [6] a two-layer dispersive system provides an efficient and accurate approach to model the propagation of waves near coastal areas and intermediate waters.

The aim of the present paper is to derive and study a hyperbolic approximation of the proposed in [10] two-layer dispersive model for shear shallow flows. The construction of such an approximation allows us to simplify and accelerate the numerical calculations of the propagation and breaking of nonlinear waves. In Section 2 we recall a two-layer long-wave model proposed in [10] that describes the flow of a homogeneous fluid with a free surface over mild slopes, considering the effects of dispersion and vorticity. We also recall the numerical method for solving these equations, used in [10] and based on integrating a hyperbolic system with an additional inversion of an elliptic operator at each time step [25]. In Section 3, we give other formulation of this model, which is suitable for applying the method of constructing a hyperbolic approximation of the considered dispersive equations. Following [24, 22], using additional ‘internal’ variables, we derive a hyperbolic approximation of the two-layer long-wave model. Further, we find the velocities of the characteristics of this first-order system and prove that all of them are real for a sufficiently large relaxation parameter. We also perform a linear analysis and obtain the dispersion relations for the original model and its hyperbolic approximation. In Section 4, we study stationary solutions and determine the asymptotic behaviour of the solution. Numerical examples for the propagation of Favre waves, undular and monotonic turbulent bores generated by obstacle as well as the breaking of a solitary wave are considered in Section 5. In particular, we show that the obtained hyperbolic system approximates solutions of the original dispersive model. Finally, we draw some conclusions.

2 A two-layer long-wave approximation of the homogeneous Euler equations

The Euler equations for two-dimensional flows can be written as

[TABLE]

Here $u$ and $w$ are the velocity components in the horizontal $Ox$ and vertical $Oz$ directions, $p$ is the fluid pressure, $\rho$ is the fluid density, $g$ is the gravity acceleration in vertical direction and $\varepsilon=H_{0}/L_{0}\ll 1$ is the dimensionless long wave parameter ($H_{0}$ and $L_{0}$ are the characteristic vertical and horizontal scales). Equations (1) admit the conservation of energy

[TABLE]

where $E=(u^{2}+\varepsilon^{2}w^{2})/2+gz$. The boundary conditions at the bottom $z=b(t,x)$ and at the free surface $z=Z(t,x)$ are

[TABLE]

Following [10], we assume that there is the internal boundary $z=b+h$ separating the lower layer of depth $h(t,x)$, where the flow is potential, and the upper turbulent layer. The kinematic condition at this interface is

[TABLE]

where the right-hand side $M$ is responsible for the mixing between layers. Let us introduce the depth-average velocity in the lower and upper layers and the specific turbulent energy $q^{2}$ for the upper layer with respect to the formulas

[TABLE]

Here $\eta(t,x)$ is the depth of the upper turbulent layer and $Z=b+h+\eta$.

We assume the pressure distribution in the lower layer is non-hydrostatic and the flow in this domain is almost potential. The upper layer is hydrostatic and the flow in this layer is turbulent. We use an additional variable $q\geq 0$ for averaged description of the shear flow in this layer. The interaction between the layers is taken into account through a kinematic condition (4) with the right-hand side $M=\sigma q$. We also assume that $u|_{z=b+h}=U$ (this condition provides compatibility between the averaged mass, momentum and energy equations for the lower layer). Under these assumptions the following depth-averaged model for a two-layer long-wave approximation of the homogeneous Euler equations with a free surface evolving over mild slope has been derived in [10]

[TABLE]

where

[TABLE]

Here and below ‘dot’ denotes the material time derivative $\dot{f}=f_{t}+Uf_{x}$. The non-negative constants $\sigma$ and $\kappa$ are the empirical parameters responsible for the mixing and energy dissipation. According to [38, 3, 24] these parameters are as follows $\sigma\approx 0.15$ and $\kappa\in[2,6]$. A mild slope approximation means that the dimensionless bottom variation is weak [32, 11]: $z=b(\varepsilon^{\gamma}t,\varepsilon^{\gamma}x)$, $\gamma>0$. Due to this fact the terms $\varepsilon^{2}\ddot{b}$ and $\varepsilon^{2}\dot{b}$ can be neglected during the derivation of system (5).

The first two equations in (5) correspond to the balance of mass and momentum in the non-hydrostatic and almost potential lower layer. They are obtained by averaging the incompressibility and horizontal momentum equations in (1) over the depth taking into account the boundary conditions (3) and (4). We also assume that $u|_{z=b+h}=U$, since this is the only case the first two equations in (5) are compatible with the energy equation for the lower layer (see [10] for details)

[TABLE]

The averaged mass and momentum balance equations for the upper hydrostatic layer are derived similarly. However, the flow in this layer is vortex and in the averaged description it is characterized by shear velocity $q$ defined above. In order to obtain the closing relation, we average the energy equation (2) by the upper layer depth. It was shown in [10] that it is convenient to use the differential consequence

[TABLE]

instead of the last balance law in (5).

The upper turbulent layer is considered within the hydrostatic approximation for shear shallow water flows [35], while the lower layer can be described by the Serre–Green–Naghdi dispersive equations [14, 32]. The interaction between these layers is taken into account through a natural mixing process. This model reveals the main mechanism of the spilling breaker development. For waves of moderate amplitude the upper thin layer is dynamically passive and the flow is governed by the dispersive model. However, for the waves of larger amplitude, the upper turbulent layer dynamics becomes crucial. The larger the wave amplitude, the greater the difference is in the average flow velocities of layers. It finally results in an intensive growth of the upper turbulent layer.

For the numerical solution of system (5) in [10] it was proposed to introduce a new variable $K$ as follows [25]

[TABLE]

and instead of the second equation in (5) to use its differential consequence

[TABLE]

Then the numerical resolution of the governing equation is divided into two successive steps: the time evolution of the conservative variables $\mathbf{u}=(h,K,\eta,\bar{u}\eta,q)$ using the method based on modification of the Godunov’s scheme; and the resolution of an ordinary differential equation (7) to obtain the values of velocity $U$ from variables $h$ and $K$.

3 Hyperbolic approximation of the two-layer dispersive equations for shear flows

At the beginning of this section, we recall the method for constructing the hyperbolic approximation of a dispersive system proposed in [24, 22]. Further, we apply this method to model (5) and derive a first-order evolutionary system that approximates this dispersive model. We also study the hyperbolicity of the obtained system and perform a standard Fourier analysis of linearised equations.

3.1 Method for constructing a hyperbolic approximation of a dispersive model

As in the example of the following one-dimensional dispersive system

[TABLE]

we briefly explain the method of constructing a hyperbolic approximation of (9) on the basis of the averaging of instantaneous variables [24, 22]. System (9) includes the equations of bubbly flow, models of the shallow water theory, etc. [12, 14, 39]. Such equations are used to describe flows with ‘internal inertia’, i.e., heterogeneous media where a certain part of the total energy transforms to the energy of small-scale motion. Numerical implementation of dispersive models involves some difficulties caused by non-hyperbolicity of the considered system.

A hyperbolic model is obtained by additional averaging of the equations and introducing new ‘internal’ variables. The scale of averaging is assumed to be rather small, which allows the values of the variables $\rho$ and $u$ to be replaced in the equations by their mean values. However, to calculate the function $P$, derivatives of the ‘instantaneous’ variable $\tilde{\rho}$ are used. This means that

[TABLE]

The relation between the averaged and ‘internal’ variables is given by expanding the functions $\tilde{\rho}(s)$ into a Taylor series along the trajectory $x=x(s,\xi)$

[TABLE]

Here $\xi$ is the fixed Lagrangian coordinate of the particle. It follows from the previous formula that the mean value $\rho(t)$ and instantaneous variable $\tilde{\rho}(t)$ are related by

[TABLE]

and, consequently,

[TABLE]

Using the main part in representation of function $\tilde{\rho}^{\prime\prime}(t)$, we can approximate system (9) as follows [22]

[TABLE]

where $\bar{P}=(\rho,\tilde{\rho},v)=P(\rho,v,\alpha(\rho-\tilde{\rho}))$. For a wide class of fluid flows, Eqs. (10) form a hyperbolic system with two sonic and two contact characteristics. As the parameter $\alpha$ increases (in this case the averaging interval $\tau$ decreases), the solutions of hyperbolic Eqs. (10) approximate the solutions of dispersive system (9).

3.2 Hyperbolic approximation of Eqs. (5)

Let us modify system (5) so that it is convenient to apply the method of constructing a hyperbolic approximation described above. Taking into account the first equation in (5) one can deduce that

[TABLE]

and, consequently, the dispersive term in the second equation of (5) reads

[TABLE]

We recall that the governing equations for the lower dispersive layer have been derived with accuracy $O(\varepsilon^{2})$. Since the mixing parameter $\sigma\approx 0.15$ is also small ($\sigma\sim\varepsilon$), in this approximation the terms of order $O(\sigma\varepsilon^{2})$ can be omitted. That is why with the same order of accuracy we can use the second equation of system (5) in the form

[TABLE]

This form of the momentum equation for the lower dispersive layer is more convenient for constructing a first-order approximative system.

From system (5) one can clearly see conservation of the total mass and momentum (for flows over a flat bottom). We note that up to the small terms of order $O(\sigma\varepsilon^{2})$ the total energy also conserves if the additional dissipation term vanishes ($\kappa=0$).

Following [22], we introduce the new instantaneous variables $\zeta$ and $V$ so that $\dot{\zeta}=V$, $\dot{V}=\alpha(h-\zeta)$ and replace the term $\ddot{h}$ in (11) by $\alpha(h-\zeta)$. In this case, the approximation of Eqs. (5) by a system of the first-order evolutionary equations is written as follows

[TABLE]

Here we represent equations for the variables $\zeta$ and $V$ in a conservative form (the last two equations in (12)) that is convenient for numerical treatment. There are non-conservative terms $\pm g\eta h_{x}$ on the right-hand side of momentum equations in (5) and (12). This is admissible because in the lower layer the flow is described by the Green–Naghdi equations (or their hyperbolic approximation) having smooth solutions.

System (12) contains parameter $\alpha$ depending on the chosen scale of averaging of the original model (5). The solutions of Eqs. (12) approximate the solutions of Eqs. (5) as $\alpha\to\infty$. They also yield the solutions of (5) in the hydrostatic case ($\varepsilon=0$) if $\alpha\to 0$. We note that model (5) in hydrostatic approximation was derived and studied in [21]. To describe such flows we eliminate the last two equations and rewrite the third and fourth equations of (12) in the conservative form

[TABLE]

where $H=h+\eta$ and $Q=Uh+\bar{u}\eta$ are the total depth and flow rate respectively.

3.3 Characteristics of Eqs. (12)

Let us find the characteristics of system (12). It follows from (12) that there are three contact characteristics $dx/dt=U$ (multiplicities two) and $dx/dt=\bar{u}$ since the variables $\zeta$, $V$ and $s=q/\eta$ satisfy the equations

[TABLE]

The momentum equations for the lower and upper layers in system (12) are equivalent to

[TABLE]

Here and below we assume that $\varepsilon=1$.

Using the previous equations, we represent system (12) in the form

[TABLE]

where $\mathbf{U}=(\zeta,V,s,h,\eta,U,\bar{u})^{\rm T}$ is the vector of unknown variables, $\mathbf{A}$ is the $7\times 7$ matrix, and $\mathbf{F}$ is the right-hand side. The eigenvalues of $\mathbf{A}(\mathbf{U})$ are determined by equation

[TABLE]

where

[TABLE]

Obviously, for a single-layer flow ($\eta=0$ or $h=0$), the system is hyperbolic and the sonic characteristics are

[TABLE]

Let us show that for a two-layer flow ($h>0$, $\eta>0$ and $h>2\zeta/3$) system (12) is hyperbolic for a sufficiently large value of $\alpha$. A geometric interpretation of the characteristics proposed in [29, 5] for two-layer hydrostatic flows can be applied here. We introduce the new variables $\Phi$ and $\Psi$ by the formulas

[TABLE]

Then equation $\hat{\chi}(\lambda)=0$ can be rewritten in the form

[TABLE]

where $r=(a_{1}a_{2})^{-1}g^{2}h\eta$. In the $(\Phi,\Psi)$–plane, Eq. (14) describes a fourth-order curve with four symmetry axes (see Fig. 1). The variables $\Phi$ and $\Psi$ by virtue of (13) are related by

[TABLE]

The number of real roots of $\hat{\chi}(\lambda)=0$ is determined by the number of intersections of the curve (14) with the straight line (15). Each point of interaction yields a sonic characteristic with the slope $\lambda=U-\Psi\sqrt{a_{1}}$. It follows from Fig. 1 that line (15) always has two points of intersection with curve (14) in the domains $\{\Phi>1,\Psi>1\}$ and $\{\Phi<-1,\Psi<-1\}$. We note that the necessary condition for the existence of four real roots is the inequality $\mu=\sqrt{1-r}>0$, which is always fulfilled. If we increase the parameter $\alpha$, the ‘radius’ $\mu$ becomes larger (it tends to unity for $\alpha\to\infty$) while an initial ordinate $(U-\bar{u})/\sqrt{a_{1}}$ of the straight line (15) tends to zero. As a result, for sufficiently large $\alpha$ there are two points of intersection of the line with curve (14) in the square $\{-1<\Phi<1,-1<\Psi<1\}$.

For all numerical examples considered below, Eqs. (12) are hyperbolic (characteristic equation $\chi(\lambda)=0$ has seven real roots). The presence of real characteristics makes it possible to apply the standard Godunov-type schemes for numerical solving system of balance laws (12). We note that the increase in parameter $\alpha$ also increases the absolute value of the maximum characteristic velocity. Time step in the Godunov-type schemes is inversely proportional to the characteristic velocity. Therefore, increasing of $\alpha$ slows down the computation.

3.4 Linear analysis

In this section, dispersion relations of the two-layer dispersive system (5) and its hyperbolic approximation (19) are obtained and studied. This analysis allows one to follow the influence of the relaxation parameter $\alpha$ on the accuracy of the hyperbolic approximation. Here we assume a flat bottom topography $b=0$ and take $\varepsilon=1$.

Consider the perturbation of a constant state $h=h_{0}$, $U_{0}$, $\eta=\eta_{0}$, $\bar{u}=U_{0}$, $q=0$ for dispersive system (5). Since the governing equations are Galilean invariant, without loss of generality, we can choose $U_{0}=0$. At first order the system reads:

[TABLE]

The search for a non-trivial solution of the linear system in the form

[TABLE]

leads to the following dispersion relation

[TABLE]

Here $k$ is the wave number and $c$ is the phase velocity. Note that for a single-layer non-hydrostatic flow ($\eta_{0}=0$) formula (17) transforms to the well-known dispersion relation of the Serre–Green–Naghdi equations.

Let us linearise hyperbolic model (12) on the same constant solution (additionally we suppose $\zeta=h_{0}$ and $V=0$). Then the linear system takes form (16) where the second equation reads

[TABLE]

and two more equations $\zeta_{t}=V$ and $V_{t}=\alpha(h-\zeta)$ should be added to the system. As before, we are looking for a solution in the form of monochromatic perturbations (with $c\neq 0$). In this case a non-trivial solution exists if the phase velocity $c$ satisfies the equation

[TABLE]

where $H_{0}=h_{0}+\eta_{0}$. The previous equation (at least for $\eta_{0}<h_{0}+(h_{0}^{2}/6+1/k^{2})\alpha/g$) has two real positive roots $(c^{\pm})^{2}$:

[TABLE]

It follows from (17) and (18) that the phase velocities $c^{\pm}$ and $c$ satisfy to inequalities

[TABLE]

and $c^{-}(\alpha)\to c$ as $\alpha\to\infty$ (see Fig. 2). We note that similar inequalities for the phase velocities were obtained in [8] for one-layer flows ($\eta_{0}=0$).

The linear analysis allows us to estimate the values of the relaxation parameter $\alpha$ to achieve the required approximation accuracy using model (12). As we can see from Fig. 2 with the given data and $\alpha=6$ the relative error $|c-c^{-}|/c$ between the phase velocities of the hyperbolic model and the original dispersive system is not more than 5 % for $k\in(0,1.5)$. With $\alpha=12$ and $\alpha=24$ accuracy increases to about 2 % and 1 %, respectively. Thus, a comparison of the dispersion relations (17) and (18) gives a criterion for choosing $\alpha$, which provides the specified accuracy of approximation between these models.

4 Stationary solutions

Continuous stationary solutions of model (12) for flows over a flat bottom satisfy the equations

[TABLE]

Here we denote the derivative with respect to $x$ by ‘prime’ and take $\varepsilon=1$. We also point out that the governing equations (12) (as well as original model (5)) are Galilean invariant if $b=0$. Therefore, for flows over a flat bottom, solutions of Eqs. (12) in the class of travelling waves are described by equations (19).

The normal form of Eqs. (19) is

[TABLE]

where

[TABLE]

The variables describing the lower dispersive layer $U$, $h$ and ‘internal’ variables $V$ and $\zeta$ are always continuous in flow. However, the variables describing the upper turbulent layer may be discontinuous when the supercritical–subcritical transition occurs. Stationary solutions of system (5) are obtained and studied in [10]. Below we show that such solutions can be constructed in the framework of model (12). Here we restrict our consideration to smooth solutions.

Further we are looking for stationary solutions of equations (20) having supercritical constant potential flow as $x\to-\infty$:

[TABLE]

Here $H_{0}$ and $U_{0}$ are the given flow depth and velocity at infinity and $F$ is the Froude number. To construct such a solution, it is necessary first to understand the asymptotic behaviour of the supercritical solution at negative infinity. Linearising equations (19), we obtain

[TABLE]

Here and below, the ‘tilde’ symbol denotes small perturbations of the corresponding variables:

[TABLE]

Following [24], we are looking for the solutions of Eqs. (22) that vanish at negative infinity in the form

[TABLE]

Substituting this representation of the solution into Eqs. (22) and expressing the unknown amplitudes (denoted with the ‘hat’ symbol) through $\hat{U}$ and $\nu$, we get

[TABLE]

Here $\alpha_{1}=\alpha H_{0}/g$ and $\hat{U}\leq 0$ is a given amplitude of the small perturbation of velocity. A non-trivial ($\hat{U}\neq 0$) solution of (22) in the form (23) exists if the parameter $\nu$ satisfies the equation

[TABLE]

If the parameter $\alpha_{1}$ is such that

[TABLE]

(or $\alpha>\alpha_{1}^{*}g/H_{0}$) then for a given Froude number $F>1$ there is a single positive root $\nu$ of the previous quadratic equation

[TABLE]

where $W=(l_{1}^{2}l_{2}^{2}\alpha_{1}+4l_{1}\alpha_{1}+4)F^{2}-4F^{4}-l_{1}^{2}l_{2}^{2}\alpha_{1}$. We use asymptotic expressions (23), (25) when the conditions are imposed at $x=x_{0}$ in the numerical treatment of stationary system (20). It should be noted that the amplitudes $\tilde{\bar{u}}$, $\tilde{q}$, $\tilde{\eta}$ and $\tilde{h}$ coincide with the obtained in [10] for asymptotic behaviour of stationary solutions in the framework of Eqs. (5), but $\nu$ is different. To solve ODE (20) numerically, we use the standard $ode45$ procedure of the MATLAB package.

An example of a stationary solution in the form of an undular bore is shown in Fig. 3. We choose the parameters $H_{0}=1$ and $U_{0}=1.2$ for fluid flow at the negative infinity and solve numerically ODE (20) with conditions at $x=0$ perturbed according to formulas (23) and (25) (we take here $\hat{U}=0.07$). The empirical constants, parameter $\alpha$ and non-dimensional gravity acceleration are as follows: $\sigma=0.15$, $\kappa=3$, $\alpha=10$ and $g=1$. In view of restriction (24) for the Froude number $F=1.2$ we have to choose $\alpha>1.32$. In this case the upper turbulent layer develops quite slowly (Fig. 3(a)) since the flow is supercritical everywhere (the determinants $\varphi$ and $\Delta$ in (20) are positive). We also note that the instantaneous depth $\zeta$ of the lower dispersive layer almost coincides with $h$. The larger $\alpha$, the better $\zeta$ approximates the depth of the lower potential layer $h$. The velocities in the layers $U$, $\bar{u}$ and the shear velocity $q$ of the turbulent layer are shown in Fig. 3(b). Note that in the vicinity of the wave crests, the variable $q$ reaches its maximum value, while the velocities $U$ and $\bar{u}$ are minimal. The smaller the Froude number $F>1$, the larger is the domain where the solution is supercritical. For considered example with $F=1.2$ the supercritical–subcritical transition happens at $x\approx 30.5$, i.e. the determinant $\Delta$ in (20) changes sign from positive to negative. This leads to a jump of the variables of the upper turbulent layer.

System (19) can be used to describe the structure of turbulent bore if the dispersion effects are negligible ($\alpha=0$). Let the flow depth and velocity be the following $h=H_{0}$, $u=U_{0}$ and $\eta=0$ for $x<x_{0}$. At the point $x=x_{0}$ the upper turbulent layer is formed ($\eta>0$ for $x>x_{0}$). The values of the functions at $x=x_{0}$ are denoted by the subscript “0”. We choose here $\zeta_{0}=H_{0}$, $V_{0}=0$ and $\alpha=0$. From the two last equations in (19) we get $\zeta=\zeta_{0}$, $V=0$. Let us assume that as $x\to x_{0}$ (or $\eta\to 0$) there exist finite limits of the functions $\bar{u}\to\bar{u}_{0}$ and $q\to q_{0}$. Assuming that the unknown functions and their derivatives are bounded, from the third, fourth and fifth equations of system (19) we have

[TABLE]

After eliminating the limiting value $\eta^{\prime}$ from these relations, we obtain

[TABLE]

We show below that Eqs. (20) with $\alpha=0$ and conditions (26) describe the stationary turbulent bore when the hydrostatic approximation can be applied [21].

5 Numerical results

In this section, we present the results of numerical simulation of the formation and evolution of turbulent bores. First of all, we show that terms of order $O(\sigma\varepsilon^{2})$ practically do not affect the numerical results in framework of model (5). Then, we clearly demonstrate on an example of Favre waves that as the $\alpha$ parameter increases, the solutions of the hyperbolic system (12) converge to the solutions of the corresponding dispersive equations. We also present the results of modelling of undular and monotone bores formed as a result of the flow past an obstacle. In the end of this section we consider the shoaling and breaking of a solitary wave.

To solve differential balance laws (12) numerically, we implement here the Nessyahu–Tadmor second-order central scheme [28]. The same method is used to solve dispersive equations (5) for the unknown functions $\mathbf{u}=(h,\eta,K,\bar{u}\eta,q)$, but at each time step we apply the Thomas algorithm (a simplified form of Gaussian elimination) to find the velocity $U$ in the lower layer from second-order equation (7).

5.1 Favre waves

With accuracy up to the terms of order $O(\sigma\varepsilon^{2})$, system (5) for the evolutionary variables $\mathbf{u}=(h,\eta,K,\bar{u}\eta,q)$ can be written as

[TABLE]

with equation (7) for the variable $U$ to be solved at each time step. Let us show that there is almost no difference in calculations on the basis of equations (27), (7) and their original version used in [10], where equation (8) stands for the third equation in system (27). To do this, we repeat the test concerning Favre waves considered in [10].

We perform calculations in dimensionless variables in the domain $x\in[0,100]$ on a uniform grid with $N=1000$ nodes. To start calculations we specify unknown variables $\mathbf{u}$ in the node points $x=x_{j}$ at $t=0$ as follows $\eta=0.01$, $h+\eta=1$, $\bar{u}=K=U_{0}$, and $q=0$. The right boundary $x=100$ is an impermeable wall. To satisfy the impermeability condition ($\bar{u}=U=0$) we set $\bar{u}_{N+1}=-\bar{u}_{N-1}$ and $K_{N+1}=-K_{N-1}$. At $x=0$ the initial data are used as boundary condition. We choose here $g=1$, $\varepsilon=1$, $\sigma=0.15$ and $\kappa=3$.

Instead of the upstream velocity $U_{0}$, it is more convenient to prescribe the wave Froude number $F_{w}$:

[TABLE]

where $D$ denotes the propagation speed of the reflected wave. A simple relation between $U_{0}$ and $\bar{F}$ can be obtained in the form (see [10] for details)

[TABLE]

In supercritical regime ($F_{w}>1$), undulations start developing at the bore front. When the wave Froude number is approximately between 1.3 and 1.4, the transition from the undular bore to the breaking bore occurs. This problem was studied in experimental works [7, 37] and numerically [36, 10].

The results of calculations for the wave Froude numbers $F_{w}=1.28$ and $F_{w}=1.40$ are shown in Fig. 4. We clearly see that terms of order $O(\sigma\varepsilon^{2})$ in (5) almost do not affect on the numerical result, especially for an undular bore. A similar coincidence is observed for the other tests considered [10] including the shoaling and breaking of a solitary wave propagating in a long channel of mild slope. Thus, for carrying out calculations, the two-layer model (5) can be used in simplified form (27), (7).

Let us show that the solutions of the proposed hyperbolic model (12) approximate the solutions of dispersive equations (27), (7). Fig. 5(a) presents free surface $z=h+\eta$ for the previous test for $F_{w}=1.28$ obtained on the basis of system (12) with $\alpha=10$ (solid curve) and in the framework of Eqs. (27), (7) (dashed curve). For the hyperbolic equations (12) at $t=0$ we add conditions for the ‘instantaneous’ variables as follows $\zeta=h$ and $V=0$. There is a good match for the leading wave, while for the secondary waves some difference is observed. An almost complete coincidence of the solutions of the dispersion model and its hyperbolic approximation can be achieved by increasing the parameter $\alpha$ (see Fig. 5(b)). However, for this, the spatial resolution should be improved the calculation by using the hyperbolic model (12). This is due to the fact that as $\alpha$ increases, the velocity of the characteristics (in absolute value) also increases and, consequently, the time step becomes smaller. As a result, it leads to a stronger effect of the numerical viscosity.

For other Froude numbers and, in particular, in the transition to monotonous bore, solutions of hyperbolic equations (12) also approximate solutions of the dispersion model.

It is important to note that the verification of the dispersive model (5) performed in [10] is based on a comparison of the maximum and minimum amplitude of the first wave with the experimental data [7] and [37]. As we can see from Fig. 5, the results of calculations for the first wave by using model (12) with $\alpha=10$ and $\alpha=20$ almost coincide, so there is no need to choose large values of $\alpha$. In addition, a finer mesh is necessary for large values of the relaxation parameter $\alpha$.

5.1.1 Dependence of the solution on the empirical parameters

Both models (27) and (12) include the empirical parameters $\sigma$ and $\kappa$, which may differ slightly from the fixed above values. In particular, the following range of these parameters were used in [24, 21, 10] for layered hydrostatic and dispersive models: $\sigma\in[0.15,0.20]$ and $\kappa\in[2,6]$. A change in the empirical parameters in the specified range does not lead to a significant change in the solution. Let us consider this in more detail.

Fig. 6 shows the results of the Favre waves calculation ($F_{w}=1.28$) obtained with the help of hyperbolic model (12) for different empirical parameters $\sigma$ and $\kappa$. As can be seen from Fig. 6(a), an increase in the parameter $\sigma$ leads to a decrease in the thickness $h$ of the lower layer, as well as to an insignificant decrease in the wavelength (dashed curves). At the same time, the fluid depth $h+\eta$, in contrast to the ‘artificial’ internal interface, almost coincides with the case $\sigma=0.15$ and $\kappa=3$ (solid curves), especially for the leading wave. In the case of an increase in the parameter $\kappa$ (see Fig. 6(b)), the thickness of the lower layer increases and the wavelength is slightly longer. However, in the vicinity of the leading wave the difference in the results is minimal. Thus, we can use fixed parameters $\sigma=0.15$ and $\kappa=3$. Moreover, there is no need to change these parameters for other tests.

5.1.2 A comparison of computational time for the hyperbolic and dispersive models

Let us find out the dependence of the computational time for model (12) on the value of the relaxation parameter $\alpha$. We choose the same data as at the beginning of this section, in particular, $F_{w}=1.28$, $N=1000$ and $t_{\max}=90$. It is obviously, the computation time depends on the velocities of the characteristics, since the time step is determined by the Courant condition

[TABLE]

where $\Delta x$ is the spatial grid spacing, $\lambda_{k}$ is the characteristic velocity and $\mathbf{U}_{j}^{n}$ is the numerical solution at $t=t^{n}$ and $x=x_{j}$. In what follows, we take ${\rm CFL}=0.475$. Note that $\alpha=0$ corresponds to a hydrostatic approximation. In this case, the calculation time is minimal, since the maximum value of the characteristic velocity increases with the relaxation parameter $\alpha$.

The calculation time $T_{\alpha}$ (in seconds) for this test on the computer being used is presented in Table 1 for different values of $\alpha$. We also indicate here the number of iterations $M_{\alpha}$ (time steps) needed to reach the final time $t_{\max}$. This parameter is determined only by the used numeric code and does not depend on the CPU power. As we can see from Table 1, the computation time and the number of iterations increase with the relaxation parameter $\alpha$. However, the increment in computation time (or number of iterations) decreases with the growth of $\alpha$. In particular, $T_{6}-T_{3}=2.5382$ ($M_{6}-M_{3}=1083$), $T_{12}-T_{9}=1.4285$ ($M_{12}-M_{9}=737$). A further increase in the parameter $\alpha$ at the given grid resolution $N$ is impractical, since an increase in the number of steps in time enhances the numerical viscosity. Therefore, it is necessary to improve the spatial resolution, as it was done in the example shown in Fig. 5(b).

For the numerical solution of the dispersive system (27), we employ the approach proposed in [25] with the inversion of the elliptic operator (7) at each time step. Although model (27) is not a hyperbolic system, we apply condition (28) to determine the time step $\Delta t$. We take the velocity of the characteristic $\max\limits_{k}|\lambda_{k}|$ of the hydrostatic model ($\varepsilon=0$) with factor $a_{d}>1$ needed to stabilize the numerical algorithm. We note that system (27) contains a term of the form $\varepsilon^{2}(h^{2}U_{x})_{x}$. Therefore, to stabilize the numerical code the velocity of the characteristic should be larger. The choice of the factor $a_{d}$ depends on the considered problem. Here we take $a_{d}=3/2$. In this case the computation time is $T=12.1951$ seconds and it takes $M=3840$ iterations to reach the final time $t_{\max}$.

Comparing with the data from Table 1, we see that the calculation time for the dispersion model (27) approximately corresponds to the hyperbolic model with $\alpha=9$, despite the greater number of iterations required for system (12). The computational time spent per iteration for equations (27) is longer than for system (12). It is clear, if the spatial resolution improves, the time spent on iteration for the dispersive model increases significantly due to the inversion of the elliptic operator (7). In this sense the hyperbolic model works faster for all reasonable $\alpha$. In most cases $\alpha\in[5,8]$ provides a fairly good approximation, and the computation time for the hyperbolic model is shorter, despite the greater number of iterations.

Let us note that $\alpha$ should be multiplied by $g/H_{0}$ if we use dimensional variables.

5.2 Flows over a local obstacle

Let us consider the flows arising from the interaction of the upstream supercritical flow with a local obstacle located in the vicinity of the outlet section of the channel. Here we perform calculations using of the hyperbolic system (12) in the domain $x\in[0,35]$ and set the constant flow $h_{0}=\zeta_{0}=0.99$, $\eta_{0}=0.01$, $U_{0}=\bar{u}_{0}=1.2$, $q_{0}=V_{0}=0$ (in this case the Froude number $F=1.2$) as the initial data at $t=0$. These data are used as a boundary condition at $x=0$. On the right boundary we pose the Neumann condition. The relaxation parameter $\alpha$ is equal to 10. As before, we take $g=1$, $\varepsilon=1$, $\sigma=0.15$ and $\kappa=3$.

The formation of bore is carried out by controlling the obstacle. At the initial time, the bottom is flat ($b=0$). During time $t_{1}=25$, the height of the smooth obstacle of width $5$ located on the right edge is increased to $b^{*}=0.2$, then, over time $t_{2}=50$ it is reduced to $b_{*}=0.1$. It leads to the formation of an undular bore shown in Fig. 7. At the initial stage of the process, the thickness of the upper turbulent layer increases substantially and the shear velocity $q$ sharply increases at the wave front (Fig. 7(a)). Then a wave train propagating upstream the flow is formed (Fig. 7(b)). The thickness of the upper turbulent layer becomes insignificant in the vicinity of the leading wave and the flow qualitatively corresponds to the stationary solution (Fig. 3(a)). Note that the speed of waves propagation depends on the height of the obstacle. It should be greater than a certain limit value for the given number Froude of the upstream flow.

Let us repeat the previous calculation, with the difference that for $t>t_{2}=50$ the bottom topography becomes even ($b_{*}=0$). In this case, a solitary wave slowly propagating upstream is formed. This wave (the free surface $z=h+\eta$) at $t=225$ is shown in Fig. 8(a) (solid curve). The dashed line corresponds to the solution of stationary Eqs. (20) with the same parameter $\alpha=10$. As can be seen from the figure, the profile of a solitary wave is quite well described by a stationary solution.

With increase in the velocity of upstream flow (and the height of the obstacle), a monotone turbulent bore is formed. For this calculation we choose $t_{1}=15$, $b^{*}=0.3$, $t_{2}=30$, $b_{*}=0.2$, $U_{0}=\bar{u}_{0}=1.4$ and the remaining parameters do not change. The result of calculations on the basis of Eqs. (12) at $t=285$ (the free surface $z=h+\eta+b$ and the interface $z=h+\eta$) is shown in Fig. 8(b) (solid curves 1 and 2, line 3 — bottom topography). The solution of stationary equations (20) with $\alpha=0$ (it corresponds to the hydrostatic approximation) and conditions (26) at $x=x_{0}=5$ is shown by the dashed lines. This solution is in a good agreement with non-stationary calculation. Note that a local area of the subcritical flow is formed above the obstacle (the variable $\Delta$ in (20) changes the sign), at the exit from the channel the flow becomes supercritical again. For this reason, the stationary solution is constructed only in the interval before the obstacle. The given example shows that Eqs. (12) with a fixed parameter $\alpha$ allow one to describe flows when the effect of dispersion is important (Fig. 8(a)), as well as flows that can be described in the hydrostatic approximation (Fig. 8(b)).

Remark 1

There is an additional advantage of the hyperbolic approximation (12), which consists in the formulation of boundary conditions. As a rule, when implementing Godunov-type numerical schemes, the formulation of boundary conditions depends on the number of incoming and outgoing characteristics at the boundary of the computational domain. Due to seven equations in system (12), it is convenient to use central schemes [28], which do not require exact or approximate solution of the Riemann problem. In this case we have to specify the boundary conditions for all unknown functions. Let the flow at the exit from the computational domain be supercritical and the disturbances propagating upstream do not reach the left boundary of the region, as in the example with the flow over an obstacle. Then the Neumann conditions at the right boundary do not affect the solution inside of the computational domain. It is not valid for the dispersive model since at each time step a non-local operation of inversion of the elliptic operator is required. As a result, the boundary conditions affect the solution.

5.3 The breaking of a solitary wave

The evolution of breaking solitary wave on a mild sloping beach in the framework of the two-layer dispersive model (5) was considered in [10], where an excellent agreement with experimental data [16] was shown. Therefore, to verify the proposed hyperbolic model (12) in this test, it suffices to compare the numerical results with calculations using of dispersive equations (5).

Here we perform the calculations in dimensional variables ($g=9.8$ m/s2) in the domain $x\in[0,200]$ using of $N=1500$ nodes for space resolution. As in the previous tests we take $\sigma=0.15$ and $\kappa=3$. The bottom topography has a slope of 1/60 (curve 1 in Fig. 9). To avoid modelling of the wave propagation over dry bottom (run-up), we introduce a shelf zone near the right boundary. Initial data represent a solitary wave having amplitude $a_{0}=0.3344$ m propagating with the velocity $C_{0}=\sqrt{g(a_{0}+H_{0})}$, where $H_{0}=2.2$ m is the undisturbed water depth. The initial thickness of the upper turbulent layer $\eta_{0}$ is equal to $0.01$ m. On the left and right walls we assume reflecting boundary conditions and take $\alpha=11g/H_{0}$ for system (12). Fig. 9 shows the results of calculations of the free surface at $t_{*}=t\sqrt{g/H_{0}}=50$ and $t_{*}=75$ with respect to models (5) and (12) (solid and dashed curves, correspondingly).

The temporal evolution of the free-surface displacements at different positions along the channel is shown in Fig. 10. Here we also compare the calculation results for models (5) and (12). The original dispersive model was verified in [10] by experimental data [16] for the time history of the free surface. Therefore, it can be said that hyperbolic model (12) also gives an excellent agreement with the experiment [16] for the problem of the breaking of a solitary wave.

It should be noted that it is impossible to model the run-up and run-down using governing equations (5) or (12) because a singularity occurs after the upper turbulent layer reaches the bottom. Modifying the systems to describe such process is an open problem.

6 Conclusion

We derived hyperbolic equations (12) approximating the two-layer long-wave model (5) proposed in [10] and taking into account dispersion and vortex effects. The construction of the hyperbolic model is based on the use of additional ‘instantaneous’ variables describing the flow of a dispersive fluid layer [24, 22]. The advantage of the hyperbolic model is simpler numerical implementation and formulation of boundary conditions. In contrast to the solution of dispersive equations of the Green–Naghdi type, in this case there is no need for a time-consuming operation for inverting an elliptic operator at each time step.

The velocities of characteristics of the proposed equations (12) are determined and hyperbolicity of this system (for sufficient large relaxation parameter) is established. Linear analysis of the derived and original models is performed and the dispersive relations are obtained. A comparison of these relations allows us to observe the effect of the relaxation parameter on the accuracy of the hyperbolic approximation. Stationary solutions of the hyperbolic system (12) are studied. In particular, we find the asymptotic behaviour of the stationary solution at infinity and construct a solution describing an undular bore with a turbulent sub-surface layer (Fig. 3). Further, we apply the hyperbolic system to model Favre waves, non-hydrostatic flows over a local obstacle and the evolution of breaking solitary wave on a sloping beach.

Firstly, we establish that terms of order $O(\sigma\varepsilon^{2})$ in model (5) have practically no effect on the results of calculations for the Favre waves. Therefore, model (5) can be used in a simpler form (27) (see Fig. 4). Then we clearly show that solutions of hyperbolic equations (12) approximate solutions of the dispersive model (27) (Fig. 5). It is found that for a good approximation of the leading wave it is not necessary to choose large values of the relaxation parameter. We show that variation of the empirical parameters of the model has little effect on the solution of the equations of motion (Fig. 6). Moreover, these parameters are not required to change for other tests. We also compare the time consuming for the dispersive and hyperbolic models. The formation and evolution of a turbulent bore in the interaction of a constant upstream flow with a local obstacle are considered. It is established that at the initial stage of wave bore formation a near-surface turbulent layer intensively develops. The thickness of the upper turbulent layer becomes smaller at least for the leading wave (Fig. 7). By controlling the obstacle in the vicinity of the output section of the channel, various wave configurations can be realized. In particular, it is possible to obtain a solitary wave or monotone bore. Such examples are presented in Fig. 8. It is shown that the structure of these waves can be described by stationary equations. We also present the results of the numerical modelling of the shoaling and breaking of a solitary wave (Fig. 10) which verified by experimental data.

The used method for constructing the hyperbolic approximation of the dispersive system can be applied to other models, in particular, to multilayer dispersive equations of stratified flows [11, 23]. The construction of such hyperbolic approximations of dispersive equations seems especially promising for the numerical simulation of spatial flows.

Acknowledgments

This work was partially supported by the Russian Foundation for Basic Research (project 19-01-00498). The authors thank V.Yu. Liapidevskii, S.L. Gavrilyuk and I.V. Stepanova for fruitful discussions. They also address special thanks to the reviewers for their helpful comments and suggestions.