A modelling of springing, whipping and slamming for ships

The slamming phenomenon is a violent impact of the hull of a ship on the free surface of the sea. This loading case is particularly difficult to modelize for several reasons: first of all, the wet surface of the hull is an unknown; then a coupling with the springing (flexibility of the ship) is very complex and finally the interaction with the waves (even if the eigenfrequencies of the structure and the one of the waves are very different) which can be at the origin of important damage mechanisms, involves pointwise effects. This paper aims at giving a simple mathematical model which enables one to simulate the full coupling between these phenomena.


1.
Introduction. Several kinds of ships are concerned by the slamming-springing effect coupled with waves. Two of them are particularly sensitive to this phenomenon. They are on the one hand flat bottom boats like barges or supertankers and on the other hand, the slender hull ships like the fregates or the racing sailing boats. The energy of the movement comes from the waves which is nevertheless a low frequency excitation. But the slamming is a non linear phenomenon which spill over energy from the low frequencies onto much higher ones. One mentioned consequence in maritime reports is an early damage of super-structures even after a restricted number of loading cycles. Recently, this damaging event has been introduced in the building rules for ships. A lot of works have been published through out the scientific journals [23], [22], thesis [19], [33] and congresses [21], [18], for tackling this particularly important problem. But the coupling between waves, slamming and springing is barely considered. As it has been mentioned previously, this decoupling can be fully justified in a linear formulation but the slamming which is highly non-linear, destroys completely such a justification. This is the reason why a first attempt for a fully coupled model is suggested in this paper. It rests on physical assumptions which can certainly be criticized. Nevertheless, it enables one to suggest a mathematical framework which leads to several meaningful results concerning the engineering point of view and to point out major difficulties arising in the analysis. There are several steps in the approach suggested in this paper. The first one, given in section 2, consists in a description of the three media which appear in the physical modelling: the flexible ship, the free surface of the sea and the sea itself. The coupling between the three media requires two relations between each couple in contact. Generally one is the stress equilibrium and the other is the velocity continuity. This is discussed in section 3. The global formulation is then explicited in section 4 and the treatment of the interface between the water and the ship is taken into account through a penalty technique in section 5. Unfortunately, this model is not standard and the already known results (see J.L.Lions [26]) doesn't apply. The reason is that the contact forces are measures in space and time variables. Therefore a regularization method based on a time integration enables to overcome this new difficulty. This is explicited in sections 3 and 4. Finally a series of comments concerning the validity of the model, the unstabilties and the numerical scheme which can be used, are presented in section 6. assumed to be smooth enough, contains two or three main parts, depending on the value of the ship velocity U which is assumed to be constant in all the paper : Γ 1 (the free surface), Γ 0 (the bottom) and Γ 3 (the latteral edge) which can be empty if U = 0 (intake and/or outflow). The set Γ 1 is divided into two portions: Γ 00 where the contact between the water and the ship is possible (see figure 1), and Γ 2 which is always free. Furthermore, let us point out that all the equations of the coupled model are written in the reference configuration which is the frame of the ship moving forwards at the velocity U . The gradient in R 3 of a function f is denoted by ∇f and the tangential component on Γ 1 is ∇ s f . Finally, the laplacian in R 3 is denoted by ∆ and on Γ 1 it is denoted by ∆ s ; the scalar product between vectors of R 3 is denoted by (., .).

2.1.
The steady flow around the ship. The water surrounding the ship is assumed to be incompressible and unviscid. Therefore the flow velocity can be modelled by a potential function -say ϕ 0 -which satisfies the following set of equations: The existence and uniqueness of a solution to (1) are classical (P.A. Raviart and J.M. Thomas [30] and G. Duvaut [14] for instance). It simply requires that the Fredholm condition should be satisfied. In this particular case, it can be written as follows: or else, using the boundary condition satified by ϕ 0 : Γ0∪Γ3 e x .ν = 0.
There is a simple way to establish this property. Let us consider the plane open set Γ x orthogonal to e z containing Γ 2 and the open set interior to Γ 0 ∪ Γ 3 ∪ Γ x denoted by Ω x . From: Γx∪Γ0∪Γ3 e x .ν = ∂Ωx e x .ν = Ωx div(e x ) = 0,

