domain: math sde

text:
Shashwat Tiwari

Boniface Nkonga
email: boniface.nkonga@unice.fr

Praveen Chandrashekar
email: praveen@math.tifrbng.res.in

Sergey Gavrilyuk
email: sergey.gavrilyuk@univ-amu.fr

Finite volume approximations of shear shallow water model on unstructured grids

We develop and implement a finite volume method based on a conservative approach for the hyperbolic non-conservative system of the shear shallow water model. This model is a generalization of the shallow water model to include the effects of shear, which are ignored in the standard shallow water model. The existing path conservative scheme is extended to unstructured grids composed of triangles or quadrilaterals with the solution stored at the cell centers. The high-resolution scheme is based on a local solution reconstruction idea for which a new set of variables based on a Cholesky decomposition of velocity fluctuation tensor is developed, which is found to increase the robustness of the schemes by maintaining positive definiteness property. The scheme is applied to some test problems like 1D dam break, square dam break, roll waves, and radial hydraulic jump.

Introduction

Shallow flows are common in many situations, like the flow in the sea and oceans and even in some channel flows. Recently [START_REF] Foglizzo | Shallow Water Analogue of the Standing Accretion Shock Instability: Experimental Demonstration and a Two-Dimensional Model[END_REF][START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF] the shallow water model has been used as an analog of the accretion shock instability during a stellar core collapse. For these situations, shallow modeling can provide a simplified and accurate model that can be used in a numerical approach to estimate the behavior of these complex flows that are usually very expensive to solve with more general models. Indeed, the flow under concern can involve effects like hydrodynamics, free surface, neutrino interactions, general relativity, etc. The domain of validity of the shallow assumptions is that of thin film flows. In this context, since the depth is much smaller than the horizontal scales, the depth variations are ignored, and only the horizontal variations and velocity are modeled. The usual shallow modeling is derived by depth averaging the incompressible Euler equations and neglecting some terms that involve shear in the depth direction. The shear shallow water model retains some of these terms and can give a more accurate representation of the flow in cases where vertical shear is not insignificant [START_REF] Teshukov | Gas-dynamic analogy for vortex free-boundary flows[END_REF]. This extended shallow water model was introduced in [START_REF] Teshukov | Gas-dynamic analogy for vortex free-boundary flows[END_REF][START_REF] Gavrilyuk | Multi-dimensional shear shallow water flows: Problems and solutions[END_REF] and numerically validated in [START_REF] Gavrilyuk | Multi-dimensional shear shallow water flows: Problems and solutions[END_REF][START_REF] Bhole | Fluctuation splitting Riemann solver for a non-conservative modeling of shear shallow water flow[END_REF][START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF][START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF][START_REF] Busto | On Thermodynamically Compatible Finite Volume Methods and Path-Conservative ADER Discontinuous Galerkin Schemes for Turbulent Shallow Water Flows[END_REF].

The dominant mathematical structure of this model can be classified as a genuinely non-conservative hyperbolic system. Therefore, there is no uniqueness of the jump conditions [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF], and the weak formulations for numerical strategies need to be appropriately defined. In our recent paper [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF], this difficulty has been investigated, and a path conservative [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] approach was used to design a generalized Rankine-Hugoniot (RH) jump conditions. The SSW model can be written in different forms based on the choice of independent variables. While all states are identical for smooth solutions, they can lead to other weak solutions. In [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF], a particular form of the PDE model was proposed for shear shallow flows which directly follows from the depth averaging process, and which is a re-writing of the models used in [START_REF] Gavrilyuk | Multi-dimensional shear shallow water flows: Problems and solutions[END_REF][START_REF] Bhole | Fluctuation splitting Riemann solver for a non-conservative modeling of shear shallow water flow[END_REF]. This model admits a convex entropy function and is very similar to the 10-moment model of rarefied gas dynamics [START_REF] Levermore | The Gaussian Moment Closure for Gas Dynamics[END_REF][START_REF] Berthon | Numerical approximations of the 10-moment Gaussian closure[END_REF].

Reference [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF] also constructed HLL, 3-wave HLLC and 5-wave HLLC solver, with the last one including all the waves in the Riemann problem. Unlike previous works, we do not split the model into several sub-systems, but instead, we construct a unified Riemann solver for the whole system. The same PDE form was also used in [START_REF] Yadav | Entropy stable schemes for the shear shallow water model equations[END_REF] where entropy stable schemes were constructed. In the present work, the approximate Riemann solvers from [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF] are used on the unstructured grids in the face normals directions to estimate the fluctuations. A higher-order version of the scheme is constructed following the method of lines approach based on Runge-Kutta time integration and an implicit treatment of the source terms, and the resulting semi-implicit system is solved exactly. While such path-conservative schemes provide a framework to construct stable numerical approximations, the theoretical analysis of such methods is not well developed.

Up to now, the numerical strategy in the published papers are mostly designed for structured Cartesian meshes. For the shock instability of the stellar accretion experience, the domain is cylindrical. In order to use structured meshes in this context, the formulation in cylindrical coordinates is applied in [START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF], and penalization strategy was used in [START_REF] Busto | On Thermodynamically Compatible Finite Volume Methods and Path-Conservative ADER Discontinuous Galerkin Schemes for Turbulent Shallow Water Flows[END_REF] to account for the hole in the middle of the Cartesian grid. It is helpful to design a method on unstructured grids to be able to solve problems involving complex domain shapes where Cartesian meshes cannot be used to represent the boundaries.

The rest of the paper is organized as follows. In Section 2, we present the SSW model in terms of the partial differential equations and its re-writing in an almost conservative model. Section 3 discusses path-conservative schemes and fluctuations. Section 4 deals with the reconstruction scheme and explains the finite volume method. Some numerical results are presented in Section 5, and the paper ends with some conclusions.

The SSW model

Figure [START_REF] Nkonga | Exact solution for riemann problems of the shear shallow water model[END_REF] shows the schematic of the flow problem. Gravity is acting along the x 3 axis and the flow happens in the x 1 x 2 plane. The top surface of the water is shown in red and the bottom surface shown in blue may have some spatial variations, which are measured with reference to some reference level. In the shallow water models, the top free surface is not explicitly modeled but is accounted for by the water depth function h. Let ∇ denote the two dimensional gradient operator

x 3 x 1 , x 2 b h x 3 = ξ(x 1 , x 2 , t) x 3 = b(x 1 , x 2 )
g Figure 1: Shallow water flow: b denotes the profile of the bottom surface and h is the depth of water column with respect to the variables (x 1 , x 2 ). The SSW model can be written as [START_REF] Gavrilyuk | Multi-dimensional shear shallow water flows: Problems and solutions[END_REF]

∂h ∂t + ∇ • (hv) = 0 (1a) ∂(hv) ∂t + ∇ • hv ⊗ v + gh 2 2 I + hP = -gh∇b -C f |v|v (1b) ∂P ∂t + v • ∇P + (∇v)P + P(∇v) = D (1c)
where h is the water depth, v is the horizontal velocity, P is a tensor that accounts for the shear and b is the bottom topography, see Figure [START_REF] Nkonga | Exact solution for riemann problems of the shear shallow water model[END_REF]. The tensor P has the character of a Reynolds which measures the fluctuation of the three dimensional velocity field u i relative to the depth averaged velocity field v i ,

P ij = 1 h ξ b (u i -v i )(u j -v j )dx 3
The momentum equation also contains a frictional source term to model bottom friction and the tensor D represents the dissipation model used to close the system of equations obtained by depth averaging process. The SSW model can be written in an almost conservative form as first done in [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF]. To do this, we define the symmetric tensors

