Philippe Destuynder

Caroline Fabre

A MODELLING OF SPRINGING, WHIPPING AND SLAMMING FOR SHIPS

Keywords: 2000 Mathematics Subject Classification. Primary: 74F10, 35L85, 35R35; Secondary: 76B15, 35L70, 35Q35 Fluid-structure interaction, unilateral contact, wave equation, variational inequation, algorithms for dynamical coupled models. 1

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 [START_REF] Korobkin | Second-order Wagner theory of wave impact[END_REF], [START_REF] Holloway | Added mass of whipping modes for ships at high Froude number by a free surface boundary element method coupled with strip theory[END_REF], thesis [START_REF] Th | Modélisation du phénomène de tossage pour les bateaux[END_REF], [START_REF] Storhaug | Experimental investigation of wave induced vibrations and their effect on the fatigue loading of ships[END_REF] and congresses [START_REF] Hackett | Whipping loads due to aft body slamming[END_REF], [START_REF] Th | Three-dimensional Wagner problem using variational inequalities[END_REF], 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 [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]) 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:

                
-∆ϕ 0 = 0 in Ω, ∂ϕ 0 ∂ν = -U ν.e x on Γ 0 ∪ Γ 3 and ∂ϕ 0 ∂ν = 0 on Γ 1 = Γ 00 ∪ Γ 2 , Ω ϕ 0 = 0.

(

) 1
The existence and uniqueness of a solution to (1) are classical (P.A. Raviart and J.M. Thomas [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF] and G. Duvaut [START_REF] Duvaut | Mécanique des milieux continus[END_REF] for instance). It simply requires that the Fredholm condition should be satisfied. In this particular case, it can be written as follows:

∂Ω ∂ϕ 0 ∂ν = 0, (2)
or else, using the boundary condition satified by ϕ 0 :

Γ0∪Γ3 e x .ν = 0. (3)
There is a simple way to establish this property. Remark 1. The solution ϕ 0 is proportional to U . Hence one can set:

ϕ 0 = U φ0 . (4)
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:

ϕ T = U φ0 + ϕ. (5)
The general equation satisfied by ϕ in the open set Ω is the following one ( [START_REF] Ph | Vibration des systèmes couplés[END_REF]) 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 :

φ + 2U ∇ φ0 .∇ φ + U 2 ∇ φ0 .∇(∇ φ0 .∇ϕ) -c 2 f ∆ϕ = 0, in Ω×]0, T [. (6)
When the fluid is assumed to be incompressible, this equation is reduced to the following one:

-∆ϕ = 0 in Ω×]0, T [. (7
) 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 [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF], 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

U = 0 Γ 3 = ∂Ω -(Γ 0 ∪ Γ 1 )
. 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 [START_REF] Ph | Compressibility of the fluid in a slamming model for ships[END_REF].

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:

F = λN. (8)
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:

λ ≤ 0. (9)
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):

-σ∆ s η + gην z + ̺ f φ + ̺ f U ∇ s φ0 .∇ s ϕ + c(t) = λ(N, ν), (10)
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:

σ ∂η ∂ν Γ1 = 0,
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 [START_REF] Ph | Vibration des systèmes couplés[END_REF]):

σ ∂η ∂ν Γ1 = ̺ f U ∂ φ0 ∂ν Γ1 ϕ.
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 [START_REF] Ph | Compressibility of the fluid in a slamming model for ships[END_REF]:

Γ1 η(s, t)ds = 0. (11)
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 )):

               L 2 0 (Γ 1 ) = {v| v ∈ L 2 (Γ 1 ), Γ1 v(s)ds = 0} H -1/2 0 (Γ 1 ) = {v| v ∈ H -1/2 (Γ 1 ), < v, 1 >= 0}, V = H 1 (Γ 1 ) ∩ L 2 0 (Ω). (12)
Let us recall that c(t) can not be prescribed a priori. It appears as the Lagrange multiplier of the relation [START_REF] Ph | Aéroélasticité et aéroacoustique[END_REF]. 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:

                               ∀t ∈]0, T [, η ∈ V is such that: ∀v ∈ V : a sf (η, v) + ̺ f U [ Γ1 (∇ s φ0 .∇ s ϕ)vds - ∂Γ1 ∂ φ0 ∂ν Γ1 v] + Γ1 ̺ f φvds = a sf (η, v) -̺ f U Γ1 [(∆ s φ0 ϕ v + ϕ (∇ s φ0 , ∇ s v)]ds + Γ1 ̺ f φvds = Γ1 λ(N, ν)vds, (13)
where the bilinear form a sf (., .) is defined by:

∀η, v ∈ V, a sf (η, v) = Γ1 [σ∇ s η.∇ s v + gν z ηv]ds. (14)
The mathematical properties of this bilinear form are classical (see J. Necas [START_REF] Necas | Les méthodes directes en théorie des équations elliptiques[END_REF] or P.A. Raviart and J.M. Thomas [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF]). 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 [START_REF] Duvaut | Mécanique des milieux continus[END_REF]. There exists a real number a 0 ≥ 0 such that: 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 [START_REF] Ciarlet | Mathematical Elasticity, Vol. III : Theory of Shells[END_REF] or Budiansky [START_REF] Ciarlet | An Introduction to Differential Geometry, with Applications to Elasticity[END_REF] shell model, it is isomorphic to [H 1 (ω)] 2 × H 2 (ω) up to rigid body modes (P.G. Ciarlet [START_REF] Ciarlet | Mathematical Elasticity, Vol. III : Theory of Shells[END_REF]). The time derivatives of the displacements denoted by q 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 velocities q 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:

∀v ∈ V, a sf (v, v) + a 0 ||v|| 2 -1/2,Γ1 ≥ ||v||
       ∀t ∈]0, T [, q(t) ∈ W such that: ∀t ∈]0, T [, ∀p ∈ W, m h (q, p) + a h (q, p) = - Γ00 λ(N, p)ds. (15)
Furthermore, initial conditions should be prescribed at t = 0: q(0) = q 0 and q(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 [START_REF] Ciarlet | Introduction to Linear Shell Theory[END_REF], [START_REF] Ciarlet | Mathematical Elasticity, Vol. III : Theory of Shells[END_REF]) ν ′ = ν + θ where θ = -∇ s η. Hence the normal fluid velocity is up to second order terms:

u ≃ (ν ′ , U ∇ s φ0 + ∇ϕ)
and keeping only the first order terms (with respect to the perturbation ϕ or η), one obtains:

u ≃ (∇ϕ, ν) -U (∇ φ0 , ∇ s η).
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 :

∂ϕ ∂ν = η + U (∇ s η, ∇ s φ0 ). (16)
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:

V = {ψ| ψ ∈ H 1 (Ω), ψ = 0 on Γ 3 if U = 0 and Γ1 ψ = 0 if U = 0}. (17)
The velocity potential ϕ has to satisfy the variational model:

       ∀t ∈]0, T [, ϕ(t) ∈ V, such that : ∀t ∈]0, T [, ∀ψ ∈ V, Ω (∇ϕ(t), ∇ψ)dx = Γ1 η(t)ψds + U Γ1 (∇ s η(t), ∇ s φ0 )ψds. (18)
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:

ϕ = Φ( η) + U Φ((∇ s φ0 , ∇ s η)), (19)
where the linear operator Φ is defined from the functional space H -1/2 (Γ 1 ) for U = 0 (and from

H -1/2 0 (Γ 1 ) if U = 0) into V by:        g ∈ H -1/2 (Γ 1 ) ( respectively H -1/2 0 if U = 0), → Φ(g) ∈ V such that: ∀ψ ∈ V, Ω (∇Φ(g), ∇ψ)dx = Γ1 gψds. (20)
The restriction of Φ(g) on Γ 1 is denoted by Φ Γ1 (g). It is an element of the space H 1/2 (Γ 1 ), (and even of H 1/2 (Γ 1 ) ∩ L 2 0 (Γ 1 ) if U = 0). One has the standard following statement [START_REF] Brezis | Analyse fonctionnelle[END_REF].

Theorem 3.1. There exists two strictly positive constants c 0 and c 1 such that (H -1/2 (Γ 1 ) is replaced by

H -1/2 0 (Γ 1 ) if U = 0): ∀g ∈ H -1/2 (Γ 1 )), c 0 ||Φ Γ1 (g)|| 1/2,Γ1 ≤ ||g|| -1/2,Γ1 ≤ c 1 ||Φ Γ1 (g)|| 1/2,Γ1 .
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:

∀g, h ∈ H -1/2 (Γ 1 ), (respectively H -1/2 0 (Γ 1 ) if U = 0) : m sf (g, h) = ̺ f Ω (∇Φ(g), ∇Φ(h))dx def = ̺ f < g, Φ Γ1 (h) > H -1/2 (Γ1)×H 1/2 (Γ1) . (21)
With the previous notations, the free surface model can be written:

∀v ∈ V, m sf (η, v) + a sf (η, v) +̺ f U Γ1 {Φ Γ1 ((∇ s φ0 , ∇ s η))v -Φ Γ1 ( η)[(∇ s φ0 , ∇ s v) + ∆ s φ0 v ]}ds -̺ f U 2 Γ1 [(∇ s φ0 , ∇ s v) + ∆ s φ0 v ]Φ Γ1 ((∇ s φ0 , ∇ s η))ds = Γ00 λ(N, ν)v. ( 22
)
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:

η(N, ν) ≤ (q, N ) + ψ on Γ 00 . (23)
A reasonable assumption is that (see figure 3):

(N, ν) ≥ n 0 > 0 on Γ 00 . ( 24
)
FIGURE 3. Illustration of the condition [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF] 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 [START_REF] Ciarlet | An Introduction to Differential Geometry, with Applications to Elasticity[END_REF]).

