Kseniya A Ivanova
email: kseniya.ivanova@univ-amu.fr

S L Gavrilyuk
email: sergey.gavrilyuk@univ-amu.fr

B Nkonga

G L Richard
email: gael.richard@math.univ-toulouse.fr

N Favrie

FORMATION AND COARSENING OF ROLL-WAVES IN SHEAR

Keywords: shear shallow water equations, non-conservative hyperbolic equations, Godunovtype scheme, shock waves, roll waves, convergent circular hydraulic jump équations d'écoulement cisaillé en eau peu profonde, équations hyperboliques non-conservatives, schéma de Godunov, ondes de choc, trains de rouleaux, ressaut hydraulique circulaire shear shallow water flows, roll-waves, coarsening, hyperbolic equations, Godunov method Non-conservative hyperbolic equations, Godunov-type methods, roll waves Convergent circular hydraulic jump, radial dam break problem, radial oscillations, instabilities, shock waves, non-conservative hyperbolic equations, Godunov-type scheme

These validations demonstrate the capability of the model and numerical method to solve challenging multi-dimensional problems of shear shallow water flows.

First and foremost I would like to express my sincere gratitude to my supervisor Prof. Sergey Gavrilyuk. This thesis would not have been possible without the help that I have received from him. His fruitful supervision, encouragement, constant support and push showed me the way whenever I faced hardships. I am heartily thankful to my co-supervisor Nicolas FAVRIE. His cooperation and advices have been essential for me to pursue this project.

I would like to extend my thanks to the members of the jury, in particular, Rémi AB-GRALL and Jean-Paul VILA, thesis reviewers, for careful reading of this manuscript and for the thoughtful comments and suggestions. My thanks also go to my examiners Gennady EL and Thierry FOGLIZZO for the interest they showed for my work. I want to thank Henri GOUIN who has honored me by presiding the jury.

Moreover, I would like to thank Boniface NKONGA for his encouragement, insightful comments, and pertinent questions.

My special thanks go to Sarah HANK for stimulating discussions and her friendship. She always helped me when I had some troubles with coding. I thank my fellow labmates in SHOCK Group: Jacques MASSONI, Olivier LE METAYER, Eric DANIEL, Kevin SCHMIDMAYER, Firas DHAOUADI and Sergey TKACHENKO for inspiring environment and discussion we had about research (and not only).

My special thanks and sincere gratitude to Pascal CAMPION for administrative help and moral support throughout three years. I would like to express my gratitude to the staff of IUSTI for their kindness and technical support. I am grateful to my friends: Auna BARLUKOVA, Olga OBRAZCOVA, Olga GORININA, Anna CHERNISHOVA and Maria KAZAKOVA for informal chatting about work in progress and future plans. I wish them a successful scientific carrier.

Especially, I would like to thank my big family: my mother Nina, my brothers Iliya and Alexandre, and my sisters Olesya and Inga.

Last but not the least, I thank Kirill PETROV for his love and support.

This work was partially supported by l'Agence Nationale de la Recherche, France (grant numbers ANR-13-BS01-0009, ANR-11-LABX-0092, and ANR-11-IDEX-0001-02).

Male and female represent the two sides of the great radical dualism. But in fact they are perpetually passing into one another. Fluid hardens to solid, solid rushes to fluid. There is no wholly masculine man, no purely feminine woman.

Margaret Fuller

Thermodynamics is a funny subject. The first time you go through it, you don't understand it at all. The second time you go through it, you think you understand it, except for one or two small points. The third time you go through it, you know you don't understand it, but by that time you are so used to it, it doesn't bother you any more.

Arnold Sommerfeld

Abstract

This thesis is devoted to the numerical modelling of multi-dimensional shear shallow water flows.

In 1D case, the corresponding equations coincide with the equations describing non-isentropic gas flows with a special equation of state. However, in the multi-D case, the system differs significantly from the gas dynamics model. This is a 2D hyperbolic non-conservative system of equations which is reminiscent of a generic Reynolds averaged model of barotropic turbulent flows. The model has three families of characteristics corresponding to the propagation of surface waves, shear waves and average flow (contact characteristics).

First, we show the ability of the one-dimensional conservative shear shallow water model to predict the formation of roll-waves from unstable initial data. The stability of roll waves is also studied.

Second, we present in 2D case a new numerical scheme based on a splitting approach for non-conservative systems of equations. Each split subsystem contains only one family of waves (either surface or shear waves) and contact characteristics. The accuracy of such an approach is tested on exact 2D solutions describing the flow where the velocity is linear with respect to the space variables, and on the solutions describing 1D roll waves. The capacity of the model to describe the full transition as commonly seen in the formation of roll waves is shown:

• from uniform flow to one-dimensional roll waves • 2D transverse "fingering" of roll wave profiles.

Finally, we model a circular hydraulic jump formed in a convergent radial flow of water.

Obtained numerical results are qualitatively similar to those observed experimentally: oscillation of the hydraulic jump and its rotation with formation of a singular point.

Résumé

Cette thèse porte sur la résolution numérique du modèle multi-dimensionnel d'écoulement cisaillé en eau peu profonde. Dans le cas d'un mouvement unidimensionnel, ces équations coïncident avec les équations de la dynamique des gaz pour un choix particulier de l'équation d'état. Dans le cas multi-dimensionnel, le système est complètement différent du modèle de la dynamique des gaz. Il s'agit d'un système EDP hyperbolique 2D non-conservatif qui rappelle un modèle de turbulence barotrope. Le modèle comporte trois types d'ondes correspondant à la propagation des ondes de surface, des ondes de cisaillement et à celle de la discontinuité de contact.

Dans un premier temps, nous montrons les capacités du modèle 1D conservatif d'écoulement cisaillé en eau peu profonde à prédire la formation des trains de rouleaux à partir des données initiales instables. Nous étudions également la stabilité des trains de rouleaux.

Dans un second temps, nous présentons dans le cas 2D un schéma numérique basé sur une nouvelle approche de "splitting" pour les systèmes d'équations non-conservatives. Chaque soussystème ne contient qu'une seule famille d'ondes: ondes de surface ou ondes de cisaillement, et discontinuité de contact. La précision d'une telle approche est testée sur des solutions exactes 2D décrivant l'écoulement lorsque la vitesse est linéaire par rapport aux variables spatiales, ainsi que sur des solutions décrivant des trains de rouleaux 1D.

La capacité du modèle à décrire le mécanisme complet intervenant dans le formation des trains de rouleaux est demontré:

• d'écoulement uniforme vers des trains de rouleaux 1D

• "digitation" du profil transversal des trains de rouleaux 2D.

INTRODUCTION

In this introductory chapter, I present general motivation, purpose, outline and contributions.

General Motivation

Fig. 1.1 -The image on the left is the Supernova remnant exploded 300 years ago in the Milky Way.

The aspherical explosion gave birth to a compact object. The fountain on the right uses a fluid analogy to demonstrate the physical process responsible for the asymmetric explosion to the public of the Science Museum "Palais de la Découverte" in Paris since 2015. These pictures were taken from http: // irfu. cea. fr/ Sap/ Phocea/ Vie_ des_ labos/ Ast/ ast_ sstechnique. php? id_ ast= 4058

Shallow water theory is commonly used to describe the propagation of surface gravity waves, and is applied to flows where the characteristic wavelength is much larger than the fluid depth.

In many situations it is a reasonable assumption. For example, the shallow water equations (also called Saint-Venant equations) are commonly used to describe

• fluid flow in the rivers, estuaries, channels, coastal areas

• oceanic and atmospheric circulations

• tsunami (also called seismic sea wave or tidal wave)

• debris flows

• volcanic eruptions

• supernovae explosion.

In particular, to explain the supernova explosion, Foglizzo et al. (CEA-Saclay, France) have performed an interesting experimental study, which demonstrate that the shallow water theory can help to understand this phenomenon (see Figure 1.1) [12], [13], [14].

The classical shallow water equations (Saint-Venant equations) over a flat bottom read:

             ∂h ∂t + div(hu) = 0, ∂hu ∂t + div hu ⊗ u + gh 2 2 I = 0, (1.1.1)
where t is the time, h is the water depth, u is the depth averaged horizontal velocity, g is the gravity acceleration constant, I is the identity matrix. The equations (1.1.1) admit the energy conservation law that plays the role of mathematical "entropy" of the system:

h |u| 2 2 + gh 2 t + div hu |u| 2 2 + gh = 0.
The derivation is based on the smallness of the parameter ε = H L , where H is a characteristic fluid depth, L is a characteristic wavelength, and the assumption that the flow is potential (the horizontal velocity does not depend on the vertical coordinate). The Saint-Venant equations form a nonlinear hyperbolic system of conservation laws. The long wave approximations allows us to reduce the dimension of the problem, and hyperbolicity property allows us to describe adequately the propagation of perturbations in the flows. The analogy between shallow water equations and gas dynamics equations plays a significant role in the mathematical analysis of the Saint-Venant equations. In particular, it allows us to interpret a hydraulic jump as a shock wave.

In the present work, we will focus on shear effects (creation of vorticity in the flow). Indeed, when the shear effects are neglected, shallow water equations are not able to provide the information related to the length of hydraulic jump, predict the formation of large scale eddies ("rollers") appearing in the hydraulic jumps near the free surface, the position of the hydraulic jump toe and, in particular, transverse jump toe perimeter profiles.

For shear flows with varying in space and time vorticity a model was recently proposed where the governing equations are obtained by depth averaging of Euler equations without assuming potential flow [32], [17], [26], [27], [28]. In Appendix A one can find the derivation of the multi-dimensional shear shallow water model.

The one-dimensional shear shallow water model (SSWM) representing a system of non-linear conservation laws of mass, momentum and energy without source terms is:

                           ∂h ∂t + ∂hU ∂x = 0, ∂hU ∂t + ∂ ∂x hU 2 + gh 2 2 + (ϕ + Φ)h 3 = 0, ∂hE ∂t + ∂ ∂x hU E + gh 2 2 + (ϕ + Φ)h 3 U = 0.
(1.1.2)

Here E = U 2 2 + 1 2 (gh + (ϕ + Φ)h 2 ).
The system (1.1.2) admits the "entropy" conservation law:

∂ ∂t (hΦ) + ∂ ∂x (hU Φ) = 0. (1.1.3)
This model can be considered as an extension of the standard shallow water model, involving some additional unknowns, the enstrophy (squared vorticity) Φ of large eddies formed in the roller, and an additional parameter, the enstrophy of small vortexes at the vicinity of the bottom, ϕ. We suppose that the small-scale enstrophy ϕ is constant. It is parameter of our system. It depends on the wall roughness, and the plane inclination. It increases if the plane slope increases. The shear shallow water model complemented by friction terms provides a more realistic description of some flow phenomena, such as the hydraulic jumps or travelling waves down inclined plane (also called "roll waves") (see Figure 1.2). This model is similar to the Euler equations of compressible flows (ρ → h, p → gh 2 2 + (ϕ + Φ)h 3 ). When one takes ϕ = 0, Φ = 0, the system is reduced to the classical Saint-Venant equations.

The multi-dimensional model of shear shallow water flows over a flat bottom without friction effects consists of the equation of mass balance, equation of average momentum, and evolution equation of the shear stress tensor [32,17,28]:

                           h t + div(hu) = 0, (hu) t + div hu ⊗ u + gh 2 2 I + hP = 0, DP Dt + ∂u ∂x P + P ∂u ∂x T = 0, D Dt = ∂ ∂t + u • ∇, (1.1.4)
where t is the time, x = (x, y) T are the Cartesian coordinates, u = (u, v) T is the depth averaged horizontal velocity, h is the fluid depth, g is the gravity acceleration constant, and P is shear stress tensor defined as follows:

P = 1 h h 0 (ũ -u) ⊗ (ũ -u)dz.
Here ũ = ũ(t, x, y, z) is the instantaneous horizontal velocity depending of the vertical coordinate z. The tensor P is the symmetric and positive definite, and measures the distortion of the instantaneous horizontal velocity. We obtain reduction to the multi-dimensional Saint-Venant equations in the limit P = 0.

The striking mathematical analogy with the Reynolds averaging equations of barotropic compressible turbulent flows allows us to call R = hP the Reynolds stress tensor, and P the reduced Reynolds stress tensor. For mathematical reasons, the choice of P is more convenient than tensor R. For simplicity, both R and P will be further refereed to as the Reynolds stress tensor. The system (1.1.4) admits the energy and "entropy" conservation laws: The variable Ψ is referred to as 'entropy' because this quantity is transported along the mean flow as the true entropy for the Euler equations of compressible fluids. Also, we will see that this quantity will increase through the shocks in analogy with the conventional entropy.

∂
Definition 1.

A system is in conservative form if it can be written as follows (in two space dimensions) ∂ ∂t U(x, y, t) + ∂ ∂x F(U(x, y, t)) + ∂ ∂y G(U(x, y, t)) = 0.

Here U : R 2 × R → R m is a m-dimensional vector of conserved quantities and F, G are flux functions. It is hyperbolic if for any real ξ and η the matrix A(ξ, η) = ξ ∂F ∂U + η ∂G ∂U has only real eigenvalues and a complete set of the corresponding left eigenvectors. A similar definition of hyperbolicity can be given in the case of non-conservative systems.

The multi-dimensional system (1.1.4) is hyperbolic but not in conservative form (see Definition 1). For six unknowns (the fluid depth, two components of the depth averaged horizontal velocity, and three independent components of the symmetric Reynolds stress tensor) the system has only five independent conservation laws. This fact can be formulated as theorem.

Theorem 1.

The only linearly independent conservation laws admitted by the governing equations (1.1.4) are those of the mass, momentum, energy and mathematical "entropy".

This theorem is proved in Appendix A of Chapter 3. The definition and computation of discontinuous solutions for non-conservative hyperbolic equations is a challenging problem. The non-conservative systems appear in compressible turbulence [6,2], multi-layer shallow water flows [25,4,23,5,1,22], multi-phase fluid flows [3,20,29,30,31,16,9], solid-fluid systems [15,24,19]).

Essentially, four approaches are used for numerical solving of non-conservative systems of equations. The most classical one is based on the Volpert path definition for non-conservative products [10]. The second one is the formulation of an augmented system of 'Rankine-Hugoniot relations' through the study of travelling wave solutions of an extended system of equations approximating a given system (formulation of kinetic relations) [34,21]. The third one is based on the relaxation technique where the studied system is approximated by a new hyperbolic system with all linearly degenerate in the sense of Lax eigenfields [7,9]. Finally, additional shock relations can be formulated from the compatibility between theoretical and experimental results [11,31]. Excepting the first approach (more formal and hence less precise), all the approaches mentioned above are not universal : they are usually specific to the model under study.

Purpose

The aims of this thesis are:

• to show the ability of the one-dimensional shear shallow water model to predict the formation of roll-waves from the unstable initial data down an inclined plane and study their stability is demonstrated (see Figure 1.4). Also, the stability of such waves was studied and the so-called wave "coarsening" phenomenon (see Definition 3) produced by coalescence of non-stationary roll-waves was demonstrated (see Figure 1.5).

•

Definition 3.

If the roll-wave train contains waves of different lengths, they begin to interact and coalesce.

The short waves transfer their energy to long waves, and finally a roll-wave train of larger wavelength appears. This physical process is called "coarsenning", or, as it was proposed by Brock, "growth by overtaking" [8].

We consider a constant slope channel. The equations are solved by a shock-capturing Godunov type scheme (MUSCL-Hancock extension) [33,18]. The scheme is based on the conservative finite volume formulation.

A simplified 2D extension of the model, where the Reynolds stress tensor is replaced by a spherical one, is also presented. Such an approximation corresponds to the isotropic turbulence, often used in applications. Within the framework of this assumption, we deal with a multi- dimensional conservative system, but this system is not able to describe any transverse structure (see Figure 1.6). In the case of the complete non-conservative model, the results will be different (the transverse structure will be observed).

In Chapter 3, we describe a complete multi-dimensional model of shallow water flows which includes effects relating to vertical shear in two horizontal space dimensions. The resulting non-conservative system of hyperbolic equations cannot be solved with standard approaches for conservative hyperbolic system of equations. Instead, it was split into subsystems, which allowed us the construction of HLLC-type approximate Riemann solvers. A new numerical method was thus developed which was applied further to a number of physical problems. It was shown that it can reproduce a physically realistic flow behaviour which cannot be covered in the classical model of shallow water flows. In this chapter we define an almost conservative formulation of the governing equations based on a splitting procedure allowing us to define the discontinuous solutions. A conventional geometric splitting of (1.1.4) is first applied consisting in solving of the governing equations first in x and then in y direction. We will do two steps more for each uni-directional subsystem referring to this as a 'physical' splitting. Roughly speaking, each uni-directional subsystem, describing two types of waves (surface and shear waves), is split into two subsystems which are hyperbolic and contain only one type of waves. Each physical subsystem admits its own energy conservation law, and its own 'entropy'. In physical terms, one can say that we perform a thermodynamically compatible splitting. An analogous splitting was proposed earlier in [15] in hyperelasticity. The subsystem in x-direction reads:

                                                                 h t + uh x + hu x = 0, u t + uu x + gh x + 1 h (hP 11 ) x = 0, v t + uv x + 1 h (hP 12 ) x = 0, P 11t + uP 11x + 2P 11 u x = 0, P 12t + uP 12x + P 12 u x + P 11 v x = 0, P 22t + uP 22x + 2P 12 v x = 0. (1.3.1)
The system is hyperbolic and admits three types of waves : a contact discontinuity propagating with the velocity u, surface gravity waves propagating with the velocity u±a, with a 2 = gh+3P 11

                                                                 h t + uh x + hu x = 0, u t + uu x + gh x + 1 h (hP 11 ) x = 0, v t + uv x = 0, P 11t + uP 11x + 2P 11 u x = 0, P 12t + uP 12x + P 12 u x = 0, P 22t + uP 22x = 0. (1.3.2)
The equations (1.3.2) admit the following conservative form :

                                                                 h t + (uh) x = 0, (hu) t + hu 2 + g h 2 2 + hP 11 x = 0, (hv) t + (huv) x = 0, h u 2 +v 2 2 + gh 2 + P 11 +P 22 2 t + hu u 2 +v 2 2 + gh 2 + P 11 +P 22 2 + gh 2 2 u + hP 11 u x = 0, P 12t + (uP 12 ) x = 0, (hP 22 ) t + (huP 22 ) x = 0. (1.3.3)
They also admit the 'entropy' equation :

D Dt P 11 P 22 -P 2 12 h 2 = 0, D Dt = ∂ ∂t + u ∂ ∂x .
One solves the 6 conservation laws (1.3.3) at this step. The 'entropy' will increase through the shock. Indeed, the ratio P 12 /h and the component P 22 are conserved through the shock, so we need only the increase of P 11 /h 2 . This is similar to the case of 1 D flows [26,27]. The subsystem for b-waves reads:

                                                                 h t = 0, (hu) t = 0, (hv) t + (hP 12 ) x = 0, P 11t = 0, P 12t + P 11 v x = 0, P 22t + 2P 12 v x = 0. (1.3.4)
An 'almost' conservative form of (1.3.4) for b-waves is :

                                                                 h t = 0, (hu) t = 0, (hv) t + (hP 12 ) x = 0, P 11t = 0, P 12t + P 11 v x = 0, h v 2 2 + P 22 2 t + (hP 12 v) x = 0. (1.3.5)
The shock relations for the equation for P h 2 = 0 is a consequence of (1.3.5).

Since the characteristic field corresponding to b-waves is linearly degenerate, the velocity of the shocks is equal to the corresponding eigenvalue. As a consequence, the conservation of the energy is equivalent to the conservation of the 'entropy' through the shock. Thus, the positive definiteness of P is guaranteed even in the presence of shocks. The study of the system in y-direction is analogous. The numerical approach consists in solving first the model in xdirection, and then in y-direction. For each direction, the two subsystems for a-and b-waves are solved successively. A Godunov-type scheme augmented by a correction step is used for this aim. Finally, the source terms are integrated. The accuracy of such an approach is tested on exact 2D solutions describing the flows where the velocity is linear with respect to the space variables, and 1D solutions. "Fingering" is a physical phenomenon, appearing at the unstable interfaces. It consists in the formation of transverse waves whose lengths depend only on the physical parameters of the problem and not on the initial perturbations.

In Chapter 4, we model the so-called SWASI ("Shallow Water Analogue of Shock Instability") experiment by Foglizzo et al. [12], [13], [14] which involves a circular hydraulic jump formed in a convergent radial flow of water. In this experiment, in contrast with the circular hydraulic jump in a kitchen sink (i.e. in a divergent flow), the hydraulic jump was formed in a convergent flow by injecting water inwards from a circular reservoir and evacuating it through a hollow vertical cylinder playing the role of obstacle. The experiments demonstrated the evolution of the radial hydraulic jump flow into a new flow pattern where a rotating singularity was formed at the free surface (see Figure 1.9). Foglizzo et al. also pointed out the analogy between this type of instability of circular hydraulic jump in the water fountain and a supernova formation.

The analogy stems from the fact that the shallow water equations are mathematically equivalent to the gas dynamics equations. Both the hydraulic jump and supernova represent shock waves that connect supercritical (supersonic) flows with subcritical (subsonic) flows. The sur- of physical parameters. We will call such a phenomenon "apple instability" by analogy with the place where the apple stem is attached and the geometrical form of the singularity. The existence of such a point allows us to easily observe pattern's period of rotation. The obtained numerical results are qualitatively similar to those observed experimentally (see Figures 1.9, 1.10). This demonstrates the capability of the model and numerical method to reproduce the shock wave dynamics appearing in multi -dimensional physical problems.

APPENDIX A

Derivation of the model

In this Appendix we develop the governing multi-dimensional shear shallow water equations from Euler's equations. Let v = (u, v, w) T be the spatial velocity field of the fluid, and be the fluid density. The Euler equations are:

           div(v) = 0, v t + (v∇)v + 1 ∇p = g, (1.4.1)
Here p is the fluid pressure, u and v are the horizontal velocities and w is the vertical velocity.

Or, in the projections on the axes: We now introduce a small parameter ε = H L , where H and L are the vertical and horizontal scale lengths, respectively. In the long wave approximation, the vertical momentum equation implies the hydrostatic law:

                                              
∂p ∂z = -g.
The boundary conditions on the free surface z = h(x, y, t) are:

p(z = h) = p 0 , w(z = h) = ∂h ∂t + u(h) ∂h ∂x + v(h) ∂h ∂y .
(1.4.3)

Assuming a flat bottom, we have w(z = 0) = 0. We find immediately p = p 0g(zh). So,

∂p ∂x = g ∂h ∂x , ∂p ∂y = g ∂h ∂y .
And (1.4.2) can be rewritten now in the following way:

                           ∂u ∂x + ∂v ∂y + ∂w ∂z = 0, ∂u ∂t + u ∂u ∂x + v ∂u ∂y + w ∂u ∂z + g ∂h ∂x = 0, ∂v ∂t + u ∂v ∂x + v ∂v ∂y + w ∂v ∂z + g ∂h ∂y = 0. (1.4.4) or,            div 2 (ũ) + ∂w ∂z = 0, ∂ ũ ∂t + (ũ∇ 2 )ũ + w ∂ ũ ∂z + g∇ 2 h = 0, (1.4.5) here ũ = (u, v) T , div 2 (ũ) = ∂u ∂x + ∂v ∂y , ∇ 2 = ∂ ∂x , ∂ ∂y T .
Decomposition of horizontal velocity.

Let us present ũ = U + u 1 , where U is the depth averaged horizontal velocity:

U(t, x, y) = 1 h h 0 ũ(t, x, y, z)dz, (1.4.6)
u 1 is the perturbation of the instantaneous horizontal velocity profile depending of the vertical coordinate z. We note U = (U, V ) T .

Continuity equation.

We integrate continuity equation in (1.4.4) with respect to z from 0 to h. As