and since
Γx e x .ν = 0, one can conclude that the Fredholm condition is satisfied.
Remark 1. The solution ϕ 0 is proportional to U . Hence one can set: In the following it will be assumed that the the function ϕ 0 is smooth enough in order to justify the computations. For instance ϕ 0 ∈ W 2,∞ (Ω) is sufficient and realistic from the general regularity results for elliptic PDE's (see P. Grisvard [20]).
2.2. The flow model due to the ship movements. Due to the ship movements the velocity in the fluid is modified. Hence, the new velocity potential can be written: The general equation satisfied by ϕ in the open set Ω is the following one ( [12]) where the upper dot stands for the partial time derivative, c f is the wave velocity in the fluid and T > 0 is the time delay : When the fluid is assumed to be incompressible, this equation is reduced to the following one: This is the case that we consider in this paper but one must notice that it is clearly an approximation. Concerning the boundary condition on Γ 0 ∪ Γ 3 that one can prescribe for ϕ, there are several possibilities. When U = 0, the most simple case is ∂ϕ ∂ν = 0. It corresponds to a boat floating in a tank with rigid walls. From Fredholm analysis [30], a necessary condition in order to define ϕ is Γ1 ∂ϕ ∂ν = 0. This implies a restriction (the average is zero) on the normal displacement to the free surface of the sea. When U = 0 the portion Γ 0 ∪ Γ 3 of the boundary ∂Ω can not be made of walls because the water gets in through one side and get out on an other one. In this case, one can adopt for instance the condition ϕ = 0 on these parts denoted by Γ 3 . Let then us set, for The mechanical interpretation is that on Γ 3 the velocity is normal to the boundary Γ 3 . These boundary conditions should be modified for the compressible case which is discussed [10].
2.3. The free surface of the sea. The whole free surface of the sea at rest is denoted by Γ 1 (figure 1). The normal displacement of the fluid to this boundary denoted by η. On the part Γ 2 , with Γ 2 ⊂ Γ 1 , the movement of the free surface is governed by the equilibrium between the gravity, the capillarity force and the pressure in the fluid. The unit normal to Γ 1 outwards Ω is denoted by ν. Because at rest, the portion Γ 2 is horizontal, one has ν = e z (hence ν z = 1) on it. But on the part of Γ 00 , with Γ 00 ⊂ Γ 1 , which can be in contact with the boat, the inclination of the surface Γ 1 should be considered. Concerning the free surface which can be in contact with the boat, one also has to take into account the pressure due to the ship. If N is the unit inwards normal to the hull, the pressure applied on the free surface of the sea by the ship, is: Conversely the one applied on the structure by the fluid is −F . Because the water doesn't stick on the hull, one has the condition: Let us point out that λ = 0 where the hull is not in contact with the sea. The equilibrium of the forces applied on the free surface Γ 1 leads to: (σ is the capillarity constant, ̺ f is the mass density of the fluid, g is the gravity and the tangential displacements are neglected because the viscosity is not taken into account): where the constant c(t) depends on the choice of the constant ϕ depends on (it is a potential function).
In order to complete the model, one should prescribe boundary conditions on ∂Γ 1 . Several possibilities occur. When U = 0 one can set for instance: which means that the water glides perfectly on the lateral boundary of Ω (far away from the ship). If U = 0 the perfect gliding is more complicated to write. One can adopt the following boundary condition (see [12]): In order to simplify the characterization of the constant c(t) one can prescribe the additional condition which enables to separate the computation of the constant c(t) and the function η in (10): Let us point out that when U = 0 this is a necessary condition due to the incompressibility of the fluid and to the Neumann homogeneous boundary condition satisfied by ϕ on Γ 0 . For sake of brevity this relation is kept for all cases in the following. More generally, the following functional spaces are used (< ., . > is the duality between the spaces H −1/2 (Γ 1 ) and H 1/2 (Γ 1 )): Let us recall that c(t) can not be prescribed a priori. It appears as the Lagrange multiplier of the relation (11). Finally the equation which traduces the equilibrium of the free surface can be written in a variational form as follows, making use of the boundary conditions satisfied by η on ∂Γ 1 and because of (11) the constant c(t) which is not zero, nevertheless disappears: where the bilinear form a sf (., .) is defined by: The mathematical properties of this bilinear form are classical (see J. Necas [28] or P.A. Raviart and J.M. Thomas [30]). But it should be worth noting that because ν z can be negative, one has to use an extension of the regular ness property. Let us summarize this in the following statement which is sometimes called Garding's coerciveness. It is a consequence of the compact embedding from H 1 (Γ 1 ) into L 2 (Γ 1 ); therefore, the one from V into H −1/2 0 (Γ 1 ) is also compact.
Theorem 2.1. Let a sf the bilinear and symmetrical form defined on the space V at (14). There exists a real number a 0 ≥ 0 such that: where ||.|| 1,Γ1 is the norm in the space H 1 (Γ 1 ) and ||.|| −1/2,Γ1 the one in H The result is obvious if ν z > 0. But it can happen, for a large class of hull, that ν z < 0 on a part of it (see figure 2 below). The ship is a flexible structure the movement of which can be modelled by a standard hyperbolic system which can be written in a variational formulation. Let us denote by q the generalized displacement field of the surface ω which represents the hull of the ship. A dynamical shell model can be used. Let us denote by m h (q, p) the bilinear form representing the inertia and by a h (q, p) the one for the stiffness. The functional space in which q and p are chosen is denoted by W . For a standard Koiter [8] or Budiansky [9] shell model, it is isomorphic to [H 1 (ω)] 2 × H 2 (ω) up to rigid body modes (P.G. Ciarlet [8]). The time derivatives of the displacements denoted byq belongs for classical shell or even 3D structural model, to a space denoted by H and which is isomorphic to [L 2 (ω)] 3 . If inertia terms are considered in a shell model, one should use for H the space [L 2 (ω)] 2 × H 1 (ω). But one can also use a three dimensional model on the structure occupying the open set S. Then, the components of q are in the space [H 1 (S)] 3 and the trace on ω are only in the space [H 1/2 (ω)] 3 . In this case, the velocitiesq are in the space [L 2 (S)] 3 . Because the ship is free, there are no prescribed kinematical conditions in W . Assuming that the boat is in equilibrium (the weight is compensated by the static Archimede force), the movement of the hull is governed by the variational model where λ is in fact the transient pressure applied by the ship on the surface of the sea: Furthermore, initial conditions should be prescribed at t = 0: q(0) = q 0 andq(0) = q 1 . In a large number of applications it is sufficient to restrict the previous model to a finite dimensional one by considering only a finite number of eigenmodes for the ship (which can be represented by a more complete model than a shell one). But one meaningful advantage of the shell model compared to the three dimensional model (even if only few eigenmodes are kept) is that the deformed configuration avoids angles ((q, N ) ∈ H 2 (ω)) which can induce new singularities in the distribution λ (the pressure).
3. The coupling relations. There are three mechanical models which interact: one is the ship assumed to be flexible, the second one is the free surface of the sea and finally, the sea itself. The coupling between the three media is derived from continuity conditions for the kinematics and for normal pressure (continuities of tangential components are omitted in our formulation as far as the fluid is assumed to be unviscid and therefore the boundary layers are cancelled). The two coupling which appear in the full model are discussed herafter but separately.
3.1. Coupling between the sea and the free surface. Let us recall that ∇ s is the tangential gradient of a function on the surface Γ 1 . The velocity of a material point on the free surface is : U ∇ sφ0 + ∇ϕ. The unit normal to the free surface in the deformed configuration can be approximate by the following expression as far as the curvature of this surface is small enough (P-G. Ciarlet [7], [8]) ν ′ = ν + θ where θ = −∇ s η. Hence the normal fluid velocity is up to second order terms: and keeping only the first order terms (with respect to the perturbation ϕ or η), one obtains: This is the normal velocity of the free surface which is also equal up to second order terms toη. Finally, one has the following boundary condition on Γ 1 : This boundary condition enables one to write the variational formulation for the fluid as follows. First of all, let us define the functional space for ϕ at an instant t by: The velocity potential ϕ has to satisfy the variational model: It clearly appears that ϕ is the summ of two contributions (both linear with respect to η): one proportional toη and the other to U ∇ s η. But when U = 0 the average of ∂ϕ ∂ν is zero on Γ 1 and therefore the Fredholm condition is automatically satisfied. When U = 0 this condition is no more necessary because ϕ = 0 on Γ 3 . Hence one can set: where the linear operator Φ is defined from the functional space gψds.
Therefore one can define a symmetrical bilinear form on the functional space H −1/2 (Γ 1 ) for U = 0 and H −1/2 0 (Γ 1 ) for U = 0, which is continuous and coercive, by: With the previous notations, the free surface model can be written: This equation must be discussed in the context of the full coupled model.