(q,N) ψ FIGURE 4
. Geometrical approximation of the gap between the boat and the free surface of the sea Hence, one defines the admissible convex set K for the admissible set of the displacements η and q, by:

K = {(v, p)| (v, p) ∈ V × W, (N, ν)v ≤ (p, N ) + ψ on Γ 00 }. ( 25
)
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 [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués[END_REF]). 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:

           i) λ ≤ 0, on Γ 00 , ii) (N, ν)η ≤ (q, N ) + ψ on Γ 00 ×]0, T [, iii) λ[(N, ν)η -(q, N ) -ψ] = 0 on Γ 00 ×]0, T [. (26)
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.

   ẍ + ω 2 x ≤ 0, x(0) = x 0 , ẋ(0) = x 1 ,
x ≤ 0 and x(ẍ

+ ω 2 x) = 0. ( 27
)
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 [START_REF] Schatzman | A class of nonlinear equations of second order time[END_REF] and L. Paoli-M. Schatzman [START_REF] Paoli | Penalty Approximation for dynamical systems submitted to multiple non-smooth constraints[END_REF] for pointwise contacts and by A. Bamberger and M. Schatzman [START_REF] Schatzman | Un problème hyperbolique du 2ème ordre avec contrainte unilatérale : la corde vibrante avec obstacle ponctuel[END_REF] [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

                       (η, q)(t) ∈ K such that ∀(v, p) ∈ K, m h (q, p -q) + a h (q, p -q) + m sf (η, v -η) + a sf (η, v -η) +̺ f U Z Γ 1 ΦΓ 1 ((∇s φ0, ∇s η))(v -η) -ΦΓ 1 ( η)[(∇s φ0, ∇s(v -η)) + ∆s φ0 (v -η) ] -̺ f U 2 Z Γ 1 [(∇s φ0, ∇s(v -η)) + ∆s φ0 (v -η) ]ΦΓ 1 ((∇s φ0, ∇sη)) ≥ 0. (28)
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).

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:

                                       ∀t ∈]0, T [, (η ε , q ε )(t) ∈ V × W such that: ∀(v, p) ∈ V × W, m h (q ε , p) + a h (q ε , p) + m sf (η ε , v) +̺ f U Γ1 {Φ Γ1 ((∇ s φ0 , ∇ s ηε ))v -Φ Γ1 ( ηε )[(∇ s φ0 , ∇ s v) + ∆ s φ0 v ]}ds +a sf (η ε , v) -̺ f U 2 Γ1 [(∇ s φ0 , ∇ s v) + ∆ s φ0 v ]Φ Γ1 ((∇ s φ0 , ∇ s η ε ))ds + 1 ε Γ00 [η ε (N, ν) -(q ε , N ) -ψ] + [v(N, ν) -(p, N )]ds = 0. ( 29
)
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 [START_REF] Ph | Compressibility of the fluid in a slamming model for ships[END_REF]. First of all, the eigenmodes {w h k , λ h k } of the stucture defined by:

           w h k ∈ W, λ h k ∈ R + , ∀p ∈ W, λ h k m h (w h k , p) = a h (w h k , p), m h (w h k , w h k ) = 1. (30)
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:

λ h 7 ≤ λ h 8 ≤ . . . ≤ λ h k ≤ λ h k+1 ≤ . . .

It is known that the family {w h

k } is an Hilbert basis of the space W [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF], [START_REF] Brezis | Analyse fonctionnelle[END_REF]. Similarly, we introduce the eigenmodes of the fluid -say {w sf k , λ sf k } -which are solution of:

           w sf k ∈ V, λ sf k ∈ R, ∀v ∈ V, λ sf k m sf (w sf k , v) = a sf (w sf k , v), m sf (w sf k , w sf k ) = 1. (31)
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 [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF], [START_REF] Brezis | Analyse fonctionnelle[END_REF]) 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:

W J = {p| p = k=1,J α h k w h k }. ( 33
)
V J = {v| v = k=1,J α sf k w sf k }. ( 34
)
The approximation model consists in finding (η J , q J ) ∈ V J × W J which satisfies initial conditions (for (q J (0), qJ (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:

                   ∀p ∈ W J , ∀v ∈ V J , m h (q J , p) + a h (q J , p) - 1 ε Γ00 [(N, ν)η J -(q J , N ) -ψ] + (N, p) = 0, m sf (η J , v) + a sf (η J , v) + 1 ε Γ00 [(N, ν)η J -(q J , N ) -ψ] + (N, ν)v = 0. (35)
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 = qJ in the approximate model (35). One obtains:

1 2 {m h ( qJ , qJ ) + m sf ( ηJ , ηJ ) + a h (q J , q J ) + a sf (η J , η J ) + 1 ε Γ00 |(η J (N, ν) -(q J , N ) -ψ) + | 2 ds} = constant/time
As a consequence the sequences qJ and ηJ are bounded in the spaces:

L ∞ (]0, T [; H) and L ∞ (]0, T [; H -1/2 0 (Γ 1 )).
The sequences q J and η J are respectively bounded in the spaces:

L ∞ (]0, T [; W ) and L ∞ (]0, T [; V ).
These estimates enable one to extract subsequences which converge weakly to a solution of the following variational model:

                   ∀p ∈ W, ∀v ∈ V, m h (q, p) + a h (q, p) - 1 ε Γ00 [(N, ν)η -(q, N ) -ψ] + (N, p)ds = 0, m sf (η, v) + a sf (η, v) + 1 ε Γ00 [(N, ν)η -(q, N ) -ψ] + (N, ν)vds = 0, (36)
and satisfying the initial conditions: η(0) = η 0 , η(0) = η 1 , q(0) = q 0 , q(0) = q 1 . Details can be found in [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF], [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] and for this particular case in [START_REF] Ph | Compressibility of the fluid in a slamming model for ships[END_REF].

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 [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]):

(v, p) ∈ V × W → |[(N, ν)v -(p, N ) -ψ] + | 2 .
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 ε.

Theorem 4.1. Let (η 0 , η 1 , q 0 , q 1 ) ∈ V × H -1/2 0 (Γ 1 ) × W × H. The unique solution (η ε , q ε ) of ( 36), satisfies the invariant property:

∀ε, ∀t ∈ [0, T ], E(t) = E(0),
where the energy E = E(t) is defined by:

E = 1 2 {m h ( qε , qε ) + m sf ( ηε , ηε ) + a h (q ε , q ε ) + a sf (η ε , η ε ) + 1 ε Γ00 |(η ε (N, ν) -(q ε , N ) -ψ) + | 2 ds}
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 [START_REF] Ph | Aéroélasticité et aéroacoustique[END_REF]. 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:

H ε (s, t) = t 0 η ε (s, ξ)dξ, on Γ 1 Q ε (s, t) = t 0 q ε (
s, ξ)dξ on the structure. (37) One can check directly the following relations:

Qε = qε , Ḧε = ηε Qε = q ε , Ḣε = η ε . ( 38
)
Therefore the path defined from the solution of the penalty model is also solution of:

                                                     (H ε , Q ε )(t) ∈ V × W such that ∀t ∈]0, T [, ∀(v, p) ∈ V × W, m h ( Qε , p) + a h (Q ε , p) + m sf ( Ḧε , v) +̺ f U Z Γ 1 {ΦΓ 1 ((∇s φ0, ∇s Ḣε ))v -ΦΓ 1 ( Ḣε )[(∇s φ0, ∇sv) + ∆s φ0 v ]}ds +a sf (H ε , v) -̺ f U 2 Z Γ 1 [(∇s φ0, ∇sv) + ∆s φ0 v ]ΦΓ 1 ((∇s φ0, ∇sH ε ))ds + 1 ε Z t 0 Z Γ 00 [ Ḣε (N, ν) -( Qε , N ) -ψ] + [v(N, ν) -(p, N )]ds = m h ( qε (0), p) + m sf ( ηε (0), v) +̺ f U Z Γ 1 {ΦΓ 1 ((∇s φ0, ∇sη(0)))v -ΦΓ 1 (η(0))[(∇s φ0, ∇sv) + ∆s φ0 v ]}ds. (39)
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:

                             X ε (t) = t 0 x ε (ξ)dξ, where x is solution of ẍε + ω 2 x ε + 1 ε ẋε = 0, x ε (0) = x 0 , ẋε (0) = x 1 ,
therefore X ε is solution of:

Ẍε + ω 2 X ε + 1 ε t 0 ( Ẋε ) + = x 1 , X ε (0) = 0, Ẋε (0) = x 0 . (40)
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:

                           Ẍε ⇀ Ẍ * in L ∞ (]0, T [) weak*, X ε → X * in C 1 ([0, T ]), 1 ε t 0 ( Ẋε ) + ⇀ Λ * in L ∞ (]0, T [) weak * and strongly in C 0 ([0, T ]) t 0 ( Ẋε ) + → 0 as ε in L ∞ (]0, T [). (41)
Furthermore, as 1 ε t 0 ( Ẋε ) + is bounded in the space W 1,1 (]0, T [), it converges weakly* in the space BV ([0, T ]) which ensures that Λ * is a measure on [0, T ] and it makes sense to the term Λ * Ẋ * .

The limit (Λ * , X * ) satisfies:

Ẍ * +ω 2 X +Λ * = x 1 , X * (0) = 0, Ẋ * (0) = x 0 , Ẋ * ≤ 0, Λ * ≥ 0, Λ * Ẋ * = 0. (42)
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. Hence 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:

                                                     (H ε , Q ε )(t) ∈ V × W such that ∀t ∈]0, T [: ∀(v, p) ∈ V × W, m h ( Qε , p) + a h (Q ε , p) + m sf ( Ḧε , v) ̺ f U Z Γ 1 {ΦΓ 1 ((∇s φ0, ∇s Ḣε ))v -ΦΓ 1 ( Ḣε )[(∇s φ0, ∇sv) + ∆s φ0 v ]}ds +a sf (H ε , v) -̺ f U 2 Z Γ 1 [(∇s φ0, ∇sv) + ∆s φ0 v ]ΦΓ 1 ((∇s φ0, ∇sH ε ))ds + 1 ε Z Γ 00 [ Ḣε (N, ν) -( Qε , N ) -ψ] + [v(N, ν) -(p, N )]ds = m h ( qε (0), p) + m sf ( ηε (0), v) +̺ f U Z Γ 1 {ΦΓ 1 ((∇s φ0, ∇sη(0)))v -ΦΓ 1 (η(0))[(∇s φ0, ∇sv) + ∆s φ0 v ]}ds. (43)
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:

Ẍε + ω 2 X ε + 1 ε ( Ẋε ) + = x 1 , X ε (0) = 0, Ẋε (0) = x 0 . (44)
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 5.2.3):

   Ẋ * ≤ 0, X * (0) = 0, Ẋ * (0) = x 0 , ∀V ≤ 0, ( Ẍ * + ω 2 X * -x 1 )(V -Ẋ * ) ≥ 0. (45)
Choosing V = 1.5 Ẋ * and V = 0.5 Ẋ * leads to (after a time integration between 0 and t):

ω 2 ( Ẋ * ) 2 (t) + ω 2 (X * ) 2 (t) = 2x 1 X * (t) + x 2 0 . (46)
Let us consider an instant t 0 such that Ẋ * (t 0 ) = 0 (contact point). One has from (46):

ω 2 (X * ) 2 (t 0 ) -2x 1 X * (t 0 ) -x 2 0 = 0, (47)
or else:

X * (t 0 ) = x 1 ± ω 2 x 2 0 + x 2 1 ω 2 , (only + is meaningful). ( 48
)
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 [START_REF] Ph | Compressibility of the fluid in a slamming model for ships[END_REF]. 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.

5.1. Analysis of the strong penalty model (39).

5.1.1.

Existence and uniqueness of a solution to (39) for U = 0. The existence of a solution is a direct consequence of the construction of the integrated penalty model. In fact, as far as ε = 0 (29) and (39) are equivalent. The point is to obtain a priori estimates on the solution (H ε , Q ε ) and to take the limit when ε → 0.

5.1.2.

The formulation of (39) when U = 0. There are two kinds of estimates. The first ones are obtained by setting (v, p) = ( Ḣε , Qε ) = (η ε , q ε ) in the integrated penalty model (39) where one has introduced the restriction U = 0. The second ones come from the initial penalty model. But both are necessary to get a convergence result. First of all, let us introduce few notations which are convenient in the mathematical proof. The functional space (H ε , Q ε ) belongs to, is denoted by V = V × W and the time derivatives ( Ḣε , Qε ) belongs to the space H = H -1/2 0 × H. Two bilinear forms are defined on V by:

           ∀X = (H, Q), Y = (K, P ) ∈ V : M c (X, Y ) = m h (Q, P ) + m sf (H, K), A(X, Y ) = a h (Q, P ) + a sf (H, K). ( 50
)
They satisfy the following properties which can be proved easily from classical results and theorems 2.1 and 3.1:

           ∃ c 0 , c 1 > 0, s.t. ∀Y ∈ V, c 0 ||Y || 2 H ≤ M(Y, Y ) ≤ c 1 ||Y || 2 H , ∃ d 1 > 0, s.t. ∀Y ∈ V, A(Y, Y ) ≤ d 1 ||Y || 2 V , ∃ d 0 > 0, τ > 0 such that: d 0 ||Y || 2 V ≤ A(Y, Y ) + τ ||Y || 2 H . (51)
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.

Theorem 5.1. Let us consider a sequence X n ∈ C 0 ([0, T ]; V) ∩ C 1 ([0, T ]; H) which satisfies the estimate:

M( Ẋn , Ẋn ) + A(X n , X n ) ≤ c 0
where c 0 is a constant which is independent on n. Then the sequence X n is bounded in the space C 0 ([0, T ]; V) ∩ C 1 ([0, T ]; H).

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. Hence if X n (t) = k=1,K α k (t)w k (x)
, 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):

k=1,K [ α2 k -ξα 2 k ] ≤ c 0 .
But one has also:

α k (t) = α k (0) + t 0 αk (s)ds,
which implies:

∀t ∈ [0, T ] : α k (t) 2 ≤ max(1, T ) 2 [α 2 k (0) + t 0 α2 k (s)ds].
A summation from 1 to K leads to the existence of two strictly positive constants E and F such that:

k=1,K α2 k (t) ≤ D + F t 0 k=1,K α2 k (s)ds.
From Gronwall's lemma (see for instance [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]), one can conclude that α k (t) is bounded on [0, T ] and this prove theorem 5.1.

The model (39) can also be written as follows (for U = 0):

                       find X ε (t) = (H ε (t), Q ε (t)) ∈ V such that: ∀Y = (K, P ) ∈ V, M( Ẍε , Y ) + A(X ε , Y ) + 1 ε Γ00 { t 0 [(N, ν) Ḣε -(N, Qε ) -ψ] + }[(N, ν)K -(N, P )] = m h ( qε (0), P ) + m sf ( ηε (0), K), (52)
and the initial conditions are the following ones:

X ε (0) = (0, 0), Ẋε (0) = (η ε (0), q ε (0)). (53)
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:

1 2 [M( Ẍε , Ẍε ) + A( Ẋε , Ẋε ) + 1 ε Γ00 {[(N, ν)η ε -(q ε , N ) -ψ] + } 2 ](t) = E(0), (54)
with (the initial condition is assumed to satisfy the non interpenetration condition):

E(0) = 1 2 [M( Ẍε , Ẍε ) + A( Ẋε , Ẋε )](0). ( 55
)
From theorem 5.1 one derives the following estimates:

∀t ∈ [0, T ], || Ẍε || 2 H (t) + || Ẋε || 2 V (t) ≤ c, (56)
from which one deduces that:

∀t ∈ [0, T ], ||X ε || 2 V (t) ≤ c ′ . ( 57
)
These a priori estimates enable one to extract subsequences of X ε denoted X ε ′ such that:

                   X ε ′ ⇀ ε ′ →0 X * in L ∞ (]0, T [; V) weak *, Ẋε ′ ⇀ ε ′ →0 Ẋ * in L ∞ (]0, T [; V) weak *, Ẍε ′ ⇀ ε ′ →0 Ẍ * in L ∞ (]0, T [; H) weak *, ||[(N, ν)η ε -(q ε , N ) -ψ] + || 2 L ∞ (]0,T [;L 2 (Γ00)) ≤ c 0 ε. ( 58
)
One can note that this implies:

(N, ν)η * -(N, q * ) -ψ ≤ 0 on Γ 00 .
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 )):

Λ ε ′ = 1 ε ′ t 0 [(N, ν)η ε ′ -(q ε ′ , N ) -ψ] + dt ∈ L ∞ (]0, T [; 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):

||Λ ε ′ || L ∞ (]0,T [;H -1/2 (Γ00)) ≤ c. (59)
This estimate enables one to extract a subsequence say Λ ε such that

Λ ε ⇀ Λ * in L ∞ (]0, T [; H -1/2 (Γ 00 )) weak*. ( 60
)
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:

||Λ ε ′ || L ∞ (]0,T [;(H 1 ) ′ (Γ00)) ≤ c. (61)
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 ): Γ00 Λε ′ (N, ν) = -̺g Γ1 ν z η ε ′ -c ε ′ (t)|Γ 1 |. (62)
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):

| t 0 c ε ′ (µ)dµ| ≤ c |Γ 2 | [G 0 + G 1 ||η ε || L ∞ (]0,T [;V ) + G 2 ||ϕ ε ′ (t)|| 0,Γ1 ≤ G 4 .
Hence, from the assumption: (N, ν) ≥ n 0 > 0, one gets (because Λε ≥ 0):

n 0 || Λε ′ || L 1 (]0,T [×Γ1) = n 0 t 0 Γ00 Λε ′ ≤ t 0 Γ00 Λε ′ (N, ν) ≤ G 5 .
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:

   Λ ε ′ ⇀ Λ * ≥ 0 in L ∞ (]0, T [; M( Γ00 )) weak*, Λε ′ ⇀ Λ * ≥ 0 in M([0, T ] × Γ00 ) weak*. (63)
From the preceding convergences one can take the limit in the model (52). Thus one obtains:

                   X * (t) = (H * (t), Q * (t)) ∈ V, (N, ν)η * -(N, q * ) -ψ ≤ 0, Λ * ≥ 0, Λ * ≥ 0 such that: ∀Y = (K, P ) ∈ V, M( Ẍ * , Y ) + A(X * , Y )+ < Λ * , [(N, ν)K -(N, P )] > = m h ( q * (0), P ) + m sf ( η * (0), K), (64)
and the initial conditions are the following ones:

X * (0) = (0, 0), Ẋ * (0) = (η * (0), q * (0)), (65)
But this is insufficient for a characterization of the limit model. It is necessary to prove for instance that:

< Λ * , [(N, ν)η * -(N, q * ) -ψ] >= 0, (66)
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 :

C 0 ([0, T ]; H 1 (Γ 1 )) ∩ C 1 ([0, T ]; L 2 (Γ 1 )) ⊂ C 0 ([0, T ] × Γ1 ) with compactness one can ensure that (N, ν)η ε ′ -(N, q ε ′ ) converges in C 0 ([0, T ] × Γ00
). 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:

< Λ * , [(N, ν)η * -(N, q * ) -ψ] > M([0,T ]× Γ00)/C([0,T ]× Γ00) = 0. (67)
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 [START_REF] Ciarlet | An Introduction to Differential Geometry, with Applications to Elasticity[END_REF]). 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: ∀Y = (K, P ) ∈ V, M( Ẍε , Y ) + A(X ε , Y ) + 1 ε Γ00 [(N, ν) Ḣε -(N, Qε ) -ψ] + }[(N, ν)K -(N, P )] = m h ( qε (0), P ) + m sf ( ηε (0), K), (68)
and the initial conditions are the following ones:

X ε (0) = (0, 0), Ẋε (0) = (η ε (0), q ε (0)), (69) 5.2

.1. Existence of solution to the weak penalty model (68).

The proof is similar to the one given in J.L. Lions [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]. It rests upon the convexity of the term:

Ẋ = ( Ḣ, Q) → 1 ε Γ00 |[((N, ν) Ḣ -(N, Q) -ψ] + | 2 .
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:

K 1 = {Y | Y = (K, P ) ∈ L ∞ (]0, T [; V), (N, ν)K -(P, N ) -ψ ≤ 0 on Γ 00 }. (70)
5.2.2. 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):

1 2 [M( Ẍε , Ẍε ) + A( Ẋε , Ẋε )](t)+ + 1 ε t 0 Γ00 χ Ẋ / ∈K1 {[(N, ν) ηε -( qε , N )] 2 } = E(0), (71)
where χ Ẋ / ∈K1 = 1 if Ẋ / ∈ K 1 and 0 else.Therefore one gets the estimates with respect to ε:

       || Ẍε || 2 H (t) + || Ẋ|| 2 V (t) ≤ c, t 0 Γ00 χ Ẋ / ∈K1 {[(N, ν) ηε -( qε , N )] 2 } ≤ cε. (72)
Let us come back to the equations (68). From the estimates (72), one deduces by a similar method as in the previous section that:

Λ ε = 1 ε [(N, ν)η ε -(N, q ε ) -ψ] + (73)
is bounded in the space L ∞ (]0, T [; H -1/2 (Γ 00 )) in the case where W ≡ [H 1 (S)] 3 and in L ∞ (]0, T [; (H 1 ) ′ (Γ 00 )) if one uses a shell model for the structure. Finally, from (73) and (72), one can deduce that there exists subsequences X ε ′ and Λ ε ′ such that:

                   X ε ′ → ε ′ →0 X * in L ∞ (]0, T [; V) weak *, Ẋε ′ → ε ′ →0 Ẋ * in L ∞ (]0, T [; K 1 ) weak *, Ẍε ′ → ε ′ →0 Ẍ * in L ∞ (]0, T [; H) weak *, Λ ε ′ → ε ′ →0 Λ * in L ∞ (]0, T [; H -1/2 (Γ 00 )) weak * for W ≡ [H 1 (S)] 3 for instance. (74)
