Edwige Godlewski
email: edwige.godlewski@sorbonne-universite.fr

Martin Parisot
email: martin.parisot@inria.fr

Jacques Sainte-Marie
email: jacques.sainte-marie@inria.fr

Fabien Wahl
email: fabien.wahl@inria.fr

Congested shallow water model: on floating body

Keywords: 2010 Mathematics Subject Classification. 35Q35, 70E15, 74F10, 76B07, 76M12 shallow water equations, wave-body interaction, congested model, coupling, entropy satisfying scheme

published or not. The documents may come

Introduction

We are interested in the modeling of a floating body in the vertical plane. In a previous work [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF], we considered a shallow water flow with an impermeable surface above it, namely the roof, which constrains the water surface. The roof is motionless or is given by an explicit time and space dependent function independent of the flow. In this work, we focus on a floating structure leading to the coupling of the equations on the fluid and on the motion of the structure. Applications are the drift of icebergs or flood debris, or the improvement of sustainable marine energy converters such as buoys.

Two main approaches are proposed in the literature to treat fluid-solid interaction. With the exception of work done in antiquity, in particular "On Floating Bodies" by Archimedes, the first one dates back to Fritz John who proposed a mathematical formulation of the problem [START_REF] John | On the motion of floating bodies[END_REF]. The fluid is described by a velocity potential with a linear model for the free surface evolution. The motion of the solid is supposed to be of small amplitude so that the interface between the water and the solid is constant in time. Finally the surface pressure is obtained using the linearized Bernoulli equation. Although the model is quite simplified, linear potential flow theory is still used in industrial context since it is not costly in CPU time, see [START_REF] Matt | A review of numerical modelling of wave energy converter arrays[END_REF]. However nonlinear terms play an important part for wave interaction. For instance the non-linearity should be taken into account in simulations of wave energy converter farms. Improvements have been made to include an interface depending on time [START_REF] Kashiwagi | Non-linear simulations of wave-induced motions of a floating body by means of the mixed Eulerian-Lagrangian method[END_REF] and nonlinear effects based on the boundary element methods [START_REF] Harris | Simulation of wave forces on a gravity based foundation by a BEM based on fully nonlinear potential flow[END_REF]. This progress needs higher CPU times and does not yet allow wave breaking.

The second approach involves Navier-Stokes CFD computations. This approach is considerably used in blood flow, for instance in aortic flow [START_REF] Quarteroni | Computational vascular fluid dynamics: problems, models and methods[END_REF][START_REF] Chabannes | High-order fluid-structure interaction in 2d and 3d application to blood flow in arteries[END_REF]. Some other works also use Navier-Stokes CFD computations for the simulation of the flow around a yacht [START_REF] Parolini | Mathematical models and numerical simulations for the America's cup[END_REF], fish [START_REF] Bergmann | Bioinspired swimming simulations[END_REF] or wave energy converters [START_REF] Agamloh | Application of fluid-structure interaction simulation of an ocean wave energy extraction device[END_REF][START_REF] Palm | CFD simulation of a moored floating wave energy converter[END_REF][START_REF] Yu | Reynolds-averaged Navier-Stokes simulation of the heave performance of a two-body floating-point absorber wave energy system[END_REF]. In these works the fluid is expressed in an Eulerian frame while the solid is expressed in a Lagrangian frame. The difficulty comes from the mesh. A moving grid with front tracking methods or a fixed mesh with fictitious fluid domain can be considered. The first one perfectly catches the interface but needs a moving mesh while the second is less accurate on the position of the interface but the mesh is fixed. However this strategy implies high CPU times and does not make these methods suitable for engineering optimization. To take advantage of the benefits from each approach, a coupling has been proposed [START_REF] Wu | The coupled finite element and boundary element analysis of nonlinear interactions between waves and bodies[END_REF]. Since the viscosity is not considered in potential flow theory, it can hardly be included in the coupling.

Recently the nonlinear floating body problem was reformulated in the framework of vertical-integrated models [START_REF] Lannes | On the dynamics of floating structures[END_REF]. A vertical movement of the body is considered so that the interface between the congested and the free surface domain is fixed in time and continuity of the unknowns is assumed at the interface. The surface pressure is then obtained by an elliptic equation. An analysis of this shallow water type model in the two dimensional case with radial symmetry is done in [START_REF] Bocchi | Floating structures in shallow water: local well-posedness in the axisymmetric case[END_REF]. In [START_REF] Bresch | Waves interacting with a partially immersed obstacle in the Boussinesq regime[END_REF], a phenomena of dispersive boundary layer is highlighted. In a more general framework the well-posedness of the model in the one dimensional case has been proven in [START_REF] Iguchi | Hyperbolic free boundary problems and applications to wavestructure interactions[END_REF], relaxing the assumption of vertical lateral walls. The description of the interface position is dealt with. A numerical resolution in one dimension is done in [START_REF] Bosi | A spectral/hp element depth-integrated model for nonlinear wave-body interaction[END_REF]. The free surface and the congested subdomain are solved using an appropriate numerical strategy, i.e. hyperbolic for the first and elliptic for the second. At the interfaces, the continuity of the pressure and the discharge is assumed. Unfortunately, the transmission condition is not clear physically, since a discontinuous solution is possible. In addition, when the position of the interface is also an unknown, the coupling strategy is tricky to handle, and not robust. In [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] we have proposed a uniform numerical resolution which eliminates the description of the interface and the transition conditions between the congested and the free surface domain. The surface pressure is seen as a Lagrange multiplier associated to a constraint. These models [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF][START_REF] Lannes | On the dynamics of floating structures[END_REF] are based on depth-averaging and thus they cannot represent wave breaking but nonlinear terms are accounted for. Simpler than the Navier-Stokes equations but still physically relevant, their resolution is less CPU consuming since the computational domain does not depend on time. For optimization processes this seems interesting, especially when large domains are considered. Eventually the interaction of energy converters with the bottom as in [START_REF] Benyo | Wave-structure interaction for long wave models with a freely moving bottom[END_REF][START_REF] Benyo | Numerical analysis of the weakly nonlinear Boussinesq system with a freely moving body on the bottom[END_REF] can also be considered.

In the present work, we give a formulation of the problem taking into account a freely floating object in the vertical plane, i.e. two translations and one rotation. A choice of discretization for the different operators is proposed in order to ensure a strong coupling between the two systems. We use a Newmark scheme for the solid and adapt the numerical strategy given in [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] for a 'fixed roof' to the coupling. Furthermore our method allows to write the energy of the coupled fluid-solid system. The energy transfer between the fluid and the solid is indeed a challenging problem and of major interest in energy production using buoys. From a mathematical point of view, the energy balance, acting as an entropy, is an argument for the existence of long time solutions.

Note that the modeling of submerged objects is not allowed by the choice of modeling for the fluid since the vertical averaged model considers only one water height. Thus in the following, we will restrict ourselves to a certain type of configuration for the floating objects. More details are given in Hypothesis 2.

In Section 2 the formulation at the continuous level is given while Section 3 is dedicated to the numerical strategy. In both cases, the fluid and solid dynamics are considered separately, then the coupling is handled. At the discrete level, first a naive non-entropy satisfying approach is presented followed by an entropy correction. We also give details on the practical 1D implementation to treat general shapes of the floating body. Finally in Section 4 we show simulations to validate our approach.

Mathematical modeling

In the current section, a description of the governing equations of the floating body problem using a vertical-integrated model is given. The physical context is illustrated in Figure 1. 2.1. Fluid dynamic. In [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] a shallow water type model with an additional congestion constraint modeling a 'roof' is proposed. Let us briefly introduce this model which is derived from the Navier-Stokes equations. Let us consider the one dimensional domain Ω x ⊂ R where x stands for the horizontal coordinate in the frame of the observer and t ∈ R + stands for the time variable. In practice the domain Ω x will be bounded but for simplicity the boundary conditions will not be detailed. A flow contained between two surfaces (in the 1D framework, these are actually two curves) respectively called roof and bottom is considered. The two surfaces can be parametrized by two given mono-valued functions R (x, t), the roof, and B (x, t), the bottom, which satisfy B (x, t) ≤ R (x, t). The surfaces must be regular enough where they are reached by the fluid, i.e. everywhere for the bottom and at least at the inferior surface of the floating body, see Hypothesis 1 below. The opening between the roof and the bottom is defined by H (x, t) = R (x, t) -B (x, t). The unknowns of the model are the water depth h (x, t), the vertical-averaged horizontal velocity u (x, t) and the surface pressure p (x, t) which satisfy (1)

   ∂ t h + ∂ x (hu) = 0 ∂ t (hu) + ∂ x hu 2 = -h∂ x g (h + B) + p ρ