R ij := hP ij , E ij := 1 2 R ij + 1 2 hv i v j , 1 ≤ i, j ≤ 2
where the quantity R ij is the total Reynolds stress tensor while E ij is a kinetic energy tensor that measures the energy in the average flow and the fluctuations. Then an elementary computation shows that the set of equations ( 1) can be written as the following system of non-conservative equations ∂U ∂t

+ ∂F 1 ∂x 1 + ∂F 2 ∂x 2 + B 1 ∂h ∂x 1 + B 2 ∂h ∂x 2 = S (2)
where

U =           h hv 1 hv 2 E 11 E 12 E 22           , F 1 =           hv 1 R 11 + hv 2 1 + 1 2 gh 2 R 12 + hv 1 v 2 (E 11 + R 11 )v 1 E 12 v 1 + 1 2 (R 11 v 2 + R 12 v 1 ) E 22 v 1 + R 12 v 2           , F 2 =           hv 2 R 12 + hv 1 v 2 R 22 + hv 2 2 + 1 2 gh 2 E 11 v 2 + R 12 v 1 E 12 v 2 + 1 2 (R 12 v 2 + R 22 v 1 ) (E 22 + R 22 )v 2           B 1 =           0 0 0 ghv 1 1 2 ghv 2 0           , B 2 =           0 0 0 0 1 2 ghv 1 ghv 2           S =          0 -gh ∂b ∂x1 -gh ∂b ∂x2 -ghv 1 ∂b ∂x1 -1 2 ghv 2 ∂b ∂x1 -1 2 ghv 1 ∂b ∂x2 -ghv 2 ∂b ∂x2          +         0 0 0 1 2 hD 11 1 2 hD 12 1 2 hD 22         +         0 -C f |v|v 1 -C f |v|v 2 -C f |v|v 2 1 -C f |v|v 1 v 2 -C f |v|v 2 2        
This is a hyperbolic system of equations in which the characteristic fields are either linearly degenerate or genuinely non-linear. The smallest and largest eigenvalues are genuinely non-linear, while the remaining four intermediate eigenvalues are linearly degenerate. Hence, the model can give rise to shock solutions even if the initial data is very smooth, which complicates the solution process both at the theoretical and numerical levels. The SSW model admits the conservation of total energy

∂E ∂t + ∂ ∂x 1 E + R 11 + 1 2 gh 2 v 1 + ∂ ∂x 2 E + R 22 + 1 2 gh 2 v 2 = 1 2 hTr(D) -C f |v| 3 where E = E 11 + E 22 + 1 2 gh 2 + ghb
is the total energy. The SSW model with gravity g = 0 is the same as the 10-moment Gaussian moment closure equations of gas dynamics [START_REF] Levermore | The Gaussian Moment Closure for Gas Dynamics[END_REF], where a convex entropy function and associated entropy condition is also shown. The function

η = η(U ) = -h log det(P) h 2
is a convex entropy function concerning the set of variables U for the model written as [START_REF] Balay | PETSc Users Manual[END_REF]. The enstrophy Ψ = log det(P)

h 2
satisfies the equation

∂(hΨ) ∂t + ∇ • (hvΨ) = hTr P -1 D
The shear shallow water modeling should dissipate the enstrophy. The enstrophy's dynamics will control the dissipation rate from large to smaller vortices. Therefore, the dissipation D modeling should mimic the physical property of the flow under consideration. By analogy with the Stokes hypothesis, the dissipation tensor is modeled as [START_REF] Gavrilyuk | Multi-dimensional shear shallow water flows: Problems and solutions[END_REF][START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF]

D = γ Tr (P) 2 I -ζP (3)
where I is the identity matrix, γ and ζ are positive parameters to be properly modeled such as to achieved the dissipation of the total energy and the dissipation of the enstrophy. where in the last equality we have use the relation: 4 det (P) -(Tr (P)) 2 < 0. This relation is valid for any positive definite 2 × 2 tensor. In [START_REF] Gavrilyuk | Multi-dimensional shear shallow water flows: Problems and solutions[END_REF][START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF], two families of dissipation are investigated for γ = 0. As the factor ζ scales as the inverse of a relaxation time, it can be evaluated using the ratio of a characteristic velocity over a characteristic distance. For the first family, scaling of P as v 2 is used to estimate the ratio as v 3 hTr(P) . Then

γ = 0 and ζ = 2 v 3 hTr (P) max 0, C r Tr (P) -h 2 ϕ Tr (P) (4)
For the second family, the ratio of characteristic velocity to distance is estimated as √ where C r is a positive constant and ϕ is a scaling parameter characterizing the enstrophy of the smallest vortexes formed at the interaction with the bottom topography. The profiling of ζ is designed to recover the model used in an earlier work [START_REF] Richard | The classical hydraulic jump in a model of shear shallowwater flows[END_REF]. Only the first family of models is considered in this work. Notice that, for the 10-moment equations in fluid dynamics, the dissipation is required for the entropy η.

Fluctuation distribution scheme

Our objective in this section is to formulate, at the interface between two states, the equivalent for our non-conservative system of what is the normal flow for conservative systems. The interfaces exchanges are expressed as fluctuations based on "path conservative methods" strategy. This strategy, developed in detail in [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF], is now briefly recalled. Let us first write the model (2) in quasi-linear form

∂U ∂t + A 1 (U ) ∂U ∂x 1 + A 2 (U ) ∂U ∂x 2 = S(U )
where the Jacobian matrices are given by

A 1 = ∂F 1 ∂U + [B 1 , 0 6×5 ], A 2 = ∂F 2 ∂U + [B 2 , 0 6×5 ]
Let n be a unit vector and define the Jacobian and fluxes in the direction of n by

A(U , n) = A 1 (U )n 1 + A 2 (U )n 2 , F (U , n) = F 1 (U )n 1 + F 2 (U )n 2 B(U , n) = B 1 (U )n 1 + B 2 (U )n 2
The eigenvalues of the Jacobian matrix A(U , n) are given by

λ 1 = v • n -c f , λ 2 = v • n -c s , λ 3 = λ 4 = v • n, λ 5 = v • n + c s , λ 6 = v • n + c f
where c s = √ n Pn and c f = gh + 3n Pn. Let U L , U R be two states separated by a discontinuity surface with unit normal vector n pointing from L to R. Consider a path connecting the two states Ψ(ξ; U L , U R ) such that

Ψ(0; U L , U R ) = U L , Ψ(1; U L , U R ) = U R
If the discontinuity is moving with speed S in the normal direction n, then the generalized Rankine-Hugoniot jump condition is given by

1 0 [A(Ψ(ξ; U L , U R ), n) -SI] dΨ dξ (ξ; U L , U R )dξ = 0 (5)
The choice of the correct path is important but it is usually not known. We will assume a linear path in the phase space of U variables between the two states so that

Ψ(ξ; U L , U R ) = U L + ξ(U R -U L ), 0 ≤ ξ ≤ 1
Then the jump condition (5) can be written as

F (U R , n) -F (U L , n) + B( Ū , n)(h R -h L ) = S(U R -U L ) (6)
where

Ū = 1 2 (U L + U R )
The main tool used to derive the numerical scheme is to distribute the total fluctuation D (associated to a given direction) into two parts D -, D + corresponding to left moving and right moving waves arising in the Riemann solution,

D(U L , U R , n) = D -(U L , U R , n) + D + (U L , U R , n)
where the fluctuation is defined as

