Norddine Zeroual

Cinderella Blin

Marine Saour

Hélène David

Safa Aouinti

Marie-Christine Picot

Pascal Colson

Philippe Gaudard

Restrictive Transfusion Strategy after Cardiac Surgery

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


O Geoffroy
email: olivier.geoffroy@neel.cnrs.fr

H Chazal

C Chillet

A Demier

T Waeckerlé

About the exothermic nanocrystallisation of FeCuNbSiB cores

Keywords: Nanocrystalline, exothermy, modelization

This article deals with the modelization of exothermic properties of crystallization of FeCuNbSiB nanocrystalline alloys. Investigations are based on DSC measurements, confronted to the evolution of crystalline fraction 𝑓 𝑐 (𝑡) in course of annealing. Formulations for dissipation are derived, allowing to simulate the evolution of heat produced for any temperature profile 𝑇(𝑡). Confrontations with experiments are successful.

Motivation

Nanocrystalline FeCuNbSiB alloys belong to the family featuring vanishing anisotropies, with lower eddy current losses than usual crystalline soft magnetic alloys. This combination of low losses and magnetostriction, and very high permeability, explains they are now widely used as cores in power Electronics (filtering…) or in ground fault current circuits breakers.

Resulting from the technology of elaboration, rough amorphous ribbons with thicknesses around 15-20µm are wound to shape the cores, and then crystallized. A critical temperature, near roughly 600°C, must not be exceeded in order to avoid the irreversible formation of borides, very detrimental for magnetic properties.

The nanocrystallization process being markedly exothermic, the cores temperature profile 𝑇(𝑡) has then to be carefully controlled to avoid over-heating. It is usually achieved by setting a plateau sufficiently long at the beginning of heat generation. However, an industrial requirement to increase competitiveness is to reduce as much as possible the total annealing time, involving process optimization. Empirical approach being time consuming, modelization is of course preferable. For this purpose, an accurate modelization of the heat generated, denoted 𝑞(𝑡) in the following, is needed, applicable to any 𝑇(𝑡) profile. Our work lies within this framework.

Investigations

We investigate 𝑞(𝑡) properties on a 𝐹𝑒 73.0 𝐶𝑢 1.1 𝑁𝑏 3.1 𝑆𝑖 15.7 𝐵 7.1 alloy through four items :

-Characterization of 𝑞(𝑡) in course of crystallization (DSC equipment), starting from amorphous ribbons.