h 0 ∂u ∂x dz = ∂ ∂x h 0 udz -u(h) ∂h ∂x = ∂hU ∂x -u(h) ∂h ∂x , h 0 ∂v ∂y dz = ∂ ∂y h 0 vdz -v(h) ∂h ∂y = ∂hV ∂y -v(h) ∂h ∂y , h 0 ∂w ∂z dz = w(h) -w(0) = ∂h ∂t + u(h) ∂h ∂x + v(h) ∂h ∂y ,
we have

∂hU ∂x + ∂hV ∂y -u(h) ∂h ∂x -v(h) ∂h ∂y + ∂h ∂t + u(h) ∂h ∂x + v(h) ∂h ∂y = 0.
Finally, we found the following equation

∂h ∂t + ∂hU ∂x + ∂hV ∂y = 0, or ∂h ∂t + div 2 (hU) = 0. (1.4.7) Momentum equation. Again, ∂ ũ ∂t + (ũ∇ 2 )ũ + w ∂ ũ ∂z + g∇ 2 h = 0.
Using the relation

(ũ∇ 2 )ũ = div 2 (ũ ⊗ ũ) -ũdiv 2 ũ
we rewrite the Euler equation in the form:

∂ ũ ∂t + div 2 (ũ ⊗ ũ) -ũdiv 2 ũ + w ∂ ũ ∂z + g∇ 2 h = 0.
Using now the continuity equation div 2 (ũ) + ∂w ∂z = 0, the last equation becomes

∂ ũ ∂t + div 2 (ũ ⊗ ũ) + ∂wũ ∂z + g∇ 2 h = 0. (1.4.8)
We integrate the last equation with respect to z from 0 to h. We have

h 0 ∂ ũ ∂t dz = ∂ ∂t h 0 ũdz -ũ(h) ∂h ∂t = ∂hU ∂t -ũ(h) ∂h ∂t , (1.4.9) h 0 div 2 (ũ ⊗ ũ)dz = div h 0 ũ ⊗ ũdz -ũ(h) [ũ(h) • ∇ 2 h] . (1.4.10)
Using the decomposition ũ = U + u 1 , one can develop

ũ ⊗ ũ = U ⊗ U + U ⊗ u 1 + u 1 ⊗ U + u 1 ⊗ u 1 . (1.4.11)
We obtain four integrals. Two of them are zero, because the depth averaged velocity does not depend on z by definition, and the depth averaged value of u 1 is zero.

For example,

h 0 U ⊗ u 1 dz = U ⊗ h 0 u 1 dz = 0. (1.4.12)
Definition of the Reynolds stress tensor. We define the Reynolds stress tensor by:

R = h 0 u 1 ⊗ u 1 dz. (1.4.13)
By definition it is symmetric and positive definite tensor which measures the distortion of the instantaneous horizontal velocity profile. The positive definitness of R comes from the Cauchy-Schwarz inequality. Finally,

h 0 div 2 (ũ ⊗ ũ)dz = div 2 (hU ⊗ U) + div 2 R -ũ(h) [ũ(h) • ∇ 2 h] . (1.4.14)
Then, taking into account the boundary conditions, one can find:

h 0 ∂wũ ∂z dz = wũ(h).
The last term is:

h 0 g∇ 2 hdz = g∇ 2 h h 0 dz = gh∇ 2 h = ∇ 2 gh 2 2 . (1.4.15)
The momentum equation becomes:

∂hU ∂t -ũ(h) ∂h ∂t + div 2 (hU ⊗ U) + div 2 R -ũ(h) [ũ(h) • ∇ 2 h] + wũ(h) + ∇ 2 gh 2 2 = 0.
It leads to:

∂hU ∂t + div hU ⊗ U + gh 2 2 I + R = 0, (1.4.16)
where I is the identity matrix.

Evolution equation for the Reynolds stress tensor. As the Reynolds stress tensor R appears in the equation of depth averaged velocity, we need to find the evolution equation of R

to close the problem. Previously, we defined it by

R = h 0 u 1 ⊗ u 1 dz.
We start again with the equation:

∂ ũ ∂t + (ũ∇ 2 )ũ + w ∂ ũ ∂z + g∇ 2 h = 0. ( 1
)ũ} + {(ũ∇ 2 )ũ} ⊗ ũ = ∇ 2 (ũ ⊗ ũ) • ũ.
Here

∇ 2 (ũ ⊗ ũ) • ũ = 2 k=1 ∂ ∂x k (ũ ⊗ ũ)u k .
The last equation can be rewritten in the following form:

∇ 2 (ũ ⊗ ũ) • ũ = div(ũ ⊗ ũ ⊗ ũ) -(ũ ⊗ ũ)div(ũ).
Here

(ũ ⊗ ũ ⊗ ũ)a = (ũ ⊗ ũ)(ũ • a), ∀a, div(ũ ⊗ ũ ⊗ ũ) = 2 k=1 ∂ ∂x k (ũ ⊗ ũu k )
Using the continuity equation div(ũ) + ∂w ∂z = 0, we find the equation

∂ ũ ⊗ ũ ∂t + div(ũ ⊗ ũ ⊗ ũ) + ũ ⊗ ũ ∂w ∂z + w ∂ ũ ⊗ ũ ∂z + ũ ⊗ ∇ 2 (gh) + ∇ 2 (gh) ⊗ ũ = 0. (1.4.18)
Finally:

∂ ũ ⊗ ũ ∂t + div(ũ ⊗ ũ ⊗ ũ) + ∂wũ ⊗ ũ ∂z + ũ ⊗ ∇ 2 (gh) + ∇ 2 (gh) ⊗ ũ = 0. (1.4.19)
We will integrate the last equation with respect to z from 0 to h.

• First term h 0 ∂ ũ ⊗ ũ ∂t dz = ∂ ∂t h 0 ũ ⊗ ũdz -(ũ ⊗ ũ)(z = h) ∂h ∂t .
Using the relation

h 0 ũ ⊗ ũdz = hU ⊗ U + R, we obtain h 0 ∂ ũ ⊗ ũ ∂t dz = ∂hU ⊗ U ∂t + ∂R ∂t -(ũ ⊗ ũ)(z = h) ∂h ∂t .
(1.4.20)

• We calculate now the second term:

h 0 div(ũ ⊗ ũ ⊗ ũ) = div h 0 ũ ⊗ ũ ⊗ ũdz -(ũ ⊗ ũ ⊗ ũ)(z = h)∇ 2 h. (1.4.21)
We have:

ũ ⊗ ũ ⊗ ũ = U ⊗ U ⊗ U + U ⊗ U ⊗ u 1 + U ⊗ u 1 ⊗ U + u 1 ⊗ U ⊗ U + U ⊗ u 1 ⊗ u 1 + +u 1 ⊗ U ⊗ u 1 + u 1 ⊗ u 1 ⊗ U + u 1 ⊗ u 1 ⊗ u 1 .
It means that h 0 ũ ⊗ ũ ⊗ ũdz is divided into eight integrals. Three of them are zero, because the depth averaged velocity does not depend on z by definition, and the depth average value of u 1 is zero. The others non-vanishing integrals are:

h 0 U ⊗ U ⊗ Udz = hU ⊗ U ⊗ U, (1.4.22) h 0 U ⊗ u 1 ⊗ u 1 dz = U ⊗ h 0 u 1 ⊗ u 1 = U ⊗ R, (1.4.23) h 0 u 1 ⊗ u 1 ⊗ Udz = h 0 u 1 ⊗ u 1 dz ⊗ U = R ⊗ U. (1.4.24)
We calculate now the divergence of the integral

h 0 u 1 ⊗ U ⊗ u 1 dz: div h 0 u 1 ⊗ U ⊗ u 1 dz .
Note that ij components of this tensor can be written in the form:

div h 0 u 1 ⊗ U ⊗ u 1 dz ij = ∂ ∂x k h 0 u 1i U j u 1k dz = = ∂ ∂x k U j h 0 u 1i u 1k dz = = ∂ ∂x k (U j R ik ) = = U j ∂R ik ∂x k + ∂U j ∂x k R ik .
(1.4.25)

We can conclude :

div h 0 u 1 ⊗ u 1 ⊗ Udz = div(R) ⊗ U + R ∂U ∂x T , (1.4.26)
where the superscript "T" means "transposed" and div(R)

= 2 j=1 ∂R ij ∂x j . The last integral h 0 u 1 ⊗ u 1 ⊗ u 1 dz
involves the third order fluctuations around the depth averaged velocity. This term is a new tensor of third order, and we need to find his own evolution equation which will involve a new tensor of order 4 and so on. So, the problem is not closed as it usually happens in turbulence theory. However, under the hypothesis of weakly sheared flows [32,17,28], one can prove that this term is negligible with respect to R. Finally, the second term is:

h 0 div(ũ ⊗ ũ ⊗ ũ)dz ≈ div(hU ⊗ U ⊗ U) + div(U ⊗ R) + div(R ⊗ U)+ +div(R) ⊗ U + R ∂U ∂x T -(ũ ⊗ ũ ⊗ ũ)(z = h)∇ 2 h
(1.4.27)

• We calculate now the third term:

h 0 ∂wũ ⊗ ũ ∂z dz = [wũ ⊗ ũ] (z = h). (1.4.28)
• We calculate now the fourth and fifth terms:

h 0 ũ ⊗ ∇ 2 (gh)dz = h 0 ũdz ⊗ ∇ 2 (gh) = hU ⊗ ∇ 2 (gh) = = U ⊗ ∇ 2 gh 2 2 .
(1.4.29)

Analogously, we calculate the fifth term :

h 0 ∇ 2 (gh) ⊗ ũdz = ∇ 2 gh 2 2 ⊗ U. (1.4.30)
• The final expression is:

∂hU ⊗ U ∂t + ∂R ∂t -(ũ ⊗ ũ)(z = h) ∂h ∂t + div (hU ⊗ U ⊗ U) + div(U ⊗ R) + div(R ⊗ U)+ +div(R) ⊗ U + R ∂U ∂x T -(ũ ⊗ ũ ⊗ ũ)(z = h)∇ 2 h+ +w(ũ ⊗ ũ)(z = h) + U ⊗ ∇ 2 gh 2 2 + ∇ 2 gh 2 2 ⊗ U = 0.
(1.4.31)

We can significantly simplify the last equation. First of all,

w(ũ ⊗ ũ)(z = h) -(ũ ⊗ ũ)(z = h) ∂h ∂t -(ũ ⊗ ũ ⊗ ũ)(z = h)∇ 2 h = 0, because (ũ ⊗ ũ ⊗ ũ)(z = h)∇ 2 h = (ũ ⊗ ũ)(z = h)(ũ(z = h) • ∇ 2 h)
and

w = ∂h ∂t + ũ(z = h) • ∇ 2 h.
Then,

div(hU ⊗ U ⊗ U) = ∇ 2 (U ⊗ U) • hU + U ⊗ Udiv(hU).
Also we can develop:

∂hU ⊗ U ∂t = (U ⊗ U) ∂h ∂t + h ∂U ⊗ U ∂t .
Using the continuity equation, we obtain:

U ⊗ U ∂h ∂t + U ⊗ Udiv(hU) = U ⊗ U ∂h ∂t + div(hU) = 0.
We have now

h ∂U ⊗ U ∂t + ∂R ∂t + h∇ 2 (U ⊗ U) • U + div(U ⊗ R) + div(R ⊗ U) + div(R) ⊗ U + R ∂U ∂x T + +U ⊗ ∇ 2 ( gh 2 2 ) + ∇ 2 ( gh 2 2 ) ⊗ U = 0.
This last equation can also be simplified, because

h ∂U ⊗ U ∂t + h∇ 2 (U ⊗ U) • U = h DU ⊗ U Dt ,
where

D Dt = ∂ ∂t + U ∂ ∂x .
Then we can develop:

h DU ⊗ U Dt = h DU Dt ⊗ U + hU ⊗ DU Dt . Because div(U ⊗ R) = U ⊗ div(R) + ∂U ∂x R, we can write h DU Dt + div(R) + ∇ 2 gh 2 2 ⊗ U + U ⊗ h DU Dt + div(R) + ∇ 2 gh 2 2 + + ∂R ∂t + ∂U ∂x T R T + div(R ⊗ U) + R ∂U ∂x T = 0.
Using the momentum equation h DU Dt + div(R) + ∇ 2 gh 2 2 = 0, we find:

∂R ∂t + ∂U ∂x T R T + div(R ⊗ U) + R ∂U ∂x T = 0. Now, using the relations div(R ⊗ U) = Rdiv(U) + ∂R ∂x U and DR Dt = ∂R ∂t + ∂R ∂x U, where ∂R ∂x U = 2 k=1 ∂R ∂x k U k ,
we finally obtain the evolution equation of the Reynolds stress tensor:

DR Dt + Rdiv(U) + ∂U ∂x R T + R ∂U ∂x T = 0. (1.4.32)
Note, that this tensor is symmetric: R T = R. Introducing the tensor P = R h , the last equation can be rewritten in a more compact form:

DP Dt + ∂U ∂x P + P ∂U ∂x T = 0. (1.4.33)
Finally, the full multi-dimensional system of shear shallow water equations over a flat bottom without friction effects can be written in the form:

                           ∂h ∂t + div(hU) = 0, ∂hU ∂t + div hU ⊗ U + gh 2 2 I + hP = 0, DP Dt + ∂U ∂x P + P ∂U ∂x T = 0, D Dt = ∂ ∂t + U • ∇.
(1.4.34)

Introduction

Uniform fluid flows with a free surface down an inclined plane are unstable if the inclination angle is larger than some critical value. The flow then transforms into a system of breaking waves usually called "roll-waves". Brock [4] measured the permanent wave profiles obtained by introducing periodic disturbances at the channel inlet and in different conditions (different slopes and wall roughness). He also studied natural roll-waves propagating in non-periodic manner. He noticed that the roll wave profiles always contained the following three essential parts: first, a steeply sloping wave front, second, a continuous zone where the depth increases progressively, and third, a slowly decreasing zone until reaching a new hydraulic jump (see Figure 1). The second part of the stationary wave profile corresponds to a flow region with a strong vorticity riding the wave front (called also "roller"). This part of the profile is absent when the classical Saint-Venant (SV) equations are used [6]. Recently, a new model of shear shallow water flows was developed in [13], [14], describing, in particular, the roller formation. With such a model, it was shown that the permanent roll-wave profiles correspond to experimental ones with good accuracy. A common feature of solutions obtained in [6] and [13] is that they form a one-parameter family of solutions parametrized by wave length. Roughly speaking, for a given flow discharge, roll-waves of different lengths (different amplitudes) may exist. However, it is quite obvious that very long waves (large amplitude waves) can not exist. So, a natural question arises whether the permanent roll-waves of all lengths are stable in this model. The answer to this question is not at all trivial, and several important studies were performed on this subject. The modulational stability of Dressler waves for SV model was studied, in particular, in [10] and [2]. They have shown that for any Froude number larger than two, there exists an interval (L min , L max ) such that the waves of the length belonging to this interval are stable. The waves of smaller or greater lengths are unstable in the sense that the corresponding modulation equations loss their hyperbolicity. In [1] the spectral stability of roll-waves for viscous SV equations was studied. Stability diagrams at the plane "Froude numberwave length", in particular, were constructed. The spectral instability of small length non-viscous roll-waves for SV equations was also proved in [19]. The new system derived in [13], [14] contains one equation more, so the stability analysis is more difficult. This is why we propose here a numerical approach to the roll-wave stability study.

We study a roll-waves formation through a non-stationary process. If the waves are formed, they are stable, if not, they are not stable. Also, different scenarios of a roll-wave train formation appearing in numerical study could give an important information for practical applications. In particular, the wave "coarsening" produced by irreversible coalescence of non-stationary roll-waves will be studied. The wave "coarsening" is a physical process of the interaction of waves of different lenghts where short waves transfer their energy to longer waves forming later a regular roll-wave train consisting of longer waves. The numerical study of roll-waves formation for the SV model was performed earlier in [20]. The wave coarsening was studied, in particular, in the framework of Shkadov' model where the capillarity and viscosity was taken into account [5]. This paper is divided into five sections. The governing equations of shear shallow-water theory are presented in Section 2. The corresponding finite-volume Godunov type discretization methods are described in Section 3. The numerical results concerning the formation of roll-waves, comparison with Brock's experiments and numerical solutions of SV equations are presented in Section 4. The study of a 2D simplified model of shear shallow water flows is presented in Section 5. In Appendix A the mesh convergence study is presented.

Governing equations

The one-dimensional governing equations for shear shallow water flows on an inclined plane, represent a system of non-linear conservation laws of mass, momentum and energy with source terms describing the balance between gravity and friction [13], [14]:

∂h ∂t + ∂hU ∂x = 0, (1)
∂hU ∂t + ∂ ∂x ( hU 2 + p ) = ĝh -CU |U |, (2)
∂ ∂t (hE) + ∂ ∂x (hU E + pU ) = (ĝh -C e U |U |) U. ( 3
)
Here h is the fluid depth, U is the average velocity, g ′ = g cos θ, ĝ = g sin θ, θ is the inclination angle, C is the Chézy coefficient (friction coefficient) corresponding to the dissipation in the momentum equation, Φ is the enstrophy (squared vorticity) of large eddies formed in the roller, C e = C + C r Φ/ (φ + Φ). The coefficient C r is the dissipation coefficient in the jump roller, and φ is the enstrophy of small eddies developed near the bottom. It is supposed to be constant. The total specific energy E, "internal energy" e, "pressure" p are defined as

E = 1 2 U 2 + e, e = 1 2
( g ′ h + (φ + Φ)h 2 ) , (4)
p = g ′ h 2 2 + (φ + Φ)h 3 . (5)
The second term in the expression of "internal energy" e is the sum of the "turbulent" kinetic energy in the roller ( 1 2 Φh 2 ), and "turbulent" kinetic energy in the boundary layer ( 1 2 φh 2 ). The system (1), (2), (3) admits the following equation for enstrophy Φ:

DΦ Dt = 2 h 3 (C -C e )|U | 3 < 0. (6)
System (1), ( 2) and ( 3) is a time-dependent system of non-linear hyperbolic partial differential equations with characteristic speeds given by U and U ± a s , where a s is the speed of surface waves. It plays the role of "sound speed" in this model :

a s = √ g ′ h + 3(φ + Φ)h 2 . ( 7
)
This system is derived in the framework of shallow water theory, where the ratio of water depth to the wave length is small, and by averaging over the fluid depth. Also an additional hypothesis about the horizontal velocity shear smallness was supposed. The equation (6) means that the enstrophy is decreasing along the trajectories, if C < C e . Since the equations are reminiscent of the Euler equations of compressible flows, the conservation laws imply standard Rankine-Hugoniot relations. At the shock front, the enstrophy is increasing analogously to the entropy increase for the Euler equations of compressible flows. The enstrophy production through shocks corresponds physically to the "roller" formation. Then the enstrophy dissipates over the length of the roller according to (6). When one takes φ = 0, Φ = 0, and C = C e , the system is reduced to the classical Saint-Venant equations.

The system (1) -(3) can be rewritten in conservative form :

U t + F(U) x = S(U), (8)
where the vectors of "conservative" variables U, fluxes F(U), and source term S(U) are:

U =   h hU hE   , F(U) =   hU hU 2 + p hU E + pU   , S(U) =   0 ĝh -C|U |U (ĝh -C e |U |U ) U   .
We will study formation of roll-waves by two methods. The first one is a direct numerical study of roll-waves formation in a long channel from a uniform unstable flow. The flow is perturbed by a wave maker at the inlet of the channel (situated for definiteness at x = 0) with a given constant flow discharge. The second method consists in the roll-wave formation in a "periodic box" : the outlet variables at x = L b (subscript "b" means "box") are the same as the inlet variables at x = 0. The initial velocity is also uniform, but the fluid depth is perturbed in such a way, that the averaged layer depth is equal to the unperturbed one. In this case we conserve the average depth, but not the average discharge.

The relation between these two approaches will be established later. The "periodic box" method is more comfortable from the numerical point of view, because it allows us to study the large time flow behaviour, but the long channel study is physically more understandable.

Initial and boundary conditions for a long channel

We have to impose for (8) initial conditions for x belonging to the interval [0, L] (L is the length of the open channel), and the boundary conditions at x = 0 and x = L. The number of boundary conditions corresponds to the number of characteristics entering the flow domain. The initial conditions are :

h(x, t = 0) = h 0 = const, Φ(x, t = 0) = Φ 0 = 0, U (x, t = 0) = U 0 = √ gh 0 sin θ/C. (9)
This is an exact solution of our system. It is linearly stable, if the generalized Froude number

F r g0 = U 0 a s0
is smaller then 2 (here a s is given by ( 7)), and linearly unstable in the opposite case [13], [12]. Hence, we consider the unstable case to study the roll-waves formation. Since the flow is supercritical (the generalized Froude number is greater than one), three characteristics enter the domain at x = 0. So, we need three boundary conditions at x = 0:

h(x = 0, t) = h 0 (1 + a sin(ωt)), Φ(x = 0, t) = 0, U (x = 0, t) = q h(x = 0, t) . ( 10
)
Here q = h 0 U 0 is the flow discharge, a is a constant perturbation amplitude, ω is a constant frequency. Since the flow is supercritical, we do not need to impose the boundary conditions for x = L (there is no characteristics entering the flow domain). The boundary condition for h allows us to model qualitatively a wave maker movement that was used by Brock at the channel inlet to accelerate the formation of roll-waves.

Initial and boundary conditions for a periodic box

We take periodic boundary conditions for the box of length L b :

U(0, t) = U(L b , t), t ≥ 0.
Consider the stationary solution (h 0 , U 0 , Φ 0 = 0) where the equilibrium velocity U 0 is given by (9). Suppose that the generalized Froude number is greater than two, i.e. the flow is linearly unstable. Assuming sinusoidal initial perturbation of amplitude a of the free surface :

h(x, t)| t=0 = h 0 ( 1 + a sin ( 2πx L b )) , (11)
and taking other variables constant :

U | t=0 = U 0 , Φ| t=0 = 0,
we allow the flow evolve in time.

Numerical scheme

We use here conservative, finite volume Godunov type scheme on a fixed grid. It requires the solution of the Riemann problem at every cell boundary at each time step [8], [9], [18], [15]. The MUSCL-Hancock extension of the Godunov method is used with the MinMod limiter for the depth, the velocity and the "pressure".

Hyperbolic step

Let us consider a fixed grid of size ∆x = x i+1/2x i-1/2 , the time increment is defined as ∆t = t n+1t n that must respect the Courant-Friederichs-Lewy's (CFL) condition. The value of the CFL number is taken to be 0.8. The number of cells is always given for a converged solution.

In particular, for such a solution every single roll-wave is represented at least by 100 points. The discrete values of the vector-function U(x, t) at (x i , t n ) will be denoted by

U n i ≡ U(x i , t n ).
The first step (hyperbolic one) consists in computing the source term-free system:

x i -1/2 t t+∆ t x ∆ x x i+1/2 i+1 i i -1 Riemann problem
U t + F x = 0, (12)
with the initial condition for the complete problem U(x, t n ) = U n . Integrating in space and time

[x i-1/2 , x i+1/2 ] × [t n , t n+1
] the conservation laws (12), one obtains a conservative finite volume Godunov scheme on a fixed grid :

Ūn+1 i = U n i - ∆t ∆x ( F * ,n i+1/2 -F * ,n i-1/2 ) , (13)
where F * ,n i+1/2 and F * ,n i-1/2 are numerical fluxes. They are constant across interfaces between cells during the time step. For computing the fluxes F * ,n i+1/2 and F * ,n i-1/2 , we solve the Riemann problems between cells i, i + 1 and i -1, i, respectively (see Figure 2). The HLLC Riemann solver is used for this aim [18].

Integrating the source terms

The last step is to integrate the differential equation

dU dt = S(U), (14)
with the initial condition U| t=0 = Ūn+1 given by (13). To integrate the source term, we use Strang splitting [16] :

U n+1 = S 1 2 ∆t M ∆t S 1 2 ∆t U n ,
where M is the MUSCL-Hancock operator, S is the 4 th order explicit Runge -Kutta numerical operator [11]. [4] are given. The flow discharge q 0 = h 0 U 0 can be calculated by using (9) and the value of h 0 .