D(U L , U R , n) := 1 0 A(Ψ(ξ; U L , U R ), n) dΨ dξ (ξ; U L , U R )dξ = F (U R , n) -F (U L , n) + B( Ū , n)(h R -h L )
The distributions D ± of the total fluctuation can be performed in many different ways, according to the waves patterns used in the Riemann solver. In this work, we will adopt the HLL-type approach presented in [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF] where 2-wave, 3-wave and 5-wave approximate Riemann solvers have been developed. The multi-state Riemann solution is obtained by satisfying the generalized jump conditions (6) to determine the intermediate states. This approach requires us to estimate the smallest and largest wave speeds S L , S R that arise in the Riemann problem. We refer the reader to [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF][START_REF] Nkonga | Exact solution for riemann problems of the shear shallow water model[END_REF] for details on the choice of these speeds.

C i C j n ij

Path-conservative finite volume method

We consider a non-overlapping conformal partition of the spatial domain Ω into control volumes C i containing their gravity center

x i Ω = i C i with C i ∩ C j = ∅ when i = j and x i ∈ C i
This paper will adopt the cell-centered finite volume strategy, which means that, inside each cell, the solution is approximated by the cell average: U i . The design of numerical schemes is in a general framework where the cells can be of different shape. Nevertheless, all cells will be either triangles or quadrangles in the applications we consider. Let ϑ(i) be the set of indices of the neighbor cells to the cell C i . The set ϑ(i) is the von Neumann neighbours of the i'th cell ϑ(i) = {j, such that C i and C j have a common boundary} Let n ij be the unit normal vector to the interface between cell i and cell j, which points into cell j, and let η ij be the length of this face.

n ij = η ij η ij where η ij = ∂Cij ndS,
with ∂C ij = ∂C i ∩ ∂C j the interface between the cells C i and C j (see figure 2). Let Q ≡ Q (U ) denote the variables used for high order reconstructions inside each cell. The reconstruction can be achieved either by conservative or by some other set of variables. Here we reconstruct some set of primitive variables which is found to be more robust and also easier to maintain positivity of the reconstructed states, than reconstructing the conserved variables. The interpolated states to the interface mid-point x ij , associated to the two cells surrounding it are given by

Q L ij = Q i + (x ij -x i ) • ∇Q i , Q R ij = Q j + (x ij -x j ) • ∇Q j where Q m ≡ Q (U m ) for any m and ∇Q m = ∂ x Q m , ∂ y Q m t
denotes some estimate of the gradient in the m'th cell. The finite volume formulation is obtained by integrating the equation

(2) over the cell C i . This procedure gives contributions that are, on the one hand, associated with the interior involving the non-conservative terms and, on the other hand, traces on the edge of the cell. In the present context where the considered system is non-conservative, the generalized jump relations between the states at the cell edge neighborhood are used to evaluate the traces at the edge in terms of fluctuations. The second order in space semi-discrete scheme then takes the following form

dU i dt + 1 a i j∈ϑ(i) η ij F U L ij , n ij + D -U L ij , U R ij , n ij + B 1 (U i )∂ x h i + B 2 (U i )∂ y h i = S(U i ) (7)
where a i is the area of the i'th cell, and (∂ x h i , ∂ y h i ) are approximations of derivatives of h at the cell center which is contained in the first component of ∇Q i and U

L/R ij = U Q L/R ij
. Under the Path-conservative framework, the trace of the integral at the interface ∂C ij is decomposed into left and right-going fluctuations. The left going fluctuation

D -U L ij , U R ij , n ij modifies the solution in cell C i and the right going fluctuation D + U L ij , U R ij , n ij modifies the solution in cell C j .
We evaluate the fluctuations from different 1D Riemann solvers proposed in our previous paper [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF]. These solvers use different wave structures of increasing richness, ranging from a two-wave structure (HLL) to a five-wave structure (HLLC5) through a three-wave structure (HLLC3). These Riemann solvers require the estimates of slowest and fastest wave speeds, in the normal direction n, denoted by S L , S R respectively. Once the approximate Riemann solution is computed, the split fluctuations in the direction n for the HLLC5 are obtained as

D ± (U L , U R , n) = S ± 1 (U L * -U L ) + S ± 2 (U L * * -U L * ) + S ± 3 (U R * * -U L * * ) + S ± 4 (U R * -U R * * ) + S ± 5 (U R -U R * )
where the input of the solver are the unit normal n, the outermost waves S 1 = S L , S 5 = S R in the direction n, the states L/R * / * * are computed in the Riemann solver process. The exponent (+) gives wave velocities S + = 0 if S is negative and S + = S if S is positive. On the contrary, the exponent (-) gives wave velocities S -= 0 if S is positive and S -= S if S is negative. In summary, the notation can be defined as S ± = S±|S| 2 . The stability of the numerical scheme also relies on estimating the fastest (S R ) and slowest (S L ) wave velocities such that the numerical dependence cone is larger than the physical cone for the given Riemann data. On the other hand, the scheme's robustness also depends on the properties of the reconstructions at the cell interfaces, which is summarized in the approximation of the averaged slopes.

Q L = Q U L and Q R = Q U R .
Remark 1. The properties of the fluctuation operators allow the following relations:

D -U L ji , U R ji , n ji = D + U L ij , U R ij , n ij
Moreover, for conservative models, we recover the flux formulation and a conservative scheme. Indeed, in that context, we use the fact that

U L ji = U R ij , n ji = -n ij and F U L ji , n ji = -F U R ij , n ij to derive the following relations -F U R ij , n ij + D + U L ij , U R ij , n ij = -F U L ij , n ij -D -U L ij , U R ij , n ij
This relation also says that

D + U L ij , U R ij , n ij + D -U L ij , U R ij , n ij = F U R ij , n ij -F U L ij , n ij

Variables used for slope reconstructions

We obtain the first-order scheme in space by setting all the average gradients to zero. However, we must define the gradient according to the given mean states for second-order schemes. Estimating average gradients is done using the total variation diminishing properties (TVD). The choice of variables for reconstruction is essential to ensure robustness in maintaining the TVD properties. Moreover, the reconstruction scheme may not keep the positive definiteness of the symmetric tensor R. We propose a set of primitive variables to be interpolated to ensure the positive definite property of the tensor R. Let us perform a Cholesky decomposition of the tensor R and obtain

R = 1 0 L 12 1 d 1 0 0 d 2 1 L 12 0 1
We have opted for an alternative Cholesky decomposition that avoids using square roots in the change of variables, which is usually called LDL Cholesky decomposition. The transformation between the variables (R 11 , R 12 , R 22 ) and (d 1 , d 2 , L 12 ) is given by

R 11 = d 1 , R 12 = d 1 L 12 and R 22 = d 1 L 2 12 + d 2 . (8)
and inversely

d 1 = R 11 , d 2 = det R R 11 and L 12 = R 12 R 11 (9)
If we reconstruct d 1 , d 2 , L 12 and the reconstructed values satisfy d 1 > 0 and d 2 > 0, which can be ensured by suitably limiting the gradient estimate ∇Q, then the reconstructed tensor R satisfies the properties

R 11 = d 1 > 0, R 22 = d 1 L 2 12 + d 2 > 0, det R = d 1 d 2 > 0
This ensures the positive definiteness of the reconstructed tensor R at the faces, which is not guaranteed if we reconstruct primitive or conserved variables. For classical shallow-water models, it is customary to consider the water level and velocity components as high-order reconstruction variables. For the shear shallow water model, we will consider for high order reconstruction the following set of variables