3.2.
Coupling between the free surface and the hull. In fact, the continuity relations have already been taken into account in the previous sections. But one is still to be discussed and it is the most important one. It concerns the non penetration of the water inside the ship. First of all, let us introduce the gap function denoted by ψ which represents the free space between the free surface of the sea and the hull of the boat in the stationary configuration. It is assumed that ψ ∈ C 0 ∩ H 1 (Γ 00 ) which is a realistic assumption. This function is defined on the boundary Γ 00 and ψ ≥ 0. One has ψ = 0 on the part of the free surface of the sea which is in contact with the hull at rest (see figure 4). It is measured along the direction of the vector N . The non-interpenetration of the water into the boat can be approximately explicited by the relation: A reasonable assumption is that (see figure 3): This condition is only an approximation which is consistant up to second order as far as the movements of the water are small enough. For instance, the couple of points on the free surface and on the hull which can be in contact are not known a priori. A first order approximation can be used by considering that they are on a same line parallel to the normal N to the hull of the ship. This result can be derived from differentiable manifold theory (see P. G. Ciarlet [9]). Hence, one defines the admissible convex set K for the admissible set of the displacements η and q, by: This is a closed and non empty convex set of V × W . It is assumed in all the paper that the initial conditions are in K (see (25)). If the couple (η, q) solution of the mechanical model is such that on a part of Γ 00 one has η(N, ν) < (q, N ) + ψ the water is not in contact with the hull of the ship and therefore on the same part one also has λ = 0. Finally the coupling condition between the hull and the free surface can be written as follows: But an important difficulty appears in this formulation (which is only an approximation!) because λ is a measure and the previous characterization is difficult to handle in a mathematical framework. Hence a new formulation is required in order to make sense to this non interpenetration condition. This the goal of the next sections.
This model is ill posed because there are an infinite number of solutions. The problem comes from the lack of reflexion condition on the obstacle (x = 0). One can see several possible solutions on figure 6. In order to select the true one it is necessary to traduce the energy conservation (or not) when the material point touches the obstacle at an instant t 0 . The energy conservation is in this caseẋ(t 0 − 0) = −ẋ(t 0 + 0). But this condition is difficult to traduce for the full coupled model. Hence a strategy based on a penalty method and a time integration technique enables one to overcome this difficulty. Nevertheless, in case of a damping phenomenon during the contact it is possible to take it into account in the ship model or in the free surface one. But this is out the scope of this paper. Important contributions on these questions are given in the works of M. Schatzman [31] and L. Paoli-M. Schatzman [29] for pointwise contacts and by A. Bamberger and M. Schatzman [32] [2] for a vibrating string.
4.1. The variational model. Let us summarize the model which has been presented previously but which is ill posed because of the lack of information on the conservation of energy before and after the contact. There are two unknown fields: the normal displacement η of the free surface and the displacement of the boat q. One has (η, q) ∈ K and it FIGURE 6. Solutions to (27); top: absorption, center: conservation, bottom: damping satisfies the following variational inequation directly deduced from (26): Unfortunately the non-uniqueness remark can be applied to this model. Hence it is necessary to select the solutions by introducing an additional information on the energy conservation (or not) during the shock with the obsatcle (water against the hull of the ship).