TEST h 0 [m] θ [rad] C r C φ [s -2 ] λ [m] ω [rad/s] T [s] F
The amplitude value a is not given in [4]. The amplitude increase accelerates the formation of a roll-wave train but does not influence its final form. The values C r and C correspond to those used in [13].

Numerical simulation of roll waves

Two sets of parameters are taken for the numerical study of roll-wave formation in a long channel (see Table 1). In Case 2 the inclination angle is greater than that in Case 1. Consequently, the generalized Froude number and wavelength are greater also in Case 2. The friction coefficient C in both cases is approximately the same. In Table 1, C r is the dissipation coefficient in the jump roller, λ is wave length, ω is wave frequency, T is period of perturbation, h 0 is initial depth, a is perturbation amplitude. They will be used for the flow modelling in both a long channel and a periodic box. The length of the channel varied between L = 40 [m] and L = 600 [m]. For the periodic box we will use L b of 1.3 m or 1.8 m, corresponding to Cases 1 and 2. We will also perform computations for a periodic box having a length multiple of L b , with the perturbation of the same type:

h(x, t = 0) = h 0 (1 + a sin(2πx/(nL b ))), n = 2, 3, ... (15)

Long channel

In this section we will focus on the evolution of a roll-wave train moving down in a rectangular inclined open channel of constant slope. For a prescribed average discharge, the Dressler theory gives a one-parameter family of solutions, parametrised by wavelength ranging from zero to infinity. It was shown analytically for the SV equations [10], [2], [1] that not all waves are stable. We will obtain the analogous effect for the new (1)-( 3) model. We perturb a uniform flow at the channel inlet (the perturbation amplitude is 5% of the uniform depth). First, we take the generalized Froude number smaller then two. As it was expected from the theory [13], [12], the perturbations move downstream the channel with decreasing amplitude until they completely disappear. Hence, the roll-wave train is not formed because the uniform flow is stable (see Figure 3). When we take the generalized Froude number larger than two, the perturbations are amplified with time and form rapidly a roll-wave train. All waves have the same amplitude and propagation speed (Figures 4,5). Each wave profile includes shock and a strongly sheared region after the shock (roller), where the depth continues to increase. The roller enstrophy increases through shocks and then rapidly decreases over the roller length (Figures 4,5).

The roll-waves are completely formed at the distance of about 15 -20 m from the channel inlet (in both cases). The length of a single roll wave is about 1.3 m (Case 1 ) and 1.8 m ( Case 2).

In Figure 6 we compare a single roll-wave obtained for model ( 1)-( 3) with Brock's measurements and numerical solutions of the SV equations. The profile of a single wave corresponds to that found in [12], [13]. It contains a rapidly varying part (a jump), a gradual monotonic increase of the wave profile, and a decreasing part. A good agreement between stationary solution of the new model and the experiments is observed. We note that the amplitude of the jump is always smaller than the double of the depth before the jump. This is in accordance with the limit depth ratio coming from the Rankine-Hugoniot relations (see [12], [13]). One can observe an important difference between the SV solutions and experimental observations. Saint-Venant equations is just a regular (non-modulated) roll-wave train (Figure 9). The wave length of a single wave forming such a roll-wave train is also about 5 m, but the amplitude is larger than in the case of shear waves. Thus, the long-time behaviour of the model ( 1) -( 3) results in a strongly modulated roll-wave train, while the Saint-Venant equations produce a regular wave train.

Coarsening

If the roll-wave train contains waves of different lengths, they begin to interact and coalesce. The short waves transfer their energy to long waves, and finally a roll-wave train of larger wavelength appears. This physical process is called "coarsening", or, as it was proposed by Brock [3], "growth by overtaking". To observe the coarsening, we perturb at the inlet the unstable uniform flow with two different frequencies ω 1 , ω 2 and amplitudes a 1 , a 2 :

h(0, t) = h 0 (1 + a 1 sin(ω 1 t) + a 2 sin(ω 2 t)).
A larger wavelength corresponds to a smaller frequency. Qualitatively, the scenario of "coarsening" does not depend too much on exact values of ω i and amplitudes a i , i = 1, 2. For definiteness, we take ω 1 = 6.06 [rad/s], ω 2 = 4.55 [rad/s]. This corresponds to the wavelengths of 1.6 m and of 3.55 m, respectively. The initial depth is h 0 = 0.01 m, the amplitudes are a 2 = 2a 1 = 0.1. Two different development stages can be observed. The first one is the transition phase where a strong interaction between waves occurs, and the wave length starts to increase (Figure 10, the time instants 50 s and 100 s). During the second stage, a roll-wave train is formed with

Periodic box

We partially repeat here the above numerical experiments for the periodic box. We take the initial sinusoidal perturbation of the free surface in the form given by (11). The amplitude of perturbation a is always 0.05. We will see how the perturbation evolves in time. As in the case of a long channel, one can show that the roll wave is not formed for the generalized Froude number smaller than two (Figure 12). For the generalized Froude number greater than two, one can see that just a single roll wave was formed (see Figures 13 for Case 3) and for the SV equations. A good agreement between model ( 1)-( 3) and Brock's experiments can also be observed. As it was expected, the SV equations exaggerate both the wave amplitude and profile.

Let us show that final numerical solution does not depend on the form of the initial perturbation. We consider two different initial perturbations:

h(x, t)| t=0 = h 0 (1 + a sin(8πx/L b )) , 0 ≤ x ≤ L b , and h(x, t)| t=0 = h 0 (1 + 2a sin(4πx/L b ) + a sin (2πx/L b )) , 0 ≤ x ≤ L b .
For definiteness, we take L b = 1.3 m corresponding to Case 1. At the time moment 100 s, we obtain the same single wave of length L b = 1.3 m (see Figure 14). In Figure 15 we present the

t = 0 s, F r g0 < 2 t = 100 s, F r g0 < 2
TEST h 0 [m] < h > [m] error CASE 1 0.00798 0.007917 1 % CASE 2 0.00533 0.0050327 5.5 % Table 2:
The average depth of a single steady wave in a roll-wave train is shown for Case 1 and Case 2 in a long channel at the time instant 100 s. Compared to the periodic box, where the average value must be exactly h 0 , the average value for a single steady wave in a roll-wave train is a little bit smaller. However, as it can be seen from the Table, this error is small: 1% for Case 1, and 5.5% for Case 2. The error is defined as (h 0 -< h >)/h 0 .

average depth and the average discharge per unit width in the periodic box (Case 1), calculated in the following way: The average depth is constant in time, but not the flow discharge. However, when a single wave is formed, the wave velocity becomes constant, so the average discharge becomes also constant. So, even if, a priori, the periodic box and long channel problem are not equivalent, they become equivalent in the long time limit. This allows us to use the periodic box as a mathematical tool for the study of qualitative properties of roll-waves. Below we use the periodic box for the instability study of roll-waves.

< h > (t) = 1 L b ∫ L b 0 h(x, t)dx, < q > (t) = 1 L b ∫ L b 0 h(x, t)U (x,

Wave instability in a periodic box

A periodic box of different lengths multiple of L b is considered. In particular, Case 1 is considered with L b = 1.3 m. The numerical results are shown at the same time instant of 1000 s. The computations show that a single stable wave is always formed until the box length 8 L b (Figure 16 on the left). The wave is morphologically stable : it has always a steeply sloping wave front, then a continuous zone where the depth increases progressively, and finally, a slowly decreasing monotonic zone until reaching again the wave front. However, its amplitude changes slowly in time in the same way as in the long channel where the modulations of the roll-wave train appear. Beginning from the length 9 L b this single wave becomes very non-stationary. It finally breaks into two waves for the box of length 14 L b . The minimal length of periodic box for which a single roll wave is stable, was not found.

Analogous bifurcation results were obtained for the SV equations (Figure 16 on the right). As in the case of a long channel the single stable waves for the SV equations are steady : the amplitude modulations are not present. The critical "bifurcation" lengths are here 11 L b (the existence of a stable single wave) and 12 L b (the bifurcation of a single wave into two waves).

8 L b , RG 11 L b , SV 14 L b , RG 12 L b , SV

2D simplified model

Multi-dimensional model of shear shallow water flows is much more complicated. Its structure is reminiscent of equations of compressible turbulent flows , [7], [17], [12]. Even if it is hyperbolic, it is not conservative. To avoid the treatment of non-conservative equations, we consider here a simplified model of shear flows where the Reynolds stress tensor is replaced by a spherical one. Such an approximation corresponds to the consideration of isotropic turbulence, often used in applications. The governing equations for the general case of a space varying topography are :

h t + div(hU) = 0, ( 16
) (hU) t + div (hU ⊗ U + pI) = -gh∇b -CU|U|, (17)
(hE) t + div (hUE + pU) = -C e |U| 3 . (18)
Here x = (x, y) T , z = b(x) is the bottom topography, the free surface is at

z = b(x) + h(t, x), h(t, x) is the fluid layer thickness, U = (U, V ) T is the average velocity, |U| = √ U 2 + V 2 ,
and Ψ is the enstrophy (squared vorticity). The total energy E and the "pressure" p are :

E = |U| 2 2 + gh 2 + gb + Ψh 2 2 , p = gh 2 2 + Ψh 3 .
As in 1D case, the system implies the equation for enstrophy:

h 3 2 DΨ Dt = (C -C e )|U| 3 , where D Dt = ∂ ∂t + U ∂ ∂x + V ∂ ∂y .
We decompose

Ψ = Φ + φ,
where φ is a given small constant (describing the intensity of the vortexes in the boundary layer near bottom), and Φ is the large scale enstrophy. The equations are written in a reference system where the direction of gravity is orthogonal to the (x, y)-plane. This a little bit different from 1D case, where x-direction corresponds to the tangent vector to the inclined plane. However, for a small inclination angle, the equations are equivalent. We take the bottom topography in the form b(x) = -x tan θ that corresponds to an inclined plane. It is supposed that θ = const > 0. The system of conservation laws is reminiscent of the Euler equations of compressible fluids with a right-hand side. Due to the following identity :

h D Dt gb = gh (b t + U b x + V b y ) = ghU b x = ghU tan θ,
the energy equation can be simplified to :

( h Ẽ) t + div { hU Ẽ + ( gh 2 2 + (Φ + φ)h 3 ) U } = -ghU∇b -C e |U| 3 , (19)
where the modified total energy is :

Ẽ = |U| 2 2 + gh 2 + (Φ + φ)h 2 2 .
We use a multi-dimensional unsplit extension of the one-dimensional conservative Godunov scheme presented earlier. Two-dimensional system ( 16), ( 17) and ( 19) can be rewritten in the following form :

U t + F(U) x + G(U) y = S(U), (20)
where the vectors of the "conservative" variables U, fluxes F(U) and G(U), and source term S(U) are:

U =     h hU hV h Ẽ     , F(U) =     hU hU 2 + p hU V hU Ẽ + pU     , G(U) =     hV hU V hV 2 + p hV Ẽ + pV     , S(U) =     0 g tan θh -CU (U 2 + V 2 ) 1/2 -CV (U 2 + V 2 ) 1/2 -C e (U 2 + V 2 ) 3/2 + g tan θhU,     , p = gh 2 2 + (Φ + φ)h 3 .
Consider a typical finite volume cell I i,j = ∆x × ∆y. A finite volume scheme for the homogeneous part of system (20) is :

U n+1 i,j = U n i,j + ∆t ∆x [ F * ,n i-1/2,j -F * ,n i+1/2,j ] + ∆t ∆y [ G * ,n i,j-1/2 -G * ,n i,j+1/2
] .

Across each intercell boundary in each direction we solve the Riemann problem to find the numerical fluxes F * ,n i±1/2,j , G * ,n i,j±1/2 . The integration of the source term is done in the same way as in 1D case. We use the two-dimensional model for a periodic box of the length L x in x-direction, and L y in y-direction, with the initial conditions in the following form:

h| t=0 = h 0 {1 + a sin(8πx/L x ) + a sin(8πy/L y )} , U | t=0 = U 0 , V | t=0 = 0, Φ| t=0 = 0. (21)
Here L x and L y are the length and the width of the periodic box, respectively. In Figure 17 the variation of the initial depth is shown for Case 1. In the direction y we use the wall boundary conditions. Single permanent waves which are formed in the time instant approximately 150 s (see Figure 18) correspond to those obtained in one-dimensional case (see Figure 13). Thus, the limit solution is not sensible to the initial transverse perturbations.

Conclusion

The numerical study of roll-waves that develop from a uniform unstable flow down an inclined rectangular channel is performed. The formation of the roll-waves is studied by two different approaches. In the first approach, the roll-waves were generated in a long channel. The free surface was perturbed by a wave maker at the channel inlet. The average discharge was fixed.

In the second approach, the roll-waves were produced in a "periodic box" with a uniform flow velocity. The average depth of a perturbed free surface was the same as in the long channel.

The generation of a periodic roll-wave train was studied for a long channel for two sets of experimental parameters (noted as Case 1 and Case 2) corresponding to Brock's experiments [3], [4]. In both cases, the free surface profile for model ( 1)-( 3) was found in good agreement with the experimental results, while the SV equations give a wrong wave form and amplitude. The model was also tested in the case where the perturbation frequency was much lower than the experimental one, so long waves were generated at the channel inlet. The amplitude and enstrophy of the corresponding roll-wave train are strongly modulated, while the SV equations give for the same set of parameters a steady roll-wave train. In the case where the waves of two different lengths were initially generated, the coarsening was observed. The coarsening was always accompanied by a modulation phenomenon.

The formation of a single wave composing a roll-wave train in a "periodic box" was studied for the same sets of experimental parameters. The free surface profile for the model ( 1)-( 3) was found in very good agreement with the experimental results. This allows us to justify the use of the "periodic box" as a simple mathematical tool for a qualitative study of roll-waves stability. In particular, we studied the stability of a single steady wave by taking its length as a multiple of L b = 1.3 m corresponding to Case 1 of Brock's experiments. It was shown that the wave becomes morphologically unstable after some critical wave length. This fact was also validated for the SV model.

Finally, we proved that a single steady wave corresponding to Case 1 is stable under multidimensional perturbations in the framework of a model which represents a simplification of a general multi-D model of shear shallow water flows. However, in the case of the full non-conservative model the result could be different.

Our future work will certainly be oriented to the development of a full multi-dimensional model capable to describe more complex multi-dimensional and non-stationary physical phenomena.

APPENDIX A

Mesh convergence study For the mesh convergence study, we construct a sequence of approximate solutions and then we calculate the error compared with the solution on a refined mesh (16000 cells). We took here a periodic box of 1.3 m long to calculate error by the formula :

error = √ Σ N i=1 Σ k j=1 [H(i)-H(k(i-1)+j)] 2 ∆x L b L b .
Here N is the number of cells for the sequence of approximate solutions (100, 200, 400, 500, 1000, 2000, 4000 and 8000), k is obtained using the formula : N max = kN , where N max = 16000 corresponds to the refined mesh, H(i) is the fluid depth in the cell i, H(k(i -1) + j) is the value of the fluid depth in the cell k(i -1) + j for the refined numerical solution. In Figure 19 one can see the convergence of MUSCL-Hancock scheme (MinMod limiter) with the line slope corresponding to approximately 1.4.

-

Introduction

The Saint-Venant equations [5] are of great importance both in hydraulic and oceanographic applications. The reason is that they are simpler compared to the n-dimensional Euler equations with a free surface. Indeed, the Saint-Venant equations describe the evolution of the fluid depth and depth averaged velocities defined in a fixed n -1 dimensional domain. They form a hyperbolic system of equations in conservative form for which standard numerical methods can be applied [23,29,[START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]. The derivation of Saint-Venant equations is based on the smallness of the parameter ε = H/L where H and L are the vertical and horizontal scale lengths, respectively, and the hypothesis that the dependence of the horizontal velocity on the vertical coordinate is very weak, i. e. the fluid flow is almost potential (the flow is not sheared). Since the shear effects are completely neglected in the model derivation, the Saint-Venant equations are not able to describe neither the formation of large scale eddies ('roller') appearing in the hydraulic jumps near the free surface nor the form of the hydraulic jump. A natural extension of the Saint-Venant equations are integro-differential Benney's equations [7] describing shear shallow water flows of uniform density. The notion of hyperbolicity of the Benney equations was introduced in [START_REF] Teshukov | On hyperbolicity of long-wave equations[END_REF][START_REF] Teshukov | On Cauchy problem for long-wave equations[END_REF]30,9]. Benney's equations admit, in particular, a linear vertical shear with constant vorticity as an exact solution. A numerical approach based on such a linear approximation was developed in [START_REF] Teshukov | Analytical and numerical solutions of the shallow water equations for 2-D rotational flows[END_REF] in 1D case. To the best of our knowledge, 2D computations of the Benney equations are absent in the literature.

For shear flows with varying in space and time vorticity an intermediate model was recently proposed where the governing equations are obtained by depth averaging of Euler equations without assuming potential flow [START_REF] Teshukov | Gas dynamic analogy for vortex free-boundary flows[END_REF][START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF][START_REF] Richard | Elaboration d'un modèle d'écoulements turbulents en faible profondeur: application au ressaut hydraulique et aux trains de rouleaux[END_REF]. The hypothesis of smallness of the horizontal vorticity (the hypothesis of weakly sheared flows) allows us to keep the second order depth averaged correlations in the governing equations but neglect the third order correlations, and thus to close the governing system. This approach was further extended to the case where the third-order correlations are taken into account [14].

This intermediate multi-dimensional model is reminiscent of the classical Reynolds averaged Euler equations for the compressible barotropic turbulent flows [22]. The model complemented by friction terms was used for the study of 1D travelling waves down inclined plane (roll waves) and hydraulic jumps. A strong physical adequacy of the model with the experimental observations was found [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF][START_REF] Richard | Elaboration d'un modèle d'écoulements turbulents en faible profondeur: application au ressaut hydraulique et aux trains de rouleaux[END_REF].

The multi-dimensional case is much more challenging. For the flows over a flat bottom without friction effects, the system can be written in the form [START_REF] Teshukov | Gas dynamic analogy for vortex free-boundary flows[END_REF][START_REF] Richard | Elaboration d'un modèle d'écoulements turbulents en faible profondeur: application au ressaut hydraulique et aux trains de rouleaux[END_REF]21] :

h t + div(hu) = 0, (1)
(hu) t + div hu ⊗ u + gh 2 2 I + hP = 0, DP Dt + ∂u ∂x P + P ∂u ∂x T = 0, D Dt = ∂ ∂t + u • ∇.
Here t is the time, x = (x, y) T are the Cartesian coordinates, u = (u, v) T is the depth averaged horizontal velocity, h is the fluid depth, g is the gravity, and P is the stress tensor which measures the distortion of the instantaneous horizontal velocity profile ũ(t, x, y, z) depending of the vertical coordinate z. The sign ⊗ means the tensor product, and I is the identity tensor. The definitions of u(t, x) and P are as follows:

u(t, x) = 1 h h 0 ũ(t, x, y, z)dz, P = 1 h h 0 (ũ -u) ⊗ (ũ -u) dz.
The tensor P is symmetric and positive definite. The positive definiteness of P is a consequence of the Cauchy-Schwarz inequality.

A striking mathematical analogy with the Reynolds averaging equations of barotropic compressible turbulent flows [32,[START_REF] Pope | Turbulent flows[END_REF][START_REF] Wilcox | Turbulence Modeling for CFD[END_REF]21] allows us to call R = hP the Reynolds stress tensor, and P the reduced Reynolds stress tensor. For mathematical reasons, the choice of P is more convenient than that of R. For simplicity, both R and P will be further refereed to as the Reynolds stress tensor. 1Equations (1) admit the energy conservation law:

∂ ∂t h 1 2 |u| 2 + e i + e T + div hu 1 2 |u| 2 + e i + e T + gh 2 2 I + hP u = 0, (2)
where

e T = 1 2 trP, e i = 1 2 gh,
and an additional conservation law:

∂hΨ ∂t + div (huΨ) = 0, Ψ = Det (P) h 2 . ( 3
)
The variable Ψ will be referred to as 'entropy' (mathematical) because this quantity is transported along the mean flow in the same way as the true entropy for the Euler equations of compressible fluids. Also, we will see that this quantity will increase across the shocks in analogy with the conventional entropy (see [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF][START_REF] Richard | Elaboration d'un modèle d'écoulements turbulents en faible profondeur: application au ressaut hydraulique et aux trains de rouleaux[END_REF] for the 1 D study of roll waves and hydraulic jumps).

The system is hyperbolic but not in conservative form. The hyperbolicity was established, for example, in [8], for the equations of compressible turbulent flows, generalizing system (1). However, the fact of the non-conservativity was just a 'feeling' which was not rigorously proved in the literature. In Appendix A we establish that system (1) admits only five conservation laws written above : conservation of mass, momentum, energy and 'entropy'. Since the number of scalar unknowns is six (h, u and three independent components of P), the system is not in conservative form. The definition and computation of discontinuous solutions for non-conservative hyperbolic equations is a challenging problem (see examples of non-conservative systems in compressible turbulence [8,2], multi-layer shallow water flows [START_REF] Ovsyannikov | Two-layer shallow water model[END_REF]4,33,6,1,31], multi-phase fluid flows [3,25,[START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF][START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF][START_REF] Saurel | Shock jump relations for multiphase mixtures with stiff mechanical relaxation[END_REF]20,15], solid-fluid systems [18,34]).