5.2.3. The limit characterization for the weak penalty model. The results recalled here are similar to those of J. L. Lions [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]. Therefore, details are omitted. Let us now define a linear form on V by:

∀Y = (K, P ) ∈ V, L(Y ) = m h ( q(0), P ) + m sf ( η(0), K)-< Λ * (t), ψ > . (75)
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:

∀Y ∈ K 1 , T 0 M( Ẍ * , Y ) + A(X * , Y ) ≥ T 0 L(Y )
and one can prove that:

T 0 M( Ẍ * , Ẋ * ) + A(X * , Ẋ * ) ≤ T 0 L( Ẋ * ).
Hence X * is solution of the following variational equation:

               find X * ∈ L ∞ (]0, T [; V), Ẋ * ∈ K 1 such that: ∀Y ∈ K 1 , T 0 M( Ẍ * , Y -Ẋ * ) + A(X * , Y -Ẋ * ) ≥ T 0 L(Y -Ẋ * ),
X * (0) = 0, Ẋ * (0) = (η(0), q(0)).

(

) 76
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:

T 0 M( Ẍ1 -Ẍ2 , Ẋ2 -Ẋ1 ) + A(X 1 -X 2 , Ẋ2 -Ẋ1 ) ≥ 0, (77)
or else, because the initial conditions are the same for X 1 and X 2 and because T is arbitrary:

∀T ≥ 0, [M( Ẋ1 -Ẋ2 , Ẋ1 -Ẋ2 ) + A(X 1 -X 2 , X 1 -X 2 )](T ) = 0,
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 [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] 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:

∀Y ∈ L ∞ (]0, T [; V), T 0 Ψ * (Y ) = - T 0 {M( Ẍ * , Y ) + A(X * , Y ) -L(Y )} = lim ε→0 1 ε T 0 Γ00 [(N, ν)η ε -(N, q ε ) -ψ] + [(N, ν)K -(N, Q) -ψ], =< Λ * , [(N, ν)K -(N, Q) -ψ] > . (78)
From (76) one can claim that Λ * ≥ 0 because:

∀Y ∈ K 1 , T 0 Ψ * (Y )(t)dt ≤ 0.
Furthermore:

T 0 Ψ * ( Ẋ * )(t)dt = 0, (79)
which implies that:

a.e. t ∈]0, T [, Ψ * ( Ẋ * )(t) = 0. ( 80
)
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:

∀(t, x) ∈ O, (N, ν)η ε -(N, q ε ) -ψ < 0 .
Therefore if we set:

Λ ε = 1 ε [(N, ν)η ε -(N, q ε ) -ψ] + ,
one first has that:

Λ ε ≡ 0 on O => Λε ≡ 0 on O,
and from the weak* convergence of

Λ ε to Λ * in the space L ∞ (]0, T [; H -1/2 (Γ 00 )), one can deduce that Λ * ≡ 0 on O. (81)
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 -say Xε -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 [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF]) 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:

q ∈ K(q, q), (82)
which is the general framework of quasivariational equations following the terminology introduced by C. Baiocchi [START_REF] Baiocchi | Problèmes à frontière libre en hydraulique : milieu non homogènes[END_REF]). The solution method can be a fixed point algorithm on the structural movements setting (k is the iteration index):

q k+1 ∈ K(q k , qk ). ( 83
)
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 [START_REF] Farhat | Two efficient staggered algorithms for the serial and parallel solution of the three dimensional nonlinear transcient aeroelastic problems[END_REF] 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, [START_REF] Morand | Interactions fluides-structures[END_REF]. Further details are given in [START_REF] Ph | Compressibility of the fluid in a slamming model for ships[END_REF].

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 Ω):

φ0 (x) = -(x 1 -x 10 ), => ϕ 0 (x) = -U (x 1 -x 10 ), (84)
which enables one to rewrite the model ( 29) as follows:

                                       (η ε , q ε )(t) ∈ V × W such that ∀t ∈]0, T [: ∀(v, p) ∈ V × W, m h (q ε , p) + a h (q ε , p) + m sf (η ε , v) -̺ f U Γ1 [Φ Γ1 ( ∂ ηε ∂x 1 ) vds -Φ Γ1 ( ηε ) ∂v ∂x 1 ds] +a sf (η ε , v) -̺ f U 2 Γ1 Φ Γ1 ( ∂η ε ∂x 1 ) ∂v ∂x 1 ds + 1 ε Γ1 [η ε -(q ε , e z ) -ψ] + [v -(p, e z )]ds = 0. (85)
The existence and uniqueness of a solution are given in [START_REF] Ph | Compressibility of the fluid in a slamming model for ships[END_REF] 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 [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF], that (the first property is due to the antisymmetry of the Corriolis term):

i) ̺ f U Γ1 [Φ Γ1 ∂ ηε ∂x 1 ηε ds -Φ Γ1 ( ηε ) ∂ ηε ∂x 1 ds] = 0, and
(dx = dx 1 dx 2 ) ii) -̺ f U 2 Γ1 ∂Φ Γ1 ∂x 1 ( ∂η ε ∂x 1 ) ∂ ηε ∂x 1 ds = d dt {- ̺ f U 2 2 Ω |∇Φ Γ1 ∂η ε ∂x 1 | 2 (x)dx}
such that one deduces the following energy estimate:

E U ε (t) = E U 0 (0), (86)
E U ε (t) = 1 2 {m sf ( ηε , ηε ) + m h ( qε , qε ) + a sf (η ε , η ε ) + a h (p ε , p ε )} + 1 2ε Γ1 |[η ε -(q ε , e z ) -ψ] + | 2 ds - ̺ f U 2 2 Ω |∇Φ Γ1 ( ∂η ε ∂x 1 )| 2 ds (87)
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:

η ∈ V → a sf (η, η) -̺ f U 2 Ω |∇Φ Γ1 ( ∂η ∂x 1 )| 2 (x)dx. ( 88
)
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
).

(89)

This eigenvalue problem satisfies the standard assumptions of a variational eigenvalue model (see P.A. Raviart and J.M. Thomas [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF]) because on the one hand, the bilinear form a sf (., .) is continuous and coercive on the space V (in our example a sf (., .) is directly Vcoercive because ν z = 1), and on the other hand, m sf (., .) is continuous and coercive on the space H -1/2 (Γ 1 ). Therefore:

∀v ∈ V, m sf ( ∂v ∂x 1 , ∂v ∂x 1 ) ≥ c m || ∂v ∂x 1 || -1/2 . ( 90
)
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 [START_REF] Necas | Les méthodes directes en théorie des équations elliptiques[END_REF]):

∃c p > 0, such that ∀v ∈ V, ||v|| -1/2 ≤ c p || ∂v ∂x 1 || -1/2 . ( 91
)
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 [START_REF] Hackett | Whipping loads due to aft body slamming[END_REF]).

The physical model (existence but non uniqueness if energy restitution is not prescribed; but it is difficult to do so) -→

Addition of a penalty term

The penalty model (existence and uniqueness)

Derivation      
Integration and adjustment of the penalty terms

The integrated limit model(s) (existence and sometimes uniqueness) ←-

ε → 0
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.

2 .FIGURE 1 .

 21 FIGURE 1. Geometrical notations (Γ c ⊂ Γ 00 is the contact area). It is time dependent.

  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.

FIGURE 2 .

 2 FIGURE 2. Case of a hull with negative values of ν z

4 .FIGURE 5 .

 45 FIGURE 5. The spring model

FIGURE 6 .

 6 FIGURE 6. Solutions to (27); top: absorption, center: conservation, bottom: damping

  Λε Ẋ * = 0 on the open set O. Finally, one can prove that: O Λ * Ẋ * = lim ε→0 O Λε Ẋ * = 0. Therefore: Λ * Ẋ * = 0 because O can be arbitrary small.

FIGURE 7 .

 7 FIGURE 7. A particular case of a flat bottom ship

  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 figure2 below).

	2 1,Γ1 ,
	where ||.|| 1,Γ1 is the norm in the space H 1 (Γ 1 ) and ||.|| -1/2,Γ1 the one in H	-1/2 0	(Γ 1 ).