Q ≡ Q (U ) = [h, v 1 , v 2 , d 1 , L 12 , d 2 ]
To estimate the gradient in a cell, we have the current cell's average value Q (U i ) and the neighboring cells' average values : Q (U j ) for j ∈ ϑ(i). Each cell is assumed convex, and its boundary is described by a set of connected edges so that there are at least three different neighbor states.

On the other hand, to compute a 2D gradient, we need data corresponding to three non-aligned points. Consequently, we have more data than we need to evaluate the average slope. We will then use least squares first to assess the mean gradient. This first evaluation has no reason to satisfy any particular property. In particular, it does not prevent, at the interfaces, have interpolated quantities that go out of the bounds defined by the cell's average data. It is usual to apply a slope limiting technique to fix this issue. In this way, the moderate gradients ∇Q i are computed so that, if initially, the quantities Q (U j ) were all positive, we are assured that the amounts interpolated will be positive too. In passing, let us note that in previous articles published on this topic, gradients reconstruction uses the components of the symmetric tensor P without guaranteeing its positivity. On the contrary, by guaranteeing here the positivity of h and R, the positivity of the P tensor is realized.

For the computation of the gradient using least squares for cell i, we consider its Moore neighborhood stencil ϕ(i) defined as ϕ(i) = {j, such that C i and C j has either a common edge or a common vertex} Let us write a truncated Taylor expansion around cell i,

Q j = Q i + (x j -x i ) • ∇Q i , j ∈ ϕ(i)
which is an over-determined system of equations for the unknown gradient ∇Q i . We use the least squares method as shown in [START_REF] Barth | Aspects of unstructured grids and finite-volume solvers for the euler and navierstokes equations[END_REF] which allows for weighting of the equations corresponding to each neighbour by some weighting coefficient of our choice, i.e., we minimize

j∈ϕ(i) w 2 ij Q j -Q i -(x j -x i ) • ∇Q i 2
with respect to ∇Q i . In this work, we weight each equation by inverse of the distance between cell i and neighbor cell j, i.e., w ij = 1/|x jx i |, which provides a better estimate of gradient for unstructured mesh, especially in case of stretched cells. The gradients obtained by the least squares procedure are limited such that min k∈{i}∪ϑ(i)

Q k ≤ Q i + (x j -x i ) • ∇Q i ≤ max k∈{i}∪ϑ(i) Q k , for all j ∈ ϑ(i)
is satisfied. This limitation is necessary to avoid the generation of Gibbs-type oscillations in high resolution schemes around regions of sharp solution variation. The limiting also ensures that the reconstructed values of Q 4 = d 1 , Q 6 = d 2 will be positive if the corresponding cell values are positive, and hence the reconstructed solution will ensure positive definiteness of R.

Time integration steps

The spatial discretization in (7) leads to a system of ordinary differential equations which can be written as

dU dt = R(U ) + S(U )
where R is the hyperbolic part and S is the source term. The source terms can lead to loss of positive definite-ness of the tensor P or R if explicit time stepping is used. Hence we would like to treat the source terms in an implicit manner by using an operator splitting approach based on Strang splitting. Let us first define the method to update the solution using the two individual operators.

Explicit scheme for hyperbolic part: U n+1 = R(U n ; ∆t) The update of the solution due to the hyperbolic part R(U ) is achieved by the two stage, explicit strong stability preserving Runge-Kutta scheme given by

U * = U n + ∆tR(U n ) U n+1 = 1 2 U n + 1 2 [U * + ∆tR(U * )]
Implicit scheme for source part: U n+1 = S(U n ; ∆t) The update of the solution due to the source term S(U ) is achieved by the following scheme

U * = U n + ∆tS(U * ) U * * = U * + ∆tS(U * * ) U n+1 = 1 2 U n + 1 2 U * *
Note that the first two steps in the source update scheme are backward Euler schemes (implicit) which are designed to be unconditionally stable and positive.

We use the above two schemes for the two operators R and S in a Strang splitting strategy. The update U n → U n+1 is done by the following three steps. First, we update the initial state U n to the state U (1) with source term by time step 1 2 ∆t. Second, we update the state U (1) to the state U (2) with with hyperbolic term by time step ∆t. Finally, we update the state U (2) to the state U n+1 with source term by time step 1 2 ∆t. (1) ; ∆t

U n → U (1) = S U n ; 1 2 ∆t → U (2) = R U
→ U n+1 = S U (2) ; 1 2 ∆t
The first and last steps of the solution update are implicit schemes, by the strategy used to update the source terms. Fortunately, these implicit schemes do not induce coupling between cells, and therefore, we only have to invert several nonlinear systems of six unknowns associated with each cell which can be done explicitly as explained in [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF].

Numerical results

The method is implemented in a C++ code based on MPI and PETSc [START_REF] Balay | PETSc Users Manual[END_REF], and can deal with triangular and quadrilateral grids. The grids are generated with GMSH [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF] and we use DMPLEX in PETSc to manage the unstructured grid. The time step is based on a CFL condition and is given by

∆t = CFL • min i      a i j∈ϑ(i) λ ij η ij      where λ ij = |v i • n ij | + gh i + 3n ij P i n ij
where a i is the volume of the cell C i , n ij is the unit normal, η ij denotes the length of the face between the cells C i and C j , and the minimum is taken over all cells in the mesh.

Exact smooth solution: second-order convergence analysis

In order to estimate the order of convergence of the proposed method, we use an analytical solution developed in [START_REF] Gavrilyuk | Multi-dimensional shear shallow water flows: Problems and solutions[END_REF]. The solution can be written as

h = h 0 1 + β 2 t 2 , v = β 1 + β 2 t 2 βtx + y -x + βty , P = 1 (1 + β 2 t 2 ) 2 λ + γβ 2 t 2 (λ -γ)βt (λ -γ)βt γ + λβ 2 t 2
The values of the parameters used in our numerical experiment are

h 0 = 1.0 (m), λ = 0.10 (m 2 /s 2 ), γ = 10 -2 (m 2 /s 2 ), β = 10 -3 (s -1 )
We compute L 1 and L We are tracking the evolution of the numerical error against mesh size estimated as (1/ √ number of cells). In log-scale, this error curve's asymptotic slope is the numerical convergence error. Whatever we use as a Riemann solver, we asymptotically recover the second-order convergence for Cartesian (Figure 3) and unstructured (Figure 4) meshes. However, the rate of convergence and the level of errors differ with the enrichment of the Riemann solver. This fact is visible in the errors associated with the components E ij of the kinetic energy tensor.

2D simulations of a 1D Dam break problem

In this section, we repeat the one-dimensional Riemann problems from [START_REF] Chandrashekar | A path conservative finite volume method for a shear shallow water model[END_REF][START_REF] Nkonga | Exact solution for riemann problems of the shear shallow water model[END_REF] on different 2D meshes. The aim is to show how non-alignment of the mesh with the direction of wave propagation can lead to transverse disturbances. The initial conditions of this test case define the following 1D Riemann problem. Indeed, we set the gravity g = 9.8, the source terms to zero, the topography is considered as flat, the initial velocity is zero everywhere, and the components of the stress tensor are P 11 = P 22 = 10 -4 m 2 /s 2 , P 12 = 0 m 2 /s 2 , while the initial depth has a discontinuity given by h(x, y) = 0.02 m x < 0.5 m 0.01 m x > 0.5 m