Essentially, four approaches are commonly used for numerical solving of non-conservative systems of equations. The most classical one is based on the definition of non-conservative products proposed by Volpert (Volpert's path) [16]. The second one is the formulation of an augmented system of 'Rankine-Hugoniot relations' through the study of travelling wave solutions of an extended system of equations approximating a given system (formulation of kinetic relations) [START_REF] Truskinovsky | Kinks versus shocks[END_REF]27]. The third one is based on the relaxation technique : the governing system of equations is approximated by a new hyperbolic system where all eigenfields are linearly degenerate in the sense of Lax [10,15]. Finally, the additional relations can be formulated from the compatibility between theoretical and experimental results [17,[START_REF] Saurel | Shock jump relations for multiphase mixtures with stiff mechanical relaxation[END_REF]. Excepting the first approach (more formal and hence less precise), all the approaches mentioned above are not universal : they are usually specific to the model under study.

In [19], a new splitting approach for the modelling of 3D isotropic hyperelastic materials was proposed. The model admits three types of sonic waves (one longitudinal and two transverse) and contact characteristics. The system was split into several subsystems each of which contained only one type of sonic waves (only longitudinal or transverse). Each subsystem was hyperbolic and admitted a weak formulation. Such an approach was further extended to the modelling of multi-solid materials [34]. The splitting procedure allowed us not only to define the non-conservative products, but also to increase the precision and robustness of the numerical method. Such a philosophy will also be developed here.

In section 2, the hyperbolicity of system (1) is established. The Rankine-Hugioniot relations compatible to the positive definiteness of the Reynolds stress tensor are proposed in section 3. Dissipation terms are introduced in section 4. The splitting procedure and its numerical realisation are presented in sections 5, 6, 7. The numerical results are presented in section 8. Technical details can be found in Appendices A and B.

Hyperbolicity study

The hyperbolicity study is analogous to that given in [8] for the compressible turbulent flows. For completeness, we briefly describe the main results. Denoting the components of P by P ij , i, j = 1, 2, one can rewrite system (1) in Cartesian coordinates in the form :

h t + uh x + vh y + hu x + hv y = 0, (4)
P 22t + uP 22x + vP 22y + 2P 12 v x + 2P 22 v y = 0.
Or, in matrix form :

∂W ∂t + A ∂W ∂x + B ∂W ∂y = 0,
where

W =                   h u v P 11 P 12 P 22                   , A =                      u h 0 0 0 gh + P 11 h u 0 1 0 P 12 h 0 u 0 1 0 2P 11 0 u 0 0 P 12 P 11 0 u 0 0 2P 12 0 0 u                      , (5)
B =                      v 0 h 0 0 0 P 12 h v 0 0 1 0 gh + P 22 h 0 v 0 0 1 0 2P 12 0 v 0 0 0 P 22 P 12 0 v 0 0 0 2P 22 0 0 v                      . (6)
The characteristic surfaces S(t, x, y) = 0 for (4) satisfy the relations :

det I ∂S ∂t + A ∂S ∂x + B ∂S ∂y = 0. It implies : χ = 0, (7)
χ = ± (∇S) T P∇S, (8)
χ = ± gh|∇S| 2 + 3 (∇S) T P∇S, (9)
where

χ = ∂S ∂t + u • ∇S, ∇S = ∂S ∂x , ∂S ∂y T .
The characteristic value ( 7) is double. The characteristic values ( 8) and ( 9) are simple and real since P is positive definite. One can prove that for the double root one has two left linearly independent eigenvectors of the corresponding characteristic matrix. Hence, the equations are hyperbolic.

To understand the structure of the eigenfields, consider the governing equations in x-direction. The eigenvalues of the matrix A are :

λ 1,2 = u, λ 3,4 = u ± b, b = P 11 , λ 5,6 = u ± a, a = gh + 3P 11 . ( 10
)
For the multiple eigenvalue λ 1,2 = u we have two linearly independent right eigenvectors r i , i = 1, 2, of A : r 1 = (0, 0, 0, 0, 0, 1)

T , ∇ W λ 1 • r 1 = 0, r 2 = (-h, 0, 0, gh + P 11 , P 12 , 0) T , ∇ W λ 2 • r 2 = 0. ( 11
)
Here and further, ∇ W means the gradient operator with respect to W. For the eigenvalue λ 3 = u + b one has:

r 3 = (0, 0, b, 0, b 2 , 2P 12 ) T , ∇ W λ 3 • r 3 = 0. ( 12
)
For the eigenvalue λ 4 = ub one has:

r 4 = (0, 0, -b, 0, b 2 , 2P 12 ) T , ∇ W λ 4 • r 4 = 0. ( 13
)
For the eigenvalue λ 5 = u + a one has :

r 5 = h, a, 2aP 12 a 2 -b 2 , 2b 2 , a 2 + b 2 a 2 -b 2 P 12 , 4P 2 12 a 2 -b 2 T , ∇ W λ 5 • r 5 = 3 2a a 2 + b 2 > 0. ( 14
)
For the eigenvalue λ 6 = ua one has:

r 6 = h, -a, - 2aP 12 a 2 -b 2 , 2b 2 , a 2 + b 2 a 2 -b 2 P 12 , 4P 2 12 a 2 -b 2 T , ∇ W λ 6 • r 6 = - 3 2a a 2 + b 2 < 0. (15)
Thus, the fields corresponding to the eigenvalues λ 1,2,3,4 are linear degenerate in the sense of Lax, while the fields λ 5,6 = u ± a are genuinely non-linear. The family λ 3,4 = u ± b (further on referred to as b-waves) is reminiscent of the shear waves in hyperelasticity, while the family λ 5,6 = u ± a (further on referred to as a-waves) is reminiscent of the longitudinal waves. The analogy with shear waves in hyperelasticity and b-waves was also noticed in [START_REF] Troshkin | On wave properties of an incompressible turbulent fluid[END_REF] in the case of incompressible Reynoldsaveraged Euler equations. An interesting analogy between the linearised equations of incompressible turbulence and Maxwell's equations was also underlined there.

Rankine-Hugoniot relations

Since the number of scalar conservation laws is only five (mass, momentum, energy and 'entropy' ), while the number of unknowns is six (h, u and three independent components of P), the system is not in conservative form (see Appendix A for a proof). For discontinuous solutions, the Rankine-Hugoniot relations coming from the mass, momentum and energy equations are :

[h(u • n -D n )] = 0, hu(u • n -D n ) + gh 2 2 n + hPn = 0, h(u • n -D n ) 1 2 |u| 2 + e i + e T + n T gh 2 2 I + hP u = 0.
Here for any function f we denote [f ] = f +f -, where f + , f -are the right and the left limit values of f at the discontinuity surface, D n is the normal velocity of the surface, and n is the normal unit vector to the surface. We denote also by s the tangent unit vector to the surface such that (n, s) form a Cartesian basis. Using

u = (u • n)n + (u • s)s,
one obtains from the momentum equation the following two scalar relations :

h(u • n -D n ) 2 + gh 2 2 + hn T P n = 0, [h(u • n -D n )(u • s) + hs T P n] = 0.
As usually, we distinguish two types of discontinuities: contact discontinuities (interfaces) where u • n -D n = 0, and shock waves where u • n -D n = 0.

Contact discontinuities

We consider first the interfaces u • n -D n = 0. The momentum equation is equivalent to :

gh 2 2 + hn T Pn = 0, [hs T Pn] = 0.
The energy equation implies :

[h(u • s)s T Pn] = 0.
We need to distinguish two types of contact discontinuities. The first type is determined by the condition that at each side of the contact discontinuity the tangential component of the stress vector hs T Pn vanishes. Then, a priori, the jump of the tangential velocity can be arbitrary:

[u • s] = 0. (16)
So, the sliding is admitted as in the case of contact discontinuities for the Euler equations of compressible fluids.

For the second type of contact discontinuity where the tangential component of the stress vector hs T Pn is continuous, but not necessarily vanishing, the sliding is forbidden :

[u • s] = 0. ( 17
)
So, the full velocity vector should be continuous. This kind of interfaces is not admitted by the Euler equations of compressible fluids.

Shocks

Consider now the discontinuity interfaces where u • n-D n = 0. They are called shocks. The Rankine-Hugoniot relations coming from the mass, momentum (two scalar relations) and energy equations are

: [h(u • n -D n )] = 0, h(u • n -D n ) 2 + gh 2 2 + hn T P n = 0, h(u • n -D n )(u • s) + hs T P n = 0, h(u • n -D n ) 1 2 |u| 2 + e i + e T + (u • n) gh 2 2 + hn T Pn + (u • s)hs T Pn = 0.
Contrary to the Euler equations of compressible fluids, the sliding along the shocks surfaces is not forbidden. This system of Rankine-Hugoniot relations is obviously not sufficient to describe the full structure of strong discontinuities. As it is proved in Appendix A all linearly independent conservation laws admitted by the system (1) of the form :

∂f ∂t + ∂g 1 ∂x + ∂g 2 ∂y = 0,
are only those of mass, momentum, energy and hΨ. Here f, g 1 , g 2 are functions depending on h, u and P. Thus, there is no hope to define a 'classical' weak solution to (1). Moreover, the shock relation coming from the equation for Ψ implies the continuity of Ψ across the shocks. However, this condition is not physically acceptable. Indeed, for 1D flows the system of Rankine-Hugoniot relations is closed. Moreover, the 1D system admits an analogue of Ψ having a clear physical meaning. This quantity is associated with the enstrophy (squared vorticity) and plays the role of 'mathematical entropy' : it increases across the shock. For example, for hydraulic jumps the enstrophy increase corresponds to the vorticity creation : a roller appears at the forward slope of the hydraulic jump [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF]. Note that in [20] a simplified model of multi-phase flows was proposed, and a quantity analogous to Ψ and characterising the micro-energy of interaction between phases was introduced. This quantity was also increasing across the shock.

The non-conservative nature of the multi-dimensional equations of shear shallow water flows represents an enormous difficulty from the mathematical and numerical point of view. Moreover, when the shocks appear, we also should be aware about positive definiteness of P : this property should be guaranteed for the weak solutions. If the solution is smooth, this property is easy to establish. Indeed, the equation for Ψ can be integrated in the Lagrangian coordinates X related to the mean flow. One can write in the Lagrangian coordinates :

Ψ(t, X) = Ψ(0, X), or λ 1 λ 2 h 2 = λ 10 λ 20 h 2 0 ,
where λ i , i = 1, 2 are eigenvalues of P, and the index 0 corresponds to the initial state. It is then clear that if P is initially positive, it will be positive for any time, if the solution is smooth. We will also establish a set of additional Rankine-Hugoniot relations that guarantees the positive definiteness of P even in the presence of shocks. Dissipative terms compatible with the positive definiteness of P will also be introduced into the governing equations.

'Entropy' increase across the shocks : 1D case study

This case was already studied in [START_REF] Teshukov | Gas dynamic analogy for vortex free-boundary flows[END_REF][START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF]. Here, we present it for completeness. The governing equations are :

∂h ∂t + ∂hu ∂x = 0, ( 18
)
∂hu ∂t + ∂ ∂x hu 2 + gh 2 2 + Φh 3 = 0, (19)
∂ ∂t h 2 u 2 + gh + Φh 2 + ∂ ∂x hu u 2 2 + gh + 3 2 Φh 2 = 0. ( 20
)
Here Φ = P 11 /h 2 is a one-dimensional analogue of Ψ (P 11 is a 'one-dimensional' determinant of P).

It conserves along trajectories :

DΦ Dt = 0. ( 21
)
The governing equations ( 18)-( 20) are reminiscent of the Euler equations of compressible fluids with the specific internal energy given by :

e = gh 2 + Φh 2 2 , (22)
and the pressure :

p = gh 2 2 + Φh 3 . ( 23
)
The enstrophy Φ increases across the shock as the conventional entropy for the Euler equations. Thus, P 11 remains positive across the shock if it was initially positive.

'Entropy' increase across the shocks : 2D case study

The previous 1D study suggests us a hypothesis about the following set of Rankine-Hugoniot relations. Let us suppose that

s T P n h = 0, (24)
[s T P s] = 0.

Then the 'entropy' Ψ of the system is increasing across the shock if and only if n T P n h 2 is increasing across the shock. Indeed, this fact follows from the inequalities :

Ψ = (n T P n)(s T P s) -(s T P n) 2 h 2 = n T P n h 2 s T P s - (hs T P n) 2 h 4 (26) = n T P n h 2 s T P s 0 - (hs T P n) 2 h 4 0 > n T P n h 2 0 s T P s 0 - (hs T P n) 2 h 4 0 = Ψ 0 .
Here the index '0' denotes the state before the shock. The justification of the jump relations (24), (25) which guarantees the inequality (26) will be done below by using a specific splitting procedure in solving the non-conservative equations.

Dissipative terms compatible with the positive definiteness of the Reynolds stress tensor

We add now dissipative terms in the model :

h t + div(hu) = 0, ( 27
) (hu) t + div(hu ⊗ u + gh 2 2 I + hP) = -C f |u|u, DP Dt + ∂u ∂x P + P ∂u ∂x T = D, D = D T .
Here C f is the Chézy coefficient, and D is the dissipation tensor to be defined. The equations (27) should also satisfy the energy conservation law :

∂ ∂t h 1 2 |u| 2 + e i + e T + div hu 1 2 |u| 2 + e i + e T + gh 2 2 I + hP u = -C f |u| 3 -Q, ( 28
)
where the dissipative source term Q should be positive. The positivity of Q is the analogue of the second law of thermodynamics. The compatibility of ( 27) and ( 28) implies the expression for Q in terms of D:

tr (D) = - 2 h Q. ( 29
)
The equations for the Reynolds stress tensor P can be written in Cartesian coordinates as follows : Or, in invariant form :

DP
h 2 D Dt detP h 2 = tr (P) tr (D) -tr (PD) .
By analogy with the Stokes hypotheses, we assume that the dissipation tensor D is an isotropic tensor function of P. Then, for two-dimensional case, D is linear in P :

D = - 2 h |u| 3 αP + δ 2 I ,
where α and δ are functions of invariants of P. The multiplier -2|u| 3 /h is for convenience only. So, α has the dimension s 2 m -2 , while δ is dimensionless. Consider the simplest case where δ = 0. This choice allows us the reduction to the Saint-Venant equations in the limit P = 0. One has finally :

D = - 2α h |u| 3 P. (30)
In particular, this implies the equation for Ψ in the form :

h 2 D Dt det (P) h 2 = - 4α h |u| 3 det (P) . (31)
Equations ( 29) and (30) imply the following relation between Q and α :

Q = αtr (P) |u| 3 . ( 32
)
We will finally choose Q as in [START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF] :

αtr (P) = max 0, C r tr(P) h 2 -ϕ tr(P) h 2 . ( 33
)
Here ϕ and C r are the model constants : ϕ is associated with the enstrophy of small vortexes in the vicinity of the bottom, and C r is the dissipation coefficient associated with the roller formation [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF][START_REF] Richard | Elaboration d'un modèle d'écoulements turbulents en faible profondeur: application au ressaut hydraulique et aux trains de rouleaux[END_REF]. We will give further the values of these constants evaluated from experimental data. The formula (33) allows us to recover the 1D case studied in [START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF]. As it follows from ( 31) and ( 33), 'entropy' Ψ is decreasing on continuous solutions, but always stays positive. This means that the dissipation law also guarantees the positive definiteness of P.

Another choice for the dissipation tensor could also be as follows : Formally, such a choice also guarantees the positive definiteness of P and thus is reasonable from the physical point of view. However, the wave profiles obtained with such a law do not not correspond to the experimental ones. Indeed, Brock [11,12,13] measured the stationary roll wave profiles in different conditions (different slopes and wall roughness). He noticed that the roll wave profiles contain always the following three essential parts: first, a sudden increase of the depth since all waves break i.e. acquire steeply sloping wave front, second, a continuous zone where the depth increases progressively, and third, a slowly decreasing zone until a new hydraulic jump (see such a profile in Figure 6). This is the reason to prefer the option (33) because it allows us to reproduce experimental profiles.

αtr (P) = C r . (34)

Splitting method in Cartesian coordinates

For convenience, we write here once again the governing equations (1) in Cartesian coordinates : As mentioned earlier, the system is hyperbolic but not in conservative form. We will define now an almost conservative formulation of the governing equations based on a splitting procedure. A conventional geometric splitting is first applied : the governing equations are solved first in x and then in y direction. We will do two steps more for each uni-directional subsystem referring to this as a 'physical' splitting. Roughly speaking, each uni-directional subsystem describing two types of waves (surface and shear waves), is split into two subsystems which are hyperbolic and contain only one type of waves (a or b waves). Each physical subsystem admits its own energy conservation law and its own 'entropy'. In physical terms, one can say that one performs a 'thermodynamically compatible' splitting. Analogous approach was proposed in [19]. Consider first the subsystem in x-direction :

h t +
h t + uh x + hu x = 0, ( 35
)
u t + uu x + gh x + 1 h (hP 11 ) x = 0, v t + uv x + 1 h (hP 12 ) x = 0, P 11t + uP 11x + 2P 11 u x = 0, P 12t + uP 12x + P 12 u x + P 11 v x = 0, P 22t + uP 22x + 2P 12 v x = 0.
As mentioned earlier, the system is hyperbolic and admits three types of waves : a contact discontinuity propagating with the velocity u, surface gravity waves propagating with the velocity u ± a, and shear waves propagating with a smaller velocity u ± b. The eigenfields corresponding to the contact discontinuity and b-waves are linearly degenerate in the sense of Lax, while the eigenfield corresponding to a-waves is genuinely non-linear in the sense of Lax. The idea is to split [START_REF] Ovsyannikov | Two-layer shallow water model[END_REF] into two subsystems treating separately a-waves and b-waves. Subsystem for a-waves is :

h t + uh x + hu x = 0, ( 36
)
u t + uu x + gh x + 1 h (hP 11 ) x = 0, v t + uv x = 0, P 11t + uP 11x + 2P 11 u x = 0, P 12t + uP 12x + P 12 u x = 0, P 22t + uP 22x = 0.
The equations (36) admit the following conservative form :

h t + (uh) x = 0, (37)
h 2 = 0, D Dt = ∂ ∂t + u ∂ ∂x ,
and the 'enstrophy' equation :

D Dt P 11 h 2 = 0.
At this step, one solves six conservation laws [START_REF] Ovsyannikov | Lectures on the Fundamentals of Gas Dynamics[END_REF]. The 'entropy' will increase across the shock. Indeed, the ratio P 12 /h and the component P 22 are conserved across the shock, so we need only the increase of P 11 /h 2 what is the case for 1D flows [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF]. The subsystem for b-waves is :

h t = 0, (38)
(hu) t = 0, (hv) t + (hP 12 ) x = 0, P 11t = 0,

P 12t + P 11 v x = 0, P 22t + 2P 12 v x = 0.
An 'almost' conservative form [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF] for b-waves is :

h t = 0, (39)
(hu) t = 0, (hv) t + (hP 12 ) x = 0, P 11t = 0,

P 12t + P 11 v x = 0, h v 2 2 + P 22 2 t + (hP 12 v) x = 0.
The shock relation for P 12 is well defined because for this subsystem P 11 is continuous across a shock. The 'entropy' conservation law :

∂ ∂t P 11 P 22 -P 2 12 h 2 = 0,
is a consequence of [START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF]. Since the characteristic field corresponding to b-waves is linearly degenerate, the shock velocity will coincide with b. As a consequence, the conservation of energy is equivalent to the conservation of 'entropy'. Thus, the positive definiteness of P is guaranteed even in the presence of shocks.

The study in y-direction is analogous. Indeed, one has :

h t + vh y + hv y = 0, (40)
P 22t + vP 22y + 2P 22 v y = 0.
The equations for a-waves and b-waves in y-direction are obtained from the corresponding equations ( 37) and ( 39) by the change of variables u → v, x → y, P 11 → P 22 , P 22 → P 11 .

6 Mathematical properties of split systems 6.1 Subsystem 1 for a-waves Subsystem (36) can be rewritten in the following form:

∂W ∂t + A(W) ∂W ∂x = 0, (41)
where the vector of unknowns W and matrix A are defined as:

W = (h, u, v, P 11 , P 12 , P 22 ) T , A =                   u h 0 0 0 0 (gh + P 11 )/h u 0 1 0 0 0 0 u 0 0 0 0 2P 11 0 u 0 0 0 P 12 0 0 u 0 0 0 0 0 0 u                   .
The eigenvalues of A are :

λ 1,2,3,4 = u, λ 5,6 = u ± a, a = gh + 3P 11 . (42)
For the multiple eigenvalue λ 1,2,3,4 = u we have four linearly independent right eigenvectors of A:

r 1 = (h, 0, 0, -(gh + P 11 ), 0, 0) T , ∇ W λ 1 • r 1 = 0, r 2 = (0, 0, 1, 0, 0, 0) T , ∇ W λ 2 • r 2 = 0, r 3 = (0, 0, 0, 0, 1, 0) T , ∇ W λ 3 • r 3 = 0, r 4 = (0, 0, 0, 0, 0, 1) T , ∇ W λ 4 • r 4 = 0.
For the eigenvalue λ = u + a one has:

r = (h, a, 0, 2b 2 , P 12 , 0) T , ∇ W λ 5 • r 5 = 3 2a a 2 + P 11 > 0. ( 43
)
For the eigenvalue λ = ua one has:

r 6 = (h, -a, 0, 2P 11 , P 12 , 0) T , ∇ W λ 6 • r 6 = - 3 2a a 2 + P 11 < 0. ( 44
)
The eigenvectors are linearly independent. Indeed,

det [r 1 , r 2 , r 3 , r 4 , r 5 , r 6 ] = -2a 3 h = 0. (45)
Hence, subsystem ( 41) is hyperbolic.

Subsystem 2 for b-waves

Subsystem [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF] can also be rewritten in the matrix form

∂W ∂t + A(W) ∂W ∂x = 0,
with the matrix A given by :

A =                   0 0 0 0 0 0 0 0 0 0 0 0 P 12 /h 0 0 0 1 0 0 0 0 0 0 0 0 0 P 11 0 0 0 0 0 2P 12 0 0 0                   . ( 46
)
It implies:

det(A -λI) = λ 4 (P 11 -λ 2 ) = 0. ( 47
)
Hence, one has 6 real eigenvalues:

λ 1,2,3,4 = 0, λ 5,6 = ±b, b = P 11 . (48)
For the multiple eigenvalue λ 1,2,3,4 = 0 one has 4 linearly independent right eigenvectors:

r 1 = 1, 0, 0, 0, - P 12 h , 0 T , ∇ W λ 1 • r 1 = 0, r 2 = (0, 1, 0, 0, 0, 0) T , ∇ W λ 2 • r 2 = 0, r 3 = (0, 0, 0, 1, 0, 0) T , ∇ W λ 3 • r 3 = 0, r 4 = (0, 0, 0, 0, 0, 1) T , ∇ W λ 4 • r 4 = 0. ( 49
)
For λ 5 = b one has :

r 5 = (0, 0, b, 0, b 2 , 2P 12 ), ∇ W λ 5 • r 5 = 0. ( 50
)
For λ 6 = -b one has :

r 6 = (0, 0, -b, 0, b 2 , 2P 12 ), ∇ W λ 6 • r 6 = 0. (51)
The eigenvectors are linearly independent :

det [r 1 , r 2 , r 3 , r 4 , r 5 , r 6 ] = -2P 11 P 11 = 0, if P 11 = 0. (52)
Hence, subsystem (38) is hyperbolic.

Numerical scheme

The numerical approach for (4) consists in solving first the model in x-direction, and then in ydirection. For each direction, the two subsystems for a-and b-waves are solved successively. A Godunov-type scheme augmented by a correction step (see below) is used for this aim. Finally, the source terms are integrated. The vector of unknowns for all subsystems is :

U = [h,
Here E = (u 2 + v 2 + gh + P 11 + P 22 )/2.

First subsystem : a-waves

The first subsystem (called 'subsystem 1') augmented with the energy conservation law for a-waves is :

                   h t + (uh) x = 0, (a) (hu) t + hu 2 + p x = 0, (b) (hv) t + (huv) x = 0, (c) (hP 11 ) t + (huP 11 ) x + 2hP 11 u x = 0, (d) P 12t + (uP 12 ) x = 0, (e) (hP 22 ) t + (huP 22 ) x = 0. (f ) (hE) t + (hEu + pu) x = 0, (g) (54)
with p = gh 2 /2 + hP 11 .

The equation (54.d) is not in conservative form, i.e. the product P 11 u x is not well defined across discontinuities. Hopefully, this system is overdetermined and the correct value of P 11 will be obtained using the energy equation (54.g). This system is solved in 3 steps which can be summarized as follows:

• Solve the Riemann problem using any Riemann solver.

• Evolve all conservative variables using the classical Godunov scheme.

• Compute P 11 from the energy equation.

Each step is detailed hereafter. For the first subsystem (54) the wave scheme is shown in Figure 1. The Rankine-Hugoniot relations read :

-

[h] S a + [hu] = 0, -[hu]S a + hu 2 + p = 0, -[hv]S a + [huv] = 0, -[P 12 ]S a + [uP 12 ] = 0, -[hP 22 ]S a + [huP 22 ] = 0, -[hE] S a + [hEu + pu] = 0.
Here S a is the discontinuity speed, and [f ] means the jump of any variable f . The jump relation for mass equations can be written in the following form :

[h(u -S a )] = 0, i.e. h(u -S a ) = m = const.

We denote m

L = h L (u L -S a L ) = h * L (u * -S a L ) m R = h R (u R -S a R ) = h * R (u * -S a R ).
Here the speeds of the left and right facing waves are obtained by using Davis' approximation :

S a R = max(u L + a L , u R + a R ), S a L = min(u L -a L , u R -a R ),
with a 2 L,R = gh L,R + 3P 11 L,R . The Rankine-Hugoniot relations imply the continuity of the following parameters in the 'star regions':

u * L = u * R = u * , v * L = v L , v * R = v R , p * L = p * R = p * , with p = gh 2 2 + hP 11 .
The momentum conservation law implies :

u * = p L -p R + m L u L -m R u R m L -m R ,
and

p * = m L m R (u R -u L ) + m L p R -p L m R m L -m R .
Finally, the conservative variables in the 'star region' are:

h * L,R = m L,R /(u * -S a L,R ), h * L,R u * L,R = h L,R u L,R (u L,R -S a L,R ) + p L,R -p * /(u * -S a L,R ), h * L,R v * L,R = h L,R v L,R (u L,R -S a L,R )/(u * -S a L,R ), P * 12 L,R = P 12 L,R u L,R -S a L,R /(u * -S a L,R ), h * l,r E * L,R = h L,R E L,R (u L,R -S a L,R ) + p L u L -p * u * /(u * -S a L,R
). These relations allow us to construct an approximate Riemann solver.

Godunov-type scheme for subsystem 1

Subsystem 1 can be rewritten in the following form :

∂U ∂t + ∂F ∂x = 0. ( 55
)
Here the vector of conservative variables U and the vector of fluxes F are :

U = (h,
For simplicity, we use here the same generic notation U for the vector of conservative variables, even if this vector does not contain the component hP 11 . Using the flux solution obtained in section 7.1.1 at the edge of each cells (F * ), the conservative variables are evolved as :

U n+1 i = U n i - ∆t ∆x F * ,n i+1/2 -F * ,n i-1/2 . ( 57
)
Here ∆x is the discretisation step in the x-direction, ∆t is the time step verifiing the Courant-Friedrichs-Lewy (CFL) condition (∆t ≤ ∆x/S max ). At this step the non-conservative equation for hP 11 is removed.

Computation of hP 11

Since the system is overdetermined, the non-conservative term hP 11 is obtained by using the total energy equation :

hP 11 = 2hE -gh 2 -hP 22 - (hu) 2 + (hv) 2 h .

Subsystem 2 : b-waves

The subsystem for b-waves (called 'subsystem 2'), augmented with the energy conservation law, reads:

                   h t = 0, (a) (hu) t = 0, (b) (hv) t + (hP 12 ) x = 0, (c) (hP 11 ) t = 0, (d) P 12t + P 11 v x = 0, (e) (hP 22 ) t + 2hP 12 v x = 0. (f ) (hE) t + (hP 12 v) x = 0. (g) (58)
Again, the system (58) is overdetermined. One can note that there are two nonconservative equations (58.e) and (58.f). The product P 11 v x in (58.e) is well defined since P 11 is conserved across the shock. It is not the case for the term 2hP 12 v x in equation (58.f). Since the jump relation is not well defined, there is no hope to obtain a good value of P 22 . In the first step, this equation will be removed. Then, the energy conservation law will be again used to compute the value of hP 22 . The numerical procedure will be similar to a-wave procedure :

• Solve the Riemann problem using any Riemann solver.

• Evolve all conservative variables using the classical Godunov scheme.

• Compute hP 22 from the energy equation.

These steps and some important remarks on the evolution of the 'entropy' will be detailed hereafter. Subsystem 2 for b-waves can be rewritten in the following form:

∂U ∂t + ∂F ∂x + K ∂v ∂x = 0, (59)
with U = (h, hu, hv, hP 11 , P 12 , hE) T , F = F(U) = (0, 0, hP 12 , 0, 0, hP 12 v) T , K = K(U) = (0, 0, 0, 0, P 11 , 0) T .

(60)

Again, for simplicity, the generic notation U for the 'conservative' vector is used, even if this vector does not contain now the component hP 22 . The wave scheme is shown in Figure 2. The set of jump relations for subsystem 2 is :

-[h]S b = 0, -[u]S b = 0, -[P 11 ]S b = 0, ( 61
)
-[hv]S b + [hP 12 ] = 0, ( 62
)
-h v 2 + P 22 2 S b + [hvP 12 ] = 0, ( 63
)
Across the contact discontinuity (S b = 0) a priori one has :

[h] = 0, [u] = 0, [P 11 ] = 0, ( 64
) but h * L P * 12 L = h * R P * 12 R = (hP 12 ) * , v * = v * L = v * R .
Across the shocks (S b = 0) one has from (61):

h L = h * L , h R = h * R , u R = u * R , P 11 R = P * 11 R . (65)
Equation (62) gives :

(hP 12 ) * = h R h L S b R P 12L -S b L P 12R -S b L S b R (v L -v R ) h R S b R -h L S b L , (66)
and

v * = h L (P 12 L -S b L v L ) -h R (P 12 R -S b R v R ) h R S b R -h L S b L .
The extreme wave speeds can be estimated by Davis' approximation :

S b L = min(- √ P 11L , - √ P 11R ), S b R = max( √ P 11L , √ P 11R ),
Since the eigenfields is linearly degenerate, the wave speed on the right (left) only depend on the right (left) state. Thus, another possibility is :

S b L = - √ P 11 L , S b R = √ P 11 R . ( 67
)
This choice, more precise, is used for all the numerical results presented here. These relations allows us to construct an approximate Riemann solver.

Godunov-type scheme for subsystem 2

The non-conservative equations for P 12 and P 22 necessitate a specific numerical treatment. For the equations in conservative form, we use the following Godunov-type scheme :

U n+1 i = U n i + ∆t ∆x (F * ,n i+1/2 -F * ,n i-1/2 ) + K n i (v * ,n i+1/2 -v * ,n i-1/2 ). ( 68
)

Computation of hP 22

The non-conservative term hP 22 is obtained by using the total energy equation :

hP 22 = 2hE -gh 2 -hP 11 - (hu) 2 + (hv) 2 h . ( 69
)

Conservation of the 'entropy'

It is striking that subsystem (58) conserves also the mathematical 'entropy'. Indeed, consider the Rankine-Hugoniot relations :

   -[v]S b + [P 12 ] = 0, -[(v 2 + P 22 )/2]S b + [P 12 v] = 0. ( 70
)
Here, for any f, [f ] = ff 0 , where the index '0' means the state before the shock. It implies :

   P 12 = (P 12 ) 0 + S b (v -v 0 ), P 22 = (P 22 ) 0 + [2P 12 v]/S b -[v 2 ], (71)
We will show now that [detP] = 0. Since h is continuous across the shock, it implies the conservation of the 'entropy'. As P 11 is continuous, we obtain : detP -detP 0 = (P 11 ) 0 (P 22 -(P 22 ) 0 ) + ((P 12 ) 0 -P 12 ) ((P 12 ) 0 + P 12 ) .

Replacing (71) into (72) one obtains :

detP -detP 0 = (P 11 ) 0 [2P 12 v] S b -[v 2 ] -S b [v] 2(P 12 ) 0 + S b (v -v 0 ) = = (P 11 ) 0 2 (P 12 ) 0 + S b (v -v 0 ) v -(P 12 ) 0 v 0 S b -[v 2 ] -S b (v -v 0 ) 2(P 12 ) 0 + S b (v -v 0 ) .
As (S b ) 2 = P 11 = (P 11 ) 0 , this can be simplified to :

detP -detP 0 = (P 11 ) 0 2(v -v 0 )v -(v 2 -v 2 0 ) -(v -v 0 ) 2 = 0.
The mathematical 'entropy' Ψ = detP h 2 is thus conserved because the corresponding eigenfields are linearly degenerate in the sense of Lax. Finally, the 'entropy' is increasing after solving the first subsystem, and does not change in solving the second subsystem.

Integration of the source terms

To add the source terms, we integrate the ordinary differential equations :

dU dt = S(U) (73)
with the full unknown vector U = (h, hu, hv, hP 11 , P 12 , hP 22 , hE)

T . The initial condition U| t=0 = Ūn+1 are obtained from the previous splitting steps. Here the source term is :

S(U) = 0, -gh∇b -C f u|u|, hD 11 , D 12 , hD 22 , -gh∇b • u -C f |u| 3 -Q T , (74)
where

|u| = u 2 + v 2 , D ij = - 2α h |u| 3 P ij , Q = αtr(P)|u| 3 , αtr(P) = max   0, C r tr(P) h 2 -ϕ tr(P) h 2    , ∇b = ∂b ∂x , ∂b ∂y T .
We added here the bottom topography z = b(x, y). The equations are written in the reference frame where the gravity is orthogonal to the (x, y)-plane. So, the case of a mild slope bottom is considered.

For applications, we will consider only the case of a constant slope bottom :

-∇b = (tanθ, 0) , θ > 0,
where θ is the inclination angle. Even if we have seven equations for six variables, the equations are compatible : the energy equation is a consequence of the equations of mass, momentum and stress . This system of ordinary differential equation is solved by the first order scheme with the Euler method. A second order Runge-Kutta method is used for the higher order extension.

Algorithm summary

First, the 2D system is split into two 1D systems along each spatial direction. Then, each 1D system is successively split into two additional 1D subsystems. All the systems are hyperbolic. The numerical algorithm can be summarized as follows:

1. Subsystem 1 in x-direction.

• Solve the Riemann Problem in x-direction for subsystem 1 at each cell boundary without source terms. The approximate HLLC solver was used for this aim.

• Evolve all flow variables with the Godunov-type method.

• Correction of the non-conservative variable (hP 11 ) using the energy conservation law.

2. Subsystem 2 in x-direction.

• Solve the Riemann Problem in x-direction for subsystem 2 at each cell boundary.

• Evolve for the conservative equations flow variables with the Godunov-type method.

• The correction of the non-conservative variable (hP 22 ) using energy conservation law.

3. The same procedure is repeated in y-direction by changing u → v, x → y, P 11 → P 22 , P 22 → P 11 .

4. Integration of the source terms.

5. Start again for the next time step.

Numerical Results

In this section, we present numerical results obtained with the splitting procedure described above.

1D shear test problem

We solve here the governing equations in the case of vanishing source term. The initial discontinuity is located at x = 0.5 m. The initial depth h is 0.01 m, the normal velocity u is zero everywhere, the components of the stress tensor are P 11 = P 22 = 10 -4 m 2 /s 2 , P 12 = 0 m 2 /s 2 , the tangential velocity v is 0.2 m/s on the left, and -0.2 m/s on the right. The first order Godunov method is used, with CFL number 0.3. The tangential velocity, tangential stress P 12 and stress P 22 are shown in Figure 3 at time instant 10 s for 500, 1000 and 10000 grid cells. The other variables do not evolve in time, they are shown in Figure 4.

-0.

1D dam-break problem

We solve here the governing equations in the case of vanishing source term. The initial discontinuity is located at x = 0.5 m. The initial normal and tangential velocities are zero everywhere, the components of the stress tensor are : P 11 = P 22 = 10 -4 m 2 /s 2 , P 12 = 0 m 2 /s 2 , the fluid depth is 0.02 m at the left and 0.01 m at the right. The MUSCL extension of Godunov method is used [START_REF] Van Leer | Towards the ultimate conservative difference scheme, V. A second order sequel to Godunov's method[END_REF]. The results are obtained by using Minmod limiter for subsystem 1, and van Leer limiter for subsystem 2. The solution for 100, 1000 and 10000 grid cells is shown in Figure 5

1D roll waves

We solve here the governing equations with source term. The bottom is inclined either in x-direction, or in y-direction. Using periodic conditions in the direction of wave propagation, and the rigid wall condition in the transverse direction, we observe the formation of roll waves as in [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF]24]. The initial conditions are taken in the same form as in the above mentioned references. For example, for the flow in x-direction, one takes : Figure 6: Formation of roll waves : the depth and normal stress are shown at time instant 25 s for flows both in x -and y -direction (on the left and on the right, respectively). The results are completely symmetric. The MUSCL extension of Godunov method was used with 1000 grid cells for the parameter set mentioned above. CFL number is 0.6.

h(x, y, 0) = h 0 1 + asin 2πx L x , u(x, y, 0) = gh 0 tanθ C f , v(x,

Comparison with a 2D analytical solution

Here we present an analytical solutions to (1). This solution is a generalisation of solutions with linear velocity profile in x an y found by [START_REF] Sedov | On the integration of the equations of one-dimensional gas motion[END_REF] and [START_REF] Ovsyannikov | A new solution of the hydrodynamic equations[END_REF] for the Euler equation (for proof, see Appendix B ) :

                             h = h 0 1 + β 2 t 2 , U = β 1 + β 2 t 2   βtx + y -x + βty   , P = 1 (1 + β 2 t 2 ) 2   λ + γβ 2 t 2 , (λ -γ)βt (λ -γ)βt, γ + λβ 2 t 2   , (75)
where h 0 > 0, β, λ > 0, γ > 0 are constant. The stress tensor P is not spherical. We take here h 0 = 1 m, λ = 0.1 m 2 /s 2 , γ = 0.01 m 2 /s 2 , β = 10 -3 s -1 , CF L = 0.5. The error was calculated at time instant 10 s in the following way :

error(P 11 ) = max x,y |P 11numerical -P 11analytical | λ , error(P 12 ) = max x,y |P 12numerical -P 12analytical | λ + γ , error(P 22 ) = max x,y |P 22numerical -P 22analytical | γ .
These errors are shown in the Figure 7 as a function of the grid cell size . The calculation region is a square of length L x = L y = 10 m. Different regular Cartesian grids were used (100 × 100, 200 × 200, 400 × 400, 800 × 800, 1600 × 1600) with first order Godunov method. This test shows the convergence to the exact solution. The slope of the convergence lines are the same for all stress components.

2D roll waves

We solve here the 2D governing equations with source term. The bottom is inclined in x-direction. In the direction of wave propagation we use periodic conditions for the vector of unknowns (U(0, y, t) = U(L x , y, t)), and in the transverse direction we use the rigid wall conditions (v(x, 0, t) = v(0, L y , t) = 0, P 12 (x, 0, t) = P 12 (x, L y , t) = 0, P 22 (x, 0, t) = P 22 (x, L y , t) = 0, and zero Neumann conditions for other variables). The initial conditions are taken in the form : Here θ = 0.05011 [rad] is the inclination angle, C f = 0.0036 is the Chézy coefficient, h 0 = 7.98 × 10 -3 m, ϕ = 22.76 s -2 , g = 9.81 ms -2 , C r = 0.00035, L x = 1.3 m, L y = 0.5 m, and m and k are numbers to be chosen. The set of physical parameters corresponds to that considered in [24] to describe the formation of 1D roll waves from a uniform flow having the same structure as in Brock's experiments [11,12,13]. A necessary condition for the formation of such waves is that the corresponding Froude number is larger than two :

h(x, y, 0) = h 0 1 + asin 2πmx L x + asin 2πky L y , u(x, y, 0) = gh 0 tanθ C f , v(x, y, 0) = 0, (76)
F g = gh0tan(θ) C f (gh0+3ϕh 2 0 /2) > 2.
For the flow parameters taken above, this value is about 3.7.

Figure 8 shows the convergent numerical solution obtained from a uniform flow perturbed both in x and y direction with m = 1 and k = 1 (see initial data (76)). The existence of a characteristic transverse wave length of the jump toe perimeter (the line where the gradient of the layer thickness h jumps) is clearly visible. One can count approximately 7 transverse waves.

The transverse structure formation scenario is rather surprising. First, a one-dimensional roll wave is forming, without any transverse structure. This 1D structure is formed in approximately 10 s and corresponds to a standard 1D experimental profile [11,12,13] (see Figure 9). The transverse structure starts to form in approximately 17 s and becomes stationary after approximately 35 s. Thus, the whole scenario is described : from uniform unstable flow to 1D roll waves, and finally to the formation of transverse waves at the jump toe perimeter (formation of 'fingers'). 'Fingering' of the bore toe perimeter was already observed, in particular, in field experiments [28]. To understand if the number of transverse waves per unit length is independent on the initial perturbations, we considered also the initial data (76) with m = 1 and k = 4. Thus, the initial transverse perturbation contains four waves, and one could expect that the number of transverse Figure 9: The one-dimensional wave structure is formed in approximately 10 s (solid line) and corresponds to the classical experilmental profile (dots) measured in [11,12,13].

waves formed during the evolution will increase. The result is rather surprising : the number of waves is almost independent on the initial perturbation, only approximately eight waves were formed instead of seven (see Figure 10).

An additional test was also performed to show that the number of transverse waves per unit length is invariant with respect to the domain width. For this, we multiplied the length L y by two (L y = 1 m), and took m = 1, k = 1. As a consequence, the number of transverse waves was also multiplied by two (see Figure 11 and12).

The form of the dissipation term defined by the coefficient α given by ( 33) is very important. Indeed, this form was chosen to obtain the same dissipation source term Q as in 1D case studied in [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF][START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF]. A simplified formula for α (34) will give us only 1D profile, without any transverse structure (see Figure 13). Moreover, the corresponding 1D profile does not correspond to the experimentally observed profiles in [11,12,13] shown in Figure 9.

Conclusion

We propose a numerical method consisting in 'physical' splitting of the hyperbolic non-conservative equations for shear shallow water flows. Each split subsystem is hyperbolic, contains only one type of waves, and admits the energy and 'entropy' conservation laws. Moreover, such a splitting allows us to naturally define a weak solution to our system which is compatible to the positive definiteness of the Reynolds stress tensor P. The dissipation is introduced in invariant form and also guarantees the positive definiteness of P. In the limit of one-dimensional flows, the roll waves solutions obtained earlier in [START_REF] Richard | A new model of roll waves: comparison with Brocks experiments[END_REF][START_REF] Richard | Elaboration d'un modèle d'écoulements turbulents en faible profondeur: application au ressaut hydraulique et aux trains de rouleaux[END_REF] are recovered.

An interesting feature of the model is the formation of transverse structures at the jump toe perimeter ('fingers') from one-dimensional initial data which are harmonically perturbed in the transverse direction. The number of waves does not depend neither on the amplitude nor on the number of transverse waves in the initial perturbation. Thus, the full transition scenario is observed in the formation of roll waves : from uniform flow to one-dimensional roll waves, and, finally, to 2D transverse 'fingering' of roll wave profiles.

The method can be extended to the general 3D Reynolds averaged models of barotropic flows, because the structure of the governing equations is exactly the same. The case of non-barotropic turbulent flows demands an additional modelling because both entropies, physical and 'mathematical', will increase across the shock. Thus it will be necessary to separate such an increase into two parts : the thermodynamic and turbulent ones [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF]. This will be the subject of our future work.

Acknowledgement The authors thank Boniface Nkonga for useful discussion. This work was partially supported by l'Agence Nationale de la Recherche, France (grant numbers ANR-13-BS01-0009-01, ANR-14-ASTR-0016-01, ANR-11-LABX-0092, and ANR-11-IDEX-0001-02). (34) was taken. The result obtained is very surprising : transverse structures are not developed. The second fact, less surprising, is that the one-dimensional profile is not that of Brock's type.

Appendix A : Conservation laws

We write down once again the governing equations for shear flows without right hand sides :

h t + uh x + vh y + hu x + hv y = 0, (77)
u t + uu x + vu y + gh x + 1 h (hP 11 ) x + 1 h (hP 12 ) y = 0, v t + uv x + vv y + gh y + 1 h (hP 12 ) x + 1 h (hP 22 ) y = 0, P 11t + uP 11x + vP 11y + 2P 11 u x + 2P 12 u y = 0, P 12t + uP 12x + vP 12y + P 12 (u x + v y ) + P 11 v x + P 22 u y = 0, P 22t + uP 22x + vP 22y + 2P 12 v x + 2P 22 v y = 0.
We are looking for conservation laws admitting by (77) :

∂f ∂t + ∂g 1 ∂x + ∂g 2 ∂y = 0, (78)
where the unknown functions f, g 1 , g 2 depend on h, u and P. Developing (78) one obtains :

f h h t + f u u t + f v v t + f P11 P 11t + f P12 P 12t + f P22 P 22t + g 1h h x + g 1u u x + g 1v v x + g 1P11 P 11x +g 1P12 P 12x + g 1P22 P 22x + g 2h h y + g 2u u y + g 2v v y + g 2P11 P 11y + g 2P12 P 12y + g 2P22 P 22y = 0.
Substituting the time derivatives of unknowns from (77) we have :

f h (-h x u -hu x -h y v -hv y ) + f u -uu x -vu y -gh x -P 11x -P 12y -P11 h h x -P12 h h y +f v -uv x -vv y -gh y -P 12x -P 22y -P12 h h x -P22 h h y + f P11 (-uP 11x -vP 11y -2P 11 u x -2P 12 u y ) +f P12 (-uP 12x -vP 12y -P 12 u x -P 22 u y -P 11 v x -P 12 v y ) + f P22 (-uP 22x -vP 22y -2P 12 v x -2P 22 v y ) +g 1h h x + g 1u u x + g 1v v x + g 1P11 P 11x + g 1P12 P 12x + g 1P22 P 22x
+g 2h h y + g 2u u y + g 2v v y + g 2P11 P 11y + +g 2P12 P 12y + g 2P22 P 22y = 0. (79) As the space derivatives of unknowns are independent, and (79) should be satisfied identically, we obtain the following overdetermined system of equations :

-uf h -gf u - P 12 h f v - P 11 h f u + g 1h = 0, (80a)
-hf h -uf u -2P 11 f P11 -P 12 f P12 + g 1u = 0, (80b) -uf v -2P 12 f P22 -P 11 f P12 + g 1v = 0, (80c)
-vf h - P 12 h f u -gf v - P 22 h f v + g 2h = 0, (80d)
-hf h -vf v -2P 22 f P22 -P 12 f P12 + g 2v = 0, (80e) -vf u -2P 12 f P11 -P 22 f P12 + g 2u = 0, (80f) -f u -uf P11 + g 1P11 = 0, (80g) -f v -uf P12 + g 1P12 = 0, (80h) -uf P22 + g 1P22 = 0, (80i) -f u -vf P12 + g 2P12 = 0, (80j) -f v -vf P22 + g 2P22 = 0, (80k) -vf P11 + g 2P11 = 0. (80l)
A useful remark is that there is no conservation laws with the density f depending only on P. The remark will allow us to neglect in determining f the integration 'constants' depending only on P. Indeed, equations (80g) -(80l) imply that in this case up to an additive constant g 1 = uf , g 2 = vf . The equations (80b), (80c), (80e) and (80f) form then a linear with respect to f overdetermined system of equations. The analysis shows that the solution is only trivial.

Eliminating the functions g 1 and g 2 from the last 6 equations of (80) we obtain :

                                   f uP11 = 0, f uP22 = 0, f uP12 = 0, f vP11 = 0, f vP12 = 0, f vP22 = 0. (81)
It implies :

f (h, u, v, P 11 , P 12 , P 22 ) = f 1 (h, u, v) + f 2 (h, P 11 , P 12 , P 22 ) ,
where f i , i = 1, 2, are arbitrary functions (in the following, several new arbitrary functions will appear). Eliminating the functions g 1 and g 2 from the first 6 equations of (80) we obtain :

h 2 f h hu + (P 12 f P12 + 2P 22 f P22 ) u -(P 22 f P12 + 2P 12 f P11 ) v = 0, (82a)
h 2 f h hv + (P 12 f P12 + 2P 11 f P11 ) v -(P 11 f P12 + 2P 12 f P22 ) u = 0, (82b)
P 12 h f u + g + P 22 h f v u -(P 22 f P12 + 2P 12 f P11 ) h = 0, (82c)
P 12 h f v + g + P 11 h f u v -(P 11 f P12 + 2P 12 f P22 ) h = 0, (82d)
P 12 h f u + g + P 22 h f v v -(P 12 f P12 + 2P 22 f P22 ) h -hf hh = 0, (82e)
P 12 h f v + g + P 11 h f u u -(P 12 f P12 + 2P 11 f P11 ) h -hf hh = 0. (82f)
In particular, (81), (82a) and (82b) imply :

f h hu = 0, f h hv = 0. (83)
Then ( 81) and (83) imply :

f (h, u, v, P 11 , P 12 , P 22 ) = hc 1 (u, v) + f 2 (h, P 11 , P 12 , P 22 ) . (84)
Taking the difference (82c) and (82d), then (82e) and (82f), and differentiating these differences with respect to u and v, one can immediately derive that c 1 (u, v) is quadratic with respect to u and v :

c 1 (u, v) = kuv + n 1 u 2 + n 2 u + m 1 v 2 + m 2 v,
where k, n 1 , n 2 , m 1 , m 2 are constants. The general form f is then :

f (h, u, v, P 11 , P 12 , P 22 ) = h(kuv + n 1 u 2 + n 2 u + m 1 v 2 + m 2 v) + f 2 (h, P 11 , P 12 , P 22 ) . (85)
The next step is thus to determine the function f 2 (h, P 11 , P 12 , P 22 ) from the overdetermined system obtained by replacing (85) into (82c), (82d), (82e) and (82f):

2P 12 n 1 + (gh + P 22 ) k -(P 22 f 2P12 + 2P 12 f 2P11 ) h = 0, (86a)
2P 12 m 1 + (gh + P 11 ) k -(P 11 f 2P12 + 2P 12 f 2P22 ) h = 0, (86b)
P 12 k + 2 (gh + P 22 ) m 1 -(P 12 f 2P12 + 2P 22 f 2P22 ) h -hf 2hh = 0, (86c)
P 12 k + 2 (gh + P 11 ) n 1 -(P 12 f 2P12 + 2P 11 f 2P11 ) h -hf 2hh = 0. ( 86d
)
Taking the difference of (86c) and (86d) and differentiating this difference with respect to h, one obtains :

P 11 (f 2hh ) P11 -P 22 (f 2hh ) P22 = g (n 1 -m 1 )
.

The general solution of this equation is :

f 2hh = φ (h, P 12 , X) + g (n 1 -m 1 ) ln P 11 , with X = P 11 P 22 . (87)
Taking now the difference of (86a) and (86b) and differentiating this difference with respect to h, one obtains :

(P 11 -P 22 ) (f 2hh ) P12 + 2P 12 (f 2hh ) P22 -(f 2hh ) P11 = 0 (88)
Substituting ( 87) into (88) we obtain :

P 2 11 -X ∂φ (h, P 12 , X) ∂P 12 + 2P 12 P 2 11 -X ∂φ (h, P 12 , X) ∂X -2P 12 g (m 1 -n 1 ) = 0.
This is a polynomial of degree two in P 11 . Since the identity should be valid for any P 11 , the coefficients of the polynomial vanish. This implies : φ P12 + 2P 12 φ X = 0 and m 1 = n 1 . Integrating two times in h, one obtains the general expression of f 2 :

f 2 = φ 1 (h, P 11 P 22 -P 2 12 ) + hφ 2 (P) + φ 3 (P).

One can always take φ 3 (P) vanishing because there is no conservation law with f depending only on P. With the condition m 1 = n 1 the difference (86c) -(86d) becomes :

P 11 φ 2P11 -P 22 φ 2P22 = (P 11 -P 22 ) n 1 .
Hence, φ 2 (P) = n 1 (P 11 + P 22 ) + ν 1 (P 12 , X) , and f 2 (h, P) = φ 1 h, P 11 P 22 -P 2 12 + n 1 h (P 11 + P 22 ) + hν 1 (P 12 , X) . We substitute the function f 2 (h, P) into the equation (86a) and obtain:

(gh + P 22 )k -(P 22 ν 1P12 + 2P 12 P 22 ν 1X ) = 0. ( 89
)
This is a linear function with respect to P 22 . It identically vanishes if and only if :

k = 0, ν 1P12 + 2P 12 ν 1X = 0.
Hence, ν 1 is a function of only one argument ∆ = P 11 P 22 -P 2 12 . One can insert this function into the general expression of f 2 . One has now :

f 2 (h, P) = ψ 1 h, P 11 P 22 -P 2 12 + n 1 h (P 11 + P 22 ) , (90)
with a new function ψ 1 of h and ∆. A simplified form of ( 86) is now :

2P 12 n 1 -(P 22 f 2P12 + 2P 12 f 2P11 ) h = 0, (91a)
2P 12 n 1 -(P 11 f 2P12 + 2P 12 f 2P22 ) h = 0, (91b) 2 (gh
+ P 22 ) n 1 -(P 12 f 2P12 + 2P 22 f 2P22 ) h -hf 2hh = 0, (91c) 2 (gh
+ P 11 ) n 1 -(P 12 f 2P12 + 2P 11 f 2P11 ) h -hf 2hh = 0. ( 91d
)
The equations (91a), (91b) are identically satisfied if we replace expression f 2 given by (90). The equations (91c), (91d) give us the same equation :

2∆ (ψ 1h ) ∆ + h (ψ 1h ) h = 2ghn 1 .
Its solution is :

ψ 1h (h, ∆) = A ∆ h 2 + 2gn 1 h,
with an arbitrary function A (s) , s = Ψ = ∆/h 2 . Integration with respect to h gives us :

ψ 1 (h, ∆) = hB ∆ h 2 + n 1 gh 2 ,
where B(s) is an arbitrary function. The final representation of the density function f is :

f (h, u, P) = hB P 11 P 22 -P 2 12 h 2 + n 1 (gh 2 + h(P 11 + P 22 ) + hu 2 + hv 2 ) + n 2 hu + m 2 hv,
where n 1 , n 2 , m 2 are arbitrary constants, and B is an arbitrary function of Ψ. The conservation of mass corresponds to a particular case where B = const.

Appendix B : Analytical solutions

The analytical solution to (1) we present here is a generalisation of the solution with linear in space velocity profile found by [START_REF] Sedov | On the integration of the equations of one-dimensional gas motion[END_REF] and [START_REF] Ovsyannikov | A new solution of the hydrodynamic equations[END_REF] for the Euler equations. Thus, we are looking for the solution of the type :

u = A(t)
x, h = h(t), P = P(t).

Here A(t) is a time dependent matrix. System (1) becomes :

           ḣ + htr(A) = 0, Ȧ + A 2 = 0, Ṗ + AP + PA T = 0. (92)
Here 'dot' means the time derivative. The solution A(t) of the second equation of (92) is given in the form :

A = A 0 (I + A 0 t) -1 , A 0 = const.
Since the corresponding matrices commute, A can also be written as :

A = (I + A 0 t) -1 A 0 .
Then P(t) verifies the equation :

Ṗ + A 0 (I + A 0 t) -1 P + P I + A T 0 t -1 A T 0 = 0 (93)
The solution P of (93) is :

P = (I + A 0 t) -1 P 0 I + A T 0 t -1 ,
with a constant symmetric positive definite matrix P 0 = P T 0 > 0. The solution is well defined for all t > 0, if the matrix I + A 0 t is invertible for any t. This is a case, for example, of antisymmetric matrix A 0 : A T 0 = -A 0 . As an example, consider the initial data :

A 0 =   0 β -β 0   , P 0 =   λ 0 0 γ   , h = h 0 ,
with constants β, λ > 0, γ > 0 and h 0 > 0. The solution is :

A(t) = β 1 + β 2 t 2   βt 1 -1 βt   , P = 1 (1 + β 2 t 2 ) 2   λ + γβ 2 t 2 (λ -γ)βt (λ -γ)βt γ + λβ 2 t 2   , h = h 0 1 + β 2 t 2 .

Introduction

Much work has been done in the study of the circular hydraulic jump [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15]. For example, when a jet of water impinging vertically on a flat surface, the flow spreads radially outwards and at some distance of the point of impingement there is a abrupt circular transformation of depth. This transformation is called the "circular hydraulic jump". It is commonly observed in a kitchen sink.

Our goal is a qualitative comparison between the numerical results using a two-dimensional mathematical model of shear shallow water flows [20], [17], and the experiments of Foglizzo et al.(CEA-Saclay, France) [1], [2], [3], where, in contrast with the above mentioned circular hydraulic jump in a kitchen sink (i.e. in a divergent flow), the hydraulic jump was formed in a convergent flow by injecting water inwards from a circular reservoir and evacuating it through a hollow vertical cylinder playing the role of obstacle. Even if the inflow of water was stationary and uniform, some transverse structures at the hydraulic jump toe appeared. In particular, the experiments demonstrated the evolution of the stationary radial hydraulic jump flow into a new flow pattern where a rotation singularity was formed at the free surface (see Figure 1). The formation of this flow pattern was already established through numerical resolution of the Saint-Venant equations with a Stokes-type friction law (i.e. for laminar flows) [1].

Foglizzo et al. also pointed out the analogy between this type of instability of circular hydraulic jump in the water fountain and a supernova formation. The analogy stems from a observation that both the hydraulic jump and a supernova represent shock waves that connect supercritical (supersonic) flows with subcritical (subsonic) flows. Also the astrophysical connection relies on the analogy between water waves in the fountain and acoustic waves in the star. Researchers use Our modeling is based on a 2D model of shear shallow water flows [20], [18], [19]), [17]. The numerical method for the 2D system is based on a splitting method that was recently developed for solving the non-conservative system of shear shallow water flows [17]. The whole system contains two types of waves (surface and shear waves), while each split subsystem contains only one type of waves. Here we used the same technique in polar coordinates. We solve only 1D Riemann problems in the coordinate directions. There are obvious advantages from a practical point of view: the 1D methods are relatively easy to develop.

The structure of the article is organized as follows. In the first two sections, the classical "shallow water" equations and "shear shallow water" equations are rewritten in polar coordinates. The dissipative terms are introduced in section 4. The numerical scheme is presented in section 5. The numerical results are shown in section 6. Technical details are in Appendices A and B.

Saint-Venant equations in polar coordinates

The Saint-Venant equations are :

∂h ∂t + ∇h • U + hdivU = 0, (1)
∂U ∂t + ∂U ∂x U + g∇h = 0. ( 2
)
Here t is the time, x = (x, y) T are the Cartesian coordinates, U = (U, V ) T is the depth averaged horizontal velocity, h is the fluid depth, g is the gravity. We want to rewrite this system in polar coordinates. Even if these equations are well known, it would be useful to present this standard technique before to pass to the full 2D non-conservative model of shear shallow water flows which will be rewritten in polar coordinates here for the first time. Let us introduce the polar coordinates (r, θ) ∈]0, +∞[×[0, 2π] related with the Cartesian coordinates (x, y) as follows :

x = r cos θ, y = r sin θ, r = √ x 2 + y 2 , θ = arctan ( y x

) . We consider the usual orthonormal curvilinear basis:

e r = (cos θ, sin θ) T , e θ = (-sin θ, cos θ) T .

We have :

U = U r e r + U θ e θ , ∂U ∂x = ∂U ∂r ⊗ ∇r + ∂U ∂θ ⊗ ∇θ, (3)
divU = tr ( ∂U ∂x ) , ∇h = ∂h ∂r ∇r + ∂h ∂θ ∇θ. (4)
Since ∇r = ( ∂r ∂x , ∂r ∂y

) T = ( x/ √ x 2 + y 2 , y/ √ x 2 + y 2 ) T = (x/r, y/r) T = e r , (5)
∇θ = ( ∂θ ∂x , ∂θ ∂y

) T = ( -y/(x 2 + y 2 ), 1/(x(1 + (y/x) 2 )) ) T = e θ r , (6)
then

∂U ∂x = ∂U ∂r ⊗ e r + 1 r ∂U ∂θ ⊗ e θ . (7)
One also has :

∂U ∂r = ∂U
( U r + ∂U θ ∂θ ) e θ ⊗ e θ . (10)
The same expression in the matrix form is :

∂U ∂x =       ∂U r ∂r 1 r ( ∂U r ∂θ -U θ ) ∂U θ ∂r 1 r ( ∂U θ ∂θ + U r )       . (11)
It implies :

divU = 1 r ∂U θ ∂θ + 1 r ∂ ∂r (rU r ) , (12)
∂U ∂x U = { U r ∂U r ∂r + ( 1 r ∂U r ∂θ - U θ r ) U θ } e r + { U r ∂U θ ∂r + ( 1 r ∂U θ ∂θ + U r r ) U θ } e θ . (13)
Also, one has :

∇h = ∂h ∂r ∇r + ∂h ∂θ ∇θ = ∂h ∂r e r + 1 r ∂h ∂θ e θ . (14)
The Saint-Venant equations in the polar coordinate system are then :

                       ∂h ∂t + U r ∂h ∂r + U θ r ∂h ∂θ + h r { ∂U θ ∂θ + ∂ ∂r (rU r ) } = 0, ∂U r ∂t + U r ∂U r ∂r + 1 r ( ∂U r ∂θ -U θ ) U θ + g ∂h ∂r = 0, ∂U θ ∂t + U r ∂U θ ∂r + 1 r ( ∂U θ ∂θ + U r ) U θ + g r ∂h ∂θ = 0. (15)
Or, in conservative form, with the energy conservation law :

                                 ∂(hr) ∂t + ∂(rhU r ) ∂r + ∂(hU θ ) ∂θ = 0, ∂(rhU r ) ∂t + ∂ ∂r { r ( hU 2 r + gh 2 2 )} + ∂(hU r U θ ) ∂θ = hU 2 θ + gh 2 2 , ∂(rhU θ ) ∂t + ∂(rhU r U θ ) ∂r + ∂ ∂θ ( hU 2 θ + gh 2 2 ) = -hU r U θ , ∂ ∂t { rh ( |U| 2 2 + gh 2
)} + ∂ ∂r { rhU r ( |U| 2 2 + gh 2
) + r gh 2 2 U r } + ∂ ∂θ { hU θ ( |U| 2 2 + gh 2
) + gh 2 2 U θ } = 0. ( 16
)
The same procedure will now be applied to the shear shallow water equations.

2D shear shallow water flows in polar coordinates

The governing 2D hyperbolic system of shear shallow water equations without right-hand side is :

                       ∂h ∂t + div (hU) = 0, ∂ (hU) ∂t + div ( hU ⊗ U + gh 2 2 I + hP ) = 0, DP Dt + ∂U ∂x P + P ( ∂U ∂x ) T = 0. ( 17
)
Here t is the time, x = (x, y) T are the Cartesian coordinates, h is the fluid depth, g is the gravity, D/Dt means the material derivative with respect to the mean motion :

D Dt = ∂ ∂t + U T • ∇ , (18)
P = P T is the stress tensor, U T = (U, V ) is the depth averaged velocity. Equations ( 17) admit the energy conservation law :

∂ ∂t { h ( 1
2 |U| 2 + e i + e T )} + div { hU ( 1
2 |U| 2 + e i + e T ) + ( gh 2 2 I + hP ) U } = 0, (19)
where

e i = gh 2 , e T = 1 2 tr(P)
and the conservation of "entropy" along trajectories:

D Dt ( det(P) h 2 ) = 0.
The system ( 17) is composed of the equations of mass balance, horizontal momentum and evolution equation for the stress tensor. The governing equations are hyperbolic, if P is positive definite.

The numerical resolution of this system involves some difficulties, because it is not in conservative form.

We take (17) in the following non-conservative form to rewrite it in polar coordinates:

                     ∂h ∂t + div (hU) = 0, ∂U ∂t + ∂U ∂x U + g∇h + 1 h div (hP) = 0, DP Dt + ∂U ∂x P + P ( ∂U ∂x ) T = 0. (20)
Again, the relation between Cartesian and polar coordinates are:

x = r cos θ, y = r sin θ, r = √ x 2 + y 2 , θ = arctan ( y x
) , e r = (cos θ, sin θ) T , e θ = (-sin θ, cos θ) T .

One has :

∂U ∂x = ∂U r ∂r e r ⊗ e r + ∂U θ ∂r e θ ⊗ e r + ( 1 r ∂U r ∂θ - U θ r ) e r ⊗ e θ + 1 r ( U r + ∂U θ ∂θ ) e θ ⊗ e θ . (21)
Since

div(hP) = (divP)h + P∇h, (22)
we need to calculate divP. In polar coordinates we have :

P = P
where k 1 = r, k 2 = θ. Let a = a r e r + a θ e θ . We obtain then :

div (Pa) =
Also, one has :

DP Dt = ∂P ∂t + ∑ i ( ∂P ∂k i ⊗ ∇k i ) u = ∂P ∂t + ∑ i ∂P ∂k i (∇k i • u). (27)
It implies :

DP Dt = ∂P ∂t + U r ∂P ∂r + 1 r U θ ∂P ∂θ , (28)
div ( gh 2 2 I ) = ∇ ( gh 2 2 ) = gh ( ∂h ∂r e r + 1 r ∂h ∂θ e θ ) , (29)
div ( gh 2 2 U + hPU ) = gh 2 2 div(U) + ghU • ∇h + div (hP) U + tr ( hP ∂U ∂x ) , (30)
     . (31)
Hence, the system (20) in polar coordinates is : )} = 0.

                                                                 ∂h ∂t + 1