with g the gravitational constant. The fluid is assumed to be homogeneous, then the water density ρ is constant in time and space and given. In addition, the congestion constraint

(2) min Hh, p -P = 0 has to be satisfied where P (x, t) denotes the given atmospheric pressure. In the following, the atmospheric pressure is taken constant in time and space and for simplicity equal to zero, i.e. P = 0. Finally the system (1)-( 2) is completed with the initial conditions h (x, 0) = h 0 (x) and u (x, 0) = u 0 (x). We recall the following energy balance law, the proof of which can be found in [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF]Proposition 2.3].

∂ t E + ∂ x G = -p∂ t R + (ρgh + p) ∂ t B
with the mechanical energy and the energy flux respectively given by

E = ρ 2 hu 2 + ρgh B + h 2 and G = ρ 2 u 2 + ρg (h + B) + p hu.
By smooth enough solution, we assume enough regularity on the solution to perform each step of the computation, in particular we assume that h > 0 and h ∈ C 0 (Ω x ), u ∈ C 0 (Ω x ) and p ∈ C 0 (Ω x ) and as usual the above conservation law hols almost everywhere in practice. Unlike the work done in [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF], if B (x, t) is still a given function, this is no longer the case for R (x, t). The interaction between the fluid and the roof is considered, i.e. R (x, t) is now an unknown of the problem and its evolution has to be determined, see Section 2.2.

If a smooth solution is unique, it is well known that for hyperbolic equations such as the shallow water equations, there may be several weak solutions, and one generally adds some entropy criterion. In the case of the shallow water model, uniqueness can be ensured by selecting the solution which ensures mechanical energy dissipation, acting as a mathematical entropy, see for instance [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF].

2.2. Solid dynamic. Let us describe the planar motion of a body subjected to the gravity field and a surface pressure force p. The reader can refer to [START_REF] Knudsen | Elements of Newtonian Mechanics[END_REF] for more details about solid dynamics. In the vertical plane (e x , e z ), the scalar product between two vectors a = (a 1 , a 2 ) and b

= (b 1 , b 2 ) is denoted a • b = a 1 b 1 + a 2 b 2 and the cross product a × b = a 1 b 2 -a 2 b 1 .
We describe a general planar movement of a homogeneous solid body by three degrees of freedom, denoted χ, ζ and θ, see Figure 1. The mass of the solid is denoted by M and its center of mass by G = (χ, ζ), where χ ∈ Ω x denotes the horizontal and ζ ∈ R the vertical component. The variable θ stands for the angle between the unit vector in the vertical direction e z and the vector GX where X is a point in the solid different from the center of mass. The moment of inertia of the solid around an axis passing through the center of mass is denoted J . We consider a solid subjected to gravity and a surface pressure of p. Additional horizontal (resp. vertical) forces will be denoted by F χ (resp. F ζ ) and their torque around the center of mass is denoted by T θ .