Wei Wang

Ting Zhang
email: ting.zhang@ccnu.edu.cn

Yiwen Zhao

Xinxin Jin

Harold Mouchere

Xinguo Yu

Improving Oracle Bone Characters Recognition via A CycleGAN-based Data Augmentation Method

Keywords: Characters recognition, Class imbalance, Data augmentation, CycleGAN

Oracle bone inscription is the earliest writing system in China which contains rich information about the history of Shang dynasty. Automatically recognizing oracle bone characters is of great significance since it could promote the research on history, philology and archaeology. The proposed solutions for oracle bone characters recognition are mainly based on machine learning or deep learning algorithms which rely on a large number of supervised training data. However, the existing dataset suffers from the problem of severe class imbalance. In this work, we propose a CycleGAN-based data augmentation method to overcome the limitation. Via learning the mapping between the glyph images data domain and the real samples data domain, CycleGAN could generate oracle character images of high-quality. The quality is evaluated using the quantitative measure. Totally, 185362 samples are generated which could serve as a complementary to the existing dataset. With these generated samples, the state of the art results of recognition task on OBC306 are improved greatly in terms of mean accuracy and total accuracy.

Introduction

Oracle bone inscriptions, usually inscribed on turtle shells and animal bones, are the oldest hieroglyphs in China. They record the development of civilization of late Shang dynasty. The study of oracle bone inscriptions could promote the research on the origin of Chinese characters, as well as the research on the history of Shang dynasty. Since firstly discovered in the city of Anyang, China, around 160000 bones or shells (see Figure 1) were unearthed gradually from which about 4500 classes of characters were found however more than 2000 characters remain not deciphered so far. Currently, the recognition and decipherment of oracle bone characters majorly rely on the manual method by domain experts with high expertise which to a certain degree blocks the development of this field. Thus, developing automatic recognition methods is in need and of great significance. The state of the art methods for oracle bone characters recognition are mainly based on deep learning algorithms. This task is more difficult than other handwritten characters recognition tasks since it exhibits several challenges. (1) Large intra-class variance. As oracle bone characters were carved by different ancients from multiple regions over several historic periods without much standardization [START_REF] Flad | Divination and power: A multi-regional view of the development of oracle bone divination in early China[END_REF], the same category could have several variants. Practically in China, characters were unified for the first time in Qin dynasty, after Shang dynasty. (2) Class imbalance. Due to the scarcity of the unearthed bones or shells, few instances were found for some characters. The distribution of the instances with respect to each category is quite imbalanced, following the long-tail distribution [START_REF] Huang | OBC306: A large-scale oracle bone character recognition dataset[END_REF]. (3) Low quality. Some inscriptions were damaged to some degree during the long burial period or the excavation process. Accordingly, there exist some noisy or broken character instances in the collected dataset [START_REF] Huang | OBC306: A large-scale oracle bone character recognition dataset[END_REF].

Several oracle bone character datasets have been built up to now for the research aim such as Oracle-20K [START_REF] Guo | Building hierarchical representations for oracle character and sketch recognition[END_REF] (20039 instances belonging to 261 categories), Oracle-AYNU [START_REF] Zhang | Oracle character recognition by nearest neighbor classification with deep metric learning[END_REF] (39062 instances belonging to 2583 categories) and OBC306 [START_REF] Huang | OBC306: A large-scale oracle bone character recognition dataset[END_REF] (309551 instances belonging to 306 categories) among which OBC306 has the largest scale and is closer to the real case as it uses the published oracle prints as the raw materials. As stated previously, the existing oracle datasets are quite imbalanced due to the scarcity of the unearthed materials. To overcome this limitation, in this work we propose a data augmentation method for oracle characters using Cycle-Consistent Adversarial Network (CycleGAN) [START_REF] Zhu | Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks[END_REF]. CycleGAN is trained to learn the mapping between the glyph images and the real samples to develop the capability of generating oracle images of high-quality.

The main contributions of this paper are as follows:

-We propose a feasible solution to generate oracle character images using CycleGAN by enabling it to learn the mapping between the glyph images collection and the real samples collection. Totally, 185362 samples are gener-ated which will be released later as a supplementary to the existing dataset for boosting the research in this field. -To our knowledge, it is the first attempt to use the GAN-based methods to generate oracle character images. -With the generated samples, the state of the art results of recognition task on OBC306 are greatly improved in terms of mean accuracy and total accuracy.

The rest of this paper is organized as follows. The related works are reviewed shortly in Section II. Section III introduces the proposed method. The experiments and results are discussed in Section IV. Section V concludes the work and proposes several perspectives for future work.

2 Related works

Oracle bone characters recognition

The research on oracle bone characters recognition is quite limited. Before the deep learning era, the proposed methods [START_REF] Li | The coding principle and method for automatic recognition of jia gu wen characters[END_REF] [3] [START_REF] Meng | Two-Stage Recognition for Oracle Bone Inscriptions[END_REF] for oracles recognition mainly follows the traditional pipeline, being pre-processing, feature extraction and classification. The accuracy of the models relies much on features extracted. Researchers usually design features representing oracles based on graph theory and topology theory. The idea proposed in [START_REF] Li | The coding principle and method for automatic recognition of jia gu wen characters[END_REF] regards an oracle character as a non-directed graph and uses the topological features for classification. The method in [START_REF] Li | Recognition of inscriptions on bones or tortoise shells based on graph isomorphism[END_REF] firstly transforms the inscriptions to labelled graphs, then extracts the quasi-adjacency matrix of the graphs as features, finally perform the recognition using an algorithm based on graph isomorphism. [START_REF] Meng | Two-Stage Recognition for Oracle Bone Inscriptions[END_REF] takes the skeleton and line points as features which later are fed into a two-stage recognition process. Even though these manually designed features were reported to work well on specific datasets, they have a low generalization capability.

Owing to the advantage of deep convolutional neural networks (DCNNs) [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], feature representations could be learned from data automatically. These feature representations are robust, and have strong generalization capability. Guo et al. [START_REF] Guo | Building hierarchical representations for oracle character and sketch recognition[END_REF] develop a novel hierarchical representation which combines a Gabor-related low-level representation and a sparse-encoder-related mid-level representation for oracle bone characters. This hierarchical representation serves as a complementary to CNN-based models. One point should be noted that [START_REF] Guo | Building hierarchical representations for oracle character and sketch recognition[END_REF] discards oracle characters with less than 25 instances which means they skip from the problem of class imbalance. To overcome the limitations, [START_REF] Zhang | Oracle character recognition by nearest neighbor classification with deep metric learning[END_REF] propose a deep metric learning based classification model which uses a CNN to map oracle images into an Euclidean space and then carries out the classification by the Nearest Neighbour (NN) rule. This method achieves the state of the art results on the public dataset Oracle-20K and Oracle-AYNU. Almost at the same time, a large-scale oracle bone character dataset was released which consists of 309551 samples. However, this dataset also suffers the problem of class imbalance by nature. To solve the imbalance problem, [START_REF] Li | Mix-Up Augmentation for Oracle Character Recognition with Imbalanced Data Distribution[END_REF] proposed a mix-up augmentation strategy to synthesize samples in the feature space by leveraging information from both majority and minority classes. Experimental results verify the effectiveness of the proposed method. However, [START_REF] Li | Mix-Up Augmentation for Oracle Character Recognition with Imbalanced Data Distribution[END_REF] synthesize samples in the feature space which makes it difficult to evaluate the quality of the synthesized samples explicitly. Different from work [START_REF] Li | Mix-Up Augmentation for Oracle Character Recognition with Imbalanced Data Distribution[END_REF], we propose to generate samples directly for minority classes with GAN-based methods.

Generative Adversarial Networks (GANs)