(32) The last system of equations ( 32) can be rewritten in the following form:

∂W ∂t + A ∂W ∂r + B ∂W ∂θ = F, (33)
where

W =                   h U r U θ P rr P rθ P θθ                   , A =                     U r h 0 0 0 0 g + P rr h U r 0 1 0 0 P θr h 0 U r 0 1 0 0 2P rr 0 U r 0 0 0 P θr P rr 0 U r 0 (35)
The energy equation ( 19) can be rewritten in the following form:

∂ ∂t { h ( 1 2 |U| 2 + E )} + div(V) = 0, (36)
where

E = gh 2 + tr(P) 2 , (37)
and for any vector

V : div(V) = 1 r ∂(rV r ) ∂r + 1 r ∂V θ ∂θ . Here V =   V r V θ   =       hU r ( 1
2 |U| 2 + E ) + gh 2 2 U r + h(P rr U r + P rθ U θ ) hU θ ( 1 2 |U| 2 + E ) + gh 2 2 U θ + h(P rθ U r + P θθ U θ )       . (38)
More comfortable form of these equations is :

                                                                                     ∂(
( ∂U θ ∂θ + U r )} = 0, ∂ ∂t { hr ( 1 2 |U| 2 + E )} + ∂ ∂r { r [ hU r ( 1 2 |U| 2 + E ) + gh 2 2 U r + h(P rr U r + P rθ U θ ) ]} + + ∂ ∂θ { hU θ ( 1 2 |U| 2 + E ) + gh 2 2 U θ + h(P rθ U r + P θθ U θ ) } = 0.
(39) The system (39) admits the 'entropy' conservation law:

D Dt ( det(P) h 2 ) = 0, D Dt = ∂ ∂t + U r ∂ ∂r + U θ r ∂ ∂θ . (40)

Splitting technique in polar coordinates

The system [START_REF] Richard | The classical hydraulic jump in a model of shear shallow-water flows[END_REF] is naturally split into two subsystems in r and θ directions. Consider first the subsystem in r -direction:

                                                                       ∂(
) = 0, ∂ ∂t { rh ( 1 2 |U| 2 + E )} + ∂ ∂r { r [ hU r ( 1 2 |U| 2 + E ) + gh 2 2 U r + h(P rr U r + P rθ U θ ) ]} = 0. ( 41
)
The system is hyperbolic with the following eigenvalues of the matrix A given by (34):

λ 1,2 = U r , λ 3,4 = U r ± √ P rr , λ 5,6 = U r ± √ gh + 3P rr . (42)
The fields corresponding to the multiple eigenvalues λ 1,2 and to the eigenvalues λ 3,4 are linear degenerate in the sense of Lax, while the fields λ 5,6 are genuinely non-linear. The family λ 3,4 is reminiscent of the shear waves in hyperelasticity, while the family λ 5,6 is reminiscent of the longitudinal waves.

The sub-system in θ-direction is :

                                                                     ∂(rh) ∂t + ∂(hU θ ) ∂θ = 0, ∂(rhU r ) ∂t + ∂ ∂θ (hU r U θ + hP rθ ) = 0, ∂ (rhU θ ) ∂t + ∂ ∂θ ( hU 2 θ + gh 2 2 + hP θθ ) = 0, ∂P rr ∂t + U θ r ∂P rr ∂θ + 2 P rθ r ∂U r ∂θ = 0, ∂P rθ ∂t + U θ r ∂P rθ ∂θ + P θθ r ∂U r ∂θ + P rθ r ∂U θ ∂θ = 0, ∂P θθ ∂t + U θ r ∂P θθ ∂θ + 2 P θθ r ∂U θ ∂θ = 0, ∂ ∂t { rh ( 1
2 |U| 2 + E )} + ∂ ∂θ { hU θ ( 1 2 |U| 2 + E ) + gh 2 2 U θ + h(P rθ U r + P θθ U θ ) } = 0. ( 43
)
The system is hyperbolic with the following eigenvalues of the matrix B given by ( 35) :

λ 1,2 = U θ r , λ 3,4 = U θ ± √ P θθ r , λ 5,6 = U θ ± √ gh + 3P θθ r . (44)
Each subsystem ( 41) and ( 43) can also be split into two sub-systems, each of which contains only one family of sound characteristics, corresponding to transverse and longitudinal waves, respectively.

Two-wave splitting for (41)

The system of equations for 'longitudinal waves' of ( 41) is :

                                                                     ∂(rh) ∂t + ∂(rhU r ) ∂r = 0, ∂(rhU r ) ∂t + ∂ ∂r ( r ( hU 2 r + gh 2 2 + hP rr )) = h(U 2 θ + P θθ ) + gh 2 2 , ∂ (rhU θ ) ∂t + ∂ ∂r (rhU r U θ ) = -h (U r U θ + P rθ ) , ∂P rr ∂t + U r ∂P rr ∂r -4P rθ U θ r + 2P rr ∂U r ∂r = 0, ∂P rθ ∂t + U r ∂P rθ ∂r + U θ r (P rr -2P θθ ) + P rθ r ∂(rU r ) ∂r = 0, ∂P θθ ∂t + U r ∂P θθ ∂r + 2P rθ U θ r + 2 P θθ r U r = 0, ∂ ∂t { rh ( 1
2 |U| 2 + E )} + ∂ ∂r { r [ hU r ( 1 2 |U| 2 + E ) + gh 2 2 U r + hP rr U r ]} = 0. ( 45
)
This system of equations ( 45) can be rewritten in the form

∂W ∂t + A 1 ∂W ∂r = S r (W), (46)
(47) It is hyperbolic with the following eigenvalues of the matrix A 1 given by ( 47) :

λ 1,2,3,4 = U r , λ 5,6 = U r ± √ gh + 3P rr . (48)
The fields corresponding to the multiple eigenvalues λ 1,2,3,4 are linear degenerate in the sense of Lax, while the fields λ 5,6 are genuinely non-linear. It admits also the 'entropy' conservation :

D Dt ( P rr P θθ -P 2 rθ h 2 ) = 0, D Dt = ∂ ∂t + U r ∂ ∂r . (49)
The system of equations for 'transverse waves' of ( 41) is :

                                                                     ∂(rh) ∂t = 0, ∂(rhU r ) ∂t = 0, ∂ (rhU θ ) ∂t + ∂ ∂r (rhP rθ ) = 0, ∂P rr ∂t = 0, ∂P rθ ∂t + P rr ∂U θ ∂r = 0, ∂P θθ ∂t + 2 ∂U θ ∂r P rθ = 0, ∂ ∂t { rh ( 1
2 U 2 θ + 1 2 P θθ )} + ∂ ∂r (rhU θ P rθ ) = 0. ( 50
)
This system of equations ( 50) can be rewritten in the form

∂W ∂t + A 2 ∂W ∂r = 0, (51)      0 0 0 0 0 0 0 0 0 0 0 0 P rθ h 0 0 0 1 0 0 0 0 0 0 0 0 0 P rr 0 0 0 0 0 2P rθ 0 0 0                    . ( 52
)
The system is hyperbolic with the following eigenvalues of the matrix A 2 given by ( 52) :

λ 1,2,3,4 = 0, λ 5,6 = ± √ P rr . (53)
All fields are linear degenerate. The system admits also the 'entropy' conservation :

∂ ∂t ( P rr P θθ -P 2 rθ h 2 ) = 0. (54)
3.3 Two-wave splitting for (43)

The system of equations for 'longitudinal waves' of ( 43) is :

                                                                     ∂(rh) ∂t + ∂(hU θ ) ∂θ = 0, ∂(rhU r ) ∂t + ∂ ∂θ (hU r U θ ) = 0, ∂ (rhU θ ) ∂t + ∂ ∂θ ( hU 2 θ + gh 2 2 + hP θθ ) = 0, ∂P rr ∂t + U θ r ∂P rr ∂θ = 0, ∂P rθ ∂t + U θ r ∂P rθ ∂θ + P rθ r ∂U θ ∂θ = 0, ∂P θθ ∂t + U θ r ∂P θθ ∂θ + 2 P θθ r ∂U θ ∂θ = 0, ∂ ∂t { rh ( 1 2 |U| 2 + E )} + ∂ ∂θ { hU θ ( 1 2 |U| 2 + E ) + gh 2 2 U θ + hP θθ U θ } = 0. ( 55
)
This system of equations (55) can be rewritten in the form

