Philippe Destuynder
email: philippe.destuynder@cnam.fr

Caroline Fabre
email: caroline.fabre@u-psud.fr

A discussion on Neumann-Kelvin model for progressive water waves

Keywords: AMS Subject Classification: Primary: 74F10, 35L85, 35R35; Secondary: 76B15, 35L70, 35Q35 water waves, unstabilities, wave equations, spectral theory, non-linear models

This paper focus on the role of the capillary in the modelling and the analysis of progressive waves at the surface of the sea. It is proved that the classical Neumann-Kelvin model is fully instable and should be modified in order to be mechanically reliable. It is shown that the easiest way to take the capillary into account is to use a formulation where the normal displacement of the free surface appears as a main unknown variable.

Introduction

Since the work of Kelvin [KEL 87], many contributions have been published for modelling progressive waves at the surface of the sea. In most of them, the authors consider that the capillary is neglectible [STO 57], [DOU 95]. This is a well founded hypothesis as far as one compares the forces in presence at the surface of the sea. Unfortunately, as it has been pointed out by several authors [CHEN 01] [IOS 96] [DIKH 99], the model is ill posed in a sense that is made precise in the following. In this paper, an extension of classical models is presented which takes into account the surface energy of the sea. A mathematical 1 analysis of it is given, and it is explained why it should be numerically more stable than the usual one. The discussion is governed by the value of the stationary flow's velocity. Comparaisons with the regular Neumann-Kelvin model enable one to understand the limits of the classical progressive wave theory. Before detailing the context of our work, let us precise that the literature on water waves is very plentiful. Let us mention three works which seem to us particularly interesting in our context: the work of B. Alvarez-Samaniego and D. Lannes [ALLA 08], based on Zakharov 's model, presents and justifies the main asymptotic models such as shallow water equations, Boussinesq's systems (and some more) and its references. Let us note that their framework neglects the capillary and the free surface of the water is unbounded (which enables one to use Fourier transform). Capillary-gravity waves on the free surface of a two dimensional inviscid fluid of infinite depth are studied in [IOS 96] and [DIKH 99] (and their references) where they prove the existence of solitary waves.

In this paper, one uses two functions in order to describe the wave phenomenon occuring at the surface of a water pool which is moving forwards at a given permanent velocity U . The fluid is assumed to be incompresssible and the ordinary waves are due to the gravity and the capillary on the free surface. The mathematical models that are analyzed in the following can be used in a very general framework as soon as the fluid is assumed to be inviscid. Furthermore it is also necessary to assume that if the boundary is not connected, the flow velocity has a zero circulation around each connected component of the boundary of the pool. For instance, there can be a moving submarine under the free surface of the water but its rotation should be zero. Then the two functions which enable to describe the movement of the water are a velocity potential function -say Φ-and the normal displacement -say η-of a geometrical point of the free surface. Let us underline that this is not the displacement of the water particle at this point. In our case, the kinematical continuity on the free surface which will be considered as a separated media from the interior of the sea, has to be written for a material point. The coupling between the two media (the surface of the sea and the interior of the sea) requires as usual two sets of relations: one for the kinematical continuity and the other for the equilibrium of the stresses. But let us notice that the two waves (gravity waves and capillary waves) are essentially surface waves. The model that is discussed is based on the Neumann-Kelvin formulation [ In the present analysis it is suggested that the formulation with respect to η is more appropriate as far as numerical simulations are concerned. Our discussion is based on mathematical arguments and more precisely on the spectral theory of the involved operators. The dominating result is that the wave model by Neumann-Kelvin is ill-posed and one can't neglect the capillary effect on the free surface. A consequence is that a numerical simulation of the classical model is meaningless unless an artificial numerical cut-off of the high frequencies is used. The capillary is the right physical method for defining this cut-off. There are other artificial methods used by engineers for instance by considering a larger mesh of the free surface or by averaging the solution on a larger area than the one used for the solution method. The theoretical results given in this paper can be extended to the hydrodynamic behavior of a ship for the modelling of the whipping and the slamming phenomenae [ Our plan is the following: in section 2, the interior fluid model is presented and the added mass operator is defined. Its mathematical properties are discussed in section 3. This operator is a corner stone in order to analyze the model of the free surface which is given in section 4, where the steady state is also studied. In section 5, the mathematical analysis of the full linearized coupled model is done when it is expressed in terms of the normal displacement of the free surface. It is proved that the model used in this paper is more stable than the one suggested by Neumann and Kelvin. Let us split our analysis into two main cases according to a small or a large value of the velocity of the flow. In section 6, the formulation based on the velocity potential is discussed for the Neumann-Kelvin's model. In section 8, few remarks on the non linear case are given and the section 9 concerns the numerical schemes.

Finally, let us mention that the main geometrical notations are explicited on figure 1 and that constant c may change from one line to another.

The interior fluid model

The pool containing the water occupies in R d (d = 2 or 3) a non empty, open, bounded and connected set denoted by Ω, the boundary of which is denoted by Γ = ∂Ω. Three parts of this boundary are distinguished and they are assumed to be with non empty interiors: one denoted by Γ 0 corresponds to the bottom of the pool and the boundary of a submarine (see figure 1); a second one -say Γ fis the free surface and the last one Γ 1 , contains the intake and the exit through which the water is entering and exhausting from the pool. It is assumed that Γ = Γ 0 ∪ Γ 1 ∪ Γ f and that Γ f is open. Let us denote by ν the unit outwards normal to the boundary Γ.

In all the paper, the following boundary conditions are used :

• On Γ 1 , the normal velocity of the steady flow is prescribed equal to U (e 1 , ν) where ( , ) denotes the scalar product of two vectors in R d . The transcient pressure is assumed to vanish on Γ 1 . This leads to an homogeneous Dirichlet boundary condition for the transcient component of the velocity potential. Let us remark that this is not the only possibility. Indeed, another one, may be more realistic, would be to assume that the free surface Γ f is a stream surface for the steady flow (or a stream line in 2D). Let us suppose that the free surface Γ f is initially located in the space x 3 = 0, whereas the coordinates (x 1 , x 2 ) are varying in the plane containing Γ f (see figure 1).

Entry

The velocity potential which describes the steady flow is denoted by Φ 0 . It is solution of the following system :

                 -∆Φ 0 = 0 in Ω, ∂Φ 0 ∂ν = 0 on Γ 0 , ∂Φ 0 ∂ν = U (e 1 , ν) on Γ 1 , ∂Φ 0 ∂ν = 0 on Γ f , Γ f Φ 0 = 0.
(1)

It is known that there is a unique solution in the space H 1 (Ω) to (1) as soon as the following Fredholm's condition (see for instance [NEC 67]) is satisfied:

Γ1 (e 1 , ν) = 0. ( 2
)
It is assumed to be valid. Let us remark that the hypothesis (2) is not restrictive and is satisfied in each physical case because it is a consequence of the incompressibility of the fluid. If Ω is a parallelepiped, the condition ( 2) is true and one can easily prove that Φ 0 (x) = U (x 1 + c) with some constant c. It satisfies ∆ s Φ 0 = 0 on Γ f where ∆ s is the inplane Laplace operator on Γ f . In other cases, one can write Φ 0 = U (x 1 + Φ 1 ) where:

                 -∆Φ 1 = 0 in Ω, ∂Φ 1 ∂ν = -(e 1 , ν) on Γ 0 , ∂Φ 1 ∂ν = 0 on Γ 1 ∪ Γ f Γ f Φ 1 = - Γ f x 1 . Let us assume that Φ 1 | Γ f ∈ C 3 (Γ f ) and thus Φ 0 | Γ f ∈ C 3 (Γ f ). This hypothesis is fullfilled if Γ0 ∩ Γf = ∅ and if the boundary of Ω is smooth enough (see [GRI 86] or [KON 06]).
As usual, the last relation on Γ f enables to fix the constant, Φ 0 depends on as far as a Neumann's boundary condition is prescribed on all the boundary of Ω.

Let us write Φ 0 = U Φ0 where Φ0 is solution of (1) with U = 1.

Let us now consider the dynamic behaviour of the water. The total velocity potential Φ is written as the summ of two contributions: one is the steady state defined previously and the other one is the dynamic state:

Φ = Φ 0 + ϕ = U Φ0 + ϕ. (3)
The unit normal to the deformed boundary Γ f is denoted by ν ′ . It is different from ν when there is a rotation of this free surface, which is clearly the case for waves and the coupling between the rotation of ν and the velocity of a geometrical point of the free surface leads to gyroscopic forces as it is shown in the following. The boundary condition on the deformed free surface is

∂Φ ∂ν ′ = ∂η ∂t (ν, ν ′ ),
where ∂η ∂t (ν, ν ′ ) is the velocity on the (deformed) free surface in the direction ν ′ The displacement of the free surface is neglected but not its rotation. In fact, the gradient of η appears in the first order term and therefore, it is kept.

Finally, the dynamic velocity potential ϕ satisfies the following set of equations:

                 -∆ϕ = 0 in Ω, ∂ϕ ∂ν = 0 on Γ 0 , ϕ = 0 on Γ 1 , ∂Φ ∂ν ′ = ∂η ∂t (ν, ν ′ ) on Γ f ( which is the velocity on Γ f in the direction ν ′ ).
(4) The last relation is a non-linear one. It is one of the main difficulty in progressive wave modelling. Let us make it more explicit. From a classical computation, one has at the first order in η:

ν ′ = ν -∇ s η
where ∇ s is the inplane gradient operator on Γ f . This boundary condition can be written as follows up to the second order (let us keep in mind that

∂Φ 0 ∂ν = 0 and (∇ s η, ν) = 0 on Γ f ): ∂ϕ ∂ν -∇ s η.∇ s Φ 0 -∇ s η.∇ s ϕ = ∂η ∂t .
The linearized version of this relation around η = 0 and ϕ = 0 is:

∂ϕ ∂ν = ∂η ∂t + ∇ s η.∇ s Φ 0 . (5)
If one considers an initial deformation of the free surface described by a function η 0 , one would obtain:

∂ϕ ∂ν = ∂η ∂t + ∇ s η.∇ s Φ 0 + ∇ s η 0 .∇ s ϕ. (6)
This is the case for instance if one considers the digging effect which is discussed later on. Several possibilities could occur for the boundary condition satisfied by ϕ on Γ 1 . The one used in (4) traduces a vanishing acoustic pressure. Another possibility would be to prescribe the total pressure (linearized expression) to zero (assuming also that the steady velocity is normal to Γ 1 ):

∂ϕ ∂t + U ∂ϕ ∂ν = 0. ( 7
)
Clearly there are many other possibilities in the definition of the boundary condition that should be satisfied by ϕ on this part of the boundary ∂Ω. Furthermore it is certainly validate to use different boundary conditions at the flow-intake and at the exit.

In many cases, the solution Φ 0 , in the neighbourhood of the free surface, is very close to U x 1 + c. Hence, the steady velocity on Γ f is almost uniform and close to U e 1 . This enables one to write approximately:

∂Φ ∂ν ′ = ∇Φ.ν ′ ≃ (U e 1 + ∇ϕ, ν ′ ) = ∂ϕ ∂ν -U ∂η ∂x 1 , (8)
or else:

∂ϕ ∂ν ≃ ∂η ∂t + U ∂η ∂x 1 , (9)
which is the formula used in most publications on progressive wave theory (see J.J. Stoker [STO 57]). But for the general cases, the right linearized boundary condition is:

∂ϕ ∂ν = ∂η ∂t + U ∇ s Φ0 .∇ s η. ( 10
)
This is the one that is considered in this paper. The difference between these two expressions occurs for instance if a submarine is close to the free surface or if the depth of the sea is small and if the bottom is not flat (shore effect).

It is worth noting that one needs an additional regularity on the term ∇ s Φ0 in order to be able to make sense to the previous boundary condition. The point is that Φ0 is a smooth function excepted occasionally near corners of the boundary ∂Ω. For instance one can assume that as far as the angles between Γ 1 and Γ f are near from π/2 then Φ0 ∈ C 1 (Ω) and therefore the restriction on Γ f of the gradient of Φ0 is continuous. But, as mentionned, we assume more regularity on Φ0 . Details concerning the regularity of the solution to elliptic partial differential equations can be found in

W 0 = {ψ ∈ H 1 (Ω), ψ = 0 on Γ 1 }. (11)
At any time, the function ϕ = ϕ(t, x) should be solution of:

       ϕ ∈ W 0 , ∀ψ ∈ W 0 , Ω ∇ϕ.∇ψ = Γ f ∂η ∂t ψ + U ∇ Φ0 .∇ s ηψ. (12)
Hence, ϕ appears as a linear function of η. More precisely, there are two contributions: one is a linear function of ∂η ∂t and the other one is a linear function of U ∇ s Φ0 .∇ s η. This operator, which maps the right hand side of (12) into ϕ is named the added mass operator and is denoted by G. It plays a crucial role in our analysis and this is why the next section gives a description and few properties of G which are useful in all what follows. Having in mind these properties related to G, we will then be able to introduce the model of equilibrium related to the free surface.

Definition and properties of G

Let us first state some mathematicals notations and let us introduce functionnal spaces that will be used in the following.

3.1 The spaces H1/2 (Γ f ) and its dual

Standard notations are used for Sobolev spaces and theirs norms; || . || s,X denotes the norm in H s (X) where X is an open set of R d (d ≥ 1). For instance for s = 0 one has: H 0 (X) = L 2 (X). Let us consider the functional space W 0 defined at (11) equipped with the norm:

||ψ|| W0 = ||∇ψ|| 0,Ω .
Let us define the following closed subspace of H 1/2 (Γ f ):

H1/2 (Γ f ) = {v ∈ L 2 (Γ f ), ∃ψ ∈ W 0 , v = ψ |Γ f }. ( 13
)
The space H1/2 (Γ f ) is equipped with its natural norm which induces on it, an Hilbert's structure :

|||v||| 1/2 = inf ψ∈W0,v=ψ |Γ f ||ψ|| W0 = min ψ∈W0,v=ψ |Γ f ||ψ|| W0 . (14) Obviously, D(Γ f ) ⊂ H1/2 (Γ f ), H1/2 (Γ f ) is continuously embedded in H 1/2 (Γ f ) and we have ∀v ∈ H1/2 (Γ f ), ||v|| 1/2,Γ f ≤ |||v||| 1/2 (15) Remark 1. The space H1/2 (Γ f ) can be H 1/2 (Γ f ) or H 1/2 00 (Γ f ) (see [LM 68
]) in particular cases depending on the geometry of the open set Ω and the conditions prescribed on Γ 1 in the space W 0 . For instance it is 1).