This problem is purely 1D, and we can compute the exact solution under the constraint of consistency with the conservation of the total energy [START_REF] Nkonga | Exact solution for riemann problems of the shear shallow water model[END_REF]. The detailed solution profile contains one rarefaction wave and one shock wave separated by a contact discontinuity. The numerical resolution is performed in a 2D domain [0, 1] × [0, 0.02]. We applied periodic boundary conditions in the y-direction, and the final simulation time is such that no wave has had time to interact with the x = 0 or x = 1 boundaries. On these edges, we apply Dirichlet boundary conditions. Figure [START_REF] Bhole | Fluctuation splitting Riemann solver for a non-conservative modeling of shear shallow water flow[END_REF] and ( 6) compares with the exact 1D solution, numerical solutions obtained, at the final time T = 0.5, with Cartesian, unstructured quadrangular, and triangular meshes. All these meshes are of approximately the same number of points: 500 × 10 = 5000 points. We use different Riemann solvers: HLL, HLLC3, and HLLC5. The pattern of these Riemann solvers contains 2-waves (HLL), 3-waves (HLLC3), and 5-waves (HLLC5). The 1D plots as a function of (x) have all the values computed on the 2D points: for the same position x, we plot the points values for different values of y in the computational domain. Numerical results (Figure 5) on the depth (h) show that, for the meshes used, we preserve at the discrete level the 1D property of the solution (no important difference between values associated with the same coordinate x). However, more deviation exists between the numerical solution on unstructured meshes and the exact solution. Using second-order accurate approximations reduced the differences. The situation is slightly different with the profiles of the P 11 component of the Reynolds stress tensor (Figure 6). In particular, we can see a large discrepancy between the exact and numerical solutions, irrespective of the mesh, the Riemann solver, and numerical accuracy. Here, the discrete solution no longer matches with the 1D profile of the exact solution. This tendency, already present in previous work, is largely accentuated here by the non-alignment of the mesh. Probably, this situation results due to the fact that the conservation of total energy is not fully satisfied at the discrete level, and future work will specifically focus on the improvement of this point.

Square dam break in a square domain

This test consists of a square dam-break in a square domain Ω = [-5, 5] × [-5, 5], with periodic boundary conditions. In this context, we set the source terms to zero and assume a flat bottom. Initially, the flow is at rest, and the depth has a 2D pattern with discontinuities: In this context, the stress tensor is initially very small and almost at the limit where there is no shear: det(P) = 10 -8 . The initial data has two lines of symmetry, horizontal (y = 0) and vertical (x = 0). Periodic boundary conditions along x and y ensure the preservation of the property of the initial solution at each time.

h =    0.02 if -1.5 ≤ x ≤ 1.
The numerical strategy uses the HLLC5 Riemann solver with a second-order scheme in time and space. Figures [START_REF]Periodic permanent roll waves[END_REF], ( 8) and ( 9) draw the numerical solutions (h and P 11 variables ) for structured and unstructured meshes. Computed solutions, using structured or unstructured meshes, preserve the symmetries numerically. The Cartesian and unstructured mesh of triangular elements give quantitative and qualitatively similar results. However, the unstructured mesh results are slightly more jagged in places without affecting the global dynamics. This test case demonstrates the numerical strategy's robustness in computing multiple interacting wave contexts. Particularly, the tensor P determinant remains positive at each time step of these numerical simulations.

Roll wave test

Uniform stationary flows on an inclined plane can, under certain circumstances, become unstable. The flow then forms a more or less periodic succession of breaking waves that roll down the slope. These waves, often large, are characterized by a turbulent roll, a swirling zone just behind the wavefront. This phenomenon is called roll waves and has been experimentally studied by Brock [START_REF] Brock | The development of roll-waves trains in open channels[END_REF][START_REF]Periodic permanent roll waves[END_REF]. Brock applied a frequent disturbance at the channel's beginning to obtain periodic rollers.

For the shear shallow water model, in [START_REF] Richard | The classical hydraulic jump in a model of shear shallowwater flows[END_REF], the conditions of the appearance of roll wave trains are established. Then, the parameters of the SSW modeling that fit with some Brock experiments were provided. In this context, we use one of these settings with periodic boundary conditions in the slope direction. According to [START_REF] Richard | The classical hydraulic jump in a model of shear shallowwater flows[END_REF], the computational domain is Ω = [0, L x ] × [0, 0.5]and the boundary conditions are periodic in both directions. The initial state is given by

h = h 0 1 + a sin 2π L x , v 1 = gh 0 C f , v 2 = 0, P 11 = P 22 = 1 2 ϕh 2 , P 12 = 0, b = -x tan(θ)
where L x (m) is the critical wave length, a is the level of transverse perturbation, θ(rad) is the bottom slope angle, h 0 (m) the still water depth, C f and C r the dissipation coefficients; the values of the various parameters are taken as

L x = 1.3, a = 5 10 -2 , θ = 0.05011, h 0 = 0.00798, C f = 0.0036, C r = 0.00035, ϕ = 22.76
The quadrilateral (Cartesian) mesh has 554400 cells, while the triangular mesh has 1464510 cells. Both meshes were generated using Gmsh with the same characteristic grid spacing of 1.3/1200. We have carried out numerical simulations with various Riemann solvers. In line with previously published results, we first observe the development of a 1D front independent of the Riemann solver.

In the second phase, using a Riemann solver with at least three waves, transverse perturbations of the one-dimensional front develop (Figure 10 (a) and (b)). These developments appear everywhere for meshes aligned with the 1D front (Cartesian quadrangular and triangular). Non-aligned meshes (triangular and quadrangular) reduce and localize the front perturbations at the interaction with boundaries (Figure 10 (c)). We can also observe in Figure [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] that the solution profiles behave similarly except at the front location. These results suggests that the stability of the 1D roll wave profile depends on the nature of the transverse perturbations generated by both the mesh and the numerical solver used. Unfortunately, no theoretical results exist on the stability of such 1D profiles for this non-conservative hyperbolic equation system. Our results are consistent with those we obtained in [START_REF] Bhole | Fluctuation splitting Riemann solver for a non-conservative modeling of shear shallow water flow[END_REF] using completely different numerical strategies. Namely, the 1D profile appears stable (at least over very long times) for non-aligned meshes (unstructured triangles and quadrangles). Nevertheless, in [START_REF] Busto | On Thermodynamically Compatible Finite Volume Methods and Path-Conservative ADER Discontinuous Galerkin Schemes for Turbulent Shallow Water Flows[END_REF], with a numerical strategy based on thermodynamical compatibility construction combined with a Path-conservative ADER-DG scheme, the 1D profile is stable even when using an aligned mesh. In the same article [START_REF] Busto | On Thermodynamically Compatible Finite Volume Methods and Path-Conservative ADER Discontinuous Galerkin Schemes for Turbulent Shallow Water Flows[END_REF], simulations of the hydraulic jump in a radial flow show transverse structures at the front surface. They are similar to those present when Brock's profile is transversely unstable (roll wave). Consequently, we can only speculate on the nature of the observed instabilities of Brock's profile: physically present in the model or a purely numerical artifact as the carbuncle phenomenon in high speed aerodynamics [START_REF] Gallice | Entropy stable and positivity preserving godunov-type schemes for multidimensional hyperbolic systems on unstructured grid[END_REF].

Hydraulic jumps in radial flows