GAN [START_REF] Goodfellow | Generative adversarial nets[END_REF] [10] was firstly proposed by Ian GoodFellow in 2014 and since then it has quickly become a hot research topic. GAN consists of two models, being a generator and a discriminator, which can be implemented by any form of differentiable system such as neural network. The generator aims to capture the distribution of real samples for new data generation while the discriminator struggles to discriminate generated data from the real data as accurately as possible [START_REF] Gui | A Review on Generative Adversarial Networks: Algorithms, Theory, and Applications[END_REF]. GANs are widely applied in the field of image processing and computer vision, such as image super-resolution, image synthesis and manipulation, etc. Many interesting variants of GANs have been proposed since 2014.

Conditional GANs (cGANs) [START_REF] Mirza | Conditional generative adversarial nets[END_REF] extend GANs to a conditional model by conditioning the generator and discriminator on some extra information. With cGANs, samples can be generated conditioning on class labels, texts, etc. In [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF], cGAN is used for image-to-image translation. The developed software is named as pix2pix. The generator aims to learn a mapping from the input image (instead of random noise) to the generated image. pix2pix successfully applied cGANs for image-to-image translation however the need for the aligned image pairs hinders its generalisation.

The limitation of pix2pix is well solved by Cycle-consistent GANs (Cycle-GAN) [START_REF] Zhu | Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks[END_REF] which could learn to translate an image from one domain to another without paired examples. By exploiting the property that translation should be cycle consistent, CycleGAN introduce a cycle consistency loss in addition to the typical adversarial loss.

Even though the limitation of pix2pix seems not to be a hinder in the context of our task and a paired image-to-image translation network can be trained to generate higher quality images, our task benefits a lot from the unpaired image-to-image training. Because being different from other image generation tasks which pursue only the high quality of images, our main goal is to improve the recognition accuracy, and it's mostly influenced by the diversity of the data. CycleGAN [START_REF] Zhu | Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks[END_REF] tries to learn to the mapping between two domain collections without aligned image pairs. This working mechanism is beneficial to the oracle character images generation task especially for minority classes which could utilize rich stroke variety information from majority classes. Therefore, we propose to generate oracle bone character images based on CycleGAN.

Method

To our knowledge, OBC306 is the existing largest dataset for oracle bone characters recognition. However, it suffers a severe problem of class imbalance due to the scarcity of the unearthed inscriptions. Our aim is to generate samples for minority classes in OBC306 using GAN-based methods. The vanilla GAN model [START_REF] Goodfellow | Generative adversarial nets[END_REF] generates samples without designated class labels. The generated sample could belong to any class. However, our target is to generate samples belonging to designated classes. One possible way is to train a GAN model for each class. Unfortunately, training one model for each class requires a large amount of samples which is in conflict with minority classes. This contradiction makes the task based on the vanilla GAN model be challenging. The propose of cGAN [START_REF] Mirza | Conditional generative adversarial nets[END_REF] makes it possible to generate samples conditioning on specific class labels. cGAN is trained to learn the mapping from latent code to the corresponding image which is especially complex for a domain like oracle characters which have complicated glyphs. Furthermore, it also requires considerable number of samples from each class. pix2pix [START_REF] Isola | Image-to-image translation with conditional adversarial networks[END_REF] is one successful application of cGAN in the image-to-image translation task. It uses aligned image pairs from two different domains (such as gray scale and color) to train the model to learn the mapping from one domain image to the other. The major limitation of pix2pix is its requirement for aligned image pairs. Even though this limitation seems not to be a hinder in the context of our task and a paired image-to-image translation network can be trained to generate higher quality images, our task benefits a lot from the unpaired image-to-image training. Because being different from other image generation tasks which pursue only the high quality of images, our main goal is to improve the recognition accuracy, and it's mostly influenced by the diversity of the data. CycleGAN [START_REF] Zhu | Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks[END_REF] tries to learn to the mapping between two domain collections without aligned image pairs. This working mechanism is beneficial to the oracle character images generation task especially for minority classes which could utilize rich stroke variety information from majority classes. Based on the analysis above, A CycleGAN-based data augmentation method for improving oracle bone characters recognition is proposed in this work. Two data domains are assumed as glyph images and real character images respectively.

Statistics of OBC306

OBC306 [START_REF] Huang | OBC306: A large-scale oracle bone character recognition dataset[END_REF] consists of 309551 instances belonging to 306 categories. There are 70 categories with more than 1000 samples, including 5 categories with more than 10000 images. 236 categories with less than 1000 samples, including 85 categories with less than 100 samples and 29 categories with only one sample. Figure 2 illustrates the distribution of data in OBC306. As can be seen, the distribution is quite imbalanced. We set 1000/category as the target of augmentation.

Extraction of the glyphs of oracle characters

We extract the glyph of oracle from the real samples. The process for extraction of the glyph is quite simple, being binarization, noise erasing and image inverse ∈ Y . The same case is for G 2 and D X . In order to constrain the learned mapping further, the consistency loss (Equation 1) is introduced in addition to the adversarial loss (Equation 2). The consistency loss term tries to make the model work as:

CycleGAN

x → G 1 (x) → G 2 (G 1 (x)) ≈ x and y → G 2 (y) → G 1 (G 2 (y)) ≈ y.
The full objective is shown as Equation 3.

L cyc (G 1 , G 2 ) = E x∼p data (x) [ G 2 (G 1 (x)) -x 1 ] + E y∼p data (y) [ G 1 (G 2 (y)) -y 1 ]
(1)

L GAN (G 1 , D Y , X, Y ) = E y∼p data (y) [log D Y (y)] + E x∼p data (x) [log (1 -D Y (G 1 (x))] (2)
L (G 1 , G 2 , D X , D Y ) = L GAN (G 1 , D Y , X, Y ) + L GAN (G 2 , D X , Y, X) + λL cyc (G 1 , G 2 ) (3)
As CycleGAN does not require paired samples, we randomly select 1000 real images and 1000 glyph images to train a CycleGAN model. Both generators and discriminators are implemented with neural networks. We use the same network architecture with literature [START_REF] Zhu | Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks[END_REF]. The generator consists of three convolution layers for extracting features from input images, nine residual blocks which could translate feature vectors from the source domain to the target domain, two deconvolution layers and one convolution layer for restoring low-level features from feature vectors and generating images following the distribution of the target domain. The architecture is presented in Figure 5. For the discriminator (Figure 6), we use five convolution layers among which the last one is used to produce a 1-dimensional output telling input images are real or fake.

Samples generation for categories with few instances

Our object is to generate instances for oracle characters with few real instances to solve the problem of class imbalance. 1000/category is the target of augmentation. However for some categories with few real instances, the corresponding glyph images are quite few also. To achieve the target, we enrich the glyph images set first via general augmentation methods such as rotation, noise generation, occlusion, etc. Then feed these glyph images to multiple trained CycleGAN models to generate oracle character images. For example, assume that we need to generate 900 instances for some category. First step is to get 30 glyph images via general augmentation methods if required and next is to feed these 30 glyph images to 30 trained CycleGAN models to yield 900 instances.

Experiments

Experiments on oracle character images generation

Settings As stated previously, we evaluate the proposed method on the OBC306 dataset. We randomly select 1000 real images and 1000 glyph images (unpaired) to train a CycleGAN model. The model is trained for 200 epochs with Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]. For the first 100 epochs, the learning rate is set to 0.0002; it linearly decays to zero over the next 100 epochs.

Metrics GAN-train and GAN-test as defined in [START_REF] Shmelkov | How good is my GAN[END_REF] are used as quantitative metrics to evaluate the quality of generated samples.

-GAN-train represents the accuracy of a classification model trained with the synthetic samples and evaluated on real samples. -GAN-test is the accuracy of a classification model trained with real data, and evaluated on synthetic data. -GAN-base denotes the accuracy of a baseline network trained and tested on real data [START_REF] Bissoto | The Six Fronts of the Generative Adversarial Networks[END_REF].

Via comparing GAN-train accuracy and GAN-base accuracy, GAN-test accuracy and GAN-base accuracy, we can check the quality of the generated samples. Ideally, both GAN-train and GAN-test should be close to GAN-base accuracy.

Results