H 1/2 (Γ f ) if Γ 1 is disconnected from Γ f and: H 1/2 00 (Γ f ) if Γ 1 is connected to Γ f and Γ0 ∩ Γf = ∅ (as on figure
The dual space of H1/2 (Γ f ) is denoted by H-1/2 (Γ f ), its norm is denoted by ||| . ||| -1/2 . and the duality bracket between H1/2 (Γ f ) and its dual is denoted by , .

For regular functions η ∈ L 2 (Γ f ) and v ∈ H1/2 (Γ f ), one has

η, v = Γ f ηv =< η, v > -1/2,1/2
where < , > -1/2,1/2 denotes the duality between H -1/2 (Γ f ) and

H 1/2 (Γ f ). The space H -1/2 (Γ f ) is continuously embedded into H-1/2 (Γ f ) and ∀η ∈ H-1/2 (Γ f ), |||η||| -1/2 ≤ ||η|| -1/2,Γ f (16)
Let us now turn to the study of the added mass operator G.

Analysis of the operator G.

The definition of the linear operator G is :

g ∈ H-1/2 (Γ f ) → G(g) = Φ |Γ f ∈ H1/2 (Γ f ) , (17)
where Φ = Φ(g) ∈ W 0 is the unique solution of :

∀ψ ∈ W 0 , Ω ∇Φ.∇ψ = g, ψ |Γ f . (18)
On the one hand Φ is uniquely as a function of W 0 and on the other hand the trace on Γ f of a function in W 0 is in the space H1/2 (Γ f ). Therefore G is well defined. Furthermore, this operator satisfies the following properties.

Theorem 2. Let G be the operator defined at (17). G is a self-adjoint positive operator from H-1/2 (Γ f ) onto H1/2 (Γ f ) and there exist two constants c 1 > 0 and c 0 > 0 such that:

∀g ∈ H-1/2 (Γ f ), c 1 |||g||| -1/2 ≤ |||G(g)||| 1/2 ≤ c 0 |||g||| -1/2 , |||G(g)||| 2 1/2 ≤ g, G(g) .
Hence, the operator G is an isomorphism between H-1/2 (Γ f ) and H1/2 (Γ f ).

Proof. For (g, h) ∈ H-1/2 (Γ f ) 2 , let us set ψ = Φ(h) in (18). One obtains

Ω ∇Φ(g).∇Φ(h) = g, G(h) = h, G(g) (19)
which proves that G is self adjoint. Furthermore, (19) implies that

||Φ(g)|| 2 1,Ω ≤ |||g||| -1/2 |||G(g)||| 1/2 ,
and from the definition of the norm in H1/2 (Γ f ),

|||G(g)||| 2 1/2 ≤ ||Φ(g)|| 2 1,Ω ≤ |||g||| -1/2 |||G(g)||| 1/2
which proves the right hand side of the first inequality of theorem 2. The operator G maps continuously H-1/2 (Γ f ) into H1/2 (Γ f ) and it is clearly one to one because G(g) = 0 implies that Φ solution of (18) satisfies ∇Φ = 0 in Ω and also Φ = 0 on Γ 1 . Thus Φ = 0 and therefore g = ∂Φ ∂ν = 0.

The onto character of G is a consequence of its closed range. Let us prove it.

Let g ∈ H-1/2 (Γ f ) and let ψ 0 ∈ H1/2 (Γ f ) be such that

|||g||| -1/2 = g, ψ 0 |||ψ 0 ||| 1/2 . Let ψ ∈ W 0 whose trace on Γ f is ψ 0 and such that ||ψ|| W0 = |||ψ 0 ||| 1/2 . One has (with (19)) |||g||| -1/2 ≤ ||Φ(g)|| W0 ||ψ|| W0 |||ψ 0 ||| 1/2 ≤ ||Φ(g)|| W0 . (20)
Let us now consider a continuous extension operator from

H 1/2 (Γ) into X = {u ∈ H 1 (Ω), ∆u ∈ L 2 (Ω)} -say R-with (see J.L. Lions et E. Magenes [LM 68]): -∆R(v) ∈ L 2 (Ω) in Ω, ∂R(v) ∂ν = 0 on Γ, R(v) = v on Γ. Let v ∈ H1/2 (Γ f ).
From the definition of this space there exists w ∈ W 0 such that w |Γ f = v and |||v||| 1/2 = ||w|| W0 . Let v 1 = w |Γ . There exists c > 0 such that for any v ∈ H1/2 (Γ f ), one gets: 20) and ( 21) lead to the existence of a constant c > 0 such that (constant c might change from a line to another) for any