∂W ∂t + B 1 ∂W ∂θ = 0, (56)
λ 1,2,3,4 = U θ r , λ 5,6 = U θ ± √ gh + 3P θθ r . (58)
The multiple eigenfields are linear degenerate, while the 'sound fields' are genuinely nonlinear. The system admits also the 'entropy' conservation :

D Dt ( P rr P θθ -P 2 rθ h 2 ) = 0, D Dt = ∂ ∂t + U θ r ∂ ∂θ . (59)
The system of equations for 'transverse waves' of ( 43) is :

                                                                     ∂(rh) ∂t = 0, ∂(rhU r ) ∂t + ∂ ∂θ (hP rθ ) = 0, ∂ (rhU θ ) ∂t = 0, ∂P rr ∂t + 2 P rθ r ∂U r ∂θ = 0, ∂P rθ ∂t + P θθ r ∂U r ∂θ = 0, ∂P θθ ∂t = 0, ∂ ∂t { rh ( 1 2 |U| 2 + E )} + ∂ ∂θ (hP rθ U r ) = 0. ( 60
)
This system of equations (60) can be rewritten in the form

∂W ∂t + B 2 ∂W ∂θ = 0, (61)
λ 1,2,3,4 = 0, λ 5,6 = ± √ P θθ r . (63)
All eigenfields are linear degenerate. The system admits also the 'entropy' conservation :

∂ ∂t ( P rr P θθ -P 2 rθ h 2 ) = 0. (64)

Numerical scheme

We solve first the model in r-direction, then in θ-direction. For each direction the two subsystems are solved separately. For each subsystem the following steps will be applied:

• Riemann problem

• Godunov type scheme

• Correction step

Each subsystem must have the following properties:

• to be hyperbolic

• to admit the energy conservation law compatible with the entropy equation

The 'conservative' variables for all subsystems are: W = [ rh, rhU r , rhU θ , rhP rr , rP rθ , rhP θθ , rh

( E + 1 2 |U| 2 )] T , (65)
with E = (gh + P rr + P θθ )/2.

Subsystem 1 in the 'r'-direction

The first subsystem (called also 'subsystem 1' for a-waves ) is:

                                                     ∂(rh) ∂t + ∂(rhU r ) ∂r = 0, (a) ∂(rhU r ) ∂t + ∂ ∂r ( r ( hU 2 r + gh 2 2 + hP rr )) = h(U 2 θ + P θθ ) + gh 2 2 , (b) ∂ (rhU θ ) ∂t + ∂ ∂r (rhU r U θ ) = -h (U r U θ + P rθ ) , (c)
+ P θθ U r ), (f ) (66)
The subsystem (66) admits the energy conservation law

∂ ∂t { rh ( 1 2 |U| 2 + E )} + ∂ ∂r { r [ hU r ( 1 2 |U| 2 + E ) + gh 2 2 U r + hP rr U r ]} = 0. ( 67
)
The equation (66.d) is not conservative. The jump relation for this equation is not well defined. The value of P rr will be corrected using the energy equation (67). Numerically, we solve first the system (67) where the equation (66.d) is replaced by (rhP rr ) t + (rhU r P rr ) r = 0.

(68)

We update then rhP rr by using (67).

Godunov type method and correction of the non-conservative term rhP rr

The subsystem 1 can be rewritten in the following form:

∂W ∂t + ∂(rF) ∂r = S r , (69)
where the vector of so-called "conservative variables" W, the vector of fluxes F(W) and the geometric source term vector S r (W) here are W = ( rh, rhU r , rhU θ , rhP rr , rP rθ , rhP θθ , rh Let us consider now a fixed grid of size ∆r = r i+1/2r i-1/2 , the time increment is defined as ∆t = t n+1t n that must respect the Courant-Friederichs-Lewy's (CFL) condition. The discrete values of the vector-function W(r, t) at (r i , t n ) will be denoted by

( E + 1 2 |U| 2 )) T , F(W) = ( hU r , hU 2 r + p, hU r U θ , hU r P rr , U r P rθ , hU r P θθ , hU r ( 1 2 |U| + E ) + pU r ) T , S r = ( 0, h ( U 2 θ + P θθ ) + gh 2 2 , -h (U r U θ + P rθ ) , 4hP rθ U θ , U θ (2P θθ -P rr ) -P rθ U r , -2h (P rθ U θ + P θθ U r ) , 0 ) (70
W n i ≡ W(r i , t n ).
Consider the initial condition for the complete problem W(r, t n ) = W n . Integrating in space and time I i × [t n , t n+1 ] the conservation laws (69) one obtains the conservative finite volume Godunov scheme on a fixed grid

Wn+1 i = W n i - ∆t ∆r ( r i+1/2 (F * ,n i+1/2 -S r,n i ) -r i-1/2 (F * ,n i-1/2 -S r,n i ) ) , (71)
where F * ,n i+1/2 and F * ,n i-1/2 are the numerical fluxes. They are constant across interfaces between cells during the time step. For computing the fluxes

F * ,n i+1/2 = F * ,n i+1/2 (W n i , W n i+1 ) and F * ,n i-1/2 = F * ,n i-1/2 (W n i-1 , W n i )
, we solve the Riemann problems between cells i, i + 1 and i -1, i respectively. Using the energy conservative law, the non-conservative term rhP rr is updated by the following way:

(rhP rr ) = 2(rhE) - g(rh) 2 r -(rhP θθ ). ( 72
)
4.2 Subsystem 2 in the 'r'-direction

The subsystem 2 (called also subsytem 2 for b) is: in equation (73.f). Since the jump relation is not well defined, there is no hope to obtain the good value of P θθ . In the following this equation will be replaced by ∂ (rhP θθ ) ∂t = 0 at the first step, and then the energy conservation law will be used for update the value of P θθ .

                                                     ∂(rh) ∂t = 0, (a) ∂(rhU r ) ∂t = 0, (b)
∂ (rhU θ ) ∂t + ∂ ∂r (rhP rθ ) = 0, (c) ∂ (rhP rr ) ∂t = 0, (d)
In particular, the estimation of P rθ at time instant (t + ∆t) is given by: (rP rθ ) (r, t + ∆t) = (rP rθ ) (r, t) + (rP rr )(r, t)

( U * θi-1/2 (t) -U * θi+1/2 (t) ) ∆t ∆r . (75)
Here U * θ it is speed at the contact discontinuity.

Discretization of equation for P rθ and correction of term hP θθ

The subsystem 2 in 'r'-direction can be rewritten in the following form:

∂W ∂t + ∂(rF) ∂r + (rK) ∂U θ ∂r = 0, (76)
with W = ( rh, rhU r , rhU θ , rhP rr , rP rθ , rhP θθ , rh

( E + 1 2 |U| 2 )) T , F = F(W) = (0, 0, hP rθ , 0, 0, 0, hU θ P rθ ) T , K = K(W) = (0, 0, 0, 0, P rr , 2hP rθ , 0) T . (77)
For the equations in the conservative form, we use the Godunov scheme. The global scheme is :

W t+∆t i = W t i + ∆t ∆r ( r i+1/2 F * ,t i+1/2 -r i-1/2 F * ,t i-1/2 ) + ∆t ∆r K t i+1/2 ( U * ,t θ i+1/2 -U * ,t θ i-1/2 ) (78)
The non-conservative equation for P rθ we discretize in the following way:

P rθ (r, t + ∆t) = P rθ (r, t) + P rr (r, t)

( U * θ i-1/2 (t) -U * θ i+1/2 (t) ) ∆t ∆r . (79)
Here U * θ is the speed of the contact discontinuity. Equation for (rhP θθ ) is also non-conservative. This equation is solved numerically first with zero flux, and then the energy conservation law was used to update the value of rhP θθ :

(rhP θθ ) = 2(rhE) - g(rh) 2 r -(rhP rr ). (80)

Subsystem 1 in the 'θ'-direction

The first subsystem in the 'θ'-direction for 'longitudinal a-waves' is :

                                                     ∂(rh) ∂t + ∂(hU θ ) ∂θ = 0, (a) ∂(rhU r ) ∂t + ∂ ∂θ (hU r U θ ) = 0, (b) ∂ (rhU θ ) ∂t + ∂ ∂θ ( hU 2 θ + gh 2 2 + hP θθ ) = 0, (c) ∂ (rhP rr ) ∂t + ∂ (hU θ P rr ) ∂θ = 0, (d) ∂(rP rθ ) ∂t + ∂(U θ P rθ ) ∂θ = 0, (e) ∂(rhP θθ ) ∂t + ∂(hU θ P θθ ) ∂θ + 2hP θθ ∂U θ ∂θ = 0. (f ) (81)
The subsystem (81) admits the energy conservation law :

∂ ∂t { rh ( 1 2 |U| 2 + E )} + ∂ ∂θ { hU θ ( 1 2 |U| 2 + E ) + gh 2 2 U θ + hP θθ U θ } = 0. ( 82
)
The system is hyperbolic with the eigenvalues given by :

λ 1,2,3,4 = U θ r , λ 5,6 = U θ ± √ gh + 3P θθ r . (83)
The equation (81.f) is not conservative. The jump relation are not well defined. The value of P θθ will be corrected using the energy equation (82).

Subsystem 2 in the 'θ'-direction

The subsystem 2 (called also subsystem 2 for b-waves) is :

                                                     ∂(rh) ∂t = 0, (a) ∂(rhU r ) ∂t + ∂ ∂θ (hP rθ ) = 0, (b) ∂ (rhU θ ) ∂t = 0, (c) ∂(rhP rr ) ∂t + 2hP rθ ∂U r ∂θ = 0, (d) ∂(rP rθ ) ∂t + P θθ ∂U r ∂θ = 0, (e) ∂(rhP θθ ) ∂t = 0. (f ) (84)
The system is hyperbolic with the eigenvalues given by :

λ 1,2,3,4 = 0, λ 5,6 = ± √ P θθ r . (85)
The subsystem (84) admits the energy conservation law ∂ ∂t to obtain the good value of P rr . In the following this equation will be replaced by ∂ (rhP rr ) ∂t = 0 and then the energy conservation law was used for update the value of P rr .

( rh ( 1
2 |U| 2 + E )) + ∂ ∂l (rhP rθ U r ) = 0. (86)
The estimation of P rθ at time instant (t + ∆t) is given by: (rP rθ ) (θ, t + ∆t) = (rP rθ ) (θ, t) + (P θθ )(θ, t)

( U * r i-1/2 (t) -U * r i+1/2 (t) ) ∆t ∆θ . (87)
Here U * r is a speed at the contact discontinuity.

Integrating the source terms

To add the source terms, we integrate the ordinary differential equations :

dW dt = S(W), (88)
with the full unknown vector W = ( rh, rhU r , rhU θ , rhP rr , rP rθ , rhP θθ , rh

( E + 1 2 |U| 2 )) T .
The initial condition W| t=0 = Wn+1 are obtained from the previous splitting steps. Here the vector of source terms is

S(W) = ( 0, -grh∇b -C f r √ U 2 r + U 2 θ U r , -rC f √ U 2 r + U 2 θ U θ , rhD rr , rD rθ , rhD θθ , -grh∇b • U r -rC f ( √ U 2 r + U 2 θ ) 3 -rQ ) T , (89)

Nonstationary exact solution

We look for the solution in the form :

h = h 0 f h (t), U r = rf u (t)
, U θ = 0, P rr = P 0 f P (t), P rθ = 0, P θθ = P 0 f P (t).

Substituting the relations in the governing shear shallow water equations and resolving the corresponding ordinary differential equations, we obtain the following analytical solution:

                                               h(r, θ, t) = h 0 (1 + t/T ) 2 , U r (r, θ, t) = r T + t , U θ (r, θ, t) = 0, P rr (r, θ, t) = P 0 (1 + t/T ) 2 , P rθ (r, θ, t) = 0, P θθ (r, θ, t) = P 0 (1 + t/T ) 2 .
(92)

Here P 0 and T are constants. In fact, this is a particular solution of a more general solution presented in Appendix B which corresponds to a symmetric matrix A 0 = I T . We choose T equal

The errors for U r and P rr were calculated in the same way. Different regular cylindrical meshes were used : (100 × 100, 200 × 200, 400 × 400, 500 × 500) with a first order Godunov method. This test shows the convergence to the exact solution. The slope of the convergence lines are almost the same for the depth, normal velocity and normal stress tensor component.

Radial Dam-break problem

We solve here the governing equations without dissipative terms (just a geometrical source terms is applied). By analogy with "shock tube problem" of gas dynamics, that illustrates the behavior of solutions to conservation laws, here we present the test case to simulate radial dam break problem. The result shows the shock wave propagating radially outwards, across which all of the state variables are discontinuous. It is followed by a contact discontinuity, across which the depth is again discontinuous. The third wave moves in the opposite direction (inwards). This is a rarefaction wave. We demonstrate here that our splitting scheme in polar coordinates is able to reproduce the correct shock wave front, contact discontinuity and rarefaction wave. The numerical schemes are implemented in the following computational domain :

D = { (r; θ)|R -= 0.005[m] ≤ r ≤ R + = 1[m], 0 ≤ θ ≤ 2π } .
The initial conditions are :

Radial hydraulic jump: comparison with SWASI experiment

Recently, Foglizzo et al. (CEA-Saclay, France) [1], [2], [3] have performed an interesting experimental study in which the circular hydraulic jump was seen to form, separating convergent shallow supercritical flow from deeper subcritical flow. This jump is similar to the circular hydraulic jump observed in a kitchen sink when a vertically falling water jet strikes a horizontal plate and then spreads radially outwards [4], [5], [6], [7], [8], [9], [10], [11], [12]. The main difference in the experiment by T. Foglizzo et al. is that the fluid is radially injected inwards and strikes a hollow cylinder by forming a hydraulic jump.

In these experiments, a hydraulic jump exhibits the appearance of radial asymmetries: after its formation, the hydraulic jump starts to oscillate randomly, then the amplitude of oscillations quickly grows and finally, the flow becomes very asymmetric and develops a angular rotation (see Figure 1). This phenomenon is an analogue of the instability occurring in the supernova core. The analogy stems from the fact that the shallow water equations are mathematically equivalent to the gas dynamics equations. The surface gravity waves represent the analogue of sound waves, and the hydraulic jumps are the analogue of shocks. The authors call their installation the SWASI: shallow water analog of a shock instability.

We use here the bottom topography different from that used in T. Foglizzo et al. (see Figure 6 ):

b(r) =        A [ (r -R --L 1 ) 2 -L 2 1 ] 2 L 4 1 , if 0 < r -R -< 2L 1 (r -R --2L 1 ) tan β, if 2L 1 + R -< r < R + , (96)
Here β is a small inclination angle. The bump models the hollow cylinder. Now, we need to impose a boundary condition at the bump. The variation of the bump amplitude allows us first to decelerate the water flow (the flow becomes sub-critical somewhere between r = R -and r = R + ) by creating the hydraulic jump, and then to accelerate the flow (the flow becomes supercritical at r = R -).

To simulate it numerically, we implemented the numerical schemes described above in a computational domain R -≤ r ≤ R + , 0 ≤ θ ≤ 2π.

The initial conditions are :

h(r, θ, t = 0) =    1.5h 0 [m]; if r -R -≤ 0.3L r = 0.3(R + -R -), ∀θ h 0 [m];
else ∀θ U r (r, θ, t = 0) = -q 0 /(rh), U θ (r, θ, t = 0) = 0, P rr = φh 2 , P rθ = 0, P θθ = ϵφh 2 .

Here ϵ = 10 -8 . The other parameters are shown in Table 1. Obviously, we respect the inequality 2L 1 < 0.3L r . First, we control that at the inner boundary r = R -in radial direction the corresponding Froude number is larger than one:

F g = U r √ gh + 3P rr > 1.
In this case, we can use Neumann boundary condition at the outlet corresponding to r = R -. In Figure 7 we show 1D space variation of the fluid depth and the Froude number at time instant 0 s and 100 s.

At the inlet (r = R + ) boundary we apply the following conditions h = h 0 , U r = -q 0 (1 + 0.01 sin(16θ))/(R + h 0 )), U θ = 0, P rr = P θθ = 0.5φh 2 0 , P rθ = 0.

So, the radial velocity is not uniform, we add some perturbation of 1% and perturb the flow in θ-direction. The reason to take the perturbation of the form sin(16θ) comes from the fact that in the above mentioned experiments 16 water pumps were used to inject water. The periodic TEST boundary condition is used in θ-direction. The numerical scheme is based on a splitting technique described in the previous sections. The first-order Godunov scheme with HLLC Riemann solver was used. The corresponding friction and dissipative terms were added using a first-order Euler method. The required values of parameters are given in Table 1. The code is parallelized using the Message Passing Interface (MPI). First, in Figure 8 one can see the formation of hydraulic jump that radially oscillates. Then, in Figures 9, 10 the hydraulic jump is destabilized into a pattern which radially oscillates but these oscillations are random in amplitude depending on θ. Remarkably, in a certain parameter domain, as the amplitude of the oscillations grows, a new free surface pattern appears (an angular point at the free surface of the hydraulic jump) which starts to rotate (see Figures 11,12,13). The sense of rotation can be changed by a small initial perturbation. Indeed, if the initial condition for the tangential velocity are taken in the form : U theta = ±ϵr (ϵ is a small positive parameter), the sense of rotation will be determined by the sign of U theta .

h 0 [m] β [rad] C r C f φ [s -2 ] q 0 [L/s] R + [m] R -[m] A L 1 [m
One can also see in Figure 13 (schlieren image) the emergence of double jump structures, where the free surface changes abruptly.

One can notice that the singular (angular) point appears only in some region of parameters φ and C r (see 2). We will call such a phenomenon "apple instability" by analogy with the place where the apple stem is attached and the geometrical form of singularity. The existence of such a point allows us to easily observe the period of rotation of pattern. In particular, in our case the period of rotation was about 6 seconds. It does not depend on the values of φ and C r . In the experiments of Foglizzo et al. the period was was about 3 s, but the geometry was different : the radius R + was 0.33 m. The main reason why we did not use the same geometrical characteristics is that we could not guarantee the supercritical flow at the outlet (at r = R -) which was necessary to assure correct boundary conditions for the governing equations.

The results obtained are qualitatively similar to those observed experimentally. It means that our mathematical model of shear shallow water flows and the corresponding numerical method are capable of simulating the jump formation and bifurcation from the radial oscillations to the appearance of the transverse oscillations.

The transverse structures (formation of the "triangles" behind the moving singular point) that one can see in our numerical simulations in the Figures 11,12,13,14 is a physical phenomenon and not a numerical instability. They could be also seen in the experiments. When we put φ = 0, we do not obtain any transverse structures. Recall, that φ is associated to the enstrophy of small vortexes at the vicinity of the bottom. In Table 2 we studied how this parameter φ and the dissipation coefficient C r have an influence on the flow behavior. One can see that for φ bigger than 10 the rotating angular point does not form anymore.

In our model, the quantity det(P) h 2 plays role of "entropy" because it is transported along the mean flow as the true entropy for the Euler equations of compressible fluids, and increase through the shocks. It turns out that this condition is sufficient to recognize precisely those discontinuities that are physically correct and specify a unique solution. The evolution of this quantity is shown in the Figures 891011121314.

A surprising fact also mentioned in the experiments that the sense of rotation of the angular point and the tangential fluid velocity in the inner region is opposite, is also described by our model (see

Conclusion

In this paper, some numerical results in polar coordinates are presented and compared to the analytical solutions and experimental observations. These validations demonstrated the capability of the model and numerical method to reproduce the shock wave dynamics appearing in multi -dimensional physical problems. In particular, in some region of parameters φ and C r , the formation of rotating singular point (angular point) at the circular hydraulic jump was found numerically.

The numerical method is based on a new splitting approach for non-conservative systems of equations.

We hope that our model of shear "shallow water" flows can be useful for studying the hydrodynamic aspects of supernova physics.

Appendix B

Here we present some analytical solutions to the system describing shear shallow water flows. In the Cartesian coordinates the solution is [17] :

U = A(t)x, dh dt = -htr(A),
A(t) = A 0 (I + A 0 t) -1 , P(t) = (I + A 0 t) -1 P 0 (I + A T 0 t) -1 . Here A 0 and P 0 are constant matrices. The solution is well defined for all t > 0, if the matrix I+A 0 t is invertible for any t. This is a case, for example, of antisymmetric matrix A 0 : A T 0 = -A 0 , or symmetric positive definite matrix : A T 0 = A 0 > 0. If A 0 is antisymmetric :

A 0 =   0 β -β 0   ,
we obtain :

A = β 1 + β 2 t 2   βt 1 -1 βt   .
Here β is a constant. Let (r, θ) be the polar coordinates. The relations between the Cartesian and polar bases are :

                                 h = 1 1 + β 2 t 2 , U r = rβ 2 t 1 + β 2 t 2 , U θ = - rβ 1 + β 2 t 2 , P = I 1 + β 2 t 2 , if P 0 = I.
where F * ,n i+1/2 is the numerical flux function, which prescribes order of accuracy, with the time average of the flux at the edge x i+1/2 of the cell. Note that taking ∆ i ≡ 0 for all i and n recovers first order Godunov's method.

Fig. 1 . 2 -

 12 Fig. 1.2 -Comparison of roll-wave profiles for model (1.1.2)(solid line), Saint-Venant model (dashed line) with Brock's experimental results (dots) for two different sets of experimental parameters.

1 • Chapter 2 :• Chapter 4 :

 124 to propose a new numerical method for solving of multi-dimensional non-conservative CHAP. Introduction 9 hyperbolic equations of shear shallow water flows • to demonstrate the capability of the model and numerical method to reproduce the shock wave dynamics appearing in physical problems by comparing numerical results with the experiments and analytical solutions Outline and contributions The thesis is devided into three main chapters. The chapters are based on articles which are published, submitted or under preparation: Formation and coarsening of roll-waves in shear shallow water flows down an inclined rectangular channel • Chapter 3: Multi-dimensional shear shallow water flows : problems and solutions Convergent circular hydraulic jump: formation of transverse structures Definition 2. "Roll waves" are periodic series of smooth profiles separated by hydraulic jumps, occurring in turbulent water flows down an inclined open channel due to a natural instability (an example is presented in the Figure 1.3).In Chapter 2, we present a numerical study of the generation and stability of roll-waves (see Definition 2) developed from uniform unsteady flow in inclined channel. This subject is of great interest in the field of channel hydraulics. Two different approaches (one is for a long channel, and the other one is for a "periodic box") are used for the system of one-dimensional SSWM. The SSWM is hyperbolic and reminiscent of the non-isentropic Euler equations. The role of entropy is played by a new variable, called entropy (squared vorticity). In particular, two types of enstrophies can be distinguished in the model: a small-scale enstrophy, generated near the wall, and a large-scale enstrophy corresponding to the flow in the roller region near the free surface. The source terms in the SSWM describes bed topography, friction forces and dissipation process in the roller. The ability of the SSWM to predict the formation of roll-waves

Fig. 1 . 3 -

 13 Fig. 1.3 -Series of roll waves moving down a spillway.

Fig. 1 . 4 -

 14 Fig. 1.4 -The roll-wave formation (at the left) and a magnified view (at the right) are shown both for the wave profile and the enstrophy (Φ + ϕ).

Fig. 1 . 5 -

 15 Fig. 1.5 -The space-time diagram showing the coarsening in a long open channel of 200 m long. The trajectories of wave crests are shown by white lines. The coarsenning corresponds to the intersection of these lines (left part of the computation domain).

Fig. 1 . 6 -

 16 Fig. 1.6 -The numerical solution of a single roll wave to the simplified 2D model in a periodic box. The final stationary solution does not depend on transverse perturbations and corresponds to that obtained in one-dimensional case. The simplified 2D model is not able to describe transverse profile.

(

  a-waves), and shear waves propagating with a smaller velocity u ± b, b 2 = P 11 (b-waves). The eigenfields corresponding to the contact discontinuity and b-waves are linearly degenerate in the sense of Lax, while the eigenfield corresponding to a-waves is genuinely non-linear in the sense of Lax. The idea is to split (1.3.1) into two subsystems treating separately a-waves and b-waves. The subsystem for a-waves reads :

Fig. 1 . 7 -

 17 Fig. 1.7 -Formation of a transverse structure at the jump toe perimeter consisting of seven waves for the computation domain 1.3 m long and 0.5 m wide.

Fig. 1 . 8 - 1 Fig. 1 . 9 -

 18119 Fig. 1.8 -Formation of a transverse structure of the jump toe perimeter consisting of fourteen waves for the computation domain 1.3 m long and 1 m wide.

Fig. 1 . 10 -

 110 Fig. 1.10 -2D simulation of the hydraulic jump. The profiles of h, U θ and det(P) h 2 are shown. A singular point rotating in positive direction is formed. Transverse triangular structures appearing at the front follow the singular point.

Figure 1 :

 1 Figure 1: Sketch of a typical profile of roll waves down an inclined plane of angle θ in Brock's experiments (1967).

Figure 2 :

 2 Figure 2: The Riemann problem is solved between neighbouring cells.

t 2 Figure 3 :Figure 4 :Figure 5 :

 2345 Figure 3: A computed solution in the long channel of 40 m for the case of F r g0 = 1.15 is shown at time instant 100 s. The number of cells N = 8000. The roll waves are not formed as it is expected from the theory. The perturbation moves downstream the channel with rapidly decreasing amplitude.

2 Figure 6 :

 26 Figure 6: Comparison of roll-wave profiles for model (1)-(3) (solid line), SV model (dashed line) with Brock's experimental results (dots) for the Cases 1 and 2.

Figure 7 :Figure 8 :

 78 Figure 7: The profile of a roll-wave train for ω = 2.0 [rad/s] at the time instant 1000 s in the channel of 600 m long is at the left, and a magnified view is at the right. The other parameters correspond to Case 1. The length of a single roll wave is about 5 m. Number of cells is 60000. The roll-wave train is composed of single roll-waves having the same length and velocity. The wave amplitude is strongly modulated.

Figure 9 :Figure 10 :

 910 Figure 9: A regular roll-wave train is formed when the Saint-Venant equations are used (at the left). The results are obtained for ω = 2.0 [rad/s] at the time instant 1000 s in the channel of 600 m long. The other parameters correspond to Case 1. A magnified view is at the right. The length of a single roll wave is about 5 m. Number of cells is 60000. The modulations are absent.

Figure 11 :

 11 Figure 11: The space-time diagram showing the coarsening in a long open channel of 200 m long. The trajectories of wave crests are shown by black lines. The coarsening corresponds to the intersection of these lines (left part of the computation domain ).

  1 and Case 2, respectively). The comparison of profiles of a single periodic roll wave with that measured by Brock is shown for the model RG (1)-(

Figure 12 :

 12 Figure 12: Evolution of the numerical solution in the periodic box of L b = 1.3 m . The generalized Froude number is F r g0 = 1.15 < 2, a = 0.05. The number of cells N = 1000. The roll-wave is not formed as it was expected.

Figure 13 :

 13 Figure 13: Evolution of the numerical solution in the periodic box of 1.3 m (at the left) and 1.8 m long (at the right) with the parameter values corresponding to Case 1 and Case 2, respectively. Number of cells N = 1000. The numerical solution becomes stationary very rapidly.

Figure 14 :

 14 Figure 14: Large-time behaviour of numerical solutions in the periodic box does not depend on initial perturbations. A single roll-wave of length 1.3 m was formed. Number of cells N = 1000, F r g0 = 3.63. The initial conditions for each case for the depths was taken in the following form: a) h(x, t = 0) = h 0 (1 + a sin(8πx/L b )) , 0 ≤ x ≤ L b ; b) h(x, t = 0) = h 0 (1 + 2a sin(4πx/L b ) + a sin(2πx/L b )) , 0 ≤ x ≤ L b .