Circular hydraulic jumps are commonly observed in a kitchen sink when a vertically falling tap water jet strikes a horizontal plate and spreads radially outwards. Radial symmetry can break at the hydraulic jump front based on the value of some physical parameters such as viscosity or capillarity [START_REF] Bush | An experimental investigation of the stability of the circular hydraulic jump[END_REF]. A similar pattern arises for convergent radial flows when a liquid is radially injected inwards and strikes a hollow cylinder in the center, through which the water flows out. Using this setup, a fundamental hydrodynamic instability, playing a determining role in the explosion of massive stars, has been explained with the analogy between acoustic waves in a gas and a thin layer of convergent radial flow. Based on the flow of a thin layer of water, it was possible to reproduce experimentally and study an analog of the instabilities causing a profound asymmetry in the explosion of a star and the birth of a neutron star. In [START_REF] Foglizzo | Shallow water analogue of the standing accretion shock instability: Experimental demonstration and a two-dimensional model[END_REF], an experiment is set up with a converging radial shallow water flow to reproduce the non-circular dynamics of the shock front during a supernova's explosion. The topography has a cylindrical symmetry for this shallow water analog of a standing accretion shock instability. The bottom radial profile is a segment of a straight line inclined towards the center of symmetry where there is a hollow cylindrical obstacle of a radius R -, symbolizing a bottomless hole, see Figure [START_REF] Foglizzo | Shallow Water Analogue of the Standing Accretion Shock Instability: Experimental Demonstration and a Two-Dimensional Model[END_REF]. The hollow cylinder, relative to the point of intersection with the inclined plane, has a height of A. In practice, we are using a modified profile as shown in Figure [START_REF] Foglizzo | Shallow water analogue of the standing accretion shock instability: Experimental demonstration and a two-dimensional model[END_REF] for numerical simulation by enlarging the thickness of the bottomless cylinder.

A similar strategy was used in [START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF] and allowed progressively accelerating the fluid to make the flow supercritical at the neighborhood of the R -position. This strategy is advantageous for imposing boundary conditions at the bottomless cylinder outlet efficiently and straightforwardly.

The bottom topography for the SWASI test case is given by b

(x, y) =      A (r-R -+L * ) 2 -L * 2 L 4 * if R --2L * ≤ r ≤ R - (r -R -) tan β else where r = x 2 + y 2
Notice that L b = (R + -R -) tan β. At the outer boundary, where r = R + , there are n = 16 injectors with an average mass flux q 0 . In practice, we use perturbed velocity to start-up the shock instability. The boundary conditions for the incoming flow, at r = R + and at any time t, are given by h = h 0 , u r = -q 0 R + h 0 1 + 10 -2 sin nθ , u θ = 0.0, P rr = P θθ = ϕh 2 0 2 , P rθ = 0 where u r = v • e r , u θ = v • e θ , P rr = P : e r ⊗ e r , P rθ = P : e r ⊗ e θ and P θθ = P : e θ ⊗ e θ with e r , e θ the unit vectors of the cylindrical coordinates. We complete the shallow water analogue to standing accretion shock instability (SWASI) by the following inputs data

h 0 [m] β[rad] C r C f ϕ[s -2 ] q 0 [m 3 /s] R + [m] R -[m] A[m] L 1 [m] 3 • 10 -3 0.07 1.0 3.6 • 10 -3 2 1.2 • 10 -3 √ 2 0.12 5 • 10 -3 6 • 10 -3
Compared to the data used in [START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF], we applied a factor of √ 2 to the outer circle radius to better match the scales used in the experiments [START_REF] Foglizzo | Shallow water analogue of the standing accretion shock instability: Experimental demonstration and a two-dimensional model[END_REF], particularly the total flux on injected mass. Moreover, the definition of R -is slightly different but without changing the dimensionless parameters of the simulations. All the other parameters used are left unchanged. The input data ensures that the incoming flow Froude number is more significant than one, a necessary condition for developing an unstable shock front.

The numerical solution we are proposing is intended to complement work carried out recently [START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF][START_REF] Busto | On Thermodynamically Compatible Finite Volume Methods and Path-Conservative ADER Discontinuous Galerkin Schemes for Turbulent Shallow Water Flows[END_REF]. We share with them a similar physical model with a few nuances. In [START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF], the model uses the specific Reynolds stress tensor as a variable and a first-order numerical strategy on structured meshes based on polar coordinates reformulation. In [START_REF] Busto | On Thermodynamically Compatible Finite Volume Methods and Path-Conservative ADER Discontinuous Galerkin Schemes for Turbulent Shallow Water Flows[END_REF], reformulation of the specific Reynolds tensor introduces a new object as a variable to guarantee some valuable model properties at the discrete level. The numerical strategy uses a high-order path-conservative ADER-DG on Cartesian meshes, and the penalization strategy enforces boundary conditions on circular immerse boundaries. In contrast, we use the unstructured mesh fitting with the physical domain boundaries and secondorder approximation methods with the kinetic energy tensor as a variable in place of the Reynolds stress tensor [START_REF] Nkonga | Exact solution for riemann problems of the shear shallow water model[END_REF].

We perform numerical simulations until the final time T = 150 s, for unstructured quadrangular and triangular meshes of 261584 and 241680 cells, respectively. In both cases, we recover the scenario observed experimentally and numerically [START_REF] Foglizzo | Shallow Water Analogue of the Standing Accretion Shock Instability: Experimental Demonstration and a Two-Dimensional Model[END_REF][START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF][START_REF] Busto | On Thermodynamically Compatible Finite Volume Methods and Path-Conservative ADER Discontinuous Galerkin Schemes for Turbulent Shallow Water Flows[END_REF], that starts with the appearance of a symmetric hydraulic jump, followed by the development of transverse structures similar to what is also observed for planar hydraulic jump (see Figure 11), and after all the destabilization of the hydraulic jump with the appearance of the sloshing-type oscillatory pattern. At this latter step, a macroscopic singular point (cusp) appears at the jump front and starts rotating with an oscillatory period of about six seconds ( T s = 6s). This specific point, visible in the figures 16 and 17, is at the junction between the main front and a protrusion facing the inner hole. We have plotted the evolution of the hydraulic jump front with the positions of the macroscopic singular point for quadrangular (Figure 16) and triangular (Figure 17) meshes. We use the plot of the Schlieren function of the depth (h) to highlight areas of relatively strong depth gradients. We define the Schlieren as exp -

∇h max( ∇h )
where ∇h is the L2-norm of the local depth gradient.

Summary and conclusions

This study extends the scope of the SSW model by exploring its use in the context of complex geometries. In this direction, we have improved the construction of high-order schemes by applying slope reconstruction on a set of adapted variables that guarantee the symmetry and positivity properties of the Reynolds stress tensor. Similarly, the source term integration step uses an implicit scheme, a more numerically stable approach, which fits into a second-order Strang decomposition. The results confirm that the proposed numerical strategy is globally second-order for a regular analytic solution. Applications to a complex interaction between multiple waves while preserving positivity demonstrate the method's robustness. Finally, the results obtained for a multi-dimensional turbulent hydraulic jump formed in convergent radial flow., qualitatively compare to the experimental observations. Nevertheless, there are still a few grey areas. In the 2D roll waves test cases, transverse instability patterns are present only when aligned meshes and multiples waves HLLC Riemann solvers are used. We can't definitively reject the purely numerical origin of these transverse patterns that, as in the Quirk test case, can come from perturbations along the propagation direction. On the contrary, hydraulic jump in radial flow, the rotation pattern experimentally observed, is reproduced numerically in a reasonable time scale only when we consider misaligned meshes. Aligned polar mesh of quadrangular elements produces radial structures, but transverse perturbations are not powerful enough to trigger experimental observations where transverse asymmetries are present. We emphasize that other existing numerical results presented in the literature use strategies that generate strong numerical perturbations. Whether through the geometric terms of the polar coordinate formulation [START_REF] Ivanova | Structure of the hydraulic jump in convergent radial flows[END_REF] or the boundary conditions on a circular immersed boundary in a Cartesian mesh [START_REF] Busto | On Thermodynamically Compatible Finite Volume Methods and Path-Conservative ADER Discontinuous Galerkin Schemes for Turbulent Shallow Water Flows[END_REF]. Beyond these grey areas in the physical interpretation of numerical results, computations over complex geometries call for a mesh that can precisely adapt to boundary contours. The primary purpose of this paper was to propose a numerical strategy that accurately fits complex boundaries through discretization with unstructured meshes. Moreover, unstructured meshes also open the way to local mesh adaptation.