||v 1 || 1/2,Γ ≤ c|||v||| 1/2 . ( 21
) Applying this result with v = Φ(g)| Γ f = G(g), one obtains Φ(g) = R(v 1 ) + Φ 1 where ∆Φ 1 = -∆R(v 1 ) in Ω, Φ 1 = 0 on Γ f ∪ Γ 1 and ∂Φ 1 ∂ν = 0 on Γ 0 . One also has ||Φ(g)|| 1,Ω ≤ ||R(v 1 )|| 1,Ω + ||Φ 1 || 1,Ω ≤ ||v 1 || 1/2,Γ + ||∆R(v 1 )|| 0,Ω thus (
g ∈ H-1/2 (Γ f ), |||g||| -1/2 ≤ c||v 1 || 1/2,Γ ≤ c||v|| 1/2,Γ f ∪Γ1 ≤ c|||v||| 1/2 .
Recalling that v = G(g), one has proved that G has a closed range in H1/2 (Γ f ).Thanks to the fact that it is a self adjoint operator, one deduces that it is an isomorphism from H-1/2 (Γ f ) onto H1/2 (Γ f ).

Finally, let us prove the second relation in theorem 2. We have

c 1 |||g||| 2 -1/2 ≤ |||G(g)||| 2 1/2 ≤ c||Φ(g)|| 2 W0 = g, G(g) hence G is positive the proof of theorem 2 is complete. ✷
Let us now notice few properties of the operator L s defined on H 1/2 (Γ f ) by

L s (v) = ∇ s Φ0 .∇ s v. ( 22
) Since Φ0 ∈ C 1 ( Γf ), one has L s (v) ∈ H -1/2 (Γ f
) and there exists c > 0, depending on Φ0 such that

||L s (v)|| -1/2,Γ f ≤ c||v|| 1/2,Γ f . ( 23
)
As a consequence, one deduces :

Corollary 3. There exists c > 0 such that the following estimates hold : for every

u, v ∈ H 1/2 (Γ f ), | L s (v), G • L s (u) | ≤ c||v|| 1/2,Γ f ||u|| 1/2,Γ f and for every u, v ∈ H1/2 (Γ f ), | L s (v), G • L s (u) | ≤ c|||v||| 1/2 |||u||| 1/2 Proof. Since ∇ s Φ0 .∇ s v ∈ H -1/2 (Γ f ), one has ∇ s Φ0 .∇ s v ∈ H-1/2 (Γ f
) and the term G(∇ s Φ0 .∇ s v) makes sense. Applying Theorem 2, ( 16) and (15), it leads to

∀u, v ∈ H 1/2 (Γ f ), | L s (v), G • L s (u) | ≤ c 0 |||L s (v)||| -1/2 |||L s (u)||| -1/2, ≤ c||L s (v)|| -1/2,Γ f ||L s (u)|| -1/2,Γ f ≤ c||v|| 1/2,Γ f ||u|| 1/2,Γ f (24)
which ends the corollary's proof since the embbeding from H1/2 (Γ f ) in H 1/2 (Γ f ) is continuous and thus the second part follows. ✷

Notice that the kernel of L s plays an important role in the following.

The set of equations ( 4) and ( 10) is the dynamic fluid model. One can express ϕ in terms of η and its derivatives as follows, but η is also an unknown. In the case of boundary conditions on Γ f stated in (10), we have

ϕ| Γ f = G( ∂η ∂t ) + U G(L s (η))
whereas with boundary conditions (6), one has

ϕ| Γ f = G( ∂η ∂t ) + U G(L s (η)) + U G(∇ s η 0 .∇ s ϕ).
In this paper, we only consider boundary condition (10) on the free surface. The case of (6) would require a different definition of operator G. This will be studied in forthcoming papers.

An equation is still lacking: it is the one which traduces the equilibrium of the free surface under the combine effects of the gravity, the water pressure and the capillary forces. By eliminating ϕ as a function of η in this equilibrium equation, one obtains the surface model which is an up-dated version of Neuman-Kelvin's model. This is explicited in the next section.

4 The equilibrium of the free surface

p = p 0 -̺ ∂ϕ ∂t - ̺ 2 |∇ϕ| 2 -̺U ∇ s Φ0 .∇ s ϕ -̺ U 2 2 |∇ s Φ0 | 2 -̺gη + σ∆ s η. ( 25
)
In this expression σ > 0 is the capillary constant, ̺ > 0 is the mass density of the fluid and g is the gravity (let us recall that the notation ∆ s is the Laplace operator restricted to the free surface). One can linearize this formula by cancelling the term |∇ϕ| 2 . But one can also introduce a stationary solution η 0 named the digging effect. It corresponds a stationary solution for the velocity potential in the fluid denoted by ϕ 0 . Finally it is possible to linearize the non linear term ̺ 2 |∇ϕ| 2 around this new stationary terms. The linearized pressure at 0 is given by:

p = p 0 -̺ ∂ϕ ∂t -̺U ∇ s Φ0 .∇ s ϕ -̺ U 2 2 |∇ s Φ0 | 2 -̺gη + σ∆ s η (26)
Let us adopt for instance homogeneous Dirichlet boundary conditions for η on ∂Γ f (but many other possibilities exist). The equilibrium of the free surface is therefore governed by (p 0 which acts on both side of Γ f is eliminated):

-σ∆ s η + ̺gη + ̺ ∂ϕ ∂t + ̺U ∇ s Φ0 .∇ s ϕ = -̺ U 2 2 |∇ s Φ0 | 2 on Γ f . (27)
with :

ϕ = G( ∂η ∂t ) + U G(L s (η)) on Γ f . (28)
The term -̺ 2 U 2 |∇ s Φ0 | 2 is a Venturi effect which induces the digging effect. It means that even for a stationary case, one hasn't exactly η = 0. The solution of the previous model where the transcient terms have been cancelled is η 0 . Considering that both G and L s are linear, system ( 27) and ( 28) is linear and one has η = η 0 + η (and ϕ = ϕ 0 + φ) where η is the solution of ( 27)-( 28) with a null right hand side. In the next subsection, an analysis of the steady state of system ( 27)-( 28) is given and several concepts which will be useful in the following are introduced. Without loss of generality, the digging effect is omitted concerning the dynamic behavior of the surface.

On the digging effect

The digging component η 0 is solution of:

     η 0 ∈ H 1 0 (Γ f ), -σ∆ s η 0 + ̺gη 0 + ̺U 2 ∇ s Φ0 .∇ s G(∇ s Φ0 .∇ s η 0 ) = - ̺U 2 2 |∇ s Φ0 | 2 , (29)
where U G(∇ s φ0 .∇ s η 0 ) is the solution ϕ of the model (4)-(10) but where the term ∂η ∂t is cancelled. For example, if Φ0 = x 1 which is a classical case, and if one restricts the preceding model to the one dimensional case (see section 7), one obtains the following model for the digging effect:

         η 0 ∈ H 1 0 (Γ f ), -σ∆ s η 0 + ̺gη 0 + ̺U 2 ∂G( ∂η 0 ∂x 1 ) ∂x 1 = - ̺U 2 2 . ( 30
)
The expression of the digging effect η 0 depends on the steady flow which is described by Φ0 . For instance it could be much more meaningful for shallow water and with different boundary condition on Γ 1 . But the analysis given here-after still works.

In order to explain the difficulties which will be encountered in the tanscient analysis, let us discuss the model (29) assuming the properties of G stated in theorem 2 and in section 3. It can be written in a variational form as follows:

                   find η 0 ∈ H 1 0 (Γ f ) such that ∀v ∈ H 1 0 (Γ f ) : Γ f σ∇ s η 0 .∇ s v + ̺gη 0 v -̺U 2 G(∇ s Φ0 .∇ s η 0 )∇ s Φ0 .∇ s v -̺U 2 Γ f ∆ s Φ0 G(∇ s Φ0 .∇ s η 0 )v = - ̺U 2 2 Γ f |∇ s Φ0 | 2 v. (31)
We first prove the Theorem 4. There exists a critical velocity -say U c > 0-such that for every 0 ≤ U < U c , the bilinear form:

(η, v) ∈ H 1 0 (Γ f ) → σ Γ f ∇ s η.∇ s v + ̺g Γ f ηv -̺U 2 Γ f G(∇ s Φ0 .∇ s η)∇ s Φ0 .∇ s v -̺U 2 Γ f ∆ s Φ0 G(∇ s Φ 0 .∇ s η)v, is H 1 0 (Γ f ) elliptic. The operator: v ∈ H 1 0 (Γ f ) → -σ∆ s v + ̺gv + ̺U 2 ∇ s Φ 0 .∇ s G(∇ s Φ 0 .∇ s v),
satisfies the maximum principle (for 0 ≤ U < U c ) and one has η 0 ≤ 0. The precise characterization of U c is given in section 5.3.

Proof. From the properties satisfied by the operator G, one has (see theorem 2):

c 0 |||∇ s Φ0 .∇ s η||| -1/2 ≤ |||G(∇ s Φ0 .∇ s η)||| 1/2 ≤ c 1 |||∇ s Φ0 .∇ s η||| -1/2 .
On the one hand, with ( 16) and ( 23), (let us recall that Φ0 ∈ C 3 ( Γf )) one gets:

|||∇ s Φ0 .∇ s η||| -1/2 ≤ ||∇ s Φ0 .∇ s η|| -1/2,Γ f ≤ c 2 ||η|| 1/2,Γ f ,
and on the other hand:

| Γ f ∆ s Φ0 G(∇ s Φ0 .∇ s η)v| ≤ c 3 ||η|| 1/2,Γ f ||v|| 1/2,Γ f .
Hence, from the natural continuous inclusion of the space

H 1 0 (Γ f ) into H 1/2 (Γ f )
one can claim that there exists a constant c such that for any δ ∈]0, 1[:

                           ∀v ∈ H 1 0 (Γ f ), Γ f [σ|∇ s v| 2 + ̺gv 2 -̺U 2 G(∇ s Φ0 .∇ s v)∇ s Φ0 .∇ s v] -̺U 2 Γ f ∆ s Φ0 G(∇ s Φ0 .∇ s v)v ≥ (δ inf(σ, ̺g) -c̺U 2 )||v|| 2 1/2,Γ f + (1 -δ) inf(σ, ̺g)||v|| 2 1,Γ f . (32)
The inequality (32) ensures that for U small enough, the bilinear form which appears in the variational equation ( 31) is elliptic on the space H 1 0 (Γ f ) and (31) has a unique solution. It is worth noting that even if the bilinear form which appears in (31) is not symmetrical, the Lax-Milgram's theorem still applies. ✷ Let us extend the existence result to any velocity U with a Fredholm's alternative. Let us first assume that the term ∆ s Φ0 can be neglected which is the case in most applications (see remark 9 below for general cases).

Theorem 5. Let us assume that ∆ s Φ0 = 0 on Γ f . There exists a countable sequence {λ p } with

0 < λ 1 ≤ λ 2 ≤ λ 3 ≤ ...... ≤ λ p ≤ λ p+1 ≤ . . .
and there exists finite dimensional subspaces S p such that 1 If λ p < U 2 < λ p+1 . Then the model (31) has a unique solution.

2 If U 2 = λ p then • If: ∀v ∈ S p , Γ f |∇ s Φ0 | 2 v = 0
there is a solution to (31) which is defined up to an arbitrary element of S p .

• If : ∃v ∈ S p such that:

Γ f |∇ s Φ0 | 2 v = 0,
there is no solution to (31).

Proof. Let us make few preliminary remarks which will also be helpful in the the following. The bilinear form:

u, v ∈ H 1 0 (Γ f ) → a(u, v) = Γ f σ∇ s u.∇ s v + ̺guv,
is symmetrical continuous and coercive on the space H 1 0 (Γ f ). From theorem 2 and its corollary, the bilinear form:

u, v ∈ H 1 0 (Γ f ) → b(u, v) = ̺ Γ f G(∇ s Φ0 .∇ s u)∇ s Φ0 .∇ s v,
is symmetrical positive and continous on the space H 1 0 (Γ f ) and defines a seminorm equivalent to

||∇ s Φ0 .∇ s v|| -1/2,Γ f on this space. But it is not always a norm excepted if v ∈ H 1 0 (Γ f ), ∇ s Φ0 .∇ s v = 0 implies that v = 0. Therefore one sets: K 0 = {v ∈ H 1 0 (Γ f ), ∀q ∈ H 1 0 (Γ f ), b(v, q) = 0} (33)
which is the same than the kernel of b. It is a closed subspace of H 1 0 (Γ f ) and a partial description of it is given in remark 7. The mapping:

v ∈ H 1 0 (Γ f ) → b(v, v), (34)
is a norm on the quotient space H 1 0 (Γ f )/K 0 which is isomorphic to the orthogonal of K 0 in H 1 0 (Γ f ). The completed space with respect to this norm is denoted by H # . It is also an Hilbert space. From (24), H # contains the space H 1/2 (Γ f ) and thus H1/2 (Γ f ). Hence the embedding from H 1 0 (Γ f )/K 0 into H # is compact. Finally, the spectral theorem (see for instance [RATO 83]) can be applied to the following variational model:

                   find λ ∈ R + * , w ∈ H 1 0 (Γ f )/K 0 such that: ∀v ∈ H 1 0 (Γ f )/K 0 , a(w, v) = λb(w, v),
with (one can use the norm on the space H 1 0 (Γ f )/K 0 ):

inf v0∈K0 b(w + v 0 , w + v 0 ) = 1. ( 35
)
The result of the spectral theorem is stated in the following:

Lemma 6.
There is a countable family of solutions (λ n , w n ) ∈ R + * × H 1 0 (Γ f )/K 0 such that the multiplicity of each eigenvalue λ n is finite, the only accumulation point is ∞, the family of the eigenvectors {w n } is an Hilbert's basis in the space H # and { w n √ λ n } is an Hilbert's basis in the space H 1 0 (Γ f )/K 0 . Finally, the sequence λ n is ordered such that:

0 < λ 1 ≤ λ 2 ≤ λ 3 ≤ ...... ≤ λ p ≤ λ p+1 ≤ . . .

and tends to infinity.

Theorem 5 is a direct consequence of Lemma 6 and of the Fredholm's alternative.

Remark 7. The description of K 0 is clearly related to a unique continuation property. In particular, it is necessary to have information on the nature of the zero set of ∇ s Φ0 where Φ0 is solution of (1). Let γ ⊂ Γ f be the set

γ = {x ∈ Γ f , ∇ s Φ0 (x) = 0}.
It is a classical unique continuation result that γ has empty interior. Let us prove that in the two dimensional case, K 0 = {0}. We have : v ∈ K 0 implies that ∇ s v is null on the supplementary of γ. Since ∂ Φ0 ∂ν = 0 on Γ f , γ is the set where the gradient of Φ0 is null. Since Φ0 is a real and harmonic function, γ has no accumulation point in Γ f . There then exists a nonempty and open interval where v is constant. Let us consider a maximal interval I with this property and let us suppose that b = sup I ∈ ∂Γ f . Since Φ0 is regular, in both case ∇ s Φ0 (b) = 0 (and then b is an isolated point in γ) or

∇ s Φ0 (b) = 0, there exists α > 0 such that ]b, b + α[∩γ = ∅. This implies that v is constant on ]b, b + α[ and since v is continuous, v = v(b) and [b, b + α[⊂ I
which contracdicts the maximal character of I. One easily get that v is constant on Γ f and since v ∈ H 1 0 (Γ f ), v is null. In three dimensionnal case, let us notice that there are several examples where K 0 = {0}. The first case which is mostly important corresponds to Φ0 = x 1 + c (parallelepiped case). In this case one has:

v ∈ K 0 ⇔ ∂v ∂x 1 = 0,
which implies that v only depends on the coordinates x 2 transverse to x 1 in the plan containing Γ f . But the only function in H 1 0 (Γ f ) which is constant with respect to the coordinate x 1 and null on ∂Γ f is v = 0 and thus K 0 = {0} . More generally, as far as one can ensure that if for instance ∂ Φ0 ∂x 1 = 0 (excepted at isolated point) the result is still true. But this would be different if there was a surface boat on Γ f that we do not consider in this paper.

Remark 8. Let us consider the case where:

∃v ∈ S p such that:

Γ f |∇ s Φ0 | 2 v = 0.
When λ p < U 2 < λ p+1 there is a unique solution to (31). When U 2 → λ p the solution tends to infinity (in norm). One has a stiff problem which can be be interpreted physically by a an explosion of the progressive wave solution of (31).

Remark 9. In the case where ∆ s Φ0 = 0 on Γ f , one must consider a bilinear form b which is no more symetrical and therefore the spectral theory is much more complicated and not so well known. One can say that if U 2 is not an eigen value then the system is still well posed. If other cases, the right hand side should be conjugate with respect to the eigenvectors of the tranposed operators in order to ensure existence of solutions.

Remark 10. If the digging effect leads to a term η 0 which can't be neglected, it is necessary to modify the definition of G as follows. Instead of (18), let us introduce Φ ∈ W 0 solution of (G is the restriction of Φ on Γ f ) :

∀ψ ∈ W 0 , Ω ∇Φ.∇ψ - Γ f ∇ s η 0 .∇ s Φ ψ = Γ f gψ. ( 36
)
This is a non-symmetrical variational model. The only point which is not straightforwards is that the bilinear form implied in the preceding model is elliptic. In fact one has:

Ω |∇Φ| 2 - Γ f ∇ s η 0 .∇ s Φ Φ = ||Φ|| 2 1,Ω - 1 2 ∂Γ f ∂η 0 ∂ν Φ 2 + Γ f ∆ s η 0 Φ 2 . ( 37
)
But Φ = 0 on ∂Γ f because Φ = 0 on Γ 1 . The bilinear form implied in (36) is therefore strictly coercive on W 0 as soon as ∆ s η 0 ≥ 0 or is small enough.

5 Formulation of the model in η without digging effect

The model

In this section, let us consider dynamical part of system ( 27) and (28) neglecting η 0 but having in mind that it should be added to the solution η that is discussed hereafter in a numerical simulation. Let us consider the functional space in which η is looked for:

V = H 1 0 (Γ f ). ( 38
)
By multiplying equation ( 27) by an arbitrary element of the space V , one obtains:

∀v ∈ V, σ Γ f ∇ s η.∇ s v+̺g Γ f ηv+̺ Γ f ∂ϕ ∂t v+̺U Γ f ∇ s Φ0 .∇ s ϕv = 0. ( 39
)
Using the map G which has been introduced in section 3, one can eliminate the velocity potential ϕ. When U = 0 it leads to an inertia term on the free surface of the water. But when U = 0 it gives new terms including gyroscopic and hydrodynamic negative stiffness contributions. Let us explains how. Let us recall that the potential ϕ is solution of:

       ϕ ∈ W 0 , ∀ψ ∈ W 0 , Ω ∇ϕ.∇ψ = Γ f ∂η ∂t ψ + U ∇ s Φ0 .∇ s η ψ, (40)
which can be written (28) as

ϕ |Γ f = G( ∂η ∂t ) + U G(∇ s η.∇ s Φ0 ). ( 41
)
By introducing this expression into the equation (39) which traduces the equilibrium of the free surface, one obtains formally:

     ∀t ∈ [0, T ], η = η(t) ∈ H 1 0 (Γ f ), ∀v ∈ H 1 0 (Γ f ), m s ( ∂ 2 η ∂t 2 , v) + 2U c( ∂η ∂t , v) + a(η, v) = 0, (42)
with the notations for η and v smooth enough :

                                   m s (η, v) = ̺ Γ f G(η)v, c(η, v) = ̺ 2 Γ f [G(∇ s Φ0 .∇ s η)v -∇ s Φ0 .∇ s v G(η) -∆ s Φ0 G(η)v], a(η, v) = Γ f [σ∇ s η.∇ s v + ̺gηv] -̺U 2 Γ f [G(∇ s Φ0 .∇ s η)∇ s Φ0 .∇ s v + ∆ s Φ0 G(∇ s Φ0 .∇ s η)v]. (43)
Of course, one must add given initial data (η 0 , η 1 ) to equation ( 42).

Let us present some properties of the bilinear forms involved by (42) which will be usefull in what follows. Let us set a = a s + a r and c = c s + c r with

                               c s (η, v) = ̺ 2 Γ f [G(∇ s Φ0 .∇ s η)v -∇ s Φ0 .∇ s v G(η)] c r (η, v) = - ̺ 2 Γ f ∆ s Φ0 G(η)v, a s (η, v) = Γ f [σ∇ s η.∇ s v + ̺gηv] -̺U 2 Γ f G(∇ s Φ0 .∇ s η)∇ s Φ0 .∇ s v a r (η, v) = -̺U 2 Γ f ∆ s Φ0 G(∇ s Φ0 .∇ s η)v.
(44) One has a = a s and c = c s when the term ∆ s Φ0 can be neglectible which is the case in many examples. The first step is to check few useful properties of the bilinear forms m s , c and a. Let us notice that they are well defined if η and v are functions lying in the space H 1 0 (Γ f ). The bilinear form m s and a s are symmetrical and the bilinear form c s is anti-symmetrical. This is due to the fact that it represents a gyroscopic coupling, in fact, ∇ s η is a rotation of the unit normal ν. The main properties are the coerciveness properties which are summarized hereafter. Nevertheless the additional terms which would appear in the expression of c and a can be easily handled as soon as one assumes that Φ0 ∈ C 2 (Γ f ), (see P. Grisvard [GRI 86] or M. Borsuk and V. Kondratiev [KON 06] for the justification) because one has in this case:

                   |c r (η, v)| ≤ c|||η||| -1/2,Γ f |||v||| -1/2,Γ f , and: |a r (η, v)| ≤ c|||v||| 1/2,Γ f ||∇ s Φ0 .∇ s η|| -1/2,Γ f ≤ c ′ ||η|| 1/2 |||v||| 1/2,Γ f (45)
Theorem 11. Let us assume that min(̺g, σ) > 0. There exists constants c i > 0 such that for every η and v in H 1 0 (Γ f ), one has (i) Concerning a and a s :

|a s (η, v)| + |a(η, v)| ≤ c 3 ||η|| 1,Γ f ||v|| 1,Γ f a(v, v) ≥ c 0 ||v|| 2 1,Γ f -c 6 U 2 ||v|| 2 1/2,Γ f (ii) Concerning m s : c 1 |||v||| 2 -1/2 ≤ m s (v, v) ≤ c 2 ||v|| 2 -1/2,Γ f ≤ c 2 |||v||| 2 -1/2 . (iii) Concerning c and c s : |c s (η, v)| + |c(η, v)| ≤ c 7 ||η|| 1/2,Γ f ||v|| 1/2,Γ f c s (v, v) = 0, ∀v ∈ H 1 0 (Γ f ).
Proof. Notice first that H 1 0 (Γ f ) ⊂ H1/2 (Γ f ) thus Theorem 9 makes sense. Theorem 2 proves that the bilinear form m s induces a scalar product on H-1/2 (Γ f ) which is equivalent to the usual one. Furthermore, the continuity of m s with respect to the H -1/2 (Γ f ) norm is a consequence of ( 16) and thus assertion (ii) is proved. Concerning a s , let us notice that:

∀η ∈ H 1 0 (Γ f ), a s (η, η) = σ Γ f |∇ s η| 2 + ̺g Γ f η 2 -̺U 2 Γ f G(∇ s Φ0 .∇ s η)∇ s Φ0 .∇ s η ≥ min(σ, ̺g)||η|| 2 1,Γ f -̺U 2 Γ f G(∇ s Φ0 .∇ s η)∇ s Φ0 .∇ s η.
But, from the definition of the operator G and using again theorem 2 and its corollary :

Γ f G(∇ s Φ0 .∇ s η)∇ s Φ0 .∇ s η ≤ c ′ 0 ||η|| 2 1/2,Γ f
and assertion (i) is proved. Assertion (ii) is clear with (i) and the second part of (45). Finally, (iii) is an easy consequence of (23). The theorem 11 is now proved. ✷ Remark 12. A natural question is to obtain the best constants which appears in the previous theorem. This is discussed in a section 5.3 in the case where ∆ s Φ0 = 0. .

Well-posedness of (42) (σ > 0) for small velocity

Let us prove the following result:

Theorem 13. There exists U c > 0 such that for every 0 ≤ U < U c and every (η 0 , η 1 ) ∈ V × H-1/2 (Γ f ), there exists a unique solution to (42) in

C 0 ([0, T ]; V ) × C 1 ([0, T ]; H-1/2 (Γ f )) with η(0) = η 0 and η(0) = η 1 .
Proof. Let us consider a basis -say {w n }of the space H 1 0 (Γ f ) and let us denote by V N the finite dimensional subspace of H 1 0 (Γ f ) spanned by the N first basis vectors w n . Furthermore η N 0 and η N 1 are approximations of η 0 = η(0) and η 1 = ∂η ∂t (0) in the space V N . They are such that :

     lim N →∞ ||η N 0 -η 0 || 1,Γ f = 0, lim N →∞ |||η N 1 -η 1 ||| -1/2,Γ f = 0. ( 46
)
Let us now consider the approximate models (for any N ):

                 find η N ∈ V N such that: ∀v ∈ V N , m s ( ∂ 2 η N ∂t 2 , v) + 2U c( ∂η N ∂t , v) + a(η N , v) = 0, η N (0) = η N 0 , ∂η N ∂t (0) = η N 1 . (47)
This is a finite dimensional linear differential equation with initial conditions and constant coefficients. Therefore there is a unique solution to (46). Furthermore the time dependence is as smooth as one wishes because one has:

η N = k=1,N α N k (t)w k ,
and the time derivatives only imply the coefficients α N k . Let us set U c = c 0 c 6 where the constants c 0 and c 6 are introduced in theorem 11. Let us assume that U < U c . The bilinear form a is coercive on the space H 1 0 (Γ f ) (non symmetrical) and one obtains an a priori estimate for the approximate solution by setting v = ηN in the variational equation (47). From the symmetry of the bilinear forms m s and a s , the anti-symmetry of c s , one obtain with

E(t) = 1 2 [m s ( ∂η ∂t , ∂η ∂t ) + a s (η, η)](t), if ∆ s Φ0 = 0, E(t) = E(0) (48)
and in others cases,

E(t) ≤ c[E(0) + t 0 (c r (η, ∂η ∂t ) + a r (η, ∂η ∂t ))]. (49)
Thanks to (45), Theorem 11 and Gronwall's Lemma, the previous estimate enables one to extract from η N , a subsequence denoted by η N ′ which converges to η * in the following topology:

       η N ′ ⇀ η * in L ∞ (]0, T [; H 1 0 (Γ f )) weak*, ∂η N ′ ∂t ⇀ ∂η * ∂t in L ∞ (]0, T [; H-1/2 (Γ f )) weak *. (50)
In order to take the limit for N ′ → ∞, one uses the equivalent variational formulation which makes sense to the second order time derivative of η * :

                 ∀Ψ = k=1,N0 β k (t)w k , with β k (t) ∈ D([0, T [), N 0 ≤ N, -m s ( ∂η N ∂t , ∂Ψ ∂t ) -2U c(η N , ∂Ψ ∂t ) + a(η N , Ψ) = -{m s (η N 1 , Ψ) + 2U c(η N 0 , Ψ)}. ( 51
)
But the result is true for any N 0 ≤ N and the approximate initial conditions converge to the one of the continuous model; hence one can claim that η * is a solution to the weak (variational) model ( 42). The last point concerns the uniqueness. The proof is slightly different from the usual one as far as there is a gyroscopic term. First of all let us point out that the result is straightforwards for smooth enough solutions. If m s ( ∂ 2 η ∂t 2 , ∂η ∂t ) and a s (η, ∂η ∂t ) can be defined, the result is derived from the energy conservation property. In fact one can introduce the difference between two solutions which satisfy the same initial conditions and therefore the energy -which is constant with respect to the time variable-is always zero. Let us now consider the general case for which the regularity condition of the weak solution is not necessarily satisfied. Let us set (following an idea given in J.L Lions [JLL 69]):

h = t 0 η (52)
One can check that:

             h ∈ L ∞ (]0, T [; H 1 0 (Γ f ) ∩ W 1,∞ (]0, T [; H 1/2 (Γ f )), ∀v ∈ H 1 0 (Γ f ), m s ( ∂ 2 h ∂t 2 , v) + 2U c( ∂h ∂t , v) + a(h, v) = m s (η 1 , v) + 2U c(η 0 , v). ( 53
)
But h has also to satisfy the initial conditions:

h(0) = 0 et ∂h ∂t = η 0 . ( 54
)
This enables one to ensure the uniqueness of h because the regularity assumption is now satisfied. The uniqueness of η is a consequence of the one of h.

It has been proved that for small velocity U < U c , the system is well-posed. ✷

The value of U c given by Theorem (11) is not optimal and we precise it in the next section.

About the static instabilities

Let us discuss the best value for the constant U c in the particular case where the term ∆ s Φ0 is neglected. One has in this case a = a s . Two possibilities have to be considered depending if there is a capillary term or not (σ > 0 and σ = 0). Let us introduce the eigenvalue model which has already been used in the analysis of the digging effect (31):

               find (λ, w) ∈ R + × H 1 0 (Γ f ) such that: ∀v ∈ H 1 0 (Γ f ), λ̺ Γ f G(∇ s Φ0 .∇ s w)∇ s φ0 .∇ s v = σ Γ f ∇ s w.∇ s v + ̺g Γ f wv. (55)
It is a non-classical formulation because the bilinear form:

(η, v) ∈ H 1 0 (Γ f ), → ̺ Γ f G(∇ s Φ0 .∇ s w)∇ s φ0 .∇ s v, (56)
is not always definite. One has from corollary 2:

       ∃c > 0, s.t. ∀η ∈ H 1 0 (Γ f ) : Γ f G(∇ s Φ 0 .∇ s η)∇ s Φ 0 .∇ s η ≥ c|||∇ s Φ0 .∇ s η||| 2 -1/2 . ( 57
)
The kernel of the bilinear form (56) has been defined at (33) as the subspace of H 1 0 (Γ f ) such that:

K 0 = {v ∈ H 1 0 (Γ f ), ∇ s Φ0 .∇ s v = 0 on Γ f }. ( 58
)
It is a closed subspace of H 1 0 (Γ f ) and few properties of it are given in the remark 7 where it has been proved that in 2D and in many cases in 3D, one has K 0 = {0}. In what follows and in order to simplify, let us assume that K 0 = {0} which is the most realistic case. Notice that our results can be easily adapted to the case where K 0 = {0} using the quotient space H 1 0 (Γ f )/K 0 as in subsection 4.2.

Since K 0 = {0}, the mapping:

η ∈ H1/2 (Γ f ) → |||η||| # = |||∇ s Φ0 .∇ s η||| -1/2,Γ f , (59)
is a norm on the space H1/2 (Γ f ) and with Theorem 2 :

       ∃c, c ′ > 0, s.t. ∀η ∈ H1/2 (Γ f ) : c|||η||| 2 # ≤ Γ f G(∇ s Φ 0 .∇ s η)∇ s Φ 0 .∇ s η ≤ c ′ |||η||| 2 # . ( 60
)
Let H # (Γ f ) be the completed space of H1/2 (Γ f ) equiped with the norm ||| . ||| # . Using ( 16) and ( 23) then (15), we get the existence of a constant c > 0 such that for every v ∈ H1/2 (Γ f ),

|||v||| # ≤ c||L s (v)|| -1/2,Γ f ≤ c||v|| 1/2,Γ f ≤ c|||v||| 1/2,Γ f , thus H1/2 (Γ f ) ⊂ H # (Γ f ).
Remark 14. The space H # (Γ f ) is larger than H1/2 (Γ f ). An interesting case (see theorem 15) corresponds to the inclusion:

H # (Γ f ) ⊂ L 2 (Γ f ).
Let us first consider the two dimensionnal case with Φ0 = x 1 + c. One has:

∇ s Φ0 .∇ s v = ∂v ∂x 1 ∈ H-1/2 (Γ f ) hence v ∈ H 1/2 (Γ f ) ⊂ L 2 (Γ f ). Therefore H # (Γ f ) ⊂ L 2 (Γ f ). In the three dimensionnal case, ∂v ∂x 1 ∈ H-1/2 (Γ f )
and there is no reason to hope the previous inclusion.

From general spectral theory of linear operators one can deduce the following theorem:

Theorem 15. Let us assume that K 0 = {0}. Two cases are considered in which the results are inverted. 1 • If σ > 0, there exists a countable set of elements denoted by (λ n , w n ) in R + * × H 1 0 (Γ f ) solution of (55). Each term of the sequence (ordered by increasing values) λ n has a finite multiplicity and +∞ is the only accumulation point. The family {w n } is an Hilbert basis of the space H # (Γ f ). Furthermore, the family { w n √ λ n } is an Hilbert basis of the space H 1 0 (Γ f ). The smallest eigenvalue denoted by λ 1 will be the best constant in theorem 11.

2 • If σ = 0 the result is different. Let us assume that H # (Γ f ) ⊂ L 2 (Γ f )
and that the inclusion is compact. There exists a countable set of solutions (λ n , w n ) ∈ R + * × H # (Γ f ) and the terms of the sequence λ n are ordered by decreasing values (the largest value is denoted by λ 1 ). The multiplicity of each term is finite and 0 is the only accumulation point. The family {w n } is an Hilbert basis of the space L 2 (Γ f ) while the family { w n √ λ n } is an Hilbert basis of the space H # (Γ f ).

Proof. Let us begin with σ > 0. The result is a direct consequence of the spectral theory for linear operator as it is presented P.A. Raviart and J.M. Thomas [RATO 83]. In fact the bilinear form :

(η, v) → σ Γ f ∇ s η.∇ s v + ̺g Γ f ηv,
is symmetrical, continuous and coercive on the space H 1 0 (Γ f ). Furthermore the bilinear and symmetrical form:

(η, v) → Γ f G(∇ s Φ 0 .∇ s η)∇ s Φ 0 .∇ s η,
is continuous and elliptic on the space H # (Γ f ). The assumptions of the spectral theory formulated in [RATO 83] are satisfied and the result is proved.

Let us consider the second case: σ = 0. Now the bilinear and symmetrical form:

(η, v) → ̺g Γj ηv,
is only equivalent to the square of the norm in the space L 2 (Γ f ). The variational eigenvalue problem is inverted. The compact embedding is from the space H # (Γ f ) into L 2 (Γ f ). This completes the proof of theorem 15. ✷ Remark 16. In the first case where σ > 0, one can characterize the constant λ 1 which appears in theorem 11 as the smallest eigenvalue of the model (55). Therefore there are no confusion in the notation. Let us underline that λ 1 which is the smallest eigenvalue solution of (55) satisfies for every v ∈ H 1 0 (Γ f ) :

λ 1 ̺ Γ f G(∇ s Φ 0 .∇ s v)∇ s Φ 0 .∇ s v ≤ σ Γ f |∇ s v| 2 + ̺g Γ f v 2 . ( 61
)
Hence:

       ∀δ ∈ [0, 1], ∀v ∈ H 1 0 (Γ f ), a s (v, v) ≥ inf(σ, ̺g)(1 -δ)||v|| 2 1,Γ f +̺(δλ 1 -U 2 ) Γ f G(∇ s Φ 0 .∇ s v)∇ s Φ 0 .∇ s v. (62)
Thus λ 1 is the best constant in theorem 11 and it enables to characterize the critical velocity U c corresponding to the instability of the linear wave model in which the term ∆ s Φ0 is neglected, by:

U c = λ 1 . (63)
If U ≥ U c the coerciveness of the stiffness bilinear form on the space H 1 0 (Γ f ) is lost. But it is still true for the complementary of the first eigenvectors in this space.

Remark 17. Let us assume that σ = 0, the coerciveness of a s is no more true on the space H 1 0 (Γ f ) but only on the space L 2 (Γ f ) which is not contained in H # (Γ f ) in most cases. In fact, the contrary is very often true which is a so-called mathematical inversion of the compactness. Furthermore, for any velocity U > 0, there exists an infinite number of instable eigenmodes. They are more and more local (the wave length is smaller and smaller). Therefore, in a numerical approximation, the more the mesh is refined, the larger is the number of instabilities. This remark which has already formulated in a different numerical framework by Xiao-bo Chen [CHEN 02], condemns the Neuman-Kelvin without capillary. In fact, the variational model (42) is fully instable excepted for a finite number of eigenmodes (see theorem 15) the eigenvalues of which -say λ n solution of (55) with σ = 0-would satisfy:

λ n ≥ U 2 . ( 64
)
One can draw a strange conclusion: if U is small enough and if the mesh size in a numerical approximation is large enough, the Neumann-Kelvin is stable. This is really disturbing because the true model is not. One has a numerical filtering of the instabilities due to mesh size which would be too large.

The conclusion of this section is that the Neumann-Kelvin model is non physical (instable) as far as the capillary is not taken into account (excepted if U = 0). This result is obviously in contradiction with the usual argument that the capillary is very small and can be neglected compared to the gravity effect. A simple analytical example is discussed in section 6.

Extension of the existence and uniqueness result when

U ≥ U c = √ λ 1 and σ > 0
First of all, let us note that the assumption σ > 0 is essential in this section where the well-posedness of the system for large velocity is studied. For sake of simplicity, it is assumed that ∆ s Φ0 = 0 is neglectible. For a given velocity U, there is only a finite number of instable eigenmodes for which the stiffness operator (the static one) is negative. All the norms on the finite dimensional space spanned by these eigenmodes are equivalent. This enables one to overcome the lack of postivity. But it would be false for σ = 0 for which the instable space is not finite dimensional. Our goal in this section is to prove an existence result for the instable wave model.

Theorem 18. For every (η 0 , η 1 ) ∈ V × H-1/2 (Γ f ), there exists a unique solution to (42

) in C 0 ([0, T ]; V ) × C 1 ([0, T ]; H-1/2 (Γ f )) with η(0) = η 0 and η(0) = η 1 .
Proof. Let κ be a strictly positive constant which is specified latter on depending on the largest eigenvalue λ n solution of the model (55) and such that:

λ n ≤ U 2 < λ n+1 .
Let us set:

η = e κt z. (65)
A simple formal calculus enables one to characterize z as a formal solution of:

                 ∀v ∈ H 1 0 (Γ f ) ∩ L 2 0 (Γ f ), m s ( ∂ 2 z ∂t 2 , v) + 2U c s ( ∂z ∂t , v) + a s (z, v) + κ 2 m s (z, v) = -2κm s ( ∂z ∂t , v) -2U κc s (z, v). (66)
In order to obtain an a priori estimate on the solution z N of the approximate model (analogous to (47)) let us choose v = ∂z N ∂t for the test function in the variational formulation. This leads to:

1 2 d dt [m s ( ∂z N ∂t , ∂z N ∂t ) + a s (z N , z N ) +κ 2 m s (z N , z N )] ≤ c||| ∂z N ∂t ||| 2 -1/2,Γ f + 2U κc s (z N , ∂z N ∂t )
Let us now examine the term c s (z N , ∂z N ∂t ). It is the summ of two terms. The first one is:

̺ Γ f G(∇ s Φ0 .∇ s z N ) ∂z N ∂t ,
which is bounded as follows: (as usual with (15), ( 16) and Theorem 2)

̺| Γ f G(∇ s Φ0 .∇ s z N ) ∂z N ∂t | ≤ c 0 |||G(∇ s Φ0 .∇ s z N )||| 1/2,Γ f ||| ∂z N ∂t ||| -1/2,Γ f ≤ c 1 |||∇ s Φ0 .∇ s z N ||| -1/2,Γ f ||| ∂z N ∂t ||| -1/2,Γ f ≤ c 1 ||| ∂z N ∂t ||| -1/2,Γ f ||∇ s Φ0 .∇ s z N || -1/2,Γ f ≤ c 1 2 [||z N || 2 1/2,Γ f + ||| ∂z N ∂t ||| 2 -1/2,Γ f ].
The second term is:

-̺ Γ f G(z N )∇ s Φ0 .∇ s ( ∂z N ∂t ).
It can written as follows:

- d dt [̺ Γ f G(z N )∇ s Φ0 .∇ s z N ] + ̺ Γ f G( ∂z N ∂t )∇ s Φ0 .∇ s z N .
Let us summarize the previous results. One has:

d dt [m s ( ∂z N ∂t , ∂z N ∂t ) + a s (z N , z N ) + κ 2 m s (z N , z N )+ +2U κ̺ Γ f G(z N )∇ s Φ0 .∇ s z N ] ≤ c 2 [||z N || 2 1,Γ f + ||| ∂z N ∂t ||| 2 -1/2,Γ f ]. (67)
But one has also for any α > 0:

2U κ̺ Γ f G(z N )∇ s Φ0 .∇ s z N ≤ 2U κ̺||G(z N )|| 1/2,Γ f ||∇ s Φ0 .∇ s z N || -1/2,Γ f ≤ ̺c 3 [ακ 2 ||z N || 2 -1/2,Γ f + U 2 α ||z N || 2 1/2,Γ f ],
which enables one to derive the estimates required for the a priori estimate on z N for κ large enough. Let us explain how. From the previous inequality, one has (the norm on H-1/2 (Γ f ) can be chosen equal to m s (. .) or equivalently to ||| ||| -1/2,Γ f ):

a s (z N , z N ) + κ 2 m s (z N , z N ) + 2U κ̺ Γ f G(z N )∇ s Φ0 .∇ s z N ≥ κ 2 (1 -αc 3 ̺)|||z N ||| 2 -1/2,Γ f + a s (z N , z N ) - ̺c 3 U 2 α ||z N || 2 1/2,Γ f
In order to lower bound the previous quantity, let us consider a splitting of the space H 1 0 (Γ f ) into two subspaces as describe in the following lemma (let us denote by (( , )) the natural scalar product on H 1/2 (Γ f )).

+ A ζ (w, v) is coercive on H 1 0 (Γ f ).
Proof of Lemma 19. Applying the spectral theorem to the eigenvalue problem

a s ( w, v) = μ(( w, v)),
there exists an Hilbert basis of eigenvectors wk (k > 0) corresponding to a non negative increasing sequence of eigenvalues μk which tends to the infinity. The bilinear form A ζ is a non negative one except on a finite dimensionnal subspace of H 1 0 (Γ f ) denoted by E ζ . On that space E ζ , norms m s and || || 1/2,Γ f are equivalent. For some c 0 > 0,

∀v ∈ E ζ , δ ζ m s (v, v) + A ζ (v, v) ≥ c 0 δ ζ ((v, v)) + A ζ (v, v) ≥ a s (v, v).
Therefore, there exists

δ ζ > 0 such that δ ζ m s (v, v) + A ζ (v, v) is coercive on H 1 0 (Γ f ) which proves Lemma 19. ✷ Let ζ > 0 be fixed. The classical spectral theorem can be applied to the eigenvalue problem δ ζ m s (w, v) + A ζ (w, v) = µ m s (w, v) which is equivalent to:    find µ ∈ R, w ∈ H 1 0 (Γ f ) such that: ∀v ∈ H 1 0 (Γ f ), µ m s (w, v) = a s (w, v) -ζ((w, v)). (68)
There exists a non decreasing sequence {µ k } in R with only the infinity as accumulation point and an Hilbert basis {w k } in H-1/2 (Γ f ) equipped with norm √ m s of eigenvectors solution of (68) which are orthogonal for A ξ . Furthermore, it has been proved in the previous lemma that the space spanned by the eigenvectors w k , k = 1, P corresponding to negative eigenvalues of the model (68) is a finite dimensional one. It is denoted by

V P . Let us choose ζ = ̺ c 3 U 2 α
. This enables one to derive the following lower bound on the space V P :

∀v ∈ V P , κ 2 (1 -αc 3 ̺)m s (v, v) + a s (v, v) - ̺c 3 U 2 α ||v|| 2 1/2,Γ f ≥ [κ 2 (1 -αc 3 ̺) + µ 1 ]m s (v, v).
But on the finite dimensional space V P all the norms are equivalent. Thus choosing

α < 1 c 3 ρ and κ 2 > - µ 1 1 -αc 3 ̺ > 0,
one can conclude that there exists a positive constant cP such that:

∀v ∈ V P , κ 2 (1 -αc 3 ̺)m s (v, v) + a s (v, v) - ̺c 3 U 2 α ||v|| 2 1/2,Γ f ≥ cP ||v|| 2 1,Γ f . (69)
On the orthogonal complementary space V ⊥ P of V P in H 1 0 (Γ f ) equipped with the scalar product a 0 , the previous inequality is also satisfied and there exists d P > 0 such that :

∀v ∈ V ⊥ P , κ 2 (1 -αc 3 ̺)m s (v, v) + a s (v, v) - ̺c 3 U 2 α ||v|| 2 1/2,Γ f ≥ d P ||v|| 2 1,Γ f . (70)
Since the basis {w k } is orthogonal for the entire bilinear form

κ 2 (1-αc 3 ̺)m s (v, v)+ a s (v, v) - ̺c 3 U 2 α ||v|| 2 1/2
,Γ f , one easily gets with c P = min(c P , d P ) > 0,

∀v ∈ V P , κ 2 (1 -αc 3 ̺)m s (v, v) + a s (v, v) - ̺c 3 U 2 α ||v|| 2 1/2,Γ f ≥ c P ||v|| 2 1,Γ f (71)
and thus there exists c P > 0 and κ > 0 such that :

∀v ∈ V P , κ 2 m s (v, v) + a s (v, v) + 2ρU Γ f G(v)∇ s Φ0 .∇ s v ≥ c P ||v|| 2 1,Γ f (72)
Finally, it has been proved that for any velocity U , there exists a constant κ ≥ 0 and another one c P such that :

m s ( ∂z N ∂t , ∂z N ∂t ) + a s (z N , z N ) + κ 2 m s (z N , z N ) +2U κ̺ Γ f G(z N )∇ s Φ0 .∇ s z N ] ≥ c P [||| żN ||| 2 -1/2,Γ F + ||z N || 2 1,Γ f ]. ( 73
)
The a priori estimate on z N is then derived from the Gronwall's lemma (see for instance [JLL 69]) on the inequality (67). It is obtained in the space

L ∞ (]0, T [; H 1 0 (Γ f )) ∩ W 1, ∞ (]0, T [; H-1/2 (Γ f )),
and one can conclude to the existence and the uniqueness of a solution to the instable model exactly as in the case where U < √ λ 1 . The theorem 18 is proved. ✷ Remark 20. One can consider several kind of boundary conditions on ∂Γ f concerning η. For instance, one can set:

∂η ∂ν + αη = 0 sur ∂Γ f , (74)
which is a Robin's boundary condition. If α = 0 there is a perfect gliding of the water, if α = ∞ there is an adhesion of the water. But one can also consider a condition as:

∂η ∂t + γ ∂η ∂ν + δη = 0, (75)
and so on. The most suitable boundary condition for η is certainly a transparency one (L Halpern [Hal 06]) in order to avoid reflection. But it would also be useful to use such a condition for ϕ on Γ 1 . This will be discussed in a future work.

The model formulated with ϕ

The elimination of η is quite easy from the theoretical point of view, but much more complicated in the practical applications as far as the capillary is kept. In fact the goal is to transfom a Lagrangian formulation into an Euler one using advective derivation:

D t (.) = ∂ ∂t (.) + U ∇ s Φ0 .∇ s (.)
Thus, one has on Γ f ×]0, T [:

                                
on the one hand:

∂ϕ ∂ν = D t η,
and on the other hand:

-σ∆ s D t η + g̺D t η = -[̺ ∂ 2 ϕ ∂t 2 + 2̺U ∇ s Φ0 .∇ s ∂ϕ ∂t + ̺U 2 ∇ s Φ0 .∇ s (∇ s Φ0 .∇ s ϕ)]. (76)
Furthermore, ϕ satisfies the relations (4) in Ω and on Γ 0 ∪ Γ 1 .

Let us begin with the usual -but physically wrong-case where σ = 0 (no capillary term). Let us consider ψ ∈ W 0 and multiplying (4) by ψ, and integrating by parts, one obtains (let us recall that ϕ = 0 on ∂Γ f and therefore boundary terms disappear in an integration by parts):

8 > > > > > > > > > > > < > > > > > > > > > > > : ∀ψ ∈ W0, Z Ω ∇ϕ.∇ψ + 1 g Z Γ f ∂ 2 ϕ ∂t 2 ψ +U Z Γ f ∇s Φ0.∇s( ∂ϕ ∂t )ψ -(∇s Φ0.∇sψ) ∂ϕ ∂t -∆s Φ0 ∂ϕ ∂t ψ -U 2 Z Γ f (∇sΦ0.∇sϕ)(∇s Φ0.∇sψ) -U 2 Z Γ f ∆s Φ0(∇s Φ0.∇sϕ)ψ = 0. ( 77
)
This variational equation which corresponds to the local equation ( 4) is the wellknown Neumann-Kelvin model [STO 57]. Let us point out again that its numerical implementation is fully instable excepted for U = 0.

Let us now go to the case σ > 0. Let Z 0 be the distribution space on Γ f defined by (it would be different for other boundary condition on ∂Γ f ):

Z 0 = H -1 (Ω). (78)
Let us define by A the operator which maps Z 0 onto H 1 0 (Γ f ) defined by:

A = [- σ ̺ ∆ s + gI d ] -1 . (79)
This is a regularizing operator which enables to define the up-graded (with capillary) Neumann-Kelvin variational model by:

8 > > > > < > > > > : ∀ψ ∈ W0, Z Ω ∇ϕ.∇ψ + 1 g Z Γ f A[ ∂ 2 ϕ ∂t 2 ]ψ +2U Z Γ f A[∇sΦ0.∇s( ∂ϕ ∂t )]ψ + U 2 Z Γ f A[(∇sΦ0.∇s(∇sΦ0.∇sϕ)]ψ = 0. (80)
Remark 21. The mathematical analysis of Neumann-Kelvin model upgraded with the capillary term is not necessary because it is equivalent to the one in η. But the numerical approximation is easier for the model in η because the order of the tangential derivatives are smaller.

A simple analytical example

                         -∆ϕ = 0 in Ω, ∂ϕ ∂ν = 0 on Γ 0 ∪ Γ 1 , -σ∆ s η + ̺gη + λ̺ ∂ϕ ∂x 1 = 0 and ∂ϕ ∂ν = ∂η ∂x 1 on Γ f , η = 0 on ∂Γ f . (81)
One can solve it using a separation of the variables x 1 et x 2 (see figure 2). Thus, one obtains a sequence of solutions depending on the integer n:

                             η n (x 1 , x 2 ) = 2 L sin( nπx 1 L ), ϕ(x 1 , x 2 ) = 2 L cos( nπx 1 L ) coth( nπ(x 2 + H) L ), λ n = [gH + n 2 π 2 σH L 2 ̺ ] tanh( nπH L ) nπH L H large ≃ g k + σk ̺ , k = nπ L . (82)
Let us introduce the minimizer

x 0 = g̺ σ L π of the function f (x) = gL πx + σπx ̺L .
Let us define the integer [x 0 ] as the entire part of x 0 . One can easily prove that

U c = min( f ([x 0 ]), f ([x 0 ] + 1) =    f ([x 0 ]) if x 0 ≤ [x 0 ] 2 + [x 0 ] f ([x 0 ] + 1) if x 0 > [x 0 ] 2 + [x 0 ] (83)
In the same context than J. Lighthill, we obtain a value close to (corresponding

to f (x 0 )) 2g k = √ 2( gσ ̺ ) 1/4
, which is a value given by J. Lighthill [LIH 78] but from a different formulation (using a Fourier's transform in the x 1 direction).

One can find the same expression which is discussed in the paper by D. Richard and de P. G. Gennes [RIC 96]. Their framework is also the one of J. Lighthill.

For instance, for g = 9.81 m/s 2 , σ = 7, 5 10 -2 N/m, ̺ = 1000 kg/m 3 , one obtains: x 0 = 361.67 and thus

U c = f (362) = 0.233 m/s, (84)
which is a well known value by engineers. It is worth to comment it. If the flow velocity is smaller than U c , there are no instability on the surface. Engineers are used to say that there is no wave on the surface in this case. But the formulation given here enables to have a complete explanation of the phenomenon and can be applied to any geometry. In particular, this enables to take into account the shape of the bottom of the sea in shallow-water or the influence of a submarine under the surface. The evolution of the eigenvalues λ n have been plotted for different values of σ with respect to the wave lenght number k = nπ L . In particular, one can observe that for σ = 0 the sequence of eigenvalues tends to zero as n → ∞. The point is that this result is true for arbitrary geometry and not only the particlar one used by Lighthill. In case of Dirichlet boundary conditions on Γ 1 , computations are much more complicated and not so explicit.

8 Remark on the non linear behaviour for the solution given at section 7

Let us consider that the velocity U is larger than the critical value U c of the system (55). There exists a single eigenmode which is instable. It is denoted by: (λ 1 , w 1 ). The model in ϕ is assumed to be still valid and it is up-graded by the addition of the quadratic term contening the gradient of ϕ in the expression of the pressure.

Let us recall that this term is |∇ϕ| 2 = ∂ϕ ∂ν 2 + |∇ s ϕ| 2 . Another non-linear term should be considered see (4) but is omitted for sake of simplicity. Thus, one obtains the new wave model in the example of section 7:

                 -σ∆ s η + ̺gη + ̺G( ∂ 2 η ∂t 2 ) + ̺U [ ∂G( ∂η ∂t ) ∂x 1 + ∂G( ∂η ∂x 1 ) ∂t ] +̺U 2 ∂G( ∂η ∂x 1 ) ∂x 1 + ̺ 2 [| ∂G( ∂η ∂t + U ∂η ∂x 1 ) ∂x 1 | 2 + | ∂η ∂t + U ∂η ∂x 1 | 2 ] = 0, (85)
where the operator G is the one defined at (17). The variational formulation of this non-linear model is (the functional space V N L should defined in order to make sense to the non-linear terms):

                     find η(t) ∈ V N L such that: ∀v ∈ V N L m s ( ∂ 2 η ∂t 2 , v) + 2U c s ( ∂η ∂t , v) + a s (η, v) + ̺ 2 Γ f [| ∂η ∂t + U ∂η ∂x 1 | 2 + | ∂G( ∂η ∂t + U ∂η ∂x 1 ) ∂x 1 | 2 ]v = 0. ( 86
)
Let us look for an approximate solution which has the simple expression:

η = ξ(t)w 1 .
+ (λ 1 -U 2 )ξ + A ξ2 + B ξξ + Cξ 2 = 0 où A > 0 et C > 0. (87)
In order to illustrate our discussion, the numerical solutions of this equation (87) are plotted on figures 4 and 5. In the first case there is an oscillation of 9 Remark on the time step numerical scheme for η

Let us just give few remarks (in the case where ∆ s Φ0 = 0 on Γ f ) for solving the progressive wave model which can generate instabilities. First of all, let us underline that the finite element method is certainly more appropriate as far as the non-linear term should be taken into account in case of an instability and the eigenvalues of several model should be computed. The difficulty is to use a time-step scheme which remains stable until the physical instability (at least). One criterion is to check the energy conservation for the numerical approximation. One trap to avoid is to transform the gyroscopic effect into a negative damping which would generate an artificial instability. Therefore, a centered scheme in time for this term is better. An up-wind sheme would delayed the apparition of the critical velocity but increase the stability of the scheme. It is obtained for instance by writting the approximate equation at the middle between n∆t and (n + 1)∆t where ∆t is the time step. Let η n be the approximation of η at time t = n∆t. A first possible scheme is:

m s ( η n+1 -2η n + η n-1 ∆t 2 , v) + 2U c s ( η n+1 -η n ∆t , v) + a s ( η n+1 + η n 2 , v) = 0.
By choosing the test function v equal to η n+1η n and making use of Cauchy-Schwarz's inequality one deduces that:

1 2 m s ( η n+1 -η n ∆t , η n+1 -η n ∆t ) + 1 2 a s (η n+1 , η n+1 ) ≤ 1 2 m s ( η n -η n-1 ∆t , η n -η n-1 ∆t ) + 1 2 a s (η n , η n ) ( 88
)
which proves that the scheme is stable as far as a s (defined in (44)) is stricly positive (the energy is decreasing which suggests that a numerical damping is introduced). This is the case if U ≤ √ λ 1 . But if an instability occurs, it is more appropriate to center the scheme around the time n∆t in order to avoid an upward glidding of the critical velocity due to the artificial damping. Therefore one can consider the following scheme which can be stable even if the critical velocity is delayed:

m s ( η n+1 -2η n + η n-1 ∆t 2 , v) + 2U c s ( η n+1 -η n-1 2∆t , v) + a s ( η n+1 + η n-1 2 , v) = 0.
By choosing for the test function v = η n+1η n-1 one obtains now an energy conservation property which ensure that no artificial damping is introduced. But the expression which ensure the stability is different from the energy. It is the next one (which contains twice the kinetical energy):

m s ( η n+1 -η n ∆t , η n+1 -η n ∆t ) + 1 2 a s (η n+1 , η n+1 ) = m s ( η n -η n-1 ∆t , η n -η n-1 ∆t ) + 1 2 a s (η n-1 , η n-1 ) ( 89
)
This result is a useful indicator which proves that the scheme doesn't introduce any artificial damping and therefore the instability can only occur when the bilinear form a s is no more strictly positive. This a strong argument for this scheme even if the physical energy of the numerical sheme is not constant with respect to the time variable.

Conclusion

In this paper a modification of the Neumann-Kelvin's wave model for the instable progressive wave is discussed and analyzed. Let us underline the main points:

1 The role of the capillary is explained in mathematical framework for the physical equations for progressive waves. It is proved that the classical Neumann-Kelvin's model is fully instable and should be modified in order to to be mechanically reliable.

2 It is shown that the best way to take into account the capillary is to use a formulation where the normal displacement η of the free surface appears as a main unknown variable. Thus one can eliminate the velocity potential ϕ making use of the added mass operator. This operator plays a basic role in the equilibrium energy between the gravity and the capillary. On simple examples the first instability due to the velocity U occurs when the energy of the capillary reach the same level as the one of the gravity (see the example treated in section 7). This critical velocity is characterized by an eigenvalue problem.

3 As far as one has to consider a coupling of the sea with a ship, the formulation in η is much more appropriate than the one in ϕ. For instance the slamming between the shell of the ship and the surface of the water can be formulated with respect to η. Furthermore the mathematical model of progressive waves formulated in η -even if it is equivalent to the one in ϕ-has better properties for the numerical implementation.

There are other physical mechanisms which should be added to the progressive wave model. For instance the vortex have been eliminated in the potential model, but the viscosity should play a role in the damping of wave and the characterization of the critical velocity U c . The compressibility of the water is certainly acceptable inside the water, but much less for the surface waves. The coupling with the atmosphere is a fundamental problem and it can change a lot the behaviour of the progressive waves, mainly if temperature gradient are meaningful.

Our goal was to focus on the validity (in fact the non-validity) of the Neumann-Kelvin model) for which a modification is suggested in order to overcome classical difficulties mentioned in many scientific contributions.

  STO 57] [MOA 81] [DOU 95] [KEL 87] [STOE 92] [CAM 02]. It is different because the capillary is taken into account. Furthermore, in the Neumann-Kelvin model, one eliminates the normal displacement η in order to obtain a model where the only unknown is Φ.

  DCF 08] [DCF 09] [HAC 02] [HOL 04] [GAZZ 05][GAZ 07].

•

  Figure 1: Main notations

4. 1

 1 The modelLet us denote by p the pressure in the fluid. From Bernoulli theorem (see for instance G. Duvaut [DUV 90] or L. Landau-E. Lipschitz [LAL 71]), one obtains on the free surface of the water -say Γ f (x 3 = 0)-by considering only the linearized expression in η and ϕ of the pressure (p 0 is the pressure in the air over the free surface):

Lemma 19 .

 19 Let ζ > 0 and let A ξ (w, v) = a s (w, v)ζ((w, v)). There exists a finite dimensionnal subspace E ζ and a real number δ ζ > 0 such that the bilinear form δ ζ m s (w, v)

Figure 2 :

 2 Figure 2: Geometry used for the analytical example

Figure 3 :

 3 Figure 3: First eigenvalue of (55) versus k for several values of σ

Figure 4 :

 4 Figure 4: Non-linear oscillations of the linear instability Choosing the test function v in the one dimensional space spanned by w 1 , one obtains a non-linear differential equation (c s is anti-symmetrical and thus c s (w 1 , w 1 ) = 0): ξ + (λ 1 -U 2 )ξ + A ξ2 + B ξξ + Cξ 2 = 0 où A > 0 et C > 0.(87)

Figure 5 :

 5 Figure 5: Instability of the non linear model

η Wave H L d (depth) U e 1

Acknowledgments The authors thank Th. Gazzola and X.B. Chen from Bureau Veritas for their judicious remarks and valuable comments during this work.