Figure 15 :

 15 Figure 15: The average depth and discharge in the periodic box of 1.3 m are shown. The flow parameters correspond to Case 1. Number of cells N = 3000, F r g0 = 3.63.

Figure 16 :

 16 Figure 16: The instability scenario in the periodic box of different lengths multiple of L b = 1.3 m. The parameters correspond to Case 1. Number of cells N = 4000 for all computations. A stable single wave exists only for the periodic box of length 8 L b (RG-model) and 11 L b (SV model). For larger lengths, the waves becomes very non-stationary and finally bifurcate into two wave configuration for the box of 14 L b (12 L b for the SV model). The Figures are shown atthe same time instant 1000 s.

Figure 17 :Figure 18 :

 1718 Figure 17: The initial data for the fluid depth are given by (21). The variation of the initial depth is shown here for Case 1.

4 Figure 19 :

 419 Figure19: The mesh convergence study of the numerical method for the periodic box of 1.3 [m] with the data corresponding to Case 1. In this Figure the logarithm of the error is shown as a function of the grid size with 100, 200, 400, 500, 1000, 2000, 4000 and 8000 mesh cells (dots). The "exact solution" was replaced by a numerical solution with 16000 points. The 1.4 convergence slope is shown with a continuous line.

Then Q = C r |u| 3

 3

Figure 1 :

 1 Figure 1: HLLC approximate Riemann solver for subsystem 1 describing a-waves. The solution in the 'star' region consists of two constant states separated by a middle wave of speed S a 0 = u * . The fastest and slowest waves denoted by S a L and S a R , respectively, are estimated by using Davis' approximation.

7. 2 . 1 HLLCFigure 2 :

 212 Figure 2: HLLC approximate Riemann solver for subsystem 2 describing b-waves. Solution in the 'star' region consists of two constant states separated by a middle wave of speed S b 0 = 0. The fastest and slowest waves S b R and S b L , respectively, are estimated by using Davis' approximation.

Figure 3 :

 3 Figure 3: Shear test problem. The initial discontinuity is located at x = 0.5 m. The tangential velocity is 0.2 m/s on the left, and -0.2 m/s on the right. The tangential velocity, shear stress P 12 and stress P 22 are shown at time instant 10 s for 500, 1000 and 10000 grid cells. CFL number is 0.3.

Figure 5 :

 5 Figure 5: Convergence test in dam-break problem: the initial discontinuity of the depth ( 0.02 m at the left and 0.01 m at the right) is located at x = 0.5 m. The fluid depth, normal stress component P 11 and normal velocity u are shown at time instant 0.5 s with 100, 1000 and 10000 grid cells. The MUSCL extension of Godunov method is used. CFL number is 0.8.

  y, 0) = 0, P 11 (x, y, 0) = P 22 (x, y, 0) = ϕh 2 (x, y, 0) 2 , P 12 = 0.Here θ = 0.05011 [rad] is the inclination angle, C f = 0.0036 is the Chézy coefficient, h 0 = 7.98 × 10 -3 m, a = 0.05, ϕ = 22.76 s -2 , g = 9.81 ms -2 , C r = 0.00035, L x = 1.3 m. For the flow in y-direction the changes in initial conditions are obvious.

Figure 7 :

 7 Figure 7: The convergence lines for the stress tensor components as functions of the grid cell size are shown. The number of grid cells was 100 × 100, 200 × 200, 400 × 400, 800 × 800, 1600 × 1600. The convergence to the exact solution is clearly visible.

P 11 (

 11 x, y, 0) = P 22 (x, y, 0) = ϕh 2 (x, y, 0) 2 , P 12 = 0.

Figure 8 :

 8 Figure 8: Formation of a transverse structure of the jump toe perimeter consisting of seven waves for the initial data (76) with m = 1 and k = 1. The convergent solution is obtained for the domain 1.3 m long and 0.5 m wide and shown at time instant 39 s. The Godunov method was used with 600 × 600 grid cells for the parameter set mentioned above. CFL number is 0.1.

Figure 10 :

 10 Figure 10: Top view for the schlieren image of the fluid depth, horizontal velocity, and normal stress components P 11 and P 22 for the initial perturbation (76) with k = 1, m = 1 (on the left) and k = 1, m = 4 (on the right) at time instant 39 s. The computational domain is 1.3 m long and 0.5 m wide. The Godunov method was used, with 600 × 600 grid cells for the parameter set mentioned above. CFL number here is 0.1. Seven transverse waves are observed on the left, and almost eight waves on the right.

Figure 11 :

 11 Figure 11: Formation of a transverse structure of the toe perimeter consisting of fourteen waves. The convergent solution is obtained for the domain 1.3 m long and 1 m wide, and shown at time instant 39 s. The Godunov method was used with 600 × 1200 grid cells for the parameter set mentioned above. CFL number is 0.1.

Figure 12 : 3 depthFigure 13 :

 12313 Figure 12: Top view for the schlieren image of the fluid depth h, horizontal velocity u and normal stress components P 11 and P 22 are shown at the time instant 39 s. The computational domain is 1.3 m long and 1 m wide. The Godunov method was used, with 600×1200 grid cells for the parameter set mentioned above. CFL number here is 0.1. Fourteen transverse waves are observed.

Figure 1 :

 1 Figure 1: The shape of the rotating hydraulic jump observed in the SWASI experiment. The rotating singular (angular) point occurs as a result of non-linear behavior of the free surfaces.

Figure 2 :

 2 Figure 2: A sketch of polar coordinates.

  rr e r ⊗ e r + P rθ e r ⊗ e θ + P θθ e θ ⊗ e θ + P θr e θ ⊗ e r(23) As P = P T , we obtain that P θr = P rθ . div (Pa) = tr ( ∂(Pa) ∂x) = div (P) a,for any a ≡ const. Then ∂Pa ∂x = ∂Pa ∂k i ⊗ ∇k i ,

) with p = gh 2 / 2 +

 22 hP rr , E = (gh + P rr + P θθ )/2.

ω = 0. 1 ,Figure 3 :

 13 Figure 3: A comparison between analytical and numerical solution is shown for the test-case "water glass on a record player" at time instant t = 50 s for the tangential velocity U θ and depth h with values of ω = 0.1 [s -1 ] ( on the left, with 100 mesh cells), and ω = 1 [s -1 ] (on the right, with 100, 200 and 500 mesh cells). HLLC Riemann solver was used. CF L = 0.8, h 0 = 0.005 [m].

  [m]; if R -≤ r ≤ 0.5L r = 0.5(R + -R -), ∀θ 0.01 [m]; if 0.5L r = 0.5(R + -R -) ≤ r ≤ R + , ∀θ U r (r, θ, t = 0) = 0, U θ (r, θ, t = 0) = 0, P = 0.At the internal (r = R -) and external (r = R + ) boundaries we apply the wall type boundary conditions U r (r = R -, θ, t) = U r (r = R + , θ, t) = 0. The periodic boundary condition is used in θdirection. This situation corresponds to the spreading of a cylindrical liquid column which was initially at rest. Figure (5) is the 2D visualization of the radial dam break problem simulation. When the dam breaks, a shock wave travels radially outwards while rarefaction wave moves inwards. Time: 0.00 s Time: 0.60 s

Figure 5 :

 5 Figure 5: Time evolution of the water height is shown with CF L = 0.8, R -= 0.005 m, R + = 1 m, the number of mesh cells 500 × 500.

Figure 6 :

 6 Figure 6: The bottom topography used in numerical simulations.

Figure 7 :

 7 Figure 7: The 1D depth profile and corresponding Froude number F r = U r √ gh + 3P 11 for 500 mesh cells at time instant 0 s (on the left) and 100 s (on the right). CFL = 0.8, C r = 1, φ = .

Figure 8 :

 8 Figure 8: 2D simulation of the hydraulic jump. The profiles of h, U θ and det(P) h 2 are shown. The first order Godunov method with HLLC Riemann solver is used with 500 × 500 mesh cells, CFL = 0.45, C r = 1, φ = 2. Here time instant is 25 s.

Figure 9 :

 9 Figure 9: 2D simulation of the hydraulic jump. The profiles of h, U θ and det(P) h 2 are shown. The first order Godunov method with HLLC Riemann solver is used with 500 × 500 mesh cells, CFL = 0.45, C r = 1, φ = 2. Here time instant is 50 s. The hydraulic jump is destabilized into a pattern which radially oscillates but these oscillation are random in amplitude depending on θ.

Figure 10 :

 10 Figure 10: 2D simulation of the hydraulic jump. The profiles of h, U θ and det(P) h 2 are shown. The first order Godunov method with HLLC Riemann solver is used with 500 × 500 mesh cells, CFL = 0.45, C r = 1, φ = 2. Here time instant is 52 s. The pattern is similar to that obtained in the time instant 50 s, and is almost a symmetric reflection of the last one.

Figure 11 :

 11 Figure 11: 2D simulation of the hydraulic jump. The profiles of h, U θ and det(P) h 2 are shown. The first order Godunov method with HLLC Riemann solver is used with 500 × 500 mesh cells, CFL = 0.45, C r = 1, φ = 2. Here time instant is 89 s. A singular point rotating in positive direction is formed. Transverse triangular structures appearing at the front follow this singular point.

Figure 12 :

 12 Figure 12: 2D simulation of the hydraulic jump. The profiles of h, U θ and det(P) h 2 are shown. The first order Godunov method with HLLC Riemann solver is used with 500 × 500 mesh cells, CFL = 0.45, C r = 1, φ = 2. Here time instant is 92 s. One can see that the singular point indeed rotates in positive direction.

Figure 13 :

 13 Figure 13: 2D simulation of the hydraulic jump. The profiles of h, U θ and det(P) h 2 are shown. The first order Godunov method with HLLC Riemann solver is used with 500 × 500 mesh cells, CFL = 0.45, C r = 1, φ = 2. Here time instant is 150 s. One can see the emergence of a double jump structures, where the free surface changes abruptly. The effective radius of the hydraulic jump essentially increased.

12 A 21 A 22  Q - 1 ,-

 1221221 cos θe x + sin θe y , e θ =sin θe x + cos θe y ,    e x = cos θe rsin θe θ , e y = sin θe r + cos θe θ , A is linear application, that is by definition : A(u, v) = u • Av. We have : A = A 11 e x ⊗ e x + A 12 e x ⊗ e y + A 21 e y ⊗ e x + A 22 e y ⊗ e y .For example, A 12 = e x • Ae y . The matrix A in the polar coordinates is  A rr A rθ A θr A θθ sin θ cos θ   is the rotation matrix, Q -1 = Q T , A ij , i, j = 1,2are the components of matrix A in the Cartesian coordinates. Also we have x = re r . Finally, our solution is :

Table 1 :

 1 Hydraulic parameters for two numerical tests corresponding to Brock's measurements

	r g0 a

  u t + uu x + vu y + gh x + P 11t + uP 11x + vP 11y + 2P 11 u x + 2P 12 u y = 0, P 12t + uP 12x + vP 12y + P 12 (u x + v y ) + P 11 v x + P 22 u y = 0,

		1 h	(hP 11 ) x +	1 h	(hP 12 ) y = 0,
	v t + uv x + vv y + gh y +	1 h	(hP 12 ) x +	1 h	(hP 22 ) y = 0,

  11 Dt + 2P 11 u x + 2P 12 u y = D 11 , DP 12 Dt + P 12 (u x + v y ) + P 11 v x + P 22 u y = D 12 , DP 22 Dt + 2P 12 v x + 2P 22 v y = D 22 .

	h 2 D Dt	12 P 11 P 22 -P 2 h 2

Here D ij , i, j = 1, 2 are the components of D. These equations imply the evolution equation for Ψ : = D 11 P 22 -2P 12 D 12 + D 22 P 11 .

  . The convergence is clearly visible.Figure 4: Shear test problem: initially, the discontinuity of the tangential velocity (0.2 m/s on the left and -0.2 m/s on the right ) is at x = 0.5 m. The 'entropy', depth and normal velocity are shown at time instant 10 s for 500, 1000 and 10000 grid cells. CFL number is 0.3.

	(P 11 P 22 -P 12 2 )/h 2 h (m)	0 0.01 -0.01 0 0.01 0.02 0.03 0.04 0.02 0.05 0.012 0.018 0.014 0.016	0	0.1 0.1 p p p	0.2 0.2	0.3 0.3	0.4 0.4	0.5 x (m) x (m) 0.5	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9 100 p 1000 p 10000 p	1
	h u (m/s)	-0.01 0 0 0.02 0.01 0.02 0.03 0.04 0.14 0.05 0.04 0.06 0.08 0.1 0.12	0 0		0.1 0.1 p p p	0.2 0.2	0.3 0.3	0.4 0.4	0.5 x (m) x (m) 0.5	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9 100 p 1000 p 10000 p	1
	u P 11 (m 2 2 ) /s	0.01 0.02 0.03 0.04 0.008 0.05 0.002 0.003 0.004 0.005 0.006 0.007		p p p								100 p 1000 p 10000 p
		0.001 0 0 -0.01 0	0	0.1 0.1	0.2 0.2	0.3 0.3	0.4 0.4	0.5 0.5 x (m) x (m)	0.6 0.6	0.7 0.7	0.8 0.8	0.9 0.9	1

  U θ + P rθ )) = -h (U r U θ + P rθ ) ,

										rh) ∂t	+	∂(rhU r ) ∂r	= 0,
	∂(rhU r ) ∂t	+	∂ ∂r	( r	(	hU 2 r +	gh 2 2	+ hP rr	))	= h(U 2 θ + P θθ ) + gh 2 2 ,
		∂ (rhU θ ) ∂t (rh (U r ∂P rr + ∂ ∂r ∂t + U r ∂P rr ∂r -4P rθ	U θ r	+ 2P rr	∂U r ∂r	= 0,
	∂P rθ ∂t	+ U r	∂P rθ ∂r	+	U θ r	(P rr -2P θθ ) + P rr	∂U θ ∂r	+	P rθ r	∂(rU r ) ∂r	= 0,
	∂P θθ ∂t	+ U r	∂P θθ ∂r	+ 2P rθ	U θ r	+ 2	( ∂U θ ∂r	P rθ +	P θθ r	U r

Table 1 :

 1 Parameters for numerical tests that we use in this chapter.

	]

Table 2 :

 2 Apple instability. In this table, 'YES' means that for corresponding values of φ and C r , the rotating singular (angular) point is formed, 'No' means that this point is invisible (at least in the schlieren image).'YES-No' means the limit behaviour.

The evolution of P (the last equation of (1)) is not governed by any specific objective derivative, i.e. the equation is not invariant under the change of variables t = t, x = O(t)x, u = O(t)u + Ȯx, P = OPO T . Here O(t) is a time-dependent orthogonal transformation : OO T = I, 'T' means transposition, and 'dot' denotes the time derivative. This is due to the fact that the tensor equation for the Reynolds stress tensor is not a geometric equation, but a physical one, representing a sort of 'micro' Newton's law derived from the Euler equations by depth averaging. Thus, this equation should be only Galilean invariant, which is obviously the case.

ACKNOWLEDGEMENTS

The present work could never be accomplished by the author herself. Herein I would like to express my gratitude to the people who were involved in this project.

Acknowledgement

K. Ivanova and S. Gavrilyuk were partially supported by ANR BoND, France. The authors thank S. Hank, N. Favrie and J. Massoni for help with domain decomposition methods.

Acknowledgment

The authors thank Boniface Nkonga, Jacques Massoni and Sarah Hank for many helpful comments. S.L. Gavrilyuk and K. Ivanova has been partially supported by the ANR project BoND (ANR-13-BS01-0009-01).

Chapter 4

CONVERGENT CIRCULAR HYDRAULIC JUMP: FORMATION OF TRANSVERSE STRUCTURES

This chapter corresponds to an article in preparation.

where

1 r ∂b ∂θ

) .

(90) Here b(r, θ) is the bottom topography. The equations are written in the reference frame where the gravity is orthogonal to the (r, θ)-plane. So, the case of a mild slop bottom is considered.

As the stress tensor P are positive definite, we have restriction P rr > 0, P θθ > 0, det(P rr P θθ -P 2 rθ ) > 0.

Even if we have 7 equations for 6 variables, the equations are compatible: the energy equation is a consequence of the equations of mass, momentum and stress tensor.

The first order Euler method is used for such integration.

Numerical Results

In this section we will present the numerical solutions obtained with numerical scheme described above. We compare numerical results with some analytical solutions and experimental observations.

Comparison with exact solutions

5.1.1 Stationary solution: "water glass on a record player"

The following exact solution of our system without friction corresponds to the profile of the water in the glass turning on a record player :

Let us remark that this is also a solution of the Saint-Venant equations. As the solution does not depend on θ, we test numerically only 1D case. The present test-case (ω = 0.

) shows that our splitting method preserves this solution.

In the code, the boundary conditions at r = R -and r = R + , and the initial conditions are taken from the analytical solution. If we take a greater value of ω (for example ω = 1 [s -1 ]), the numerical results will be worse.

CF L number is equal to 0.8. 100 mesh cells is used in the radial direction for the case of ω = 0.1 [s -1 ]. One can see a good agreement between the analytical and numerical solutions. 100, 200 and 500 mesh cells are used in the radial direction for the case of ω = 1 [s -1 ]. The convergence is clearly visible (see Figure 3). The radial velocity U r and the components of P oscillate about zero for all t > 0. to 1 s. The initial conditions are :

The boundary conditions are :

(94) The error was calculated at time instant 10 s in the following way: 7 Appendix A

Tensor product

The following gives some formulas that are useful for calculations of tensor product.

Q.E.D. where

) 2θ )

) .

Or, in matrix form :

CONCLUSIONS AND PERSPECTIVES

This thesis focuses on the study of a multi-dimensional model of weakly sheared flows of long waves (turbulent shallow water flows).

We propose a numerical method consisting in 'physical' splitting of the corresponding equations. Each split subsystem is hyperbolic, contains only one type of waves, and admits the energy and 'entropy' conservation laws. Moreover, such a splitting allows us to naturally define a weak solution to our system which is compatible to the positive definiteness of the Reynolds stress tensor P. The dissipation is introduced in invariant form which also guarantees the positive-definiteness of P.

An interesting feature of the model is the formation of transverse structures at the jump toe perimeter ("fingers") from one-dimensional initial data which are harmonically perturbed in the transverse direction. The number of waves does not depend neither on the amplitude nor on the number of transverse waves in the initial perturbation. Thus, the full transition scenario is observed in the formation of roll waves : from uniform flow to one-dimensional roll waves, and, finally, to 2D transverse "fingering" of roll wave profiles.

Some numerical results in polar coordinates are presented and compared to the analytical solutions and experimental observations. In particular, in some region of physical parameters, the formation of rotating singular (angular) point at the circular hydraulic jump was found numerically. Obtained results are qualitatively similar to those observed experimentally.

The method can obviously be extended to the general 3D Reynolds averaged model of barotropic flows, because the mathematical structure of the governing equations is exactly the same. The case of non-barotropic turbulent flows demands an additional modelling because the both entropies, physical and 'mathematical', will increase at the shock. Thus it will be necessary to separate such an increase into two parts : the thermodynamic and turbulent ones.

This will be the subject of our future work.

APPENDIX A

In this Appendix we will describe the higher-order extension of Godunov's Method -MUSCL("Monotonic Upstream-centred Scheme for Conservation Laws") -Hancock method. The basic idea is to generalize Godunov's method by replacing the piecewise constant representation of the solution by a more accurate piecewise linear reconstruction. This method involves three distinct steps to obtain a higher-order accurate scheme.

The first one is a reconstruction of data. Recall that U n i represents an integral average in

U(x, t n )dx . Now we replace the constant states U n i by piecewise linear functions U i (x) :

where ∆ i ∆x is a slope of U i (x) in cell I i (see Figure 4.1). The boundary extrapolated values are given by

As to the choice of slopes ∆ i we define here ∆ i = 1 2 MinMod(∆U i-1/2 , ∆U i+1/2 ), where

and MinMod limiter is defined as :

The second one is a prediction step, when we evaluate The third step is the resolution of the Riemann problem for the system :

with the initial condition :

To compute Godunov fluxes F * ,n i+1/2 = F( ŪR i , ŪL i+1 ) and F * ,n i-1/2 = F( ŪR i-1 , ŪL i ), we solve the Riemann problems between i, i + 1 and i -1, i cells correspondingly. For that we use the HLLC Riemann solver. We use finally the Godunov scheme :