Hereafter, we give the mathematical description of the buoy. Note that in practice, it needs a numerical implementation given in Section 3.4.1. For every θ ∈ [-π, π[, the geometry of the solid body in the frame centered on G is the set of coordinates (x, z) satisfying S (x, z, θ) = 0, with S (x, z, θ) a given function. In the context of a floating body, from now on called the buoy, the surface pressure only acts on its inferior surface. This inferior surface elevation is given by

z = S (x, θ) = min { z ∈ R | S (x, z, θ) = 0}
that leads to the following expression of R, the inferior surface elevation in the frame of the observer

(3) R (x, χ, ζ, θ) = S (x -χ, θ) + ζ.
Note that R (x, χ, ζ, θ) is not necessarily defined for all (x, χ, ζ, θ) Let us give some examples of the previous mathematical objects in particular cases:

∈ Ω x × R × R × [-π, π[.
• In the case of a discal solid with radius r, since the geometry is independent of the rotation, it reads

S dis (x, z, θ) = x 2 + z 2 -r 2 and S dis (x, θ) = -r 2 -x 2 .
The domain of definition is given by

Ω dis R (χ, θ) = {x ∈ Ω x | χ -r < x < χ + r} .
• In the case of an elliptical solid with semi-minor and semi-major axis respectively equal to a and b, it reads

S ell (x, z, θ) = (x cos θ + z sin θ) 2 a 2 + (x sin θ -z cos θ) 2 b 2 -1
and S ell (x, θ) is the smaller root of the second order polynomial P x,θ (z) = S ell (x, z, θ). The domain of definition reads

Ω ell R (χ, θ) = x ∈ Ω x | χ + min ϑ∈[-π,π[ f θ (ϑ) < x < χ + max ϑ∈[-π,π[ f θ (ϑ)
with f θ (ϑ) = a cos (θ) cos (ϑ) + b sin (θ) sin (ϑ).

Since there is not always an analytical expression for R, we will indicate in Section 3.4 how to handle the problem from a practical point of view.

We are now able to write more precisely the assumption on the surface pressure.

Hypothesis 1. The pressure varies only on the inferior surface of the solid, i.e. there exists a function p (x, t) such that for any 0 ≤ t ≤ T the surface pressure reads p (x, z, t) = p (x, t) if x ∈ ΩR (χ (t) , θ (t)) and z = S (x, θ (t)) 0 else.

In addition, on the inferior surface ΩR (χ (t) , θ (t)), the roof R is regular. In particular, we assume that no vertical part of the floating object is in contact with water.

Note that the pressure can still vanish on the inferior surface of the solid when p (x, t) = 0. Proposition 1. Assume that Hypothesis 1 is satisfied. Then the movement of the solid body is described by the system

(4) C 1 Λ = F Λ, Λ, p, t with Λ =   χ ζ θ   , C 1 =   M 0 0 0 M 0 0 0 J   and F = F E p, t + F P where F E p, t =         - Ω R p ∂ x R dx + F χ (t) Ω R p dx + F ζ (t) Ω R p ∂ θ R dx + T θ (t)        
and

F P =   0 -M g 0   completed with the given initial conditions χ (0) = χ 0 , χ (0) = χ0 , ζ (0) = ζ 0 , ζ (0) = ζ0 , θ (0) = θ 0 and θ (0) = θ0 .
Proof. A general planar motion of a rigid body can be separated into a translational motion of a point in the body and a rotational motion around an axis through that point. It is generally convenient to choose this point as the center of mass G, see [START_REF] Knudsen | Elements of Newtonian Mechanics[END_REF]. The translational motion can be described by Newton's second law of motion.

Taking into account the forces F χ , F ζ , the weight together with the pressure applied on the solid from the water gives the following equations ( 5)

       M χ = Ω R p n • e x 1 + (∂ x R) 2 dx + F χ M ζ = -M g + Ω R p n • e z 1 + (∂ x R) 2 dx + F ζ
with n the interior unit normal to the solid surface. The rotational movement is described by the angular momentum theorem, i.e.

(

) ∂ t L G = Ω R p (GX × n) 1 + (∂ x R) 2 dx + T θ 6
where L G denotes the angular momentum around G, X a point at the surface with coordinates (X x , R (X x , χ, ζ, θ)). By definition of the kinetic momentum in the plane,

we have ∂ t L G = J θ. Replacing the inward normal n = 1 1 + (∂ x R) 2 -∂ x R 1 in
(5) leads to the first two equations of ( 4) and for (6) it yields ( 7)

J θ = - Ω R p (GX • e z ∂ x R + GX • e x ) dx + T θ .
Let us denote by X x = X xχ the horizontal coordinate in the frame centered on G. In this frame, the movement is a simple rotation of angle θ and as a consequence the coordinates X x , S X x , θ of the point X verify

X x S X x , θ = M θ-θ 0 X 0 x S X 0 x , θ 0 with X 0 x , S X 0
x , θ 0 the coordinates of the point X with an angle θ 0 and the matrix of rotation ( 8)

M ϑ = cos ϑ sin ϑ -sin ϑ cos ϑ .
The previous relation implies

(9) S = ∂ θ X x .
Moreover the trajectory of the points of the solid are concentric circles, i.e. ∂ θ GX = 0. Thus

∂ θ X 2 x + S X x , θ 2 = 0 and it follows X x ∂ θ X x + S∂ θ X x ∂ x S + S∂ θ S = 0.
Using now the relation ( 9), we get (10)

X x + S∂ x S + ∂ θ S = 0.
Finally introducing this relation in [START_REF] Bocchi | Floating structures in shallow water: local well-posedness in the axisymmetric case[END_REF] and noticing that GX • e z = S and GX • e x = X x gives the last equation in (4).

The system (4) requires to know the pressure p on the body. This pressure will be estimated from the fluid system (1), see Section 2.3. The derivatives of the inferior surface elevation ∂ x R and ∂ θ R are not well-defined at the boundary of its domain of definition Ω R . However, thanks to Hypothesis 1, the support of the pressure p does not reach the boundary of the domain where the inferior surface is defined Ω R . Assuming that the inferior surface elevation R is regular enough on the support of the pressure p, we ensure that the integrals in (4) are well-defined.

The solid system admits the following energy balance law.

Lemma 2.2. The solution of (4) satisfies the following energy law

∂ t E = Ω R p dx + F ζ ζ - Ω R p∂ x R dx -F χ χ + Ω R p∂ θ R dx + T θ θ
with the mechanical energy of the buoy

E χ, ζ, θ, ζ = M 2 ζ2 + χ2 + J 2 θ2 + M gζ.
Proof. Multiplying the first equation of ( 4) by χ, the second equation by ζ, the third equation by θ and summing gives the result.

2.3. Full system. The two models (1) and ( 4) have to be coupled. More precisely, the roof R is an external unknown of the model (1) but can be determined from the position of the buoy computed in the model ( 4). Similarly the pressure p is an external unknown of the model ( 4) but is estimated by the model [START_REF] Agamloh | Application of fluid-structure interaction simulation of an ocean wave energy extraction device[END_REF].

According to Newton's third law, i.e. the action-reaction principle, the action of the buoy on the water is opposite to the reaction of the water on the buoy. In the mathematical description introduced previously, it implies that for any 0 ≤ t ≤ T and x ∈ Ω R (χ (t) , θ (t)) we have [START_REF] Bresch | Waves interacting with a partially immersed obstacle in the Boussinesq regime[END_REF] p (x, t) = p (x, t) .

Let us design by supp (f ) the support of the function f . Hypothesis 1 and ( 11) together leads to the following assumption.

Hypothesis 2. We assume that the water does not reach the extremities of the buoy, i.e. for any

0 ≤ t ≤ T supp (p (•, t)) ⊂ ΩR (χ (t) , θ (t)) .
From now on, we will only use the unknown p and we assume that Hypothesis 2 is satisfied.

In the part of the domain where there is no buoy, the roof is defined high enough not to touch the fluid surface and otherwise the roof is given by R. More precisely for Π ∈ R large enough (in practice we set Π = 10 40 ) we define

(12) R (x, t) = R (x, χ (t) , ζ (t) , θ (t)) if x ∈ Ω R (χ (t) , θ (t)) Π
elsewhere.

An energy law for the coupled fluid-solid system (1)-( 2)-( 4)-( 11)-( 12) is obtained.

Proposition 2. Considering a no-flux boundary on Ω x , any smooth solution of (1)-( 2)-( 4)-( 11)-( 12) satisfies the following energy balance law

∂ t E + E = Ωx (ρgh + p) ∂ t B dx + F ζ ζ + F χ χ + T θ θ with E = Ωx E dx.
Proof. It follows from ( 3) and ( 12) that ( 13)

∂ t R (x, t) = χ∂ χ R (x, χ, ζ, θ) + ζ∂ ζ R (x, χ, ζ, θ) + θ∂ θ R (x, χ, ζ, θ) = -χ∂ x R (x, χ, ζ, θ) + ζ + θ∂ θ R (x, χ, ζ, θ) .
Now integrating the energy balance law from Lemma 2.1 over Ω x and summing with the result from Lemma 2.2 concludes the proof.

Numerical approximation

This section is devoted to the numerical approximation of a solution of ( 1)-( 2)-( 4)-( 11)- [START_REF] Bristeau | An energy-consistent depth-averaged Euler system: derivation and properties[END_REF]. We first describe separately the numerical methods chosen for each model before considering their coupling.

We consider a 1D mesh built from an increasing sequence x k+ 1 2 ∈ R, where k belongs to a set T ⊂ N of indices. A control volume is an interval

V k = x k-1 2 , x k+ 1 2
and we denote by k = x k+ 1 2x k-1 2 its length and by

x k = x k-1 2 +x k+ 1 2 2
its center. In addition, the time is discretized by t n+1 = t n + δ n+1 t , with t 0 = 0 and the time step δ n+1 t will be adapted to the flow through a CFL condition see [START_REF] Guerber | A fully nonlinear implicit model for wave interactions with submerged structures in forced or free motion[END_REF] below. The mean value of the bottom level in the control volume V k , at time t n , is denoted by B n k . We also introduce the following notation

∂ n+1 t ψ = ψ n+1 -ψ n δ n+1 t and (ψ) ± = |ψ| ± ψ 2 .
3.1. Discretization for the fluid dynamic. Let us firstly recall the numerical strategy proposed in [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] for the congested shallow water model ( 1)-( 2). A pseudocompressibility method is used, with λ the relaxation parameter, to take into account the congestion constraint. Namely the following relaxed model is considered ( 14)

     ∂ t h λ + ∂ x (h λ u λ ) = 0 ∂ t (h λ u λ ) + ∂ x h λ u 2 λ = -gh λ ∂ x h λ + B + h λ -H + λ 2
where we recall that the opening is defined by

H = R -B.
The relaxed model ( 14) satisfies the energy conservation law written in Lemma 2.1, except that the relaxed pressure is given by

p λ = ρg (hλ-H) + λ 2
and the relaxed mechanical energy by

E λ = ρ 2 h λ u 2 λ + ρgh λ B + h λ 2 + g 2λ 2 (h λ -R λ + B) 2 + .
The flux of energy G and the source term (right hand-side) are not affected by the relaxation. In [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF], we show that the solution of ( 14) tends formally to the solution of (1) when λ goes to zero. In particular, (h -H) + = O λ 2 which ensures that the energy E λ tends to E and the constraint (2) is satisfied when λ goes to zero. Although the relaxed model ( 14) is hyperbolic, the eigenvalues go to infinity when λ goes to zero, that is why a low-Froude scheme is used. The CPR scheme [START_REF] Parisot | Centered-Potential Regularization for the advection upstream splitting method[END_REF] was chosen for its strong robustness at low-Froude regime. Following [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF], one particularity of the proposed scheme is that the numerical unknowns are φ n k and u n k which are respectively the approximation of the average of the potential defined by ( 15)

φ λ = g h λ + B + h λ -H + λ 2
and the velocity u λ in the control volume V k at time t n . The potential can be used as a parametrization of the water depth and the surface pressure, thus we set ( 16)

H (φ, B, R) =        φ g -B, if φ ≤ gR R -B + λ 2 φ g -B 1 + λ 2 , else and P (φ, B, R) = ρ (φ -g (H (φ, B, R) + B)) .
For readability, we set (17)

h n k = H (φ n k , B n k , R n k ) and p n k = P (φ n k , B n k , R n k ) . where R n
k is an approximation of the mean value of the roof level in the control volume V k at time t n . Note that in the present work, R n k depends on the discrete values of χ, ζ and θ in a similar way to the continuous case [START_REF] Audusse | A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation[END_REF], see more detail in Section 3.2. The conservation of the water volume reads [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] H

φ n+1 k , B n+1 k , R n+1 k = h n+1 k = h n k - δ n+1 t k F n+1 k+ 1 2 -F n+1 k-1 2 with F n+1 k+ 1 2
an approximation of the mean mass flux hu through the interface at

x k+ 1 2 between the times t n and t n+1 . More precisely, the numerical flux reads

F n+1 k+ 1 2 = F             φ n+1 k B n+1 k R n+1 k u n k k       ,       φ n+1 k+1 B n+1 k+1 R n+1 k+1 u n k+1 k+1       , δ n+1 t       with F             φ l B l R l u l l       ,       φ r B r R r u r r       , δ t       = H (φ l , B l , R l ) u l + H (φ r , B r , R r ) u r 2 -γ δ t 2 H (φ l , B l , R l ) l + H (φ r , B r , R r ) r φ r -φ l 2 .
The given regularization parameter γ controls the numerical diffusion. More precisely, the larger γ is, the more important the numerical diffusion is, but γ has to be larger than 1 to make the computation stable, see Lemma 3.1. The discrete conservation of the water volume ( 18) is implicit and non-linear. A fixed point is used to compute the potential φ n+1 k from the implicit and non-linear discrete conservation of the water volume [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF], see Section 3.4.2 for details on the fixed-point. It is clear from [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] that the volume of water is exactly preserved in the domain of computation.

Once the potential φ n+1 k is estimated, the water depth h n+1 k and the pressure p n+1 k are computed using [START_REF] Fernández-Nieto | A hierarchy of dispersive layer-averaged approximations of Euler equations for free surface flows[END_REF] and the velocity is computed explicitly by ( 19)

h n+1 k u n+1 k = h n k u n k -δ n+1 t h n+1 k φ n+1 k+1 -φ n+1 k-1 2 k - δ n+1 t k u n k F n+1 k+ 1 2 + -u n k+1 F n+1 k+ 1 2 - - δ n+1 t k u n k F n+1 k-1 2 - -u n k-1 F n+1 k-1 2 + .
The numerical scheme requires boundary conditions which depend on the regime of the flow. In the current work, we do not deal with the theoretical issues of the boundary conditions. More precisely, we assume that the flow in not congested at the boundary, i.e. h n k < R n k -B n k . It follows that classical boundary conditions of the shallow water model can be used, see [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]Section 21.8].

The following energy dissipation law has been proven in [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF]Proposition 2.3].

Lemma 3.1. Let γ ≥ 1 and assume that for any k ∈ T and any 0 ≤ n < N the following CFL condition is satisfied

(20)   u n k + u n k+1 2 + γ 2 φ n+1 k+1 -φ n+1 k 2   δ n+1 t ≤ min h n+1 k , h n+1 k+1 2 h n+1 k + h n+1 k+1 min ( k , k+1 ) .
Then the scheme (18) admits the following energy dissipation law

(21) ∂ n+1 t E k + 1 k G n+1 k+ 1 2 -G n+1 k-1 2 ≤ -p n+1 k ∂ n+1 t R k + ρgh n k + p n+1 k ∂ n+1 t B k
with the discrete mechanical energy

E n k = ρ 1 2 h n k (u n k ) 2 + gh n k B n k + h n k 2 + g 2λ 2 (h n k -R n k + B n k ) 2 +
and the numerical energy flux

G n k+ 1 2 = ρ |u n k | 2 2 F n+1 k+ 1 2 + - u n k+1 2 2 F n+1 k+ 1 2 - + φ n+1 k + φ n+1 k+1 2 F n+1 k+ 1 2 - h n+1 k+1 u n k+1 -h n+1 k u n k 2 φ n+1 k+1 -φ n+1 k 2 +δ n+1 t h n+1 k+1 k+1 - h n+1 k k φ n+1 k+1 -φ n+1 k 2 2 .
3.2. Discretization for the solid dynamic. The current section is devoted to the numerical approximation of a solution of the solid dynamics (4). The approximation of the three degrees of freedom χ, ζ and θ of the solid at time t n are respectively denoted by χ n , ζ n and θ n . Due to its simple implementation, a Newmark scheme is chosen. For (α, β) ∈ [0, 1] 2 , the Newmark scheme writes

(22)          Λ n+1 = Λ n + δ n+1 t Λn + δ n+1 t 2 2 β Λn+1 + (1 -β) Λn Λn+1 = Λn + δ n+1 t α Λn+1 + (1 -α) Λn C 1 Λn+1 = F P + F n+1 E .
In our case, we have Λ = (χ, ζ, θ)

T while C 1 is given in (4). The forcing term is discretized by

F n+1 E = F E p n+1 k , ∂ k x R n+1 , ∂ k θ R n+1 , t n+1 with (23) F E p, ∂ k x R, ∂ k θ R, t =         - k∈T k p k ∂ k x R + F χ (t) k∈T k p k + F ζ (t) k∈T k p k ∂ k θ R + T θ (t)         where ∂ k x R n (resp. ∂ k θ R n
) are discretizations of the derivative of the roof with respect to x (resp. θ) at time t n in the control volume V k . These terms will be discussed below in order to obtain the energy stability for the coupled fluid-solid model. In the following, we take α = β = 1. This choice is motivated by the fact that the surface pressure source term in ( 21) is implicit and it is required to get the discrete energy stability, see Proposition 4 below. With this choice, the Newmark scheme is completely implicit and expected to be first order accurate.

At the discrete level, the Hypothesis 1 and Hypothesis 2 are not required for the stability of the scheme but must be considered for the relevance of the simulation result. It is actually similar to the case of discontinuous bottom in a classical shallow water model.

Let us first state the energy stability for the Newmark scheme. Lemma 3.2. For any 0 ≤ n < N the scheme ( 22)-( 23) satisfies the following energy law

(24) ∂ n+1 t E = F n+1 E • ∂ n+1 t Λ with E n = E χn , ζn , θn , ζ n defined in Lemma 2.2.
Proof. The proof for the general ODE case is detailed in [26, Section 3]. Using a Newmark scheme for the linear equation of motion on a vector of unknowns Λ

C 1 Λn + C 2 Λn + C 3 Λ n = F n
where C 1 , C 2 and C 3 are matrices and F n a given vector, the following energy equation is proven

∂ n+1 t 1 2 Λ • C 1 + δ n+1 t 2 C 2 Λ + 1 2 Λ • C 3 Λ = ∂ n+1 t Λ • F n+1 - 1 2 ∂ n+1 t Λ • C 2 + δ n+1 t C 3 ∂ n+1 t Λ - 1 2 Λn+1 + Λn 2 • C 2 Λn+1 + Λn 2 .
With the definition given in Proposition 1, the following estimate holds

∂ n+1 t 1 2 Λ • C 1 Λ = ∂ n+1 t Λ • F P + F n+1 E = -∂ n+1 t (M gζ) + ∂ n+1 t Λ • F n+1 E .
The time evolution of the mechanical energy reads

∂ n+1 t E = ∂ n+1 t M 2 ζ2 + χ2 + J 2 θ2 + M gζ = ∂ n+1 t 1 2 Λ • C 1 Λ +∂ n+1 t (M gζ) .
Summing the two previous equalities leads to [START_REF] Kashiwagi | Non-linear simulations of wave-induced motions of a floating body by means of the mixed Eulerian-Lagrangian method[END_REF].

3.3.

Coupling strategy between the fluid and the solid. Let us now concentrate on the discrete coupling. We will discuss the form of 23) and analyze the energy of the coupled system. First an estimate of the discrete entropy law for a general discretization of the operators ∂ x R n and ∂ θ R n is given. In a second paragraph an adapted choice for ∂ k x R n and ∂ k θ R n is proposed to ensure the entropy stability of the numerical scheme.

∂ k x R n and ∂ k θ R n in (

3.3.1.

A first entropy estimate. In this section, we consider the following straightforward discretization of

∂ x R n and ∂ θ R n (25) ∂ k xC R n+1 = R n k+1 -R n k-1 2 k and ∂ k θC R n+1 = -(x k -χ n ) -(R n k -ζ n ) ∂ k xC R n+1 . The discrete operator ∂ k xC R n is the centered discretization of ∂ x R and the dis- cretization of ∂ k θC R n comes from (10). More precisely ∂ k θC R n is a discretization of -(X x -χ) -(R -ζ) ∂ x R.
The shift of the time step is used to make the computation easier. A discrete energy inequality for the coupled fluid-solid system can be written.

Proposition 3. Consider a no-flux boundary on Ω x . Under the CFL condition (20) the scheme (18)-( 22)-( 23) with any discretized operators ∂ k x R n and ∂ k θ R n , admits for any 0 ≤ n < N the following energy inequality 21) over the computational domain and adding [START_REF] Kashiwagi | Non-linear simulations of wave-induced motions of a floating body by means of the mixed Eulerian-Lagrangian method[END_REF] gives the result.

∂ n+1 t E + E ≤ k∈T k ρgh n k + p n+1 k ∂ n+1 t B k +F n+1 ζ ∂ n+1 t ζ -F n+1 χ ∂ n+1 t χ + T n+1 θ ∂ n+1 t θ + k∈T k p n+1 k -∂ n+1 t R k -∂ n+1 t χ∂ k x R n+1 + ∂ n+1 t ζ + ∂ n+1 t θ∂ k θ R n+1 with E n = k∈T k E n k . Proof. Summing (
The estimate of Proposition 3 is general in the sense that it holds for any discrete operators ∂ k x R n and ∂ k θ R n . In addition, it is consistent with the continuous energy estimate given in Proposition 2 since using classical arguments and thanks to [START_REF] Cancès | Numerical analysis of a robust free energy diminishing finite volume scheme for parabolic equations with gradient structure[END_REF] the last term vanishes when the space step and the time step goes to zero. However the last term does not vanish and its sign cannot be determined in general. In particular, in the case of the discrete operators [START_REF] Knudsen | Elements of Newtonian Mechanics[END_REF], the sign of the last term is not defined, which can lead to the artificial creation of mechanical energy.

Numerical discretization ensuring entropy stability.

Let us now discuss a discretization for ∂ x R n and ∂ θ R n such that the energy stability is ensured. We set (26)

∂ k xE R n+1 =    - R n+1 k -R x k , χ n , ζ n+1 , θ n+1 χ n+1 -χ n , if χ n+1 = χ n ∂ k xC R n+1 , else and ∂ k θE R n+1 =    R x k , χ n , ζ n , θ n+1 -R n k θ n+1 -θ n , if θ n+1 = θ n ∂ k θC R n+1 , else.
With the discrete operator [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF], the discrete counterpart of the energy dissipation law Proposition 2 can be proven. Proposition 4. Consider a no-flux boundary on Ω x . Under the CFL condition (20) the scheme (18)-( 22)- [START_REF] John | On the motion of floating bodies[END_REF]

with ∂ k x R n = ∂ k xE R n and ∂ k θ R n = ∂ k θE R n
given by (26) admits for any 0 ≤ n < N the following dissipation law

∂ n+1 t E + E ≤ k∈T k ρgh n k + p n+1 k ∂ n+1 t B k +F n+1 ζ ∂ n+1 t ζ -F n+1 χ ∂ n+1 t χ + T n+1 θ ∂ n+1 t θ.
Proof. Assume first the unknowns χ and θ at two successive time steps are not the same, i.e. χ n+1 = χ n and θ n+1 = θ n . With the discrete operators (26) the last term of the estimate of Proposition 3 reads

δ n+1 t ∂ n+1 t R k + ∂ n+1 t χ∂ k xE R n -∂ n+1 t ζ -∂ n+1 t θ∂ k θE R n = R n+1 k -R n k -R n+1 k -R x k , χ n , ζ n+1 , θ n+1 -ζ n+1 -ζ n -R x k , χ n , ζ n , θ n+1 -R n k = R x k , χ n , ζ n+1 , θ n+1 -R x k , χ n , ζ n , θ n+1 -ζ n+1 -ζ n .
Using (3) we conclude that this term vanishes. Otherwise, if χ n+1 = χ n (resp. if θ n+1 = θ n ) the evolution of the mechanical energy does not depend on ∂ k x R n+1 (resp. ∂ k θ R n+1 ). So any well defined and consistent discretization can be used in this case without affecting the mechanical energy estimate. We conclude that the using the discretization [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF], the reminder of Proposition 3 and the discrete energy dissipation is ensured.

3.4. Practical details. Some details about the practical implementation are given below.

3.4.1. Buoy implementation. The discretization of the roof R (x, χ, ζ, θ) is defined from the buoy surface written under the form S (x, z, θ) = 0, see [START_REF] Audusse | A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation[END_REF]. However in practice, the geometry of a realistic buoy is not in general parametrized in this form. Let us explain in the following the numerical implementation of the discrete roof R δ (x, χ, ζ, θ) the geometry of a buoy with a general shape in a one-dimensional framework.

In practice the buoy is represented thanks to a finite sequence S δ 0 of N P points (X i , Z i ) t 0≤i≤N P -1 . The points describe the vertices of a polygon approaching the geometry of the buoy in the frame centered on the center of mass. For a set of variables Λ = (χ, ζ, θ), the buoy is described by the sequence of points

S δ (χ, ζ, θ) = M θ X Z + χ ζ | X Z ∈ S δ 0
where the rotation matrix M θ is defined in [START_REF] Bocchi | On the return to equilibrium problem for axisymmetric floating structures in shallow water[END_REF]. Let i (resp. i) be the index of the point of the buoy at the left (resp. at the right) extremity, i.e.

i (χ, ζ, θ) = max 0≤i≤N P -1 {∀j, 0 ≤ j ≤ N P -1, X i ≤ X j } and i (χ, ζ, θ) = min 0≤i≤N P -1 {∀j, 0 ≤ j ≤ N P -1, X i ≥ X j } .
The inferior surface elevation in the frame of the buoy is defined by the subsequence

S δ (χ, ζ, θ) = X Ii(χ,ζ,θ) , Z Ii(χ,ζ,θ) t ∈ S δ (χ, ζ, θ) | 0 ≤ i ≤ J (χ, ζ, θ) with I i (χ, ζ, θ) = i + i (χ, ζ, θ) mod N P and J (χ, ζ, θ) = i (χ, ζ, θ) -i (χ, ζ, θ) mod N P .
For any position x not under the buoy, i.e. x < X i(χ,ζ,θ) or x > X i(χ,ζ,θ) , we set R δ (x, χ, ζ, θ) = Π, see [START_REF] Bristeau | An energy-consistent depth-averaged Euler system: derivation and properties[END_REF]. In the part of the domain where the buoy is, i.e.

X i(χ,ζ,θ) ≤ x ≤ X i(χ,ζ,θ)
, the roof is deduced by linear interpolation from the points in S δ (χ, ζ, θ). More precisely, let

i x (χ, ζ, θ) = max i∈I(χ,ζ,θ) X Ii(χ,ζ,θ) ≤ x and i x (χ, ζ, θ) = min i∈I(χ,ζ,θ) X Ii(χ,ζ,θ) ≥ x k .
In other words,

X I i x (χ,ζ,θ) , Z I i x (χ,ζ,θ) t (resp. X I ix(χ,ζ,θ) , Z I ix (χ,ζ,θ) t
) is among the sequence of points S δ (χ, ζ, θ) the nearest vertex from the left (respectively from

φ n , u n Λ n , Λn φ n,q Λ n,q
Λn,q

φ n,q+1 Λ n,q+1
Λn,q+1

φ n+1 , u n+1 Λ n+1 , Λn+1 δ n,q+1 t ⇐ (27) φ n,q+1 ⇐ (28) h n,q+1 p n,q+1 ⇐ (29) Λ n,q+1
Λn,q+1 ⇐ (30) the right) to x. To get an approximation of the roof a linear interpolation is done, i.e.

u n+1 ⇐ (19) φ n,q+1 ← φ n Λ n,q+1 ← Λ n Λn,q+1 ← Λn ξ n,q+1 ← ωξ n,q ξ n,q+1 ← ξ n,q convergence ? i = imax ? q ← 0 i ← 0 no no yes i ← 0 i ← i + 1 q ← q + 1 yes
R δ (x, χ, ζ, θ) = x -X I i x (χ,ζ,θ) Z I ix (χ,ζ,θ) + X I ix (χ,ζ,θ) -x Z I i x (χ,ζ,θ) X I ix (χ,ζ,θ) -X I i x (χ,ζ,θ) .
In the following test cases, the discrete roof is estimated using this strategy, i.e.

R n k = R δ (x k , χ n , ζ n , θ n )
, and in the discrete operators [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF], we use R δ (x, χ, ζ, θ) instead of R (x, χ, ζ, θ).

3.4.2.

Fixed point approach. The scheme ( 18) is implicit and non-linear, see [START_REF] Parisot | Centered-Potential Regularization for the advection upstream splitting method[END_REF]. Because of the derivative of the roof in the dynamic equation of the buoy, the computation of the Jacobian required to apply a Newton method is not trivial. In order to easily introduce new forces on the buoy, we propose a weak coupling iterative process where the buoy and the water dynamics are solved successively until convergence. In addition, the CFL condition [START_REF] Guerber | A fully nonlinear implicit model for wave interactions with submerged structures in forced or free motion[END_REF] does not take into account the buoy dynamics. A time step adaptation is added to reduce the time step when the convergence of the iterative process is too slow. The reader can find more details on time step adaptation in [START_REF] Cancès | Numerical analysis of a robust free energy diminishing finite volume scheme for parabolic equations with gradient structure[END_REF] for example. The iterative process is precisely described in the following. In practice the fixed point algorithm converges in a few iterations (less than 20) with reasonable time steps, but this observation depends on the test case. The rigorous proof of convergence for the fixed point is out of the scope of this work.

Let us describe more precisely the iterative process, illustrated in Figure 2. We consider the state (h n k , u n k ) k∈T , Λ n , Λn given with Λ = (χ, ζ, θ) T . Let the exponent q be related to the iteration number of the iterative process. The iterative process is initialized by φ n,0 k computed from h n k using (15), Λ n,0 = Λ n and Λn,0 = Λn . To obtain the next approximation q + 1 from the approximation q, an approximation of the time step δ n,q+1 t is first computed such that the CFL condition ( 20) is satisfied at convergence, i.e. [START_REF] Lannes | On the dynamics of floating structures[END_REF] δ n,q+1 t = ξ n,q min h n,q k , h n,q k+1 2 h n,q k + h n+1 k,q min ( k , k+1 )

u n k +u n k+1 2 + γ 2 φ n,q k+1 -φ n,q k 2
where 0 < ξ n,q ≤ 1 is an adaptive CFL parameter, initialized by ξ n,0 = 1.

Then we compute the potential φ n,q+1 k considering the roof R n,q k obtained from the buoy position Λ n,q , see Section 3.4.1. Since ( 18) is implicit and non-linear, we use a linearization to compute the next approximation. More precisely the potential is computed by solving the linear system (28)

J n,q φ (φ n,q ) δ n,q+1 φ = Π n,q (φ n,q ) φ n,q+1 = φ n,qδ n,q+1 φ where J n,q φ (φ) is the Jacobian of the residual Π n,q (φ) = (Π n,q k (φ)) k∈T of ( 18), i.e.

Π n,q k (φ) = H φ, B n+1 k , R n,q k -h n k + δ n,q+1 t k ×       F             φ k B n+1 k R n,q k u n k k       ,       φ k+1 B n+1 k+1 R n,q k+1 u n k+1 k+1       , δ n,q+1 t       -F             φ k-1 B n+1 k-1 R n,q k-1 u n k-1 k-1       ,       φ k B n+1 k R n,q k u n k k       , δ n,q+1 t             .
Since this formula comes from a linearization of the mass flux in the finite volume formula [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] it is clear that the volume of water is exactly preserved up to the resolution error of the linear system [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF]. From φ n,q+1 k , we compute the water depth and the surface pressure using (16), i.e.

(29) h n,q+1 k = H φ n,q+1 k , B n+1 k , R n,q k and p n,q+1 k = P φ n,q+1 k , B n+1 k , R n,q k .
This pressure is then used to compute the external force, i.e. ( 30)

Λ n+1 = Λ n + δ n+1 t Λn + δ n+1 t 2 2 C -1 1 F P + F E p n,q+1 k , ∂ k x R n,q , ∂ k θ R n,q , t n+1 Λn+1 = Λn + δ n+1 t C -1 1 F P + F E p n,q+1 k , ∂ k x R n,q , ∂ k θ R n,q
, t n+1 with the discretization of the external force given by ( 23) while C 1 is given in (4). Note that even if the implicit version of the Newmark scheme is considered, i.e. α = β = 1, it can be computed explicitly thanks to the iterative process.

At this step, we get an approximation of the state variable φ, Λ, Λ at time t n+1 , that does not necessarily satisfy the constraint (2). To impose the constraint on each time step, the previous operations are repeated until numerical convergence. Two stopping criteria are used to increase the reliability of the result, based on the L ∞ -norm of the residual of the non-linear system (18) (already used in the case without structure see [START_REF] Parisot | Centered-Potential Regularization for the advection upstream splitting method[END_REF]) and the L ∞ -norm of the difference between two iterations of the buoy variable Λ. When the two errors estimates are inferior to a tolerance ε N ew , we consider the fixed point as converged. When the vertical velocity of the buoy is not too large, the convergence of the iterative process is acceptable (less than 20 iterations). However, when the buoy vertical velocity is large and especially when it reaches the water surface for the first time, convergence can be hard. To avoid two-way wiring in the fixed point, the time step is limited if too many iterations are required. More precisely, we set an acceptable maximum number of iterations i max in the fixed point. If the Newton iterative process reaches the number of iterations i max , the adaptive CFL parameter is reduced, i.e. ξ n,q+1 = ωξ n,q with a given reduction factor 0 < ω < 1, and we start from the previous time step, i.e. φ n,q+1 , Λ n,q+1 , Λn,q+1 = φ n , Λ n , Λn . Finally, the velocity u n+1 k is estimated using [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF].

Simulations

We illustrate the validity of our coupled numerical model by some one dimensional computations. In the following we take g = 9.81 and ρ = 1. The mesh is a regular grid with the space step k = δ x . The relaxation parameter is chosen as λ = √ δ x , see [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] for more details about this choice. The regularization parameter γ is set to 2. It is slightly increased compared to the optimal value necessary to get the entropy estimate see Lemma 3.1. This has a regularizing effect (see [START_REF] Parisot | Centered-Potential Regularization for the advection upstream splitting method[END_REF]) and prevents oscillations which can be observed with the CPR scheme in the vicinity of a shock.

All the test cases are computed with the different strategies presented above, i.e. ( 25) and [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]. When the difference between the results of the two strategies are not significative, only the result with the entropy-satisfying scheme (26) is illustrated. The numerical parameters (see Section 3.4) are set to ε N ew = 10 -10 , i max = 20 and ω = 0.5.

Let us give a brief motivation of the following test cases. In Section 4.1 we compare the solution computed with the numerical strategy presented above to an analytical solution established in [START_REF] Lannes | On the dynamics of floating structures[END_REF]. This test case emphasizes the convergence of the proposed numerical strategy. However, the two different discretizations [START_REF] Knudsen | Elements of Newtonian Mechanics[END_REF] and ( 26) cannot be compared in this context since the horizontal displacement and the rotation of the buoy are frozen, and the two schemes are equivalent in this context. We compare the two different discretizations [START_REF] Knudsen | Elements of Newtonian Mechanics[END_REF] and [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF] in Section 4.2, by performing a test case where all the degrees of freedom of the buoy are free. In addition, this test case emphasizes the robustness of the scheme to detect the regions where the buoy reaches the water. Eventually, in Section 4.3, we propose a simple but realistic configuration of a wave energy converter to illustrate possible applications of this work. 4.1. Return to equilibrium. To illustrate the consistency of the scheme, we compare the computed solution to an analytical solution. In [27, Corollary 1], the analytical solution of a return to equilibrium is established. It corresponds to the vertical movement of a rectangular buoy, of horizontal length L, freely floating above a water initially at rest and with a flat bottom B = 0. A similar analytical solution in axisymmetric framework is established in [START_REF] Bocchi | On the return to equilibrium problem for axisymmetric floating structures in shallow water[END_REF]. Note this configuration satisfies neither Hypothesis 1, Hypothesis 2 nor the assumption of the derivation of the fluid equations ( 1)- [START_REF] Aissiouene | A combined finite volume -finite element scheme for a dispersive shallow water system[END_REF], see [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF][START_REF] Lannes | On the dynamics of floating structures[END_REF]. However, for what concerns the numerical scheme, it can be written without difficulty and this test case validates the accuracy of the method in a dynamic context.

The equilibrium of the system is defined by u = 0, a fixed h eq ∈ R + and ζ = ζ eq such that [START_REF] Matt | A review of numerical modelling of wave energy converter arrays[END_REF] M g = Ωx p dx = ρ Ωx (h eq -(S (xχ eq , θ eq ) + ζ eq )) + dx where χ eq ∈ R and θ eq ∈ [-π, π] are fixed. For simplicity we set χ eq = θ eq = 0. In other words the equilibrium is the situation where the surface pressure balances out the weight of the buoy. Only a vertical movement is considered for this test case. More precisely we suppose that (32)

F χ = Ωx p∂ x R dx, T θ = - Ωx p∂ θ R dx and F ζ = 0
which expresses the fact that the horizontal and angular degrees of freedom are fixed.

Initially the water is at rest so that h 0 = min h eq , S (x, 0) + ζ 0 and u 0 = 0. The length of the roof is L = 0.8 and the roof is defined so that it is symmetric with respect to its center of mass. For the above defined configuration, the position of the center of mass is determined by

(33) M ρ + 2 3ζ L 2 3 ζ = -gL (ζ -ζ eq ) -ν ζ + µ (ζ) ζ2
with the coefficients

ν ζ = gL h eq -τ L 4 √ g ζ 2 and µ (ζ) = 4 3ζ 2 L 2 3 .
The coefficient τ is given by τ (r) = 1 .

In this test case, we consider a buoy of mass M = 0.5 with the equilibrium elevation ζ eq = 1 and the water depth h eq computed to satisfy [START_REF] Matt | A review of numerical modelling of wave energy converter arrays[END_REF]. Initially the state of the buoy is ζ 0 = 1.2 and ζ0 = 0. For the computation, we set Ω x = [-20 , 20] with wall condition at boundaries.

The simulation is compared to the solution of (33) in Figure 3, where δ G = ζ -ζ eq . The agreement between the two solutions is good and the solution tends to the equilibrium position as expected. The convergence rate in norm L 2 in time until the end of the simulation (here fixed to 4) of δ G is shown in Figure 4. A convergence rate slightly better than one is observed. We recall that a first order scheme is expected since the full implicit Newmark scheme and the fluid scheme are first order.

Remark 1. Because of the symmetry of the test case, the external forces actually read F χ = T θ = 0. However, this equilibrium is unstable since small perturbations, such as those introduce by the machine errors, break the symmetry. For this reason, even with an axisymmetric buoy, the external forces (32) are imposed. 4.2. Throwing. To illustrate the robustness of the scheme, we propose the simulation of a buoy thrown into the water. This simulation can be compared to a block of ice falling in water from an iceberg, except that the parameters such as the mass and the angular moment are not physically relevant.

We consider a closed domain Ω x = [0, 2] with wall boundary conditions and a flat bottom B = 0. The water is initially at rest, i.e. h 0 = 1 and u 0 = 0. The buoy is defined following the strategy explained in Section 3.4.1 with

S δ 0 = X Z i N P 2π , 0 ≤ i ≤ N P -1 , X Z = 0.
3 cos (ϑ) + 0.03 cos (4ϑ) 0.2 sin (ϑ) -0.03 cos (4ϑ) -0.06 , its mass is M = 10 -2 and its angular moment is J = 10 -2 . We do not consider any additional force, i.e. F χ = F ζ = T θ = 0. The initial state of the buoy is χ 0 = 0, ζ 0 = 1.5, θ 0 = π 2 χ0 = 1, ζ0 = 1 and θ0 = 0.

In Figure 5, the buoy position, the water height and the surface pressure at different times are shown with the numerical parameters δ x = 10 -3 and N P = 360. At the beginning the fluid is at rest and the buoy is not in contact with the fluid. At the impact an energy transfer from the buoy to the water is observed and waves are created. It is easy to verify that the non-penetration constraint (2) is qualitatively met, more precisely at the order λ, see [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF]. Let us consider the mechanical energy E

n + E n without the potential energy of the water at rest E 0 and normalized by the energy of the buoy at the initial state E 0 , i.e.

e n = E n + E n -E 0 E 0 .
At the initial state, the main part of the energy of the system is the potential energy of the water. By removing the potential energy of the water at rest, the direction of energy variation is unchanged but the amplitude of the variation is significantly increased.

The evolution of the energy ratio e n is shown in Figure 6 (first line) for both discretizations [START_REF] Knudsen | Elements of Newtonian Mechanics[END_REF] and [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF] and for two different space steps. To give an idea of the cost of the simulation, let us indicate some CPU times for computations done on a personal laptop to reach time t = 3: for the discretization [START_REF] Knudsen | Elements of Newtonian Mechanics[END_REF], the CPU times are 0, 474s for δ x = 10 -2 and 10.7s for δ x = 10 -3, while for the discretization [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF] the CPU times are 0.533s for δ x = 10 -2 and 12.7s for δ x = 10 -3.

Accordingly to the theoretical result of Proposition 4, the mechanical energy is decreasing for the discretization [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF], see Figure 6 (first line). The dissipation of mechanical energy is low (few ‰) even with a coarse resolution. Note that even at convergence, we do not expect that the mechanical energy is preserved. In particular, when the buoy reaches the water surface for the first time, the surface pressure is discontinuous and it yields a strong dissipation of energy, see time t = 0.35. In practice the centered discretization (25) leads to satisfactory results except for coarse resolution. With δ x = 10 -2 , the mechanical energy is increasing at time t = 1.5 and the results is significantly affected, in particular for the horizontal position, see Figure 6 (second line).

4.3.

Wave energy converter. To illustrate the ability of the strategy to optimize a wave energy converter, a restoring and a damping force are introduced, see Figure 7. The structure is allowed to move vertically and to rotate around its center of mass and it is fixed in the horizontal direction, i.e. The planar movement of the energy wave converter is described by the system (34)

       M ζ = -M g + Ωx p dx -c ζ -κ ζ -ζ J θ = Ωx p∂ θ R dx.
The following energy estimates can be obtained Lemma 4.1. The solution of (34) satisfies the following energy law

∂ t E = Ωx p∂ θ R dx θ + Ωx p dx ζ -c ζ2 with E (ζ, θ) = J 2 θ2 + M 2 ζ2 + M gζ + κ 2 ζ -ζ 2 .
Proposition 5. Considering a no-flux boundary on Ω x , any smooth solution of (1)-( 2)-( 11)-( 12)-( 34) satisfies the following energy balance law

∂ t E + E = Ωx (p + ρgh) ∂ t B dx -P c with P c = c ζ2 .
The term P c is a dissipation term due to the damping of the spring. A part of this energy could be recovered for energy production.

Using again the Newmark scheme [START_REF] Iguchi | Hyperbolic free boundary problems and applications to wavestructure interactions[END_REF]

       M ζn = -M g + k∈T k p n k -c ζn -κ ζ n -ζ J θn = k∈T k p n k ∂ k θ R n-1 .
The following discrete estimates can be obtained Lemma 4.2. The scheme (22)-( 35) satisfies for any 0 ≤ n < N the following energy law

∂ n+1 t E + δ n+1 t 4 ζ 2 = k∈T k p n+1 k ∂ k θ R n ∂ n+1 t θ + k∈T k p n k ∂ n+1 t ζ -P n+1 c -κ δ n+1 t 2 ∂ n+1 t ζ 2 with E n = E (ζ n , θ n ) and P n+1 c = c 2   ∂ n+1 t ζ 2 + ζn+1 + ζn 2 2   .
Proof. The proof is similar to the proof of Lemma 3.2. The coefficients are (with obvious notations)

f n =     -M g + κζ + k∈T k p n k k∈T k p n k ∂ k θ R n-1     , C 1 = M 0 0 J , C 2 = c 0 0 0 and C 3 = κ 0 0 0 . Proposition 6.
Consider a no-flux boundary on Ω x . Under the CFL condition (20) the scheme (18)-( 19)-( 22)-( 35) with ∂ k θ R n = ∂ k θE R n given by (26) admits for any 0 ≤ n < N the following dissipation law

∂ n+1 t E + E + δ n+1 t 4 ζ 2 ≤ k∈T k ρgh n k + p n+1 k ∂ n+1 t B k -P n+1 c -κ δ n+1 t 2 ∂ n+1 t ζ 2 .
The optimization of wave energy converters is a key point in the development of such devices. To show the feasibility of an optimization procedure, a naive estimate of the recoverable energy is done. The recovered energy, depending on the technology used to convert wave energy, is a part of the recoverable energy defined by

E c = N n=1 δ n t P n c
where N is the number of time iterations and P n c is defined in Lemma 4.2.

In the following a simple case of a wave energy converter attached to a spring over a flat bottom B = 0 is simulated. An elliptical shaped buoy with semi-minor and semi-major axis respectively equal to 0.2 and 0.4 is considered, see Section 2.2 for the description of the buoy geometry. The mass is M = 10 -2 and the moment of inertia J = 5 • 10 -3 . The stiffness coefficient is κ = 10 -2 and the characteristic elevation is ζ = 1.4. Initially, the water and the buoy are at equilibrium. More precisely, the initial state of the buoy reads χ = 1.5, ζ = 1.4, θ 0 = π 2 , ζ = 0 and θ0 = 0, and the water state reads u 0 = 0 and h 0 = min h eq , S xχ, θ 0 + ζ 0 with h eq so that the buoy is at equilibrium, see [START_REF] Matt | A review of numerical modelling of wave energy converter arrays[END_REF]. The computational domain is Ω x = [0, 3]. At the left, a fixed water depth boundary condition is considered, see [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]Section 21.8], with h (0, t) = h eq + 0.05 sin (4πt). At the right boundary a wall boundary condition is considered.

In Figure 8 the ratio between the recoverable energy E c and the total energy introduced E tot is shown at time t = 1 with the space step δ x = 3 • 10 -3 and for several values of c. The total energy introduced is computed by considering a similar simulation without the buoy. More precisely it is defined by E tot = E

N -E 0 .

The recoverable energy E c is about 10% of the total energy E tot that is not negligible considering that many parameters can also be optimized (shape, mass and moment of inertia of the buoy). In this case, the best choice seems to be obtained around c = 0.3. Physically a small stiffness allows the spring to move as much as possible. Furthermore a small damping coefficient would not be able to absorb any In contrary if the damping coefficient is too high, the buoy will be slowed down due to the damping and the recoverable energy would not be optimal.

Conclusion

A strong coupling between a congested shallow water type model and Newton's second law of motion is presented for the modeling of floating structures. We have taken a particular care to the energy transfer between the solid and the water. An entropy correction is made at the discrete level in order to ensure an entropy dissipation law and the numerical scheme proposed in [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] is adapted to take into account a freely floating body. Finally numerical simulations are proposed to validate the approach and to show the feasibility of physically relevant cases. Now real life applications need an extension to the two dimensional case. Since the computations become CPU costly in higher dimensions, some scientific computing issues should be resolved beforehand. One can think of linking our method with a less CPU consuming scheme for free surface flow far from the congested area.

The physical description of water waves is important when considering wave energy converters. More complex models such as dispersive models [START_REF] Bristeau | An energy-consistent depth-averaged Euler system: derivation and properties[END_REF][START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF], layerwise vertical discretization [START_REF] Audusse | A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation[END_REF] or the combination of both [START_REF] Fernández-Nieto | A hierarchy of dispersive layer-averaged approximations of Euler equations for free surface flows[END_REF] could be used to represent the flow. Dispersive models do not neglect the vertical acceleration and are therefore more adapted to characterize heave effets. In any numerical strategies where the dispersive terms are introduce as a correction of a shallow water scheme [START_REF] Aissiouene | A combined finite volume -finite element scheme for a dispersive shallow water system[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Favrie | A rapid numerical method for solving serre-green-naghdi equations describing long free surface gravity waves[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF][START_REF] Parisot | Entropy-satisfying scheme for a hierarchy of dispersive reduced models of free surface flow[END_REF], the strategy presented in the current work can be use instead of the classical shallow water scheme to obtain floating body modeling with a dispersive model. In addition, a layerwise discretization can better reproduce wind effects and friction on the body by considering a vertical profile of the velocity. We expect that the improvement of our results with a layerwise discretization not raise particular difficulties since the velocity is computed separately to the constraint in the current numerical strategy. Last but not least, our method is directly adaptable to take into account several bodies, as long as no collision between them appear. It is in particular interesting for the simulation and optimization of farms of wave energy converters [START_REF] Penalba | A numerical study on the hydrodynamic impact of device slenderness and array size in wave energy farms in realistic wave climates[END_REF]. A logical follow-up in another direction is to consider the collisions between floating bodies. It should be essential to study the drift of floating fragments during flooding. Since these phenomena take place on a large space scale, the proposed model seems well adapted. Eventually, a challenging objective is the handling of submerged body [START_REF] Ducassou | A fictitious domain approach based on a viscosity penalty method to simulate wave/structure interaction[END_REF][START_REF] Guerber | A fully nonlinear implicit model for wave interactions with submerged structures in forced or free motion[END_REF]. This issue seems not at hand at the moment since a more complex description of the flow is required.

Figure 1 .

 1 Figure 1. Floating body configuration and illustration of the unknowns.

Lemma 2 . 1 .

 21 Any smooth enough solution of the congested shallow water model (1)-(2) satisfies the following energy balance law

  For a fixed (χ, θ), we denote by Ω R (χ, θ) ⊂ Ω x the domain of definition of R (x, χ, ζ, θ). The domain of definition Ω R actually does not depend on the buoy elevation ζ since it is a vertical shift.

Figure 2 .

 2 Figure 2. Scheme of the weak coupling iterative process. The notation a ← b is the numerical affectation, i.e. a gets the value of b. The notation a ⇐ (X) means that a is computed using the equation (X).

3 h

 3 eq + C (r) +

Figure 3 .

 3 Figure 3. Return to equilibrium: time evolution of the distance to the equilibrium position δ G = ζζ eq

Figure 4 .

 4 Figure 4. Return to equilibrium: Convergence rate in L 2 -norm for δ G .

  F ζ = -c ζκ ζζ , T θ = 0 and F χ = Ωx p∂ x R dx.The damping (resp. stiffness) coefficients are denoted by c ∈ R +

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. Throwing: Buoy (black dotted line), water height (blue solid line), velocity (filled surface in the water -blue< 0 <red) and surface pressure (green filled surface in the buoy) at several times.

Figure 8 .

 8 Figure 8. Spring: energy ratio for different values of the damping coefficient

Acknowledgments

The authors wish to express their warm thanks to Cindy Guichard for many fruitful discussions and kind advices. This research is part of a collaborative program with the Anr-Fem (France Energies Marines) project Hyflo-Eflu through the project call ite emr 2015.