First order, HLLC5

Second order, HLLC5

Figure 2 :

 2 Figure 2: Interface between two triangular cells (C i and C j ) and the associated normal (n ij ).

5 and - 1 . 5 ≤ y ≤ 1 . 5 0

 1515 .01 elsewhere The remaining quantities are constant in space and are given by v 1 (m/s) v 2 (m/s) P 11 (m 2 /s 2 ) P 12 (m 2 /s 2 ) P 22 (m 2 /s 2 )

Figure 5 : 1 -Figure 6 : 1 -

 5161 Figure 5: 1-D dam break problem solved on 2-D mesh. Profile of the depth (h) at the final time T = 0.5. The meshes are Cartesian (quadrangular), unstructured quadrilateral and unstructured triangular.

Figure 7 :

 7 Figure 7: Square dam break: contours of h. Top row: t = 12.5, middle row: t = 20.0, bottom row: t = 100.0. Left column: Cartesian mesh, right column: triangular mesh.

Figure 8 :

 8 Figure 8: Square dam break: evolution of the contours plot of P 11 in log scale. Top row: t = 12.5, middle row: t = 20.0, bottom row: t = 100.0. Left column: Cartesian mesh, right column: triangular unstructured mesh.

Figure 9 :Figure 10 :

 910 Figure 9: Comparison between Cartesian and unstructured triangular meshes : evolution of h (left column) and P 11 (right column) profiles along diagonal. Top row: t = 12.5, middle row: t = 20.0, bottom row: t = 100.0. Profiles of P 11 are in log-scale.

Figure 11 :Figure 12 :

 1112 Figure 11: Roll waves problem (2D) with second order scheme: 1D profiles of the water depth for shifted Cartesian quadrangular (blue), Cartesian triangular (red), and shifted unstructured triangular (green) meshes.

Figure 13 :

 13 Figure 13: Zoom around the bottomless hole : augmented radial bottom profile (blue line) used for numerical simulations.

Figure 15 :

 15 Figure 15: Zoomed view near inner boundary of unstructured meshes used to simulate the SWASI experiment : R -= 0.12 and L * = 0.02. (a) mesh of quadrilateral cells, (b) mesh of triangular cells. The circle r = R -is visible inside the meshes.

Figure 16 :

 16 Figure 16: SWASI test on an unstructured quadrilateral mesh (261584 cells): Schlieren of the fluid depth (h) starting at t = 144.5 seconds and plot at time intervals of 0.5 seconds, top-left to bottom-right. The cusp rotates counter clockwise with a period of about 6 seconds.

Figure 17 :

 17 Figure 17: SWASI test on an unstructured triangular mesh (241680 cells): Schlieren of the fluid depth (h) starting at t = 144.5 seconds and plot at time intervals of 0.5 seconds, top-left to bottomright. The cusp rotates clockwise with a period of about 6 seconds.

  2 norm of the numerical solution errors at time t = 50 seconds on the domain of size [0, 10] × [0, 10] and use Dirichlet boundary conditions by setting ghost cell values to the analytical solution. To perform convergence analysis, we use different levels of mesh refinement. The following table shows the number of cells used at each level of mesh refinement. We can not straightforwardly generate Cartesian and unstructured mesh of the same number of cells. Therefore, we perform comparisons with slightly different cell numbers.

	Mesh Cartesian	10 × 10 20 × 20 40 × 40 80 × 80 160 × 160 100 400 1600 6400 25600
	Unstructured	104	374	1578	6396	25350
	Number of cells for each level of refinement.

Acknowledgments

B. Nkonga's work was partially supported both by the INRIA Associate team program (AMFoDUC) and the CNRS program "Dispositif de Soutien aux Collaborations avec l'Afrique subsaharienne" (AMoFlowCiGUC). The work of P. Chandrashekar was supported by Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.01-0520. PC also thanks CNRS for support to visit Univ. of Côte d'Azur (UCA), Nice, where part of this work was performed.


F Moisy
email: moisy@fast.u-psud.fr

C Morize

M Rabaud Fluides

THE DECAY LAW OF GRID TURBULENCE IN A ROTATING TANK

The energy decay of grid-generated turbulence in a rotating tank is experimentally investigated by means of particle image velocimetry. For times smaller than the Ekman timescale, a range of approximate self-similar decay is found, in the form u 2 (t) ∝ t -n , with the exponent n decreasing from 2 to values close to 1 as the rotation rate is increased. This observation is interpreted in the frame of a phenomenological model based on the exponent of the energy spectrum, in which both the effects of the rotation and the confinement are taken into account.

INTRODUCTION

Assuming self-similarity, the decay of high Reynolds number homogeneous and isotropic turbulence is usually described by a power law [START_REF] Comte-Bellot | The use of a contraction to improve the isotropy of gridgenerated turbulence[END_REF][START_REF] Saffman | Large scale structure of homogeneous turbulence generated at initial instant by distribution of random impulsive forces[END_REF][START_REF] Mohamed | The decay power law in grid-generated turbulence[END_REF][START_REF] Skrbek | On the decay of homogeneous isotropic turbulence[END_REF],

u 2 ∝ (t -t * ) -n , ( 1
)
where u is the velocity variance, n is the decay exponent and t * is a virtual origin. The value of the exponent n depends on whether the size of the energy-containing eddies is free to grow (n = 6/5) or is is bounded by the experiment size (n = 2), with a possible changeover between these two regimes [START_REF] Skrbek | On the decay of homogeneous isotropic turbulence[END_REF].

In the presence of rotation, in addition to the turnover time /u (where is the size of the energy-containing eddies), two other timescales are present in the problem, which have opposite effects on the turbulence decay: the rotation timescale, Ω -1 , and, for bounded systems, the Ekman timescale, t E = h(νΩ) -1/2 (where h is the size of the experiment along the rotation axis). The ratio of the turbulent time scale and the rotation rate is the Rossby number, Ro = u/2Ω . While the primary effect of the rotation is to reduce the energy dissipation [START_REF] Jacquin | Homogeneous turbulence in the presence of rotation[END_REF], the effect of the Ekman friction is to accelerate the decay at large times, shortening the temporal range for the turbulent decay even at large Reynolds numbers.

Based on the assumption that the energy transfers time scale is governed by Ω -1 , Squires et al. [START_REF] Squires | The asymptotic state of rotating homogeneous turbulence at high Reynolds numbers[END_REF] have proposed a self-similar decay with an exponent half of that without rotation, u 2 ∝ t -3/5 , in the limit of zero Rossby number.