4.2.
The penalty strategy (energy conservation). In order to select a rebounding strategy let us imagine that the non interpenetration method is replaced by the action of a very stiff spring, the stiffness of which being denoted by 1 ε , with ε << 1. The variational inequation is then replaced by a non linear model but which has the huge advantage to be governed by an equation. It is explicited hereafter: The existence and the uniqueness of a solution to this penalized model can be derived from a Galerkin approximation method. For sake of brevity, let us just give the result for U = 0.
The case U = 0 is discussed in subsection 6.2 for a particular situation and more generally in [10]. First of all, the eigenmodes {w h k , λ h k } of the stucture defined by: The six first eigenvalues are zero because of the rigid body motions of the structure. For sake of commodity, the eigenmodes are ordered such that: . . It is known that the family {w h k } is an Hilbert basis of the space W [30], [3]. Similarly, we introduce the eigenmodes of the fluid -say {w sf k , λ sf k } -which are solution of: Here again the classical spectral theory can be applied (using the shift of the eigenvalues due to the Garding coerciveness (see theorem 2.1). The family {w sf k } is an Hilbert basis of the space V (see [30], [3]) and the eigenmodes are also ordered by increasing values of the eigenvalues.
Let us define by V J and W J the finite dimensional spaces included respectively in W and V such that: The approximation model consists in finding (η J , q J ) ∈ V J × W J which satisfies initial conditions (for (q J (0),q J (0))) and (η J (0),η J (0)), which are defined as the projection (with respect to m sf and m h ) of the true one (those of the continuous model) on the finite dimensional spaces W J and V J and such that: This is a second order ordinary differential equation (not inequation) satisfying the regular Lipschitz assumption. Therefore, for any ε one can ensure that there is a unique solution for given initial data. The point is that one can derive an energy estimate on the sequences (η J , q J ). It is necessary to assume that the initial condition which satisfy the non-interpenetration condition (in the set K defined at (25)), are in the space of finite energy. This is obtained by setting v =η J and p =q J in the approximate model (35). One obtains: These estimates enable one to extract subsequences which converge weakly to a solution of the following variational model: and satisfying the initial conditions: η(0) = η 0 ,η(0) = η 1 , q(0) = q 0 ,q(0) = q 1 . Details can be found in [30], [26] and for this particular case in [10]. The uniqueness (for ε > 0!) can be obtained from the same energy estimate by considering two solutions satisfying the same initial conditions and because of the monotonicity of the term (see J.L. Lions [26]): In order to study the limit of this model when the small parameter ε tends to zero, it is convenient to use a change of variables based on a trajectory method. It means that the unknown are the integrals of the displacements with respect to time. But first of all let us point out that the solution of (29) satisfies an energy invariant property when U = 0 (the general case is briefly discussed at subsection 6.2). The following result proves that the energy of the initial conditions is kept by the model for any ε.
The unique solution (η ε , q ε ) of (36), satisfies the invariant property: Proof. The proof is quite classical and obtained by choosing v =η ε and p =q ε in (36).
A classical consequence of this conservation property is that the penalty contact-model doesn't damp the energy defined in theorem 4.1. In fact, the result is obtained for U = 0 and, for U = 0, it can't be obtained because the system is open since kinetical energy enters the open set Ω and some is getting out. But nothing guarantees that this is a balanced phenomenon. This characteristic can be at the origin of instabilities which are discussed in subsection 6.2 and in [11]. Therefore there is not global energy conservation in this case.

The time integration technique (strong penalty model).
In order to analyze the asymptotic behaviour of this penalty model it is much more convenient to use a path method instead of a state variable one (the path is the time integration of a state variable). The basic reason is that a measure can appear at each contact point and the estimate are not obvious at all as far as the state variable are kept. Let us set formally: One can check directly the following relations: Therefore the path defined from the solution of the penalty model is also solution of: In order to explain how this model can be interpreted from a mahematical point of view, let us come back to the simple spring model. The integrated penalty model consists in setting: therefore X ε is solution of: One can derive estimates onẌ εẊ ε , X ε and on Λ ε (t) = 1 ε t 0 (Ẋ ε ) + (ξ)dξ in the space L ∞ (]0, T [) as soon as the initial condition satisfiesẊ ε (0) ≤ 0 which enables to state a convergence result to a limit model. Therefore, up to a subsequence:  The limit (Λ * , X * ) satisfies: In (42) the difficulty is to prove the last relation and to give a pointwise interpretation to it. It is obtained by observing thatẊ ε converges in C 0 ([0, T ]) toẊ * when ε → 0 becauseẌ ε is bounded in the space L ∞ (]0, T [). Let us consider an instant t 0 such thatẊ * (t 0 ) < 0 (the result is true forẊ * (t 0 ) = 0). For ε small enough, there exists a neighbourhood of t 0 such thatẊ * (t) < 0 on this neighbourhood. ButẊ ε converges uniformly toẊ * when ε → 0. Finally, one can define ε 0 and a neighbourhood -say O -of t 0 on which ∀ε < ε 0 ,Ẋ ε < 0. Therefore Λ ε is constant on this neighbourhood, and its time derivative is zero. Hencė Therefore:Λ * Ẋ * = 0 because O can be arbitrary small.
In fact, this implies that the time derivative of X * is a solution of the initial model. And from the energy conservation property of the penalty model, the limit is a solution which keeps the energy. Unfortunately this is not always sufficient to ensure the uniqueness as far as there are several degrees of freedom. The advantage of the integrated formulation is that the fonctions are more regular and this enables a priori estimates and compactness results that are not true with the initial model of the full coupled model.

A weak penalty model (energy absorption).
The penalty term which appears in the integrated model, is very stiff because it takes into account all the history of the solution.
One can weaken this penalty term by replacing it by an instantaneous one. Thus a new integrated model which is no more equivalent to (36) -but which is mechanically motivated -is defined by: In order to explain the modifications induced by this change in the penalty term, let us come back to the simple spring model. The new weak penalty (integrated) model consists in solving:Ẍ Thus, the penalty term appears as a viscous damping on X ε whenẊ ε > 0 (first order time derivative). One can check that when ε → 0 the limit of X ε denoted by X * , is the unique solution of the following model (see the proof for the full coupled model at subsection Choosing V = 1.5Ẋ * and V = 0.5Ẋ * leads to (after a time integration between 0 and t): Let us consider an instant t 0 such thatẊ * (t 0 ) = 0 (contact point). One has from (46): or else: Therefore X * (t n ) takes always the same value for each contact instant t n after t 0 . Setting this solution in (46), one obtains thatẊ * (t n ) = 0.
From the uniqueness result, one can claim that X * (t) = 0 ∀t ≥ t 0 . Hence the solution sticks exactly on the obstacle after the first contact. One has a perfect absorption. This is not surprising because the damping coefficient is 1 ε and it tends to the infinity when ε → 0.

5.
The mathematical analysis. In order to simplify the notations, the following analysis is restricted to the case where U = 0. The extension to a more general situation (U = 0) can be obtained in a similar way as far as no instability appears. Such phenomenae seem to be similar to breaking waves problems. This is discussed in [10]. The plan of the analysis given in this section is the following one: first of all the existence and uniqueness of a solution to (39) or (43) is checked (in fact, for U = 0 the model (39) is equivalent to (36)).
Then several a priori estimates are derived by energy methods for both models (strong and weak integrated penalty ones). This enables one to extract subsequences which converge in a weak sense to a solution of a limit model in each case. Finally, the uniqueness of a solution to the limit model is proved for the so-called weak penalty model. The time derivative of each integrated limit model is a solution of the initial one.
They satisfy the following properties which can be proved easily from classical results and theorems 2.1 and 3.1: Let us point out how these properties can be used in the following analysis. First of all, the difficulty is that one has not the coerciveness of the bilinear form A(., .). But one has the Garding's inequality above which is sufficient as shown in the next result.
where c 0 is a constant which is independent on n. Then the sequence X n is bounded in the Proof. From Helmotz's theorem (spectral decomposition of V using the eigenvectors of A with respect to M), one can restrict the proof of theorem 5.1 to a finite dimensional space of V for which A is negative.
, the estimate satisfied by X n is (using the orthogonality with respect to the bilinear and symmetrical forms M and A and ξ (ξ > 0) being the opposite of the smallest negative eigenvalue of A): But one has also: which implies: A summation from 1 to K leads to the existence of two strictly positive constants E and F such that: From Gronwall's lemma (see for instance [24]), one can conclude that α k (t) is bounded on [0, T ] and this prove theorem 5.1.

5.1.3.
Estimates onẌ ε ,Ẋ ε and X ε versus ε when U = 0 for the strong penalty model. Let us consider the initial penalty model the solution of which is X ε . The energy estimate recalled at theorem 4.1 can also be interpreted in terms of X ε (for instance (η ε ,q ε ) = ∂ 3 X ε ∂t 3 ) such that: with (the initial condition is assumed to satisfy the non interpenetration condition): From theorem 5.1 one derives the following estimates: from which one deduces that: These a priori estimates enable one to extract subsequences of X ε denoted X ε ′ such that: One can note that this implies: From the equation (39) setting (let us consider for instance that the functional space for the structure is W = [H 1 (S)] 3 , S being the open set occupied by the structure, and therefore the trace (N, q) is onto in the functional space H 1/2 (Γ 00 )): and using (52) one has also the following estimate (the function (N, ν)v − (N, p) describes the space H 1/2 (Γ 00 ) because (N, p), p ∈ W , can be an arbitrary function in this space): This estimate enables one to extract a subsequence say Λ ε such that If a shell model is used for the structure, the function (N, q) ∈ H 2 (ω) and the leading term in the estimate is (N, ν)v which belongs to the space H 1 (Γ 00 ). Hence the estimate (59) would be replaced by: But additional estimates on the time derivative of Λ ε ′ are necessary in order to characterize the limit. Let us consider the equation of the free surface with the penalty term. By integrating on Γ 1 , one obtains (the support of Λ ε ′ is Γ 00 ): In order to obtain an estimate on the constant c ε ′ (t) which doesn't depend on space coordinates, one can localize the equation (52) on the boundary Γ 2 where Λ ε = 0. Again, from (52) and with test functions v the support of which are restricted to Γ 2 , one obtains (G j , j = 1, 2, 3, 4, 5 are five constants): Hence, from the assumption: (N, ν) ≥ n 0 > 0, one gets (becauseΛ ε ≥ 0): As a consequence, Λ ε ′ is also bounded in the space W 1,1 (]0, T [; L 1 (Γ 1 )). With another respect, a sequence bounded in L 1 (]0, T [×Γ 00 ) is also bounded in the space of measures denoted M([0, T ] ×Γ 00 ) which is also the dual of C([0, T ] ×Γ 00 ). Therefore, one can extract subsequences such that: From the preceding convergences one can take the limit in the model (52). Thus one obtains: and the initial conditions are the following ones: But this is insufficient for a characterization of the limit model. It is necessary to prove for instance that: where the duality has to be precised. Unfortunately this results is still out of our reach for 3-D models and furthermore, it is not obvious to make sense to it. This is possible for 2-D models (ie. Γ 1 is a one dimensional manifold) if the functions (N, q), with q ∈ W , are for instance in the space H 1 (Γ 00 ). This is the case for instance if the structure is modelled by a bending shell because one has then: (N, q) ∈ C 0 ([0, T ]; H 2 (ω)) and (N,q) ∈ C 0 ([0, T ]; L 2 (ω)). A similar property could be drived from a finite dimensional approximation for W in case of a three dimensional modelling using smooth enough eigenvectors of the ship-structure for instance. From the property : ). This is due to the compact embedding of H 1 (Γ 1 ) into C 0 (Γ 1 ). Therefore the restriction to [0, T ] ×Γ 00 of the sequences (η ε ′ , q ε ′ ) converges strongly to (η * , q * ) in the space C 0 ([0, T ] ×Γ 00 ).
The consequence of the preceding convergences and the meaning of (66), is such that: These equations (64)-(67), enable one to state that the time derivative of the limit model (for a subsequence!) is a solution of the initial model. Furthermore, the global energy is kept in the elementary case of the spring because of the C 0 convergence of Λ ε towards Λ * . In this simple example (0 − D) the conservation of energy ensures the uniqueness of a solution. But the uniqueness is not proved and certainly false, for the more general situation which is discussed in this paper.
In 3-D the hypothesis of this analysis is satisfied by the components of the displacements of the structure if a bending shell model is used. But an additional regularity on the deflection η of the free surface is required. It could be obtain by a more accurate capillarity model including the change of the curvature which require a fourth order operator (see for instance P.G. Ciarlet [9]).

5.2.
Mathematical analysis of (43) for U = 0. There are three steps: 1. Existence and uniqueness of a solution to (43), 2. A priori estimate on this solution with respect to ε, 3. Derivation of a limit model. Let us make explicit the model for U = 0 using the notations introduced previously: find X ε (t) = (H ε (t), Q ε (t)) ∈ V such that: and the initial conditions are the following ones:

Existence of solution to the weak penalty model (68).
The proof is similar to the one given in J.L. Lions [26]. It rests upon the convexity of the term: Therefore its weak derivative is a monotone operator. Furthermore, the uniqueness is also a consequence of this property. The proof can be explicited using an eigenvector basis as explained briefly in subsection 4.2. The following closed convex set is useful in the next sections: The a priori estimates with respect to ε for (68). An estimate is needed for constructing the limit of the penalty term. It can be obtained by considering the time derivative of (43) and choosing Y =Ẍ as a test function. This leads to the following upper bound where E(0) is the energy at time t = 0 of the initial model recalled at (55): where χẊ / ∈K1 = 1 ifẊ / ∈ K 1 and 0 else.Therefore one gets the estimates with respect to ε: Let us come back to the equations (68). From the estimates (72), one deduces by a similar method as in the previous section that: (72), one can deduce that there exists subsequences X ε ′ and Λ ε ′ such that: 5.2.3. The limit characterization for the weak penalty model. The results recalled here are similar to those of J. L. Lions [26]. Therefore, details are omitted. Let us now define a linear form on V by: From (58) one can claim thatẊ * ∈ K 1 . Let us take the limit for ε ′ → 0 in (68). From the weak convergence of the sequence X ε ′ given at (74), one deduces that: and one can prove that: Hence X * is solution of the following variational equation: The next step consists in proving that (76) has a unique solution. Let us consider two solutions of (76) denoted respectively by X 1 and X 2 . Let us set respectively Y =Ẋ 2 and Y =Ẋ 1 in the inequations satisfied respectively by X 1 and X 2 , one obtains after adding the two expressions: or else, because the initial conditions are the same for X 1 and X 2 and because T is arbitrary: which proves the uniqueness of a solution. A classical consequence of the uniqueness of the weak limit is that the whole sequence converges to this limit.

Remark 2.
The proof for uniqueness rests on an assumed regularity of the solution. But one can get rid of it by using the trick introduced by J.L. Lions [26] and which consists in using another time integration in the test functions.
A natural question here again, is to prove that the unique solution of the weak penalty model is a solution of the initial model. First of all, let us set: From (76) one can claim that Λ * ≥ 0 because: which implies that: a.e. t ∈]0, T [, Ψ * (Ẋ * )(t) = 0.
But this relation is difficult to interprete. Let us consider the two dimensional case. In this case one has:Ẋ ε ∈ [C 0 ([0, T ]; C 0 (Γ 1 ))] 2 . This is due to the inclusion: H 1 (Γ 1 ) ⊂ C 0 (Γ 1 ) and because one hasẊ ε →Ẋ * in the space [C 0 ([0, T ] ×Γ 1 )] 2 . Let us consider a point S 0 ∈ Γ 00 and an instant t 0 ∈]0, T [ such that:Ẋ * (t 0 , S 0 ) < 0. The convergence ofẊ ε toẊ * < 0 is uniform in time and space in [C 0 ([0, T ] ×Γ 1 )] 2 . Hence one can find a neighbourhood O of t 0 × S 0 in ]0, T [×Γ 00 such thatẊ ε = (η ε , q ε ) satisfies the condition: Therefore if we set: one first has that: and from the weak* convergence of Λ ε to Λ * in the space L ∞ (]0, T [; H −1/2 (Γ 00 )), one can deduce that This is a way to interprete the formula (80). In the three dimensional case, a similar result could be obtained by considering lines drawn on Γ 00 and defining an approximate value -sayX ε -forẊ ε on this line by taking (for instance) the average along a tranverse segment L 0 to this line. Then the convergence in C 0 ([0, T ] × L 0 ) enables to extract some interpretation of (80) similar to the previous one.
Let us come back to the local equations satisfied by the time derivatives of X * . From (80) and (81), one can conclude that it is a solution of the initial model.

5.3.
Remarks on the the limit model for the weak penalty one. The main feature of this weak penalty model is that it can swallow the energy during the contact as it has been shown on the spring model. In fact, this mechanical behaviour is certainly the nearest to the physical reality. The rebounding of water on the hull is a complex phenomenon (see the breaking wave problem [4] [21]) but which would require to model the sparkling effect (breaking the surface into small drops) which are mainly governed by the capillarity. In our model this aspect is not taken into account. Hence the full damping of the kinetical energy of the water seems to be a consistant choice. Nevertheless, one can adjust a restitution coefficient e during the contact just by combining the two model which have been described in this paper. For instance one can consider a convex combination of the two penalty terms but with different penalty parameters. This will be discussed in a forthcoming paper.
6. Few extensions of this study.

6.1.
The quasi-variational model in case of added mass strategy. Let us discuss in this section the numerical approximation of the coupled fluid-structure model. In fact, the operator Φ (see (20)) which is used in the bilinear form m(., .) couples all the points of the free surface Γ 1 . Therefore the inertia matrix associated to m(..), is a wide surface one. This is due to the elliptic operator Φ which couples all the points of the free surface. Because the movement of the global ship is obviously slower than the one of the surface of the sea (as far as fundamental mode are concerned and mainly for the six rigid body modes), it is a great temptation to solve the free surface model for a given movement of the boat. In this case, the force due to the sea and applied on the ship appears as a function of its movement. The constraint concerning the non penetration of the water into the hull can be written as: which is the general framework of quasivariational equations following the terminology introduced by C. Baiocchi [1]). The solution method can be a fixed point algorithm on the structural movements setting (k is the iteration index): One advantage is that one can use two distinct time steps: one for the free surface movement and the other which will be larger, for the structure (as far as rigid body modes and low frequencies are concerned). This is the idea developed by C. Fahrat and M. Lesoine [15] called staggered algorithm. When the contact non linearity can be omited this strategy leads to the so -called added mass method (see J.P. Morand and R. Ohayon, [27]. Further details are given in [10]. 6.2. Influence of the ship velocity. Let us come back to the case where U = 0. In order to make simple the points that we would like to emphasize, let us consider the two dimensional case of a flat bottom barge floating on the free surface of the sea as shown on figure 7. In this situation the solutionφ 0 (stationary flow) is obtained analytically and one has (x 10 is the abscissa of the middle of the rectangle Ω): which enables one to rewrite the model (29) as follows: The existence and uniqueness of a solution are given in [10] following the same method as for U = 0 explained at section 4.2. Setting v =η ε and p =q ε , one can observe first of all from the definition of Φ at (20), that (the first property is due to the antisymmetry of the Corriolis term): such that one deduces the following energy estimate: The last term which is proportional to U 2 , represents the energy stored in the waves of the fluid and due to the rotation of the unit normal to the free surface (or the inclination). From this expression one can observe that a static instability (negative eigenvalue of the stationary bilinear form) can occur if U is large enough. It is a so called breaking wave instability. More precisely, let us consider the static symmetrical bilinear form defined on the space V by: Let us introduce the following eigenvalue model for the free surface (gravity waves and capillarity riddles): find (λ, η) ∈ R × V such that: ∀v ∈ V, a sf (η, v) = λm sf ( ∂η ∂x 1 , ∂v ∂x 1 ).
But any function -say vin the space V satisfies Γ1 v = 0. Hence one has for example, the following generalization of Poincaré's inequality (see J. Necas [28]): ∃c p > 0, such that ∀v ∈ V, ||v|| −1/2 ≤ c p || ∂v The natural injection from V into H −1/2 (Γ 1 ) is compact. Hence, there exists a family of eigenvectors which is a basis of V . Let us us denote by λ 1 the smallest eigenvalue solution of (89). For U larger than √ λ 0 , the bilinear form (88) can become negative, and therefore a so-called static instability (breaking wave in our terminology) can occur.
7. Comments on this study and prospect works. The method which has been used can be summarized on figure 6. The physical model is written in terms of continuity (for the displacements and for the normal stresses) or non-penetration (of the water inside the ship). But if these conditions are sufficient in statics, they are not as far as dynamical models are concerned. A constitutive relationship is necessary in order to manage the contact between the water and the structure. One way consists in ensuring the energy conservation -but globally -between the ship and the sea. This is done through a penalty approximation of the unilateral contact condition. Nevertheless, the limit of this penalty model is difficult to handle for ε → 0 (non uniqueness). This difficulty is due to the fact that the acceleration in the wave movement can contains mathematical measures (and therefore the contact force between the free surface of the water and the hull of the ship also contains measures). In particular this distribution can be supported by the line (in 3-D) separating the wet surface of the hull and the dry one. They are clearly connected to the jet flow observed in the slamming phenomenon. But many additional phenomena should be taken into account in order to get clother the physical reality. For instance, the wave is a nonlinear phenomenon (terms like ̺ f 2 |∇ϕ| 2 or the computation of the true points of the free surface which can be in contact with the hull), the velocity U of the boat can reach an instability value inducing a breaking wave phenomenon, a trapping of the air which can be compressed between the hull and the sea are also important aspects in this challenge [21]). The integrated penalty model(s) (existence and uniqueness) Figure 6. Functional scheme of the models introduced 8. Conclusion. A quite simple model for the interaction between slamming, springing and whipping phenomenon for ships, has been discussed in this paper. There are several restrictive assumptions which enable one to derive a mathematical model. For instance, the fluid is unviscid and incompressible, flow-jets are omitted, the points of the free surface of the sea and those of the hull which can be in contact are defined a priori, everything is linearized.... Let us underline two results which have a physical importance. First of all, the absorption of the kinetical energy of the sea, for which the spring effect is very small, can be modelled from penalty path models. Secondly, if the velocity of the ship is large enough, an instability can appear on the free surface of the sea. The model that has been suggested is mathematically well founded and reliable solution methods exist. But, the restrictions concerning the linearization and the definition a priori of points which can be in contact, are certainly too restrictive and future investigations should be carried out for a non linear coupling model from both the mathematical and the numerical point of view.