We carry out samples augmentation for 236 categories with less than 1000 instances. Totally, there are 185362 samples generated. Figure 7 compares the generated samples with the real samples of the same category. From the perceptual point of view, there is still a discrepancy between the real ones and the generated ones. We also take an automatic quantitative measure on the generated samples using GAN-train and GAN-test. As CycleGAN models are trained with samples randomly selected, in other words, with unbalanced data samples, we question that if the quality of generated samples is different or not among categories with variable numbers of real instances. Thus, experiments are carried out respectively on categories with different numbers of real instances. Specifically, we select 5 categories which have around 500 real instances to conduct Experiment A and 5 categories which have 10-100 real instances to conduct Experiment B. From Table 1, we can see in Experiment A, GAN-train is relatively close to GAN-base (87.4% VS 83.2%) which tells that CycleGAN models generate for the selected 5 categories high-quality images with similar diversity to the real data. In Experiment B, GAN-train is 9.5% less than GAN-base (81.6% VS 72.1%) which means the generated samples for the selected categories lost some patterns compared to the real data. For both cases, GAN-test is not far away from GAN-base and is not larger than GAN-base telling that the trained CycleGAN models do not have the problem of over-fitting. In conclusion, CycleGAN tends to generate samples of higher quality for categories with a lot of real instances.

Experiments on oracle characters recognition

Settings We adopt AlexNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], ResNet50 [START_REF] He | Deep residual learning for image recognition[END_REF] and Inception-v4 [START_REF] Szegedy | Inception-v4, Inception-ResNet and the Impact of Residual Connections on Learning[END_REF] which have shown good performance on OBC306 as the classification models. For each category, 80% of the samples are used for training, and the remaining 20% are used for test. Classifiers are trained with Adam optimizer also.

To evaluate the contribution of generated samples to recognition task, we set three experiments. The three experiments use the same test data. In terms of training data, Experiment I uses samples from the original OBC306 dataset only.

Experiment II uses samples from the original OBC306 dataset + samples generated with traditional methods. Experiment III uses samples from the original OBC306 dataset + samples generated via CycleGAN. The numbers of training data for Experiment II & III are the same, larger than the number of training data in Experiment I .

Metrics To account for long-tail effects, we measure the performance of the model with regard to multiple metrics as same as [START_REF] Li | Mix-Up Augmentation for Oracle Character Recognition with Imbalanced Data Distribution[END_REF].

-Mean accuracy denotes the accuracy that is calculated by averaging the accuracy of per class. -Total accuracy refers to the overall accuracy of all the test samples on all classes. ). [START_REF] Li | Mix-Up Augmentation for Oracle Character Recognition with Imbalanced Data Distribution[END_REF] proposed a mix-up augmentation strategy to synthesize samples by fusing features linearly from two examples. This method synthesizes samples in the feature space which makes it difficult to evaluate the quality of the synthesized samples explicitly. Furthermore, our method outperforms the mix-up strategy in terms of both mean and total accuracy. While their total accuracy is improved by around 2% (89.81% → 91.56%). our total accuracy is improved by around 6% (86.41% → 92.02%).

Here is one point we need to clarify. For 29 classes with only one sample, we put them in the test set. Thus, there is no data from these classes to train CycleGAN models. But, we feed glyph images of these classes into trained Cy-cleGANs to generate samples (1000 per class). These operations maybe influence the mean accuracy a lot. However, the total accuracy will not be affected much as there are around 60000 samples in the test set. Furthermore, the total accuracy is improved by around 6% which verifies the effectiveness of the proposed method.

To further analyse the effects of generated samples on the recognition accuracy of classes with variable numbers of real instances, we design a new recognition experiment. (1) Firstly randomly select 5 categories which have 1000-1100 real instances. (2) Secondly remove a part of the real instances to make these 5 categories be minority classes. (3) Use the traditional augmentation method and the CycleGAN-based method respectively to recover the original quantity. (4) ResNet50 is adopted as the classifier. Figure 8 shows the recognition results under different settings. As can be seen, when 90% real data removed, CycleGAN-based method exhibits significantly better performance than the traditional method. When 70% real data removed, CycleGAN-based method slightly outperforms the traditional method. When 50% real data removed, there is no obvious discrepancy between the two methods. Form this experiment, we could tell that the CycleGAN-based augmentation method has a good performance for minority classes. For classes with around 500 real instances, only a slight improvement is observed. The result of the case using 50% real data is close to the case using 100% real data which means the classifier learns the patterns of oracle characters with around 500 samples each category.

Conclusion

In this work, we propose a CycleGAN-based data augmentation method to solve the problem of class imbalance in the existing oracle characters dataset. Totally, 185362 samples are generated which could be a supplementary to the existing dataset. With these generated samples, the SOTA results of recognition task on OBC306 are improved greatly in terms of mean accuracy and total accuracy. However, some limitations exist in our work. Currently, we extract the glyph from real samples only and take some simple transformations if required. In fact, during the training process of CycleGAN, the features learned from real samples also guide the structure of glyph. Thus, extracting the glyph from generated samples could be an interesting attempt in the future. Another point is there is no attempt to solve the problem of large intra-class variance which will be left to the future.

Fig. 1 :

 1 Fig. 1: Examples of unearthed oracle bone inscriptions.

Fig. 2 :

 2 Fig. 2: Illustration of data distribution in OBC306.

Fig. 3 :

 3 Fig. 3: Some examples of the extracted glyph.

Fig. 4 :Fig. 5 :

 45 Fig. 4: The framework of CycleGAN.

Fig. 6 :

 6 Fig. 6: The network architecture of the discriminator used in this work.

Fig. 7 :

 7 Fig. 7: Comparison of the real samples (a) and the generated samples (b).

Fig. 8 :

 8 Fig. 8: Results of oracle characters recognition on specific classes.

Fig. 9 :

 9 Fig. 9: Illustration of samples generated via CycleGAN. The samples in a row belong to the same category.

Table 1 :

 1 Results of GANtrain and GANtest experiments.A: experiment on categories with around 500 real instances; B: experiment on categories with 10-100 real instances

	Experiment GAN-base GAN-train GAN-test
	A	87.4%	83.2%	84.3%
	B	81.6%	72.1%	76.8%

Table 2 :

 2 Results of oracle characters recognition on the OBC306 dataset. The results of these two rows are extracted directly from the published literature[START_REF] Li | Mix-Up Augmentation for Oracle Character Recognition with Imbalanced Data Distribution[END_REF]. They remove 29 classes which have only one instance and carry out the experiments on the remaining 277 classes. The rest of the results in this table are produced by us.

	Experiment	Model	Mean accuracy Top 1 Top 5 Top 10	Total accuracy
		AlexNet [8]	65.24% 74.22% 79.15%	-
	I: Original OBC306	ResNet50 [8] Inception-v4 [8] 69.51% 78.11% 81.16% 86.41% 68.77% 78.67% 81.32% -
		Inception-v4 (277 classes) * 75.62%	-	-	89.81%
		AlexNet	66.41% 75.24% 80.10%	-
	II: Traditional augmentation	ResNet50	71.89% 80.27% 81.74%	-
		Inception-v4	74.61% 82.92% 84.23% 89.92%
	III: CycleGAN-based augmentation	AlexNet ResNet50 Inception-v4 83.44% 91.39% 93.13% 92.02% 70.17% 77.22% 82.71% -77.63% 87.15% 91.83% -
	IV: Mix-up augmentation (277 classes) *	Inception-v4 [19] 79.02%	-	-	91.56%
	Results Table 2 presents recognition accuracies in different experiments. From
	the statistics, we can see the samples generated by CycleGAN could improve
	the recognition accuracy dramatically: mean accuracy by approximately 14%
	(69.51% → 83.44%), total accuracy by approximately 6% (86.41% → 92.02%).
	However, the samples generated by traditional methods only improve the mean
	accuracy by around 5% (69.51% → 74.61%) and total accuracy by around 3%
	(86.41% → 89.92%

*

This work is supported by National Natural Science Foundation of China (No. 62007014), China Post doctoral Science Foundation (No. 2019M652678) and the Fundamental Research Funds for the Central Universities(No. CCNU20ZT019).