In experiments, where both confinement and finite Rossby numbers are to be considered, the situation is more complex. In small experiments in which h O((ν/Ω)

Decay without rotation

In the absence of rotation, the decay law of an unbounded homogeneous turbulence can be derived assuming a two-range model for the energy spectrum [START_REF] Comte-Bellot | The use of a contraction to improve the isotropy of gridgenerated turbulence[END_REF][START_REF] Saffman | Large scale structure of homogeneous turbulence generated at initial instant by distribution of random impulsive forces[END_REF]: A "permanent" part E(k) Ak 2 for small k, followed by the Kolmogorov spectrum

E(k) C 2/3 k -5/3 at larger k. Solving for = -d(u 2 )/dt, where u 2 = ∞ 0 E(k)
dk is (twice) the total kinetic energy, yields u 2 ∝ (t + t) -6/5 , i.e. a decay exponent n = 6/5. In this expression, the crossover time t > 0 does not necessarily correspond to the virtual origin t * introduced in Eq. ( 1), even though their order of magnitude should be both given by the timescale of the initial large eddies. During the decay, the energy-containing scale k e (t) -1 (where k e is the the crossover between the two spectral ranges) grows without bound. In a real experiment, a bounding size L is present, defining a minimum wavenumber k 0 L -1 towards which k e (t) saturates after a given time t s . Setting E(k) = 0 for k < k 0 to account for the confinement yields, for t t s , a sharper decay, with an exponent n = 2 (Skrbek & Stalp [START_REF] Skrbek | On the decay of homogeneous isotropic turbulence[END_REF]) (see Table 1).

Decay with rotation

A very crude way to account for the effect of rotation is to modify the scaling of the high-wavenumber part of the spectrum, but without including the anisotropy. Assuming that E(k) depends on , Ω and k, dimensional analysis yields

E(k) = C p Ω (3p-5)/2 (t) (3-p)/2 k -p . ( 2
)
where the exponent p can take any value between 1 and 3, and C p is a non dimensional constant .

For Ω = 0, the Kolmogorov exponent p = 5/3 is recovered. For = 0, the Kraichnan spectrum in the enstrophy cascade regime (p = 3) is recovered. Finally, the intermediate case p = 2 is the spectrum proposed by Zhou [START_REF] Zhou | A phenomenological treatment of rotating turbulence[END_REF], on the assumption of an energy transfer timescale given by Ω -1 . Solving for = -d(u 2 )/dt, and assuming that the spectral exponent p in (2) remains constant during the decay, the decay exponent for unbounded rotating turbulence is found,

n = 3 5 3 -p p -1 . ( 3
)
The value n = 6/5 is recovered for the Kolmogorov spectrum p = 5/3, and the limiting case p = 3 of the enstrophy cascade regime yields n = 0 (i.e. conservation of energy). The exponent p = 2 of Ref. [START_REF] Zhou | A phenomenological treatment of rotating turbulence[END_REF] yields n = 3/5, which is twice smaller than the exponent without rotation, as first noticed by Squires et al. [START_REF] Squires | The asymptotic state of rotating homogeneous turbulence at high Reynolds numbers[END_REF] on dimensional grounds.

Effects of the rotation and confinement -If the confinement is now taken into account, the same procedure as above may be used, by setting E(k) = 0 for the low wavenumber range, which gives the decay exponent

n = 3 -p p -1 . ( 4
)
As for the non rotating case, the exponent with confinement is larger by a factor of 5/3 than that without confinement. One now obtains n = 2 for p = 5/3 and n = 1 for p = 2, with again the factor 2 between the zero rotation and the fast rotation regimes.

EXPERIMENTAL RESULTS

In order to characterize the effects of the background rotation and the confinement, a series of experiments of decaying grid-generated turbulence have been carried out in a rotating tank.

The experimental set-up, described in Morize et al. [START_REF] Morize | Decaying grid-generated turbulence in a rotating tank[END_REF], consists of a water filled square tank, L = 35 cm in side, mounted on a rotating turntable, whose angular velocity Ω has been varied between 0.13 and 4.34 rad s -1 . Turbulence is generated by rapidly towing a square grid of mesh M = 39 mm at a mean velocity V g = 0.65 m s -1 from the bottom to the top of the tank. The grid Reynolds number is Re g = M V g /ν = 2.5 × 10 4 , a value larger than that of most conventional wind-tunnel experiments, and the grid Rossby number Ro g = V g /2ΩM is kept sufficiently large, in the range 2 -65, in order to minimize the effect of the rotation on the turbulence production mechanism in the near wake of the grid. Instantaneous velocity fields in the horizontal plane at mid-height of the tank are measured by particle image velocimetry (PIV), using a corotating high resolution camera. Ensemble averages over 50 realizations of the decay are performed to achieve statistical convergence.

The energy decay is shown for three rotation rates in Fig. 1. For low rotation rate (Ω = 0.12 rad s -1 ), after a period of approximately constant energy, a well defined power law with a decay exponent 10 0 n 2.02 ± 0.10 is found. As the rotation rate is increased, the friction from the Ekman layers arises earlier (full arrows in Fig. 1), shortening the range of power law. Although the quality of the power law is poorer for large rotation rates, a significant power law range is found, with an exponent that decreases continuously from n 2 to values close to 1, as shown in Fig. 2. However it is not clear from this figure whether a saturation towards n = 1 is present or not.

Although the assumption of a constant spectral exponent during the decay is not realistic [START_REF] Morize | Decaying grid-generated turbulence in a rotating tank[END_REF], the trend for n is found to be consistent with the model presented above, when confinement effects are included (see Table 1). The fact that the confinement plays an important role from the early time of the decay originates from the rapid growth of the vertical integral scale through the propagation of inertial waves, which quickly saturates to the tank height h [START_REF] Morize | On the energy decay of rotating turbulence in confined geometry[END_REF]. In this aspect, decaying rotating turbulence in confined geometry strongly differs from that in a wind tunnel with a rotating honeycomb [START_REF] Jacquin | Homogeneous turbulence in the presence of rotation[END_REF], which may have lateral confinement effects but no confinement along the rotation axis. In spite of this difference, it is remarkable that the factor 2 between the exponents at high and low rotation rates (Squires et al. [START_REF] Squires | The asymptotic state of rotating homogeneous turbulence at high Reynolds numbers[END_REF]) is recovered here.

CONCLUSION

To summarize, a clear evidence of the reduction of the energy decay by the rotation has been observed for times smaller than the Ekman timescale. The most important result is that, in addition to the dissipation in the Ekman layers, the confinement plays a central role in the decay law of rotating turbulence. By making the growth of the integral scale along the rotation axis to quickly saturate to the experiment size even at very low rotation rate, the confinement leads to a sharper decay than for unbounded tur-bulence, although compatible with the reduced dissipation induced by the rotation.

When the range between the rotation-induced saturation time of the vertical lengthscale and the Ekman timescale is large enough, a significant self-similar energy decay is observed, characterized by a decay exponent n decreasing from 2 to values close to 1 as the rotation rate is increased. These exponents are found to be in qualitative agreement with a phenomenological model based on the exponent of the energy spectrum, in which both the effects of the rotation and the confinement are taken into account. These observations may be of primary importance for the modelling of turbulence in rotating containers, as the confinement is shown to significantly influence the energy decay even for very weak rotation rate and large experiment size.

BIBLIOGRAPHY

Fig. 1 .

 1 Fig. 1. Energy decay for (•) Ω = 0.12 rad/s, ( ) Ω = 0.40 rad/s, and ( ) Ω = 1.50 rad/s (2nd and 3rd curves divided by 10 and 100 for visibility). The full arrows indicate the Ekman time.

Fig. 2 .

 2 Fig. 2. Semi-logarithmic plot of the decay exponent n as a function of the rotation rates Ω.

Table 1

 1 Predictions for the decay exponent n.

		not confined	confined
	no rotation	6/5 (Ref. [1])	2 (Ref. [4])
	rotation	3/5 (Ref. [6])	1