-Confrontation of 𝑞(𝑡) values obtained for constant rate dT/dt to the evolution of 𝑓 𝑐 (𝑡), simulated following [1]- [START_REF] Heddad | Impact of the magnitude of the magnetization on the induced anisotropy in transverse field annealed nanocrystalline cores[END_REF] (we'll detail the approach, based on a careful identification of such nanocrystallisation kinetics).

-Consistently with a previous 𝑓 𝑐 (𝑡 𝑒𝑞 ) model [START_REF] Boust | Use of superparamagnetic temperature transition measurement in nanocrystalline alloys to determine low crystalline fractions by modeling of the weakcoupling behavior[END_REF],

elaboration of a universal representation 𝑞(𝑡 𝑒𝑞 ), where 𝑡 𝑒𝑞 denotes an equivalent duration related to the actual duration 𝑡 0 by the following rules of Arrhenius type

𝑡 𝑒𝑞 [𝑡 0 ] = ∫ 𝑑𝜆[𝑡] 𝑑𝑡 (𝑡 0 -𝑡) 𝑑𝑡 𝑡 0 0 (1) 𝜆[𝑡] = a exp[-𝑏/𝑇[𝑡]] (𝑇 𝑖𝑛 𝐾) (2)
-Simulation of 𝑞(𝑡) for profiles 𝑇(𝑡) different from the one used for identification, to be compared to experimental measurements, as pictured on Fig. 1. [START_REF] Lecaude | Comparison of nanocrystallization processes in 𝐹𝑒 73.5 𝐶𝑢 1 𝑁𝑏 3 𝑆𝑖 22.5-𝑥 𝐵 𝑥 for 𝑥 = 7 and 𝑥 = 9[END_REF]

Figure 1 :

 1 Figure 1 : Confrontation of the simulation of heat produced during isothermal annealing (in black) with experimental values obtained for a almost identical composition [3]


J Garres-Díaz

C Escalante
email: escalante@uma.es

T Morales De Luna

M J Castro Díaz

A general vertical decomposition of Euler equations: multilayer-moment models

Keywords: Multilayer-moment models, vertical discretizations, shallow water flows, free surface Euler equations

In this work, we present a general framework for vertical discretizations of Euler equations. It generalizes the usual moment and multilayer models and allows to obtain a family of multilayer-moment models. It considers a multilayer-type discretization where the layerwise velocity is a polynomial of arbitrary degree N on the vertical variable. The contribution of this work is twofold. First, we compare the multilayer and moment models in their usual formulation, pointing out some advantages/disadvantages of each approach. Second, a family of multilayer-moment models is proposed. As particular interesting case we shall consider a multilayer-moment model with layerwise linear horizontal velocity. Several numerical tests are presented, devoted to the comparison of multilayer and moment methods, and also showing that the proposed method with layerwise linear velocity allows us to obtain second order accuracy in the vertical direction. We show as well that the proposed approach allows to correctly represent the vertical structure of the solutions of the hydrostatic Euler equations. Moreover, the measured efficiency shows that the proposed multilayer-moment model needs just a few layers to improve the results of the usual multilayer model with a high number of vertical layers.

Introduction

The modeling and numerical simulation of geophysical flows is currently a very active topic of research. In particular, the study, analysis and numerical simulation of the well-known Navier-Stokes system [START_REF] Chorin | A Mathematical Introduction to Fluid Mechanics[END_REF] is of interest. It describes the motion of fluids and also can be employed for other complex geophysical flows as granular or debris flows.

Nevertheless, its complexity and specially the high computational cost associated to the resolution of this system makes that some other simpler approaches are usually proposed. One to the most widely used is the shallow-water (or Saint-Venant) system [START_REF] De Saint-Venant | Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF]. In the shallow water approach the characteristic length of the domain is assumed to be much longer than the vertical one, and then the vertical effects are assumed to be small when compared to the horizontal ones. Therefore, a vertical integration is done and the averaged horizontal velocity is considered. In [START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF] a rigorous derivation of this system, based on a asymptotic analysis, is introduced. This approach has been used in many recent works (see e.g. [START_REF] Bouchut | Gravity driven shallow water models for arbitrary topography[END_REF][START_REF] Gray | A depth-averaged µ(I)-rheology for shallow granular free-surface flows[END_REF]2,[START_REF] Parés | On the well-balance property of Roe's method for nonconservative hyperbolic systems. applications to shallow-water systems[END_REF]).

Although shallow water type systems are quite popular, mainly due to the low computational cost required for their simulation, they present a main drawback: the vertical information is neglected. Indeed, considering the averaged horizontal velocity and neglecting the vertical deviations from this average leads us to lose important information for many practical applications. Several ways to overcome this limitation have been proposed in the literature. In next paragraphs we describe those that are relevant for our work.

A possible approach to capture vertical effects is to use the so-called σ-coordinate framework (see [START_REF] Phillips | A coordinate system having some spetial advantages for numerical forecasting[END_REF][START_REF] Casulli | Numerical simulation of 3D quasi-hydrostatic, free-surface flows[END_REF]). The idea is first to do a change of variables so that the vertical coordinates varies between 0, the bottom, and 1, the free surface. Then, the free surface is defined as a function of the horizontal coordinates which is derived from the standard kinematic free surface condition. The interval [0, 1] is then discretized, giving a vertical approximation of the solution of the model. This approach has been used in [START_REF] Haidvogel | A primitive equation ocean circulation model using vertical sigma and orthogonal curvilinear horizontal coordinates[END_REF][START_REF] Casulli | A semi-implicit finite difference method for non-hydrostatic free-surface flows[END_REF][START_REF] Casulli | Semi-implicit numerical modeling of nonhydrostatic free-surface flows for environmental problems[END_REF][START_REF] Ma | Shock-capturing non-hydrostatic model for fully dispersive surface wave processes[END_REF] among others.

Another way to recover vertical effects in shallow flows is to use the so-called multilayer approach [1, [START_REF] Audusse | A multilayer Saint-Venant system with mass exchanges for shallow water flows. derivation and numerical validation[END_REF][START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF]. Multilayer systems consider a vertical division of the fluid into M given layers, where this division is not physical, but it may be considered as a virtual vertical discretization. Within each layer, a depth-averaged process is performed leading to a layerwise constant solution for the horizontal velocity.

In [1] the vertical layers are immiscible, while [START_REF] Audusse | A multilayer Saint-Venant system with mass exchanges for shallow water flows. derivation and numerical validation[END_REF][START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF] consider a mass and momentum transfer exchange between the internal layers. In recent years several works have been proposed based on this procedure to model fluids or complex geophysical flows such as bedload and suspension sediment transport, variable density flows or granular and debris flows (see [START_REF] Fernández-Nieto | A multilayer shallow water system for polydisperse sedimentation[END_REF][START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I)-rheology: application to granular collapse on erodible beds[END_REF][START_REF] Morales De Luna | Derivation of a multilayer approach to model suspended sediment transport: Application to hyperpycnal and hypopycnal plumes[END_REF][START_REF] Bristeau | Layer-averaged euler and Navier-Stokes equations[END_REF][START_REF] Bonaventura | Multilayer shallow water models with locally variable number of layers and semi-implicit time discretization[END_REF][START_REF] Guerrero Fernández | A second-order well-balanced finite volume scheme for the multilayer shallow water model with variable density[END_REF][START_REF] Garres-Díaz | Multilayer models for shallow two-phase debris flows with dilatancy effects[END_REF][START_REF] Garres-Díaz | Flexible and efficient discretizations of multilayer models with variable density[END_REF] among others).

The last approach we describe here corresponds to the shallow water moment models [START_REF] Steffler | Depth averaged and moment equations for moderately shallow free surface flow[END_REF][START_REF] Khan | Vertically averaged and moment equations model for flow over curved beds[END_REF]. In the moment approach, rather than dividing the vertical direction into layers, one single layer is considered, but now the variables are supposed to follow a polynomial profile, with an arbitrary degree of approximation, in the vertical direction. Again, this allows us to obtain vertical effects as now the horizontal velocity is not just a constant for the whole column of the fluid, but a polynomial on the vertical variable. We refer the reader to [START_REF] Kowalski | Moment approximations and model cascades for shallow flow[END_REF][START_REF] Koellermeier | Derivation and numerical solution of hyperbolic moment equations for rarefied gas flows[END_REF][START_REF] Koellermeier | Numerical solution of hyperbolic moment models for the Boltzmann equation[END_REF][START_REF] Cantero-Chinchilla | Vertically-averaged and moment equations for flow and sediment transport[END_REF][START_REF] Garres-Díaz | Shallow water moment models for bedload transport problems[END_REF] among others for a better description and applications of this technique. An advantage of this approach with respect to the multilayer one is its lower computational cost. Typically, the vertical structure of solution could be capture more accurately with moment models with a low degree polynomial, while the number of layers required to obtain a similar solution with the multilayer approach demands more computational effort as it will be shown later on.

The contribution of this work is twofold. The first goal is to propose a common framework for vertical discretization of Euler equations, as well as a family of improved multilayer-moment methods. The idea is to avoid the usual difficulties and computational cost when solving the full 3D system, while keeping the information on the vertical direction with a low or moderate complexity and computational cost. The ideas proposed here not only can be seen as a generalization of the multilayer, σ-coordinate or moment approaches, but it will also allows us to generalize and propose new enriched models. In particular, we propose here a family of multilayer-moment method, based on this framework, to improve the accuracy of both methods, taking advantage of the strengths of each one.

The second objective will be to analyze the efficiency and accuracy of the different techniques. For instance, one may conjecture that the computational cost of moment method is lower than the multilayer one, specially when the number of layers M is large, since it increases in a quasi-linear way with M . However, moment method may not be enough to reproduce complex vertical profiles for the velocity. We aim to give answer to this conjecture and to do a thorough study of the efficiency, advantages and disadvantages of the different proposed approximations.

This paper is organized as follows: Section 2 is devoted to present the more general system considered 2 here and its representation in σ-coordinates, which will be vertically discretized. In Section 3 we focus on the case of Euler system and we develop a vertical discretization to obtain the family of multilayer-moment models, and also analyze interesting particular cases. Section 4 is devoted to numerical test cases, and some conclusions are presented in Section 5. Finally, in appendixes Appendix A and Appendix B we give detail of how viscosity effects are included, as well as the particular cases of the classical multilayer and moment models.

Governing system

Let us start by establishing the system we are going to deal with. For the sake of simplicity, we consider here a two-dimensional flow, although everything might be generalized to the 3D case. Let (x, z) be the horizontal and normal directions, respectively, and u = (u, w) the velocity vector. Let us remark that, in this section, we do not necessarily fix a Cartesian reference system, but rather proceed in a general way. Therefore, the coordinates may be considered in a general reference system, for instance performing a rotation of Cartesian reference and using tilted coordinates, which may be useful in some applications such as gravity driven flows. The incompressible Navier-Stokes system for a fluid with constant density ρ ∈ R is given by

       ∂ x u + ∂ z w = 0, ρ ∂ t u + u ∂ x u + w ∂ z u + ∂ x p = ρf x + ∂ x τ xx + ∂ z τ xz , ρ ∂ t w + u ∂ x w + w ∂ z w + ∂ z p = ρf z + ∂ x τ zx + ∂ z τ zz , (1)
where τ is the deviatoric stress tensor, p the pressure, and f = (f x , f z ) is the external force, usually representing the gravity force. For instance, if (x, z) are the usual Cartesian coordinates then f = (0, -g). Other option is to consider a tilted reference system over a fixed inclined plane with slope θ > 0, then f = (g sin θ, -g cos θ).

We shall consider for now the decomposition of the pressure into hydrostatic and non-hydrostatic components:

p = ρ (p H + p N H ) , with ∂ z p H = f z .
We denote by h(t, x) the fluid depth and we define b(x) the bottom topography. We shall denote by η(t, x) = h(t, x) + b(x) the free surface. As boundary conditions, we take the usual kinematic conditions at the bottom (z = b) and free surface (z = η):

u | z=b ∂ x b = w | z=b , ∂ t η + u |z=η ∂ x η = w |z=η .
In addition, we assume the usual normal stress balance at the free surface and a Darcy friction law at the bottom (see e.g. [START_REF] Garres-Díaz | Shallow water moment models for bedload transport problems[END_REF]):

τ xz| z=b+h = 0 and τ xz| z=b = ρ u | z=b u | z=b ,
where is a dimensionless constant friction coefficient. Now, we write system (1) in the so-called σ-coordinates. To this end, we consider the change of variable given by

ξ = z -b h , where ξ ∈ [0, 1] for z ∈ [b, η].
Notice that ξ = 0 (respectively ξ = 1) corresponds to the bottom z = b (respectively free surface z = η) level. Then, denoting by ψ(t, x, ξ) = ψ(t, x, ξh + b), the differential operators read

∂ ξ ψ = h∂ z ψ and h∂ s ψ = ∂ s h ψ -∂ ξ ∂ s (ξh + b) ψ , for s ∈ {t, x}.
Taking into account the new variables, the pressure is written as

p = ρ ( p H + p N H ) , with ∂ ξ p H = hf z ,
where it is clear that, assuming that the pressure is null at the surface, we obtain the expression for the hydrostatic pressure p H = hf z (ξ -1) .

Thus, the Navier-Stokes system in σ-coordinates reads

           ∂ x (h u) + ∂ ξ ( w -u∂ x (ξh + b)) = 0, ∂ t (h u) + ∂ x (h u 2 + h p N H ) -hf z ∂ x (z b + h) -∂ ξ (ω u -p N H ∂ x (hξ + z b )) = ∂ x (h τ xx ) + ∂ ξ ( τ xz -τ xx ∂ x (hξ + b)) ∂ t (h w) + ∂ x (h u w) -∂ ξ (ω w + p N H ) = ∂ x (h τ zx ) + ∂ ξ ( τ zz -τ zx ∂ x (hξ + b)) , (2a)
where by simplicity we consider here τ = ν (∇u + ∇u ), with ν = µ/ρ a constant kinematic viscosity.

In addition, z b = b + x f x /f z include all gravitational contributions, and ω corresponds to the vertical transference term written as

ω(ξ) = ∂ t (hξ + b) + u(ξ)∂ x (hξ + b) -w(ξ). (2b)
Notice that we have used here u 2 = u 2 and uw = u w, and in the third equation of (2a) the hydrostatic pressure has been canceled with the external force.

As a simple case, we shall consider τ xz = ν∂ z u, ant then the boundary conditions are rewritten as

ω(0) = u | ξ=0 ∂ x b -w | ξ=0 = 0, ω(1) = ∂ t η + u | ξ=1 ∂ x η -w | ξ=1 = 0, and
ν h (∂ ξ u) | ξ=0 = u | ξ=0 u | ξ=0 , ν h (∂ ξ u) | ξ=1 = 0. ( 2c
)
In the following, for the sake of simplicity, we focus on the hydrostatic Euler system, and therefore the vertical momentum equation is replaced by the hydrostatic assumption and viscous effects are neglected. However, in the numerical test we will be interested on showing flows with vertical structure, which is produced by friction forces and vertical viscosity effects mainly. To do so, we eventually consider the vertical term ∂ z τ xz in the model in usual multilayer and moment discretization, which has been already introduced in previous works and we detail it in appendixes Appendix A and Appendix B.

A general vertical decomposition of Euler equations

In this section, we perform a vertical discretization for system (2). As we said in Section 1, the first main goal is to propose a general framework that generalizes the moment [START_REF] Kowalski | Moment approximations and model cascades for shallow flow[END_REF] and multilayer [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF] approaches for shallow flows. Moreover, a family of models enhancing the vertical structure of the fluid is proposed.

For the sake of simplicity, let us assume for now that (x, z) are the Cartesian coordinates and therefore f = (0, -g) t . Recall that η = h + b denotes the free surface and we assume the pressure to be hydrostatic, p = gh(1 -ξ). This means in practice that the equation for w is no longer considered, as w may be computed from the divergence free condition. Moreover, we shall neglect viscous terms for now, ν = 0:

       ∂ x (h u) + ∂ ξ ( w -u∂ x (ξh + b)) = 0, ∂ t (h u) + ∂ x (h u 2 ) + gh∂ x η -∂ ξ ( uω) = 0, ω = ∂ t (hξ + b) + u∂ x (hξ + b) -w, (3)
together with the kinematic conditions

ω | ξ=0 = ω | ξ=1 = 0.
In what follows we shall omit the notation (•) and assume directly the dependence on (t, x, ξ). Moreover, using the definition of ω, the first equation in ( 3) is rewritten as

∂ t h + ∂ x (hu) = ∂ ξ ω, (4)
which is integrated for ξ ∈ [0, 1] to get

∂ t h + ∂ x (hu) = 0, where u = 1 0 udξ. (5)
In addition, using previous equations ( 4) and ( 5), we find an explicit expression, in terms of the spatial derivative, for the vertical transference term ω as

ω | ξ = ξ 0 (∂ t h + ∂ x (hu)) ds = -ξ∂ x (hu) + ∂ x h ξ 0 u(s)ds = ∂ x h ξ 0 u | ξ=s ds -ξu .
Notice that ω is continuous as long as u is an integrable function.

Thus, the σ-coordinate Euler equations (3) may be equivalently written as

           ∂ t h + ∂ x (hu) -∂ ξ ω = 0, ∂ t (hu) + ∂ x (hu 2 ) + gh∂ x η -∂ ξ (uω) = 0, ω = ∂ x h ξ 0 u | ξ=s ds -ξu . (6)
We consider now a vertical partition of the interval [0, 1] into M layers, that is, we consider

0 = L 1/2 < • • • < L α+1/2 < • • • < L M +1/2 = 1,
and define l α = L α+1/2 -L α-1/2 ∈ [0, 1], α = 1, . . . , M . Note that it holds then that

L α+1/2 = α β=0 l β , α = 1, 2, . . . , M, and
L M +1/2 = M α=1 l α = 1.
We shall define the layer Ω α as

Ω α = {(x, ξ) ∈ R × [0, 1] : L α-1/2 < ξ < L α+1/2 }, for α = 1, . . . , M.
Let U be a given general variable and consider its polynomial approximation within each layer:

U (t, x, ξ) = M α=1 U α (t, x, ξ)1 Ωα (ξ), with U α (t, x, •) ∈ P N [ξ],
where 1 Ωα (ξ) is the characteristic function in the subdomain Ω α . Consider now an orthogonal basis of

P N [ξ] in the interval [0, 1], consisting of the polynomials φ j : [0, 1] → R, j = 0, 1, . . . , N,
and rescale them in the interval [L α-1/2 , L α+1/2 ], defining

φ α,j (ξ) = φ j ξ -L α-1/2 l α , ξ ∈ [L α-1/2 , L α+1/2 ].
We shall consider here the particular case of the Legendre polynomial basis, although any other choice could be made. In particular, the first polynomials of the basis are

φ 0 (ξ) = 1, φ 1 (ξ) = 1 -2ξ, φ 2 (ξ) = 1 -6ξ + 6ξ 2 , φ 3 (ξ) = 1 -12ξ + 30ξ 2 -20ξ 3 .
Notice that these polynomials satisfy φ i (0) = 1, and

1 0 φ i (ξ)φ j (ξ) dξ = µ i δ ij ,
where δ ij is the Kronecker symbol and µ i is the square norm of φ i in P N [ξ], which is not necessarily 1 as the basis is assumed only to be orthogonal.

For the sake of generality, let us rewrite the system (6) as

             ∂ t h + ∂ x (hu) -∂ ξ ω = 0, ∂ t U + ∂ x F (h, U ) + ∂ ξ H(h, U ) = S(h, U )∂ x η, ω = ∂ x h ξ 0 u | ξ=s ds -ξu , (7a) (7b) (7c)
with

U = hu, F (h, U ) = u U, H(h, U, ω) = - 1 h ω U, S(h, U ) = -gh.
Remark 1. Although in our case it is quite simple, since U = hu, writing the system in the form [START_REF] Boulanger | Analytical solutions for the free surface hydrostatic Euler equations[END_REF] allows us to consider any other generic variable. For instance, in the case of non-hydrostatic pressure it would be interesting to consider U = (hu, hw) t or in the presence of a passive scalar c(t, x, ξ) transported by the fluid, we would use U = (hu, hc). Therefore, the derivation preformed in what follows will be as general as possible.

Let us denote by U α,j the coordinates of U α (t, x, •) in the basis { φ α,j } j=0,...,N :

U α (t, x, ξ) = N j=0 U α,j (t, x) φ α,j (ξ), for ξ ∈ [L α-1/2 , L α+1/2 ].
More explicitly, in our case we are considering U α,j = h u α,j .

We recall that, since the basis is orthogonal, we have

L α+1/2 L α-1/2 U α φ α,j dξ = l α µ j U α,j . (8)
Let us introduce the following definitions, as they will be useful afterwards

A ijk = 1 0 φ i φ j φ k dξ, B ijk = 1 0 φ i (ξ) ξ 0 φ j (s)ds φ k (ξ)dξ, C ijk = 1 0 φ i φ j φ k dξ, µ j = A 0jj
Remark 2. We have that the constant values A ijk coincide for any permutation of the indices i, j, k, , and the following relations hold:

A ij0 = µ j δ ij , B 0jk = 0, B ij0 = φ i (1)µ 0 δ j0 -µ j δ ij , C ijk = C ikj , C 0jk = 0.
Moreover, it is easy to show that

L α+1/2 L α-1/2 φ i φ j φ k dξ = L α+1/2 L α-1/2 (φ i φ j φ k ) ξ -L α-1/2 l α dξ = 1 0 (φ i φ j φ k )(s) l α ds = l α A ijk , L α+1/2 L α-1/2 φ i (ξ) ξ L α-1/2 φ j (s)ds φ k (ξ)dξ = L α+1/2 L α-1/2 1 l α φ i ξ -L α-1/2 l α   l α ξ-L α-1/2 lα 0 φ j (s)ds   φ k ξ -L α-1/2 l α dξ = l α 1 0 φ i (s) s 0 φ j (z) dz φ k (s)ds = l α B ijk , and
L α+1/2 L α-1/2 φ i φ j φ k dξ = L α+1/2 L α-1/2 φ i ξ -L α-1/2 l α 1 l α (φ j φ k ) ξ -L α-1/2 l α dξ = 1 0 (φ i φ j φ k )(s)ds = C ijk .
In the case of Legendre polynomials, we also have

φ j (L α-1/2 ) = φ j (0) = 1, φ j (L α+1/2 ) = φ j (1) = (-1) j .
In order to derive the final model, we first compute the mass conservation equation by integrating (7a). Afterwards we shall multiply the momentum equation (7b) by the basis functions φ α,i and perform a vertical integration in the layer. The integrated mass equation reads

l α ∂ t h + l α ∂ x (h u α,0 ) = ω α+1/2 -ω α-1/2 , for α = 1, . . . , M, (9)
where we have used that 1 0 φ j dξ = δ j0 . Summing the previous equations for α = 1, . . . , M , or equivalently integrating the mass equation along the vertical domain ξ ∈ [0, 1], we recover equation ( 5)

∂ t h + ∂ x (hu) = 0, where u = 1 0 u dξ = M β=1 l γ u β,0 . (10)
Moreover, using the expansion of U , the vertical transference term ω (7c), for ξ ∈ [L α-1/2 , L α+1/2 ], can be expressed as

ω = α-1 β=1 ∂ x h L β+1/2 L β-1/2 udξ + ∂ x h ξ L α-1/2 udξ -ξ∂ x (hu) = ω α-1/2 + L α-1/2 -ξ ∂ x (hu) + ∂ x h ξ L α-1/2 udξ (11)
with

ω α-1/2 = α-1 β=1 l β ∂ x (h ( u β,0 -u)) .
Let us focus now on the momentum equation (7b), which is multiplied by the basis function φ α,i and then integrated within each layer. We shall compute each term separately.

For the time derivative, from (8) we have

L α+1/2 L α-1/2 ∂ t U φ α,i dξ = l α µ i ∂ t U α,i . (12)
The derivative of the convective flux yields

L α+1/2 L α-1/2 ∂ x F (h, U ) φ α,i dξ = L α+1/2 L α-1/2 ∂ x (u U ) φ α,i dξ = ∂ x L α+1/2 L α-1/2 N j,k=0 u α,j U α,k φ α,j φ α,k φ α,i dξ = l α N j,k=0 A ijk ∂ x u α,j U α,k . (13)
The source term gives

L α+1/2 L α-1/2 S(h, U )∂ x η φ α,i dξ = -gh δ i0 ∂ x η. (14)
Finally, the vertical transference term gives:

L α+1/2 L α-1/2 ∂ ξ H(h, U, ω) φ α,i dξ = φ α,i (1)H | ξ=L - α+1/2 -φ α,i (0) H | ξ=L + α-1/2 - L α+1/2 L α-1/2 H(h, U, ω) φ α,i dξ. ( 15
)
In the particular case of the Legendre basis, we have φ α,i (1) = (-1) i and φ α,i (0) = 1. Let us further develop the last integral in [START_REF] Chorin | A Mathematical Introduction to Fluid Mechanics[END_REF] taking into account [START_REF] Castro Díaz | A class of computationally fast first order finite volume solvers: PVM methods[END_REF]. Then,

L α+1/2 L α-1/2 H(h, U, ω) φ α,i dξ = - L α+1/2 L α-1/2 1 h ωU φ α,i dξ = -ω α-1/2 N j=0 1 h U α,j L α+1/2 L α-1/2 φ α,i φ α,j dξ -∂ x (hu) N j=0 1 h U α,j L α+1/2 L α-1/2 L α-1/2 -ξ φ α,i φ α,j dξ - N j,k=0 1 h U α,k ∂ x (h u α,j ) L α+1/2 L α-1/2 φ α,i ξ L α-1/2
φ α,j (s)ds φ α,k dξ.

Now, taking into account that

L α+1/2 L α-1/2 L α-1/2 -ξ φ α,i φ α,j dξ = -l α 1 0 s φ i (s)φ j (s) ds = -l α B i0j ,
we have

L α+1/2 L α-1/2 H(h, U, ω) φ α,i dξ = -ω α-1/2 N j=0 1 h U α,j C ij0 + l α ∂ x (hu) N j=0 1 h U α,j B i0j -l α N j,k=0 1 h U α,k ∂ x (h u α,j )B ijk . (16)
Combining ( 12), ( 13), ( 14), [START_REF] Chorin | A Mathematical Introduction to Fluid Mechanics[END_REF], and ( 16), for α = 1, 2, . . . , M and i = 0, 1, . . . , N , together with the mass equation ( 9), we get the final multilayer-moment system

                         l α ∂ t h + l α ∂ x (h u α,0 ) = ω α+1/2 -ω α-1/2 l α µ i ∂ t U α,i + l α N j,k=0 A ijk ∂ x u α,j U α,k + gh δ i0 e 1 ∂ x η = (-1) i 1 h (ωU ) | ξ=L - α+1/2 - 1 h (ωU ) | ξ=L + α-1/2 -ω α-1/2 N j=0 1 h U α,j C ij0 + l α ∂ x (hu) N j=0 1 h U α,j B i0j -l α N j,k=0 1 h U α,k ∂ x (h u α,j )B ijk , (17)
for α = 1, . . . , M , where we recall that ω 1/2 = ω M +1/2 = 0 and

ω = ω α-1/2 + L α-1/2 -ξ ∂ x (hu) + ∂ x h ξ L α-1/2 udξ with ω α-1/2 = α-1 β=1 l β ∂ x (h ( u β,0 -u)) .
Notice that the mass equation in previous system can be immediately replaced by [START_REF] Castro | Well-Balanced Schemes and Path-Conservative Numerical Methods[END_REF].

Finally, a definition of (ωU

) | ξ=L ± α+1/2
is needed to close the system. Thanks to the continuity of ω we have (ωU

) | ξ=L ± α+1/2 = ω α+1/2 U | ξ=L ± α+1/2 = ω α+1/2 U ± α+1/2
where here we must understand the limit values at the interfaces as

U - α+1/2 = lim ξ→L α+1/2 ξ<L α+1/2 U | Ωα = N j=0 U α,j φ α,j (L α+1/2 ) = N j=0 (-1) j U α,j , U + α+1/2 = lim ξ→L α+1/2 ξ>L α+1/2 U | Ω α+1 = N j=0 U α+1,j φ α+1,j (L α+1/2 ) = N j=0 U α+1,j .
Let us focus now on the particular case U = hu. We may now use a continuity argument at the interface for some particular combinations of the variable. This can be seen as an analogous procedure to what is done in [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF], where jump conditions associated to the weak formulation of the original system are used. The detailed procedures is given in Appendix Appendix A. Then, after some straightforward computations, the term

u - α+1/2 ω α+1/2 -u + α-1/2 ω α-1/2 in (15) is replaced by u α+1/2 ω α+1/2 -u α-1/2 ω α-1/2 (18)
where here u α+1/2 = u + α+1/2 + u - α+1/2 /2 is a centered approximation unless otherwise is specified (see Appendix Appendix A for details and in particular equation (A.4)).

Therefore, the final multilayer-moment system reads

                         l α ∂ t h + l α ∂ x (h u α,0 ) = ω α+1/2 -ω α-1/2 l α µ i ∂ t (h u α,i ) + l α N j,k=0 A ijk ∂ x (h u α,j u α,k ) + gh δ i0 e 1 ∂ x η = (-1) i u α+1/2 ω α+1/2 -u α-1/2 ω α-1/2 -ω α-1/2 N j=0 u α,j C ij0 + l α ∂ x (hu) N j=0 u α,j B i0j -l α N j,k=0 u α,k ∂ x (h u α,j )B ijk , for α = 1, . . . , M . Remark 3. Written in compact form, system (17) is equivalent to L α+1/2 L α-1/2 ∂ t U φ α,i dξ + L α+1/2 L α-1/2 ∂ x F (h, U ) φ α,i dξ - L α+1/2 L α-1/2 S(h, U )∂ x η φ α,i dξ + φ α,i (1)H | ξ=L - α+1/2 -φ α,i (0) H | ξ=L + α-1/2 - L α+1/2 L α-1/2 H(h, U, ω) φ α,i dξ = 0.
This may be seen as a semi-discrete vertical discretization of [START_REF] Boulanger | Analytical solutions for the free surface hydrostatic Euler equations[END_REF] in the framework of the Discontinuous Galerkin (DG) methods. Compared to other related works, here we are using a very particular modal basis: the Legendre polynomials, in contrast with different methodologies that employ any other modal or nodal basis (see [START_REF] Dumbser | ADER schemes on unstructured meshes for non-conservative hyperbolic systems: Applications to geophysical flows[END_REF][START_REF] Escalante | An efficient hyperbolic relaxation system for dispersive non-hydrostatic water waves and its solution with high order discontinuous Galerkin schemes[END_REF][START_REF] Guerrero Fernández | Well-balanced high-order discontinuous galerkin methods for systems of balance laws[END_REF][START_REF] Kopriva | On the quadrature and weak form choices in collocation type discontinuous galerkin spectral element methods[END_REF][START_REF] Zhong | A simple weighted essentially nonoscillatory limiter for rungekutta discontinuous galerkin methods[END_REF]). It is worth mentioning that, while the regular integrals are numerically approximated in other works (see for instance [START_REF] Escalante | An efficient hyperbolic relaxation system for dispersive non-hydrostatic water waves and its solution with high order discontinuous Galerkin schemes[END_REF][START_REF] Guerrero Fernández | Well-balanced high-order discontinuous galerkin methods for systems of balance laws[END_REF]), here we are computing those integrals exactly.

Particular cases

Let us consider now some particular cases of interest for [START_REF] Dumbser | ADER schemes on unstructured meshes for non-conservative hyperbolic systems: Applications to geophysical flows[END_REF]. Concretely, we shall see that the proposed discretization gives the multilayer [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF] and moment [START_REF] Kowalski | Moment approximations and model cascades for shallow flow[END_REF] models for particular choices of M and N , including the classical shallow water system. Finally, a new model is proposed, where the horizontal velocity is assumed to be linear within each vertical layer.

Shallow water system

The easiest case is when we only have a single layer (M = 1) and we use as polynomial approximation the constants (N = 0). Then, from system (17) we get

∂ t h + ∂ x (h u 1,0 ) = ω | ξ=1 -ω | ξ=0 , and
∂ t U 1,0 + ∂ x u 1,0 U 1,0 + gh∂ x η = 1 h (ωU ) | ξ=1 - 1 h (ωU ) | ξ=0 -ω | ξ=0 1 h U 1,0 C 000 + ∂ x (hu) 1 h U 1,0 B 000 - 1 h U 1,0 ∂ x (h u 1,0 )B 000 . Now using the fact that B 0jk = 0, C 0jk = 0, ω | ξ=1 = ω | ξ=0 = 0 and U 1,0 = 1 0 U dξ = U = hu,
then [START_REF] Dumbser | ADER schemes on unstructured meshes for non-conservative hyperbolic systems: Applications to geophysical flows[END_REF] gives us the classical shallow water system

∂ t h + ∂ x (hu) = 0, ∂ t (hu) + ∂ x (hu 2 ) + gh∂ x η = 0.
Multilayer zero-moment case: Multilayer Constant system We may consider now the case of M layers and constant polynomial approximation N = 0. In this case:

           l α ∂ t h + l α ∂ x (h u α,0 ) = ω α+1/2 -ω α-1/2 l α ∂ t U α,0 + l α ∂ x u α,0 U α,0 + gh ∂ x η = 1 h (ωU ) | ξ=L - α+1/2 - 1 h (ωU ) | ξ=L + α-1/2 -ω α-1/2 1 h U α,0 C 000 + l α ∂ x (hu) 1 h U α,0 B 000 -l α 1 h U α,0 ∂ x (h u α,0 )B 000 .
Taking into account now that B 0jk = C 0jk = 0, denoting the average within the layer by U α,0 = h u α,0 = hu α , and noticing that the definition of ω α±1/2 coincides with the transference term through the layer interfaces in [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF], i.e.,

ω α+1/2 = α β=1 ∂ x (h (u β -u)) , with u = M γ=1 l γ u γ ,
with ω 1/2 = ω M +1/2 = 0, system [START_REF] Dumbser | ADER schemes on unstructured meshes for non-conservative hyperbolic systems: Applications to geophysical flows[END_REF] reduces to the multilayer system in [START_REF] Audusse | A multilayer Saint-Venant system with mass exchanges for shallow water flows. derivation and numerical validation[END_REF][START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF]

l α (∂ t h + ∂ x (hu α )) = ω α+1/2 -ω α-1/2 , l α (∂ t (hu α ) + ∂ x (hu 2 α ) + gh∂ x (b + h)) = u α+1/2 ω α+1/2 -u α-1/2 ω α-1/2 , (19)
for α = 1, . . . , M , where we consider here

u α+1/2 = u - α+1/2 + u + α+1/2 2 = u α + u α+1 2 unless otherwise stated (see Appendix Appendix A).
It is remarkable the fact that the derivation of multilayer models in this framework is simpler than the usual one, since the intermediate interfaces L α+1/2 depend no more on x, which simplifies the integration procedure.

Although the focus here is on Euler equations, it will be of interest including viscous or friction effects for some simulations. Therefore, for the sake of completeness, Appendix Appendix A.2 gives the main ideas on how these terms are included for the multilayer system. Remark 4. In view of what has been said in Remark 3, the multilayer system [START_REF] Audusse | A multilayer Saint-Venant system with mass exchanges for shallow water flows. derivation and numerical validation[END_REF][START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF] corresponds to a semi-discrete vertical discretization in the DG framework, where the polynomial space corresponds to constant functions,

u | Ωα ∈ P 0 [ξ] .

Shallow water moment model

We now consider the case of one layer (M = 1) and arbitrary N , we get

                         ∂ t h + ∂ x (h u 1,0 ) = ω | ξ=1 -ω | ξ=0 , µ i ∂ t U 1,i + N j,k=0 A ijk ∂ x u 1,j U 1,k + gh δ i0 ∂ x η = (-1) i 1 h (ωU ) | ξ=1 - 1 h (ωU ) | ξ=0 -ω | ξ=0 N j=0 1 h U 1,j C ij0 + ∂ x (hu) N j=0 1 h U 1,j B i0j - N j,k=0 1 h U 1,k ∂ x (h u 1,j )B ijk ,
Now we use again that B 0jk = 0, C 0jk = 0, ω| ξ=0 = ω| ξ=1 = 0, the identity

A 0jk = µ j δ jk = δ jk 2j + 1
, and

U 1,0 = 1 0 U dξ = h u 1,0 = hu.
Let us slightly simplify the notation by removing the subscript of the single layer, setting then

U 1,i = h u 1,i = h u i ,
the mass equation and the momentum for i = 0 takes the form

∂ t h + ∂ x (hu) = 0 ∂ t (h u 0 ) + N j=0 ∂ x u 2 j 2j + 1 + gh ∂ x η = 0.
For i > 0, we get

µ i ∂ t (h u i ) + µ i ∂ x (2h u 0 u i ) + N j,k=1 A ijk ∂ x (h u j u k ) = ∂ x (h u 0 ) N j=0 u j B i0j - N j,k=0 u k ∂ x (h u j )B ijk .
We look now the term on the right hand side,

∂ x (h u 0 ) N j=0 u j B i0j - N j,k=0 u k ∂ x (h u j )B ijk = ∂ x (h u 0 ) N j=0 u j B i0j - N k=0 u k ∂ x (h u 0 )B i0k - N j=1 u 0 ∂ x (h u j )B ij0 - N j,k=1 u k ∂ x (h u j )B ijk = - N j=1 u 0 ∂ x (h u j )B ij0 - N j,k=1 u k ∂ x (h u j )B ijk = µ i u 0 ∂ x (h u i ) - N j,k=1 u k ∂ x (h u j )B ijk ,
where the identity B ij0 = φ i (1)µ 0 δ j0 -µ j δ ij for i > 0 has been used. Thus, we obtain the usual writing of the inviscid shallow water moment model [START_REF] Kowalski | Moment approximations and model cascades for shallow flow[END_REF] with i = 1, . . . , N :

                 ∂ t h + ∂ x (h u 0 ) = 0, ∂ t (h u 0 ) + ∂ x h u 2 0 + h N j=1 u 2 j 2j + 1 + gh ∂ x η = 0, ∂ t (h u i ) + ∂ x 2h u 0 u i + h 1 µ i N j,k=1 A ijk u j u k = u 0 ∂ x (h u i ) - 1 µ i N j,k=1 B ijk u k ∂ x (h u j ). ( 20
)
Remark 5. In view of what has been said in Remark 3, the moment model [START_REF] Kowalski | Moment approximations and model cascades for shallow flow[END_REF] corresponds to a semi-discrete vertical discretization with a single element on the vertical direction, and the polynomial space corresponds to the Legendre polynomial basis of degree

N , u ∈ P N [ξ].
As in the multilayer case, we are going to perform some numerical tests involving viscous effects. Therefore, we give the details on how these terms are included in the moment model in Appendix Appendix B.

Multilayer first-moment case: Multilayer Linear system

From a practical point of view, and due to the complexity of the general system (17) from the numerical and implementation points of view, we will limit ourselves to the case M layers and N = 1, that is, a first order expansion of U α (t, x, •) inside of each vertical layer. The resulting model is a multilayer model where the velocity is linear on ξ at each layer:

U (t, x, ξ) = N α=1 U α (t, x, ξ)1 Ωα (ξ), and U α (t, x, ξ) = U α,0 (t, x) + U α,1 (t, x) φ α,1 (ξ) ,
In this case, the coefficients appearing in the system are

A 0jk = µ j δ jk , A 100 = A 111 = 0, A 110 = A 101 = µ 1 = 1 3 , C 0ij = B 0ij = 0, C 100 = -2, C 110 = 0, B 100 = -1, B 101 = µ 1 , B 110 = -µ 1 , B 111 = 0,
and the moment equation for i = 1 is

l α µ 1 ∂ t U α,1 + l α µ 1 ∂ x u α,0 U α,1 + u α,1 U α,0 = - 1 h (ωU ) | ξ=L - α+1/2 - 1 h (ωU ) | ξ=L + α-1/2 -ω α-1/2 1 h U α,0 (-2) + l α ∂ x (hu) 1 h U α,0 (-1) + U α,1 µ 1 -l α 1 h U α,0 (-1) + U α,1 µ 1 ∂ x (h u α,0 ) -l α 1 h U α,0 (-µ 1 ) ∂ x (h u α,1 ).
Let us rewrite the term involving coefficient B 1ij , C 1ij as

l α µ 1 1 h U α,0 ∂ x (h u α,1 ) + 2ω α-1/2 1 h U α,0 + 1 h U α,0 -µ 1 U α,1 ω α+1/2 -ω α-1/2 = l α µ 1 1 h U α,0 ∂ x (h u α,1 ) + 1 h U α,0 -µ 1 U α,1 ω α+1/2 + 1 h U α,0 + µ 1 U α,1 ω α-1/2 .
Then, the final system, for α = 1, . . . , M , reads

                       l α ∂ t h + l α ∂ x (h u α,0 ) = ω α+1/2 -ω α-1/2 , l α ∂ t U α,0 + l α ∂ x u α,0 U α,0 + µ 1 u α,1 U α,1 + gh∂ x η = 1 h (ωU ) | ξ=L - α+1/2 - 1 h (ωU ) | ξ=L + α-1/2 , l α µ 1 ∂ t U α,1 + l α µ 1 ∂ x u α,0 U α,1 + u α,1 U α,0 = - 1 h (ωU ) | ξ=L - α+1/2 - 1 h (ωU ) | ξ=L + α-1/2 +l α µ 1 1 h U α,0 ∂ x (h u α,1 ) + 1 h U α,0 -µ 1 U α,1 ω α+1/2 + 1 h U α,0 + µ 1 U α,1 ω α-1/2 ,
Finally, let us set the system involving just the horizontal velocity, and the mass equation obtained when summing all the layers

                     ∂ t h + l α ∂ x h M α=1 u α,0 = 0, l α ∂ t (h u α,0 ) + l α ∂ x h u 2 α,0 + h u 2 α,1 3 + gh∂ x η = ω α+1/2 u α+1/2 -ω α-1/2 u α-1/2 , l α µ 1 ∂ t (h u α,1 ) + l α µ 1 ∂ x (2h u α,0 u α,1 ) = -ω α+1/2 u α+1/2 -ω α-1/2 u α-1/2 +l α µ 1 u α,0 ∂ x (h u α,1 ) + ( u α,0 -µ 1 u α,1 ) ω α+1/2 + ( u α,0 + µ 1 u α,1 ) ω α-1/2 , (21)
for α = 1, . . . , M , where

u α+1/2 = u - α+1/2 + u + α+1/2 2 = u α,0 + u α+1,0 2 + u α+1,1 -u α,1 2
unless otherwise stated (see Appendix Appendix A).

Numerical results

We present in this section some numerical tests showing an interesting comparison between multilayer and moment models. Firstly we perform a simple test including viscosity and friction force to have diversity in the velocity profiles. Next, some more tests in the inviscid case are performed to include also the results of the first order multilayer-moment model derived in Subsection 3.1.

In the next subsections, we perform two families of tests. First, subsections 4.1 and 4.2 aim to give a comparison between usual multilayer (A.5) and moment models (B.1) commonly used, in particular in the viscous case to obtain a natural vertical structure. Concerning the numerical scheme, we have carried out the simulation using a standard first order path-conservative HLL finite volume scheme as described in [START_REF] Castro Díaz | A class of computationally fast first order finite volume solvers: PVM methods[END_REF].

Second, subsections 4.3 and 4.4 show the convergence of the multilayer models with constant (zero-moment) [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I)-rheology: application to granular collapse on erodible beds[END_REF] and linear (first-moment) (21) horizontal velocity to stationary solutions of hydrostatic Euler equations, as well as an accuracy test. For those tests, we consider a second order path-conservative HLL scheme based on a MUSCL reconstruction [START_REF] Castro | Well-Balanced Schemes and Path-Conservative Numerical Methods[END_REF].

Multilayer and moment model in subcritical viscous flow

Here we compare the results of the standard multilayer and moment models in a viscous subcritical flow. The viscous and friction forces generate complex and different velocity profiles along the domain, from almost constant to parabolic profiles and recirculations in the flow. This test aims to show how multilayer and moment models reproduce these profiles and explore each approach's limitations and advantages/disadvantages. To this aim, we consider as a reference solution the one obtained with a very fine discretization computed by the multilayer model with 640 vertical layers. We consider a classical subcritical test in a domain x ∈ [-25 , 25], where the initial free surface and discharge are η(x) = 5 m and q(x) = 4.42 m 2 /s respectively, and the bottom is b(x) = 0.05 -0.001x + 2 cos 2 πx 10 m if |x| < 5, 0.05 -0.001x m otherwise.

We use 200 point in the horizontal discretization, with CFL = 0.8, and subcritical boundary conditions (inlet discharge q in = 4.42 m 2 /s and outlet free surface η out = 5 m). The final time is t f = 1000s. We show the results for a low (ν = 0.001 m 2 /s) and high kinematic viscosity (ν = 1 m 2 /s) case, with strong friction = 100ν. For each case we show the approximation of the free surface and the errors made by the moment and the multilayer approaches. To this aim, we shall measure relative errors as follows:

Consider an horizontal discretization of the domain I = [a, b] by means of N x volumes of constant size ∆x = b -a N x ,

I i = [x i-1/2 , x i+1/2 ], for i = 1, 2, . . . , N x ,
where x 1/2 = a, and x i+1/2 = x i-1/2 +∆x for i = 1, 2, . . . , N x . We will denote as well x i = x i-1/2 + x i+1/2 2 the center of the cell I i .

Consider as well a vertical partition of the interval [0, 1] into M layers

0 = L 1/2 < • • • < L α+1/2 < • • • < L M +1/2 = 1.
For any generic variable v denote by v α,i the approximation obtained for v in the cell i and layer α. Then we shall denote by

v 2 = ∆x Nx i=1 h i M M α=1 v 2 α,i , and v ∞ = max α=1,...,M i=1,...,Nx |v α,i | .
Remark that in the particular case that v does not depend on the vertical variable ξ, then we simply consider v α,i ≡ v i for α = 1, . . . , M . We then compute the relative errors by

Err p [v] = v -v ref p v ref p , for p = 2, ∞,
where v ref is a given reference solution.

We will be interested as well in computing the errors in the vertical direction for a given cell I i , so that we denote by

v I i ,2 = h i M M α=1 v 2 α,i , and v I i ,∞ = max α=1,...,M |v α,i | ,
and we compute the relative errors

Err z p [v, x i ] = v -v ref I i ,p v ref I i ,p , for p = 2, ∞.
Here, the solution of the multilayer model, obtained with 700 nodes and 640 vertical layers, will be considered as the reference solution. Then, when computing the errors, v ref should be understood as the projection of the reference solution onto the space of the local solution v. For the case of moment models, the vertical velocity profile is computed, and its value at the centre of the corresponding layer is obtained for the errors.

In Figure 1 we see the steady state solution reached by both approaches for the free surface at the final time for both considered values of the viscosity, together with vertical profiles of the velocity at some particular locations. On the one hand, in the low viscosity case (upper figure) and for the first half of the domain, the velocity profile is almost a straight profile except near the bottom, where the velocity is smaller due to the friction force. However, the velocity profile is more complex after the Err

2 [η] Err 2 [u] Err 2 [u] Err ∞ [η] Err ∞ [u] Err ∞ [u] M
Multilayer (M layers) 5 8.57×10 -3 1.04×10 -2 2.63×10 -1 1.39×10 -2 1.94×10 -2 5.24×10 -1 10 4.58×10 -3 5.65×10 -3 1.81×10 -1 7.61×10 -3 1.04×10 -2 4.40×10 -1 20 2.44×10 -3 3.05×10 -3 1.07×10 -1 4.01×10 -3 5.51×10 -3 3.42×10 -1 40 1.13×10 -3 1.46×10 -3 5.27×10 -2 1.86×10 -3 2.58×10 -3 2.03×10 -1 80 5.17×10 -4 6.89×10 -4 2.45×10 -2 8.47×10 -4 1.20×10 -3 1.08×10 -1 160 2.22×10 -4 3.07×10 -4 1.06×10 -2 3.63×10 -4 5.34×10 -4 5.60×10 -2 320 7.44×10 -5 1.06×10 -4 3.54×10 -3 1.21×10 -4 1.83×10 -4 2.27×10 -2 N Moment (order N)

1 2.65×10 -3 2.93×10 -3 2.09×10 -1 5.88×10 -3 3.63×10 -3 4.38×10 -1 2 2.10×10 -3 2.86×10 -3 1.62×10 -1 6.90×10 -3 7.33×10 -3 4.14×10 -1 3 1.04×10 -3 1.66×10 -3 1.18×10 -1 3.66×10 -3 4.84×10 -3 3.64×10 -1 4 1.15×10 -3 1.45×10 -3 8.11×10 -2 2.09×10 -3 2.60×10 -3 2.77×10 -1 5 9.63×10 -4 1.23×10 -3 5.42×10 -2 1.61×10 -3 2.19×10 -3 1.94×10 -1
Table 1: Relative errors in the low viscosity case made on the free surface (η), averaged velocity (u) and the velocity field (u = u(x, z)) by the multilayer models (M layers) and the moment models (order N ).

bump, and we find recirculation areas (change on the sign of the velocity along the vertical direction) near the bottom. On the other hand, for the high viscosity case (lower figure), the velocity follows a parabolic profile along the whole domain.

Figures 2 and3 show the free surface computed with both the multilayer and moment approach, and Tables 1 and2 show the relative errors made on the free surface, averaged velocity and the full velocity Figure 3: Free surface approximation with moment and multilayer models for the high viscosity case..

Err 2 [η] Err 2 [u] Err 2 [u] Err ∞ [η] Err ∞ [u] Err ∞ [u] M
Multilayer (M layers) 5 1.83×10 -2 3.79×10 -2 1.23×10 -1 3.19×10 -2 4.63×10 -2 1.85×10 -1 10 8.13×10 -3 1.81×10 -2 5.62×10 -2 1.42×10 -2 2.19×10 -2 1.30×10 -1 20 3.82×10 -3 8.50×10 -3 2.65×10 -2 6.63×10 -3 1.04×10 -2 8.79×10 -2 40 1.81×10 -3 3.95×10 -3 1.25×10 -2 3.14×10 -3 4.83×10 -3 5.16×10 -2 80 8.40×10 -4 1.79×10 -3 5.79×10 -3 1.45×10 -3 2.20×10 -3 2.68×10 -2 160 3.58×10 -4 7.56×10 -4 2.47×10 -3 6.18×10 -4 9.32×10 -4 1.22×10 -2 320 1.19×10 -4 2.50×10 -4 8.21×10 -4 2.06×10 -4 3.08×10 -4 4.18×10 -3 N Moment (order N)

1 1.49×10 -2 3.98×10 -2 1.49×10 -1 2.76×10 -2 4.72×10 -2 2.81×10 -1 2 1.05×10 -3 9.67×10 -3 2.30×10 -2 2.57×10 -3 1.31×10 -2 1.14×10 -1 3 4.27×10 -4 4.26×10 -3 8.28×10 -3 1.16×10 -3 5.62×10 -3 7.77×10 -2 4 2.54×10 -4 1.79×10 -3 3.98×10 -3 6.24×10 -4 2.34×10 -3 4.87×10 -2 5 1.79×10 -4 6.65×10 -4 2.12×10 -3 3.74×10 -4 8.60×10 -4 2.83×10 -2
Table 2: Relative errors in the high viscosity case made on the free surface (η), averaged velocity (u) and the velocity field (u = u(x, z)) by the multilayer models (M layers) and the moment models (order N ).

field, measured with respect to the reference solution. In Figure 2 we see the convergence to the reference solution in the multilayer case. In the case of moment models, although the convergence towards the reference solution is not fully obtained for N = 5, especially in the left region. That is due to the fact that the vertical profiles of velocity are straight lines almost everywhere in the vertical direction, and moment models do not properly reproduce these profiles, as we discuss later. Nevertheless, we M=5 M=10 M=20 M=40 M=80 M=160 M=320 t multi(M ) /t moment(5th) see that the errors made by the moment models are of the same order of magnitude as the multilayer models. Actually, we see that the errors made by the multilayer model with 20 layers and the 2nd order moment model are similar, while the computational cost for the former is much bigger than for the latter. Moreover, to improve the 5th order moment model results, we need to use 80 layers in the multilayer system. The computational effort for both approaches is shown in Figure 4, where we see that moment models are cheaper than multilayer models from the computational point of view. Figure 3 and Table 2 show the results for the high viscosity case, where the vertical profile of velocity is approximately a parabolic profile along the whole domain. Looking at Figure 3, we observe that there is convergence to the reference solution for free surface. It can be stated that the errors for the moment models appearing in Table 2 are smaller than in the previous case, in particular for the velocity field. That is due to the better approximation of the velocity profile along the domain since they are close to a parabolic profile, which is well approximated by polynomials and, therefore, by the moment approach. Note that in this case, the 3rd order moment model improves the results of the multilayer model with 40 layers, and the 5th order moment model produces smaller errors than the multilayer model with 160 layers.

That is an interesting result since the computational time increases fast with the number of layers. This can be seen in Table 3, where we show the ratios between the elapsed times for the multilayer models and the 3rd, 4th and 5th order moment models, obtained when simulating up to the final time t f = 1000 s. We see for instance that the 5th order moment model is 6 times faster when compared to the multilayer model with 10 layers and 25 times faster when using 40 vertical layers. This ratio increases to up to 62 if we double to 80 the number of layers.

In Figure 4 we see the comparison of the wall-clock time with the L 2 relative errors on the velocity. We choose the error on the velocity field because it is the variable for which the errors are the largest.

20

In the low viscosity case, we see that the 5th order moment model is comparable to the multilayer with 40 layers and better than the multilayer with 160 layers, in the sense that these errors are lower. Let us now focus on the different profiles of velocity that we observe in these cases and how multilayer and moment models approximate them. We obtain three kinds of profiles depending on the magnitude of vertical effects. First, in the low viscosity case, we obtain a straight profile for a wide zone of the vertical domain (left-hand side plots in Figure 5), and the velocity decreases close to the bottom because of the friction force. The low order moment models badly approximate these kind of profiles as it is difficult to fit them with a first or second order polynomial, due to the lack of sufficient vertical resolution. If we

Err z 2 [u, -5] Err z 2 [u, 5] Err z 2 [u, 0] Err z ∞ [u, -5] Err z ∞ [u, 5] Err z ∞ [u, 0] (ν = 10 -3 ) (ν = 10 -3 ) (ν = 1) (ν = 10 -3 ) (ν = 10 -3 ) (ν = 1) M Multilayer (M layers) 5 2.22×10 -1 3.51×10 -1 1.25×10 -1 4.40×10 -1 3.80×10 -1 1.83×10 -1 10 1.33×10 -1 2.25×10 -1 5.66×10 -2 3.83×10 -1 2.67×10 -1 9.40×10 -2 20 6.70×10 -2 1.29×10 -1 2.67×10 -2 2.61×10 -1 1.56×10 -1 4.72×10 -2 40 3.24×10 -2 6.33×10 -2 1.26×10 -2 1.60×10 -1 7.81×10 -2 2.31×10 -2 80 1.49×10 -2 2.94×10 -2 5.84×10 -3 8.69×10 -2 3.66×10 -2 1.08×10 -2 160 6.26×10 -3 1.27×10 -2 2.49×10 -3 4.02×10 -2 1.59×10 -2 4.65×10 -3 320 2.06×10 -3 4.28×10 -3 8.28×10 -4 1.39×10 -2 5.37×10 -3 1.55×10 -3 N Moment (order N) 1 2.51×10 -1 2.31×10 -1 1.57×10 -1 4.39×10 -1 3.84×10 -1 2.77×10 -1 2 1.89×10 -1 1.48×10 -1 1.03×10 -2 4.06×10 -1 2.97×10 -1 1.34×10 -2 3 1.22×10 -1 9.66×10 -2 5.00×10 -3 3.15×10 -1 1.42×10 -1 5.98×10 -3 4 7.49×10 -2 6.95×10 -2 2.28×10 -3 2.14×10 -1 1.04×10 -1 2.69×10 -3 5 4.34×10 -2 4.79×10 -2 1.25×10 -3 1.27×10 -1 7.26×10 -2 1.50×10 -3
Table 4: Relative errors made on the vertical profile of velocity for three different profiles: quasi-straight (x = -5, ν = 10 -3 m 2 /s), profile with recirculation (x = 5, ν = 10 -3 m 2 /s) and parabolic (x = 0, ν = 1 m 2 /s) profile. Results are shown for the multilayer (M layers) and the moment models (order N ).

increase the order of the moment model, then we observe the oscillating intrinsic nature of polynomials. This fact is in particular observed for the fifth order moment model in the left bottom picture in Figure 5. By increasing the order of the moment model N , the profiles are nearer to the reference solution.

Nevertheless, the oscillations are still present. Secondly, we see a profile of velocity corresponding to recirculation of the flow (central plots in Figure 5). These profiles are well reproduced by high order moment models (N = 4, 5), whereas low order models (N < 4) fail. Remark that the multilayer model needs more than 10 layers to reproduce this effect. Finally, we observe parabolic profiles (right hand side plots in Figure 5). As expected, these profiles are efficiently and properly reproduced by moment models for order N > 1. We also see the convergence to the reference solution in all these profiles. Table 4 shows the errors with respect the reference solution for these three profiles. We see that for the parabolic profile (x = 0), the moment model with N > 1 gives better results than multilayer models. The errors are larger for the other profiles (x = -5 or x = 5), but higher moments give good results, only comparable to using many layers in the multilayer case.

In conclusion, if the vertical profile of velocity is close to a polynomial, the moment approach allows us to obtain good results with less computational effort than the multilayer model, which needs many layers. However, if one needs to reproduce cases with large areas where a straight profile is found, the moment approach will oscillate due to its polynomial nature. That is eventually a limitation of the usual formulation of moment models. In this case, the multilayer model may be more appropriate, although a sufficiently large number of layers would still be needed.

Comparison of multilayer and moment models with Navier-Stokes

Let us now consider a numerical test included in [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF] (Section 5.2) where a qualitative comparison with a Navier-Stokes solver in a simplified configuration is performed. Concretely, the test is described in [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF] as follows. They carried out a simulation with the multilayer model, and the obtained stationary free surface is imposed as the fluid domain in the FreeFem++ solver [START_REF] Hecht | New development in FreeFem++[END_REF] to solve the pressure and velocity fields. The authors checked that the simulation was in the hydrostatic framework, and they compared the obtained velocity profiles with both the multilayer method and FreeFem++. A reasonable agreement between the obtained profiles was found, where the differences are most undoubtedly due to how boundary conditions are imposed in the finite element solver and the multilayer model with a finite volume method.

To compare with the multilayer and moment models, we consider a fluid domain [0, 3] discretized with 200 nodes and CF L = 0.8. In these simulations, a non-slip condition was used in the FreeFem++ code, so we consider here a large enough friction coefficient to reproduce this boundary condition. In practice, after a convergence analysis in , we use = 10 6 . For the viscosity, two values are considered, ν = 5 × 10 -4 and ν = 10 -3 m 2 /s. We write here the initial and boundary conditions: the flow is initially at rest, and the initial bottom and free surface are given by b

0 (x) = 0 if x ≤ 1, 0.5 (x -1) if 1 ≤ x ≤ 2, 0.5 if x ≥ 2, b 0 (x) + h 0 (x) = 1 m at t = 0 s,
and subcritical boundary conditions are imposed:

q(0, t) = 2z h m 2 /s at x = 0 m; h(3, t) + b(3) = 1.2 m at x = 3 m.
Figure 6 shows the fluid domain with the stationary free surface, once it has been reached, where we focus on the velocity profiles at locations x = 1.26 and 1.74 m. We obtain the velocity profiles in these two points and compare them with the corresponding finite element solution. First, we show the convergence of the velocity profiles for both the multilayer and moment approach to a reference solution computed with a large number of layers (see Figure 7), in the case ν = 10 -3 . We see that the velocity profiles in these points exhibit a parabolic shape near the bottom due to the friction force, and then they become a linear profile. As expected, both approaches reproduce well these profiles, although we see the oscillatory character of the moment approach, especially in the straight part of the vertical profile. The same behavior is observed when comparing with the FreeFem++ solutions for ν = 10 -3 , 5 × 10 -4 m 2 /s. We see that the obtained profiles are in good agreement with the finite volume solution, in concordance with [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF]. Similarly to what was said in that work, we conclude that the observed differences are most surely due to the way of imposing boundary conditions in both solvers, and also the rest of the characteristics of both solvers (mainly, non-hydrostatic pressure and stationary with the fixed free surface approach of the performed FreeFem++ simulation in [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF]).

Numerical convergence for stationary analytical solutions of the free-surface hydrostatic Euler equations

We focus now on comparing the usual multilayer model [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I)-rheology: application to granular collapse on erodible beds[END_REF] and the one introduced here with a layerwise linear velocity [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF]. We consider two discretizations for the velocity at the vertical interfaces (see Appendix Appendix A), a centred (C) definition

u α+1/2 = u + α+1/2 + u - α+1/2 2 ,
and an upwind (U) definition

u α+1/2 = u + α+1/2 + u - α+1/2 2 - sign -ω α+1/2 2 u + α+1/2 -u - α+1/2 .
From now on, we shall denote by MCC and MCU the Multilayer Constant Central and the Multilayer Constant Upwind models, that is, the usual multilayer models with either a central or upwind definition of the velocity at the interface. Similarly, MLC and MLU will stand for Multilayer Linear Central and Multilayer Linear Upwind models with layerwise linear velocities.

We consider a set of continuous stationary solutions for the hydrostatic Euler system (6) derived in [START_REF] Boulanger | Analytical solutions for the free surface hydrostatic Euler equations[END_REF]. More precisely, we consider the general continuous stationary solution given by

       z b (x) = zb -h(x) - α 2 β 2 2g sin 2 (βh(x)) , u(x, ξ) = αβ sin(βh(x)) cos(ξβh(x)), (22)
where α, β and zb are three real constants, and h(x) is a non-negative continuous function such that sin(βh(x)) = 0 for all x in the horizontal fluid domain. This test problem aims to inspect the abilities of the proposed vertical discretizations to preserve stationary solutions of the free-surface hydrostatic Euler equations and numerically state the expected order of accuracy.

To do that we consider x ∈ [-5, 5], α = β = -1, zb = 0, and

h(x) = 1.4 -0.8e -x 2 , (23)
for the total water depth (see Figure 9). Once the height is fixed, the bottom is given by [START_REF] Garres-Díaz | Flexible and efficient discretizations of multilayer models with variable density[END_REF]. Then, given a vertical discretization into M layers, we consider as initial condition for the horizontal velocity

u α (x, ξ) = u α,0 (x) φ α,0 (ξ) + u α,1 (x) φ α,1 (ξ), α ∈ {1, . . . , M },
the piecewise linear projection of the function u in [START_REF] Garres-Díaz | Flexible and efficient discretizations of multilayer models with variable density[END_REF] as given in [START_REF] Bristeau | Layer-averaged euler and Navier-Stokes equations[END_REF], i.e.,

u α,j (x) = 1 l α µ j L α+1/2 L α-1/2
u(x, ξ) φα,j (ξ) dξ, for j = 0, 1 and α ∈ {1, . . . , M }.

(

) 24
We then compute different numerical approximations for a final time t = 10 s with a set of refined meshes and other values of the number of layers M. Dirichlet boundary conditions have been used by imposing the values of the stationary solution at both ends of the domain. We consider in this case the L 1 absolute errors. Such errors are computed against the initial condition, that coincides with the L 1 projection of the exact solution into the space of linear polynomials in ξ. More precisely, we consider in this case the errors given by

h -h ref 1 = ∆x Nx i=1 h i -h ref i , u -u ref 1 = ∆x Nx,M i,α=1 L α+1/2 L α-1/2 u α,i (ξ) -u ref α,i (ξ) dξ, (25)
where the last integral will be computed by using the Simpson's rule. As it is said above, in this test the reference solution is given by the initial state.

Table 5 shows the L 1 errors and the numerical convergence rates. It is essential to make some remarks. First, note that a usual formal order analysis has been performed, where ∆x ≈ ∆z, and an increasing number of cells has been considered for computing the expected order of accuracy of the vertical and horizontal numerical approximation methods of free-surface Euler equations. Note that the desired order of accuracy is achieved for the MLU and MLC methods since they provide second order of accuracy in the vertical direction. Concerning the MC methods, note that MCC provides a second order of accuracy on the vertical, whereas the MCU is expected to provide first order of accuracy. Finally, note the efficiency of ML methods over MC. To get a L 1 error equal or less to 1 × 10 -2 for the velocity u, the MCU and MCC need more than M = 80 and M = 20 layers, respectively, whereas the MLU and MLC approximation reach errors less than 1 × 10 -2 with less than M = 30 and M = 20 layers. A deeper analysis of the efficiency, including elapsed times, will be done in the following subsection.

Perturbed stationary analytical solutions for the free-surface hydrostatic Euler equations

We consider now a perturbation of the stationary solution ( 22)- [START_REF] Garres-Díaz | Multilayer models for shallow two-phase debris flows with dilatancy effects[END_REF]. We take the same initial conditions except for the velocity, for which we assume u(x, ξ) = αβ sin(β (h(x) -0.1)) cos(ξβ (h(x) -0.1)).

The computational domain is now [-10, 10], and the parameters are set to α = β = -1, and zb = 0, as previously. In this case, periodic boundary conditions are imposed. We then compute different numerical approximations until a final time t = 1 s with a set of refined meshes in the vertical direction with different values of the number of layers M . Table 6 shows the L 1 errors and the numerical convergence rates. Errors are computed as in ( 25), but considering a reference solution that has been computed with the Multilayer Constant Central method with N x = 1600 elements and M = 160 layers. As it can be seen, the proposed methods provide an accurate approximation as expected, and similar conclusions as in the previous test case may be reached. It is important to emphasize the outstanding performance of ML methods compared to MC ones.

M N x h -h ref 1 Order h u -u ref 1 Order u h -h ref 1 Order h u -u ref 1 Order u MCU MCC 10 100 7.24×10 -3 - 1.76×10 -1 - 6.51×10 -3 - 5.19×10 -2 - 20 200
Table 7 shows the computational time needed by both models for several layers, as well as its ratio. In terms of the computational effort, the model MLU is approximately 1.5 times more expensive than the MCU for a given number of layers M . That is natural since M extra equations, corresponding to u α,1 , are solved in the layerwise linear velocity case. However, the gain of the accuracy is worth it by far. It can be concluded from the Tables 5 and6, and also from Figure 10, where we show the computational times against the error made on the velocity field, for the cases M = 5, 10, 20, 40 layers. In particular, we see that the MLU model with M = 5 layers improves the results of the MCU model with M = 40 layers, with a lower computational effort.

M N x u -u ref 1 order u u -u ref 1 order u u -u ref 1 order u u -u ref 1 order u MCU MCC MLU MLC 5 1600 3.42×10 -2 - 9.03×10 -3 - 7.91×10 -4 - 1.47×10 -3 - 10 1600 1.71×10 -2
1.00 2.29×10 ) and MLU (t MLU(M) ) models to reach t f = 10 s, and ratios, with M = 5, 10, 20, 40 layers and N x = 1600 cells.

Conclusions

A general framework for vertical discretizations of Euler equations has been proposed. It is based on considering a change to σ-coordinates and a partition on M layers along the vertical direction of the interval [0, 1]. Then a polynomial approximation of the variables, with arbitrary degree N , is assumed. The derivation procedure performed is general enough to extend it to other cases (e.g. non-hydrostatic pressure regimes, passive scalar transport, among many others). This model may be seen as a semi-discrete vertical discretization of the Euler system written in σ-coordinates in the framework of the Discontinuous Galerkin methods.

As particular cases of the discretization introduced here, we immediately recover, as simplified configurations, the usual multilayer and moment models appearing in the literature. Therefore, both can be seen as particular vertical discretizations under the same framework. In addition, an enriched multilayer-moment model has been proposed, particularly a multilayer model with a layerwise linear horizontal velocity.

In the numerical test cases, we have shown comparisons of usual moment a multilayer approaches, pointing out some advantages/limitations of both methods. We have shown that, in general, high order moment models allows us to recover the solution of multilayer models with a much smaller computational effort, so then they are much more efficient. It is especially relevant when the velocity field exhibits a parabolic profile along the whole domain. However, multilayer models may be more appropriate when the velocity profile has straight zones since moment models introduce oscillations due to the polynomial approximation. Finally, we have shown that the multilayer model with layerwise linear velocity proposed here preserves a fully 2D steady solution of hydrostatic Euler equations, different from water at rest, up to second order accuracy. Moreover, when introducing a perturbed initial condition, we recover the desired order of accuracy of both approaches: multilayer model with layerwise constant velocity is first order accurate. In contrast, second order is obtained with the proposed model considering a layerwise linear velocity. We have also measured the efficiency of the MLU model for the MCU model. Despite being 1.5 times more expensive in terms of computational efforts for a fixed number of layers, we have shown that it makes sense to use it in terms of accuracy. Actually, with the proposed layerwise linear velocity model, we need just a few layers to improve the results of the usual multilayer model with a high number of vertical layers. In the future, it would be interesting to extend the proposed model to the non-hydrostatic case.

Recall that, from the definition of ω, we get ω + α+1/2 = ω - α+1/2 . Then, integrating the second equation for s ∈ [0, ξ] we have

L α+1/2 L α-1/2 ∂ ξ (ωu -σ)dξ = ω α+1/2 u - α+1/2 -σ - α+1/2 -ω α-1/2 u + α-1/2 + σ + α-1/2 .
where σ - α+1/2 , σ + α-1/2 should be correctly defined. In [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF], this is done by means of a centered approximation σ α+1/2 at the interface:

σ α+1/2 = σ + α+1/2 + σ - α+1/2
2 .

Nevertheless, one could do other possible choices as we will see.

In general, a given approximation at the interfaces may be done by considering a vertical numerical flux G depending on the upper and lower limit values at the interface. For the sake of generality, let us consider

σ α+1/2 = G σ - α+1/2 ,

Figure 1 :

 1 Figure 1: Steady state solution for low (ν = 0.001 m 2 /s, upper figure) and high (ν = 1 m 2 /s, lower figure) viscosity cases. Blue lines correspond to the vertical profiles of velocity at locations x = -22.5, -15, -5, -2.5, 0, 2.5, 5, 15, 22.5 m.

Figure 2 :

 2 Figure 2: Free surface approximation with moment and multilayer models for the low viscosity case.

18

 18

Figure 4 :

 4 Figure 4: Wall-clock time compared to L 2 errors on the velocity field, in log scale, for multilayer and moment models.

Figure 5 :

 5 Figure 5: Vertical profiles of velocity computed with the multilayer (upper figures) and moment (lower figures) models. Left hand side and center figures correspond to the low viscosity case whereas right hand side figures corresponds to the high viscosity case.

Figure 6 :

 6 Figure 6: Fluid domain with velocity profiles at x = 1.26, 1.74 m for the case ν = 10 -3 m 2 /s.

Figure 7 :

 7 Figure 7: Convergence of vertical profiles of velocity computed with the multilayer (left figures) and moment (right figures) models at locations x = 1.26, 1.74 m.

Figure 8 :

 8 Figure 8: Vertical profiles of velocity computed with the multilayer (M = 80) and moment (N = 5) models (lines) and the FreeFem++ solution (symbols) at locations x = 1.26, 1.74, for two values of the viscosity ν = 10 -3 (left figures) and ν = 5 × 10 -4 m 2 /s (right figures).

Figure 9 :

 9 Figure9: Initial free-surface and vertical profiles of velocity u at locations x = -5, -3.5, -2, -0.5, 1, 2.5, and 4 m. for the stationary solution given by (23)-[START_REF] Garres-Díaz | Shallow water moment models for bedload transport problems[END_REF].

Figure 10 :

 10 Figure 10: Computational time and errors on the velocity field, in log scale, for the MCU (black circles) and MLU (blue squares) models at t f = 10 s, with M = 5, 10, 20, 40 layers and N x = 1600 cells.

(Remark 6 .

 6 uω + τ xz -∂ x (hξ + b)τ xx ) | ξ -(τ xz ) | ξ =0 -∂ x (b)(τ xx ) | ξ = ∂ t h ξ 0 u(s)ds + ∂ x h ξ 0 (u 2 -τ xx )(s) ds + ghξ∂ x ηNotice that in particular this means that (uω+ τ xz -∂ x (hξ + b)τ xx ) is continuous at the interfaces L α+1/2 . In particular, [uω + τ xz -∂ x (hξ + b)τ xx ] α+1/2 = 0, or equivalently [-τ xz + ∂ x (hξ + b)τ xx ] α+1/2 = [u] α+1/2 ω α+1/2 . (A.1) Using a similar framework as the one introduced in [21], continuity of ω and uω + τ xz -∂ x (hξ + b)τ xx at the interfaces may be interpreted as the jump conditions for the mass and the horizontal momentum conservation equations verified by a weak solution of the Navier-Stokes system. Denote by σ = -τ xz + ∂ x (hξ + b)τ xx and remark that

Table 3 :

 3

Ratio between the wall-clock time needed by the multilayer models with M layers (t multi(M ) ) and the 3rd, 4th and 5th order moment models (t moment(3th) ,t moment(4th) ,t moment(5th) ).

Table 5 :

 5 Numerical convergence results for the initial condition (23)-[START_REF] Garres-Díaz | Shallow water moment models for bedload transport problems[END_REF]. The L 1 errors refer to the variables h and u at a final time t = 10 s, measured for different vertical discretizations of M layers and using uniform Cartesian meshes of N x elements.

Table 6 :

 6 Numerical convergence results for the initial condition[START_REF] Gray | A depth-averaged µ(I)-rheology for shallow granular free-surface flows[END_REF]. The L 1 errors refer to the variable u at a final time t = 1 s, measured for different vertical discretizations of N layers and using uniform Cartesian meshes of N x = 1600 elements.

Table 7 :

 7 Wall clock times of MCU (t MCU(M)

Acknowledgements

This research has been partially supported by the Spanish Government and FEDER through the research projects RTI2018-096064-B-C2(1/2) and PID2020-114688RB-I00, the Junta de Andalucía research project P18-RT-3163, the Junta de Andalucia-FEDER-University of Málaga research project UMA18-FEDERJA-16 and the University of Málaga.

Appendix A. Viscous multilayer-moment systems and continuity conditions at the interfaces This appendix is devoted to describe, for the sake of completeness, how the viscous terms are included in the derivation of the models. Although this work focuses mainly in Euler equation, they are relevant for some of the simulations and they allow to properly define the velocities at the interfaces u | L ± α+1/2 . We first detail the general procedure and then describe the viscous system in the particular case of layerwise constant velocities (Multilayer Constant method).

Appendix A.1. Multilayer-moment system

Let us first recall that for any variable f ≡ f (t, x, ξ), we shall denote by

the upper or lower limit at the interface L α+1/2 . When the variable is continuous at the interface, we simply denote it by f α+1/2 . We shall denote by

Let us focus on (2a) in the hydrostatic case, but keeping now the viscous terms:

where (A.1) has been used. Therefore,

where

In a similar way one obtains

Finally, ( 18) is recovered from (A.4) when going to the limit of ν tending to zero to recover Euler equations.

Remark 7. There is a compatibility condition between definitions (A.2) and (A.3) for τ and u, in the sense that the numerical flux we consider to define τ α+1/2 , in view of (A.1), is the flux we recover for the definition of the velocity at the interface u α+1/2 .

For instance, in [START_REF] Fernández-Nieto | A multilayer method for the hydrostatic navier-stokes equations: A particular weak solution[END_REF][START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I)-rheology: application to granular collapse on erodible beds[END_REF][START_REF] Garres-Díaz | Multilayer models for shallow two-phase debris flows with dilatancy effects[END_REF][START_REF] Bonaventura | Multilayer shallow water models with locally variable number of layers and semi-implicit time discretization[END_REF] a central flux (Λ = 0) in considered, whereas an upwind numerical flux (Λ = sign (-ω)) is taken in [START_REF] Audusse | A multilayer Saint-Venant system with mass exchanges for shallow water flows. derivation and numerical validation[END_REF][START_REF] Audusse | Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model: kinetic interpretation and numerical solution[END_REF][START_REF] Guerrero Fernández | A second-order well-balanced finite volume scheme for the multilayer shallow water model with variable density[END_REF][START_REF] Garres-Díaz | Flexible and efficient discretizations of multilayer models with variable density[END_REF], mainly for stability reasons, where the sign of -ω matches with the sign of the vertical velocity (see (2b)).

Appendix A.2. Viscous Multilayer Constant method

We consider now the viscous Multilayer Constant (MC) case, where we recall that the velocity is layerwise constant. With the goal of keeping the friction force and vertical viscosity effects, we consider the leading order approximation σ = -

which is easily justified from an asymptotic analysis. Then, from (A.4) we get

for a constant kinematic viscosity ν. By simplicity, consider here Λ α+1/2 = 0 in (A.2), and therefore

In particular, we need to approximate ∂ ξ u α+1/2

. Concretely, we choose

.

Therefore, the viscous version of system [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I)-rheology: application to granular collapse on erodible beds[END_REF], accounting for boundary conditions (2c), reads

and the boundary conditions (2c) for the stress tensor K M +1/2 = 0, K 1/2 = |u 1 | u 1 . Finally, we have

Appendix B. Viscous shallow water moment model

Including the contribution of the viscous term ∂ ξ τ xz in moment model ( 20) is a simple task. Proceeding as in [START_REF] Garres-Díaz | Shallow water moment models for bedload transport problems[END_REF], we compute the following integral for i = 0, . . . , N (see (2a)) taking into account the boundary conditions (2c):

where


Arnaud Cariou
email: arnaud.cariou.ac@gmail.com

Afef Kedous-Lebouc

Christian Chillet

Laurent Albert

Sergi Riba

New investigation in Loss Surface model to consider dynamic magnetization's history

Keywords: hysteresis modeling, Loss Surface model, dynamic magnetic history, experimental validation

This paper presents a new calculation of the dynamic part of the Loss Surface (LS) model to consider the magnetization history. It is to be combined with a static model of hysteresis to estimate iron loss. The analysis of the behavior of various hysteresis loops highlighted the history dependency of the dynamic part of the LS model. An innovative formulation is then proposed to consider the dynamic effects inertia and to construct correctly the practical hysteresis for any applied waveform. The improved model is validated through experimental measurements.

Introduction

The continuously increase of electric power in airplane urges to improve electrical machines design and, among others, iron loss calculations. Hysteresis models are characterization costly but also more precise. The LS model has been improved recently [1] and can be used properly for iron loss estimation. This paper aims to improve further its accuracy by considering dynamic magnetization history.

Dynamic effects and magnetization history

The LS model is based on the magnetic field separation into a static part 𝐻 𝑠𝑡𝑎𝑡 , depending on the flux density 𝐵 and the magnetic history 𝐻𝑖𝑠𝑡(𝐵), and a dynamic part 𝐻 𝑑𝑦𝑛 related to 𝐵 and its time derivative 𝐷 𝑡 𝐵:

𝐻 = 𝐻 𝑠𝑡𝑎𝑡 (𝐵, 𝐻𝑖𝑠𝑡(𝐵)) + 𝐻 𝑑𝑦𝑛 (𝐵, 𝐷 𝑡 𝐵) (1)
𝐻 𝑑𝑦𝑛 is determined from a single surface 𝐻 𝑠𝑢𝑟𝑓 measured under triangular 𝐵 of 𝐵 𝑚𝑎𝑥 amplitude and variable frequency. In [2], a little effect between sine and triangular waveforms is observed. Yet, these damping effects are well known and can modify the material behavior as illustrated in [3]. Here, a set of triangular and trapezoidal B(t) has been applied to compare 𝐻 𝑑𝑦𝑛 (𝐵, 𝐷 𝑡 𝐵). A large difference between 𝐻 𝑑𝑦𝑛 (𝑡𝑟𝑖) and 𝐻 𝑑𝑦𝑛 (𝑡𝑟𝑎𝑝) is observed, especially at reversal points.

New time-dependent dynamic hysteresis

A new corrected dynamic surface (𝐻 𝑠𝑢𝑟𝑓,𝑐𝑜𝑟 ) is proposed:

𝐻 𝑠𝑢𝑟𝑓,𝑐𝑜𝑟 (𝑡) = 𝜏 • 𝑑𝐻 𝑠𝑢𝑟𝑓 (𝑡) 𝑑𝑡 + 𝐻 𝑠𝑢𝑟𝑓 (𝑡) (2)
With 𝜏 = 𝜎•𝜇 𝑠 •𝑒 2
12 , 𝜎: the electric conductivity, 𝜇 𝑠 : the static magnetic permeability, 𝑒: the material thickness. Fig. 1 compares 𝐵(𝐻 𝑑𝑦𝑛 ) cycles before and after corrections. It shows that 𝐻 𝑑𝑦𝑛,𝑐𝑜𝑟𝑟 (𝑡𝑟𝑖) and 𝐻 𝑑𝑦𝑛,𝑐𝑜𝑟𝑟 (𝑡𝑟𝑎𝑝) overlap and are very similar for both 0.6 and 1.2 T. This result is also confirmed by a theorical approach involving classical eddy currents.

So, for any B(t) waveform, the time-dependent 𝐻 𝑑𝑦𝑛 is calculated from 𝐻 𝑠𝑢𝑟𝑓,𝑐𝑜𝑟 , with 𝜏′ being correlated to the static permeability 𝜇 𝑠 ′ linked to B(t) signal:

𝐻 𝑑𝑦𝑛 (𝑡) = 𝐻 𝑠𝑢𝑟𝑓,𝑐𝑜𝑟 (𝑡) • 𝛿𝑡 𝜏′ + 𝐻 𝑑𝑦𝑛 (𝑡 -1) 1 + 𝛿𝑡 𝜏′ (3)
This new formulation shows an overall improvement of the hysteresis estimation compared to the previous one. The final paper will illustrate the relevance of this correction for several excitation and its impact on the loss computation while combined with a static Play model [4].

Figure 1 :

 1 Figure 1 : Measured and corrected 𝐻 𝑑𝑦𝑛 loops at 0.6 and 1.2 T for triangular and trapezoidal B, for 𝐷 𝑡 𝐵 = 980 𝑇/𝑠


Anh-Tuan Vo

Marylin Fassenet

Afef Kedous-Lebouc

New Formulation of Dynamic LS Model for Accurate Iron Loss Modeling at High B and dB/dt

Keywords: Loss Surface model, LS model, hysteresis modeling, iron loss, high dynamic effects, high B and dB/dt

In this study, the LS model is reformulated to reliably extend the dynamic surface at extreme conditions of (B, dB/dt). A comparison between measurements and FE simulations of a highly dynamic NiFe alloy will be presented, showing that at very high values of B or dB/dt, the contribution of eddy currents to the dynamic field Hdyn is dominant and Hdyn becomes a linear function of dB/dt regardless of the value of B. An empirical function is then developed to comply with the above behavioral basis while still fitting the experimental dynamic surfaces. For each material, this function requires 14 parameters identified by a fitting program. Its robustness and accuracy are analyzed on a large set of materials. In all cases, the average deviation of iron losses is lower than 5%.

Introduction

The LS model has been used in a posteriori mode of finite element simulations for the calculation of magnetic losses in electrical machines since 20 years. Several studies have shown its ability to predict satisfactorily the experimental losses [START_REF] Frias | Rotor and stator shape optimization of a synchronous machine to reduce iron losses and acoustic noise[END_REF]. However, the particular interest in the search for high-efficiency machines, in the integration of soft magnetic materials under extreme conditions, motivates to revise this model, analyze its limits, and extend its validity range while further improving its accuracy.

Methodology

Experimental measurements and finite element simulations on the iron losses of various electrical steels are coupled to study the pattern of the dynamic surfaces. Once the general pattern is found, potential functions are tested and the most relevant one is analyzed in this paper.

Principle results

Slicing the dynamic surfaces at constant B or dB/dt, one obtain the curves Hdyn(dB/dt) or Hdyn(B), with: Hdyn(dB/dt) curves featuring an increasing function of dB/dt and Hdyn(B) curves having a piecewise form that saturates or approaches asymptotes when B is around ±Jsat. These properties of the surfaces Hdyn(B, dB/dt) can be modeled by a sigmoidal function comprising of 14 parameters. It has been tested with various electrical steels and has shown great accuracy. In all cases, the hysteresis cycles with complicated B(t) waveform (Figure 1b) or sinusoidal B(t) are reconstructed with an average iron loss deviation of less than 5% (Table 1). In addition, the unique formula can extrapolate very well the dynamic surface at high frequency and high induction (Figure 1a).

Conclusions

The dynamic LS model is accurately formulated by a modified sigmoidal function. Thanks to its extrapolation capacity, the identification of the LS model is facilitated in terms of measurements and computer tools. As a rule of thumb, the model identification requires five hysteresis cycles measured under triangle B(t) with the same induction amplitude of around 80% of saturation and at five frequencies varying from the quasi-static to the highest measureable.

Figure 1 (

 1 Figure 1 (a) Extrapolation capability of the new formula at very high dB/dt; (b) Reconstruction of a complicated hysteresis cycle.

Table 1

 1 Average iron loss deviation in the case of sinusoidal B(t).

	Material	B (T)	f (Hz)	Mean deviation (%)
	M330-35A	0.1-1.6	50-4000	3.99
	NO20	0.1-1.5	50-4000	3.04
	FeNi Supra 50 0.1-1.5	50-4000	4.65
	FeCo AFK1	0.1-2.1	50-4000	4.80

Hdyn (A/m) B (T)

Model

FE simulations

FeNi NO20

(a) (b)


Amaury Carrard-Draghi

Aktham Asfour
email: aktham.asfour@g2elab.grenoble-inp.fr

Nicolas Galopin

Sébastien Flury

Jean-Paul Yonnet

A setup for the measurement of the GMI response of amorphous wires under tensile stress

Keywords: Giant Magneto-Impedance (GMI), tensile stress, measurement setup

This paper presents the design and development of an experimental setup that allows measurement of the Giant Magneto-Impedance (GMI) response of Co-rich amorphous micro-wires under tensile stress. A full description of the setup and the first performed measurements are given.

Introduction

Giant Magneto-Impedance (GMI) is a significant change of the AC impedance of some soft magnetic materials (amorphous wires…) when they are supplied by a high frequency current and submitted to an external magnetic field. This effect is used to develop magnetic sensors based on the impedance measurement of the sensitive element. Mechanical stress like tensile, if applied to the sensitive element, considerably modifies its GMI response [1]. This is an undesirable effect for a field sensor. Nevertheless, it seems relevant to make good use of this effect and to evaluate its potential to develop a new category of strains sensors or strain gauges [2]. The well understanding of the involved mechanisms and the measurement of these effects is an issue of great importance. This paper presents the design and development of a setup that allows measurement of the GMI response under tensile stress.

Measurement Setup

The developed setup is shown in Figure 1. It allows measuring the impedance of a GMI wire, under tensile stress, as a function of the magnetic field. From electrical point of view, a signal generator and a voltage to-current, (resistor Rg) supply the high frequency current iac to the wire. The voltage vac across the wire is amplitudemodulated by the external field H produced by Helmholtz coils. This voltage is demodulated using a home-made demodulator. The output of this demodulator, which is proportional to the modulus of the impedance |Z|, is digitized by a data acquisition board for processing and visualization purposes on a PC. A straightforward calibration allows obtaining the impedance from the acquired voltage. The applied magnetic field can be automatically modified by controlling the DC power supply of the Helmholtz coils. The impedance curve could be displayed on a Graphical User Interface (GUI).

The mechanical arrangement allows applying a tensile stress on the wire in vertical position. The extremities of this wire were soldered to supports for mechanical fastening. The free extremity is also fastened to a holder ("basket"). The application of a tensile stress is ensured by placing a calibrated mass in this holder.

Results

First results were obtained using an amorphous wire (CoFeSiB) with a 100 µm and 80 mm long (Figure 2). These curves show a dependence of the GMI response on the applied tensile stress. The effect is particularly obvious around the field of the maximum impedance (inset of Figure 2). This field is directly related to the anisotropy field of the magnetic material. Further improvements of the setup should include arrangements for the measurement under torsion and bending stresses of the wire.

Figure 1 :

 1 Figure 1: Schematic of the experimental setup.

Figure 2 :

 2 Figure 2: GMI responses under tensile stresses.