Jean Gabriel Houot

Jorge San Martin

Marius Tucsnak
email: marius.tucsnak@iecn.u-nancy.fr

Existence of solutions for the equations modeling the motion of rigid bodies in an ideal fluid

In this paper, we study the motion of rigid bodies in a perfect incompressible fluid. The rigid-fluid system fils a bounded domain in R 3 . Adapting the strategy from Bourguignon and Brezis [1], we use the stream lines of the fluid and we eliminate the pressure by solving a Neumann problem. In this way, the system is reduced to an ordinary differential equation on a closed infinite dimensional manifold. Using this formulation, we prove the local in time existence and uniqueness of strong solutions.

Notation. Throughout this paper Ω denotes an open bounded and connected subset of R 3 and S 0 is a closed set with nonempty interior and with smooth boundary such that S 0 ⊂ Ω. We denote as usual by SO 3 (R) the special orthogonal group on R 3 . We will often use functions defined from a time interval to R 3 or to SO 3 (R). these functions will be denoted using bold characters, such as h : [0, T ] → R 3 or R : [0, T ] → SO 3 (R). The same kind of notation will be used for three other time dependent vector fields k, ω, η and ξ which will be defined in the sequel. The five time dependent fields mentioned above will define the state z of the fluid-solid system. A vector from R 3 or a matrix from SO 3 (R) will be denoted by h or by R, respectively. The transposed of a matrix will be denoted by * so that the column vector of components a and b is denoted either a b or by (a, b) * . Differentiation with respect to time is often denoted a dot.

The vector, respectively the inner, product of v, w ∈ R 3 will be denoted by v ∧ w and v • w, respectively . The Jacobian matrix of a vector field y → f (y) defined on an open subset of R 3 will be denoted by D y f or simply by Df .

Introduction

The interaction of rigid bodies and of ideal fluids is a topic which was probably first tackled by d'Alembert, Kelvin and Kirchhoff, who considered the case of a potential fluid (sometimes called inviscid fluid), with the solid-fluid system filling the whole space. In this case the governing equations can be reduced to system of ODE's on a finite dimensional manifold. We refer to the book of Lamb [9, chapter 6] for a detailed presentation of this early contributions and to Kanso, Marsden, Rowley and Melli-Huber [START_REF] Kanso | Locomotion of articulated bodies in a perfect fluid[END_REF] for the application of the above theory to self-propelled motions of solids in an inviscid fluid. Recently Houot and Munnier in [START_REF] Houot | On the motion and collisions of rigid bodies in an ideal fluid[END_REF] used shape sensitivity analysis techniques to deal with either bounded or unbounded domains. They also tackled the special case of a cylinder in a half space. They showed in particular that, unlike the case of a viscous fluid (see San Martin, Starovoitov and Tucsnak [START_REF] San Martín | Global weak solutions for the two-dimensional motion of several rigid bodies in an incompressible viscous fluid[END_REF], Hillairet [START_REF] Hillairet | Lack of collision between solid bodies in a 2D incompressible viscous flow[END_REF], Hesla [START_REF] Hesla | Collisions of smooth bodies in viscous fluids, a mathematical investigation[END_REF]), the cylinder can touch the wall in finite time with non zero velocity. The damping effect of the wall on the cylinder is also studied.

In the general case the system is genuinely infinite dimensional, so it cannot be reduced to ODE's on finite dimensional manifolds. As usual in fluid-solid interaction problems, a major difficulty comes from the fact that the equations for the fluid (Euler's equations in our case) hold in a time dependent domain, so that we have a free boundary value problem. As far as we know, the first papers tackling the case of a non potential flow are Ortega, Rosier and Takahashi [START_REF] Ortega | Classical solutions for the equations modelling the motion of a ball in a bidimensional incompressible perfect fluid[END_REF] and [START_REF]On the motion of a rigid body immersed in a bidimensional incompressible perfect fluid[END_REF]. The main result in these works asserts the existence and uniqueness of classical solutions in two space dimensions and with the rigid-fluid system filling the whole space. More recently, Rosier and Rosier in [START_REF] Rosier | Smooth solutions for the motion of a ball in an incompressible perfect fluid[END_REF] proved the existence of strong solutions in the case in which the solid is a ball, with the fluid-rigid system filling R n , with n 2.

The aim of the present work is to prove the existence an uniqueness of strong solutions in three space dimensions, with a bounded fluid-rigid domain and with the possibility of considering more than one solid. An idea which seems attractive, since it yields a transformed problem written in a fixed domain, is the use of groups of diffeomorphisms as proposed in Ebin and Marsden [START_REF] Ebin | Groups of diffeomorphisms and the notion of an incompressible fluid[END_REF]. Our approach, based on this idea, follows more closely Bourguignon and Brezis [START_REF] Bourguignon | Remarks on the Euler equation[END_REF]. The first new difficulty we need to tackle is that, the fluid domain being variable and the normal velocity of the fluid being different from zero on the fluid-solid interface we are not able to apply the Leray projector. Therefore, in order to eliminate the pressure we need to solve non-homogeneous Neumann problems for the Laplacian. The second difficulty consists in the fact that we need to compare solutions of these Neumann problems in different domains and to show that they depend smoothly on some geometric parameters.

To be more precise, the motion of the fluid is described by the classical Euler equations, whereas the motion of the rigid bodies is governed by the balance equations for linear and angular momentum (Newton's laws). For the sake of simplicity we state and prove our results in the case of a single rigid body, but our methods can be easily be adapted to the case of several rigid bodies. Assume that the system fluid-rigid body fills the domain Ω in R 3 and that at t = 0 the solid is located at S 0 (see the paragraph on notation from the beginning of the paper for the properties of Ω and S 0 ). The position of the solid at instant t 0 is denoted by S(t). We assume that the solid is surrounded by a perfect homogeneous incompressible fluid filling, for each t 0, the domain F (t) = Ω \ S(t). In this work we study the following initial and boundary value problem:

ρ F ∂u ∂t + ρ F (u • ∇)u + ∇p = 0 x ∈ F (t), t 0, (1.1a) div u = 0 x ∈ F (t), t 0, (1.1b) u • n = 0 x ∈ ∂Ω, t 0, (1.1c) u • n = ( ḣ + ω ∧ (x -h)) • n x ∈ ∂S(t), t 0, (1.1d) m s ḧ = ∂S(t)
pn dx, t 0, (1.1e)

d dt (Jω) = ∂S(t) p(x -h) ∧ n dx, t 0, (1.1f) Ṙ(t) = A(ω(t))R(t), t 0, (1.1g) u(0, x) = u 0 (x)
x ∈ F (0) , (1.1h)

h(0) = h 0 , ḣ(0) = k 0 , R(0) = Id M 3 , ω(0) = ω 0 , (1.1i)
where the unknowns are u (the Eulerian velocity field of the fluid), p (the pressure of the fluid), h (the trajectory of the mass center of the rigid body), R (the time variation of the orthogonal matrix giving the orientation of the solid) and ω (the time variation of the angular velocity of the rigid body). The density of the fluid ρ F is supposed to be a constant. The fluid occupies, at t = 0, the domain F 0 = Ω \ S 0 . The domain S(t) is defined by

S(t) = {h(t) + R(t)(y -h 0 ) | y ∈ S 0 , t 0} .
The domain occupied by the fluid at instant t is F (t) = Ω \ S(t). The skewsymmetric matrix A(ω) is given by

A(ω) =   0 -ω 3 ω 2 ω 3 0 -ω 1 -ω 2 ω 1 0   (ω ∈ R 3 ). (1.2)
The notation m s stands for the mass of the solid and J(t) designs its inertia matrix defined by

J i,j (t) = ρ s S(t) [(x -h(t)) ∧ e i ] • [(x -h(t)) ∧ e j ] dx (i, j ∈ {1, 2, 3}), (1.3)
where the constant ρ s stands for the density of the solid and (e k ) k=1,2,3 is the canonical basis in R 3 . It is easy to check that J(t) = R(t)J 0 R * (t) for every t 0, where J 0 is the matrix defined by

(J 0 ) i,j = ρ S S 0 [(y -h 0 ) ∧ e i ] • [(y -h 0 ) ∧ e j ] dy, (1.4)
for every i, j ∈ {1, 2, 3}. Notice that the matrix J 0 does not depend on the position of the solid and that the last formula easily implies that

d dt (Jω) = ω ∧ (Jω) + J ω. (1.5)
Moreover, we have denoted by ∂S(t) the boundary of the rigid body at instant t and by n(t, x) the unit normal to ∂S(t) at the point x directed to the interior of the rigid body.

Throughout this paper we assume that the considered boundaries are smooth in the sense that there exist the functions δ

0 , δ 1 ∈ C ∞ (R 3 , R) such that ∂Ω = x ∈ R 3 | δ 0 (x) = 0 , ∂S(0) = x ∈ R 3 | δ 1 (x) = 0 , (1.6) n(x) = -∇δ 0 (x), x ∈ ∂Ω, n(x) = -∇δ 1 (x), x ∈ ∂S(0). ( 1.7)
An important role in this work will be played by the set P (Ω, S 0 ), defined as follows:

Definition 1.1. The set of all admissible solid configurations from the solid position S 0 , denoted P (Ω, S 0 ), is the set of all pairs

h 1 R 1 ∈ R 3 × SO 3 (R) such that there exist functions h ∈ C([0, 1]; R 3 ), R ∈ C([0, 1]; SO 3 (R)), with h(0) = h 0 , h(1) = h 1 , R(0) = Id 3 , R(1) = R 1 , h(t) + R(t)(y -h 0 ) ∈ Ω (t ∈ [0, 1], y ∈ S 0 ).
Remark 1.2. For each t 0 the position of the solid and the domain filled by the fluid are completely described by the pair (h(t), R(t)) * ∈ P (Ω, S 0 ). Therefore, the evolution of the domains F (t) and S(t) is totally described by the function q ∈ C 2 ([0, T ], P(Ω, S 0 )) defined by q(t) = (h(t), R(t)) * . Consequently, in the remaining part of this work, we use the notation F q(t) and S q(t) instead of F (t) and S(t). We also denote q 0 = (h 0 , Id M 3 ) * = q(0). More generally, for every q = (h, R) * ∈ P (Ω, S 0 ) we denote

S q = {h + R(y -h 0 ) | y ∈ S 0 }, F q = Ω \ S q . (1.8)
In order to give a precise statement of our main result we first introduce some spaces. For an open set O ⊂ R 3 we denote

N m (O) = q ∈ H m (O) | O q(x)dx = 0 .
(1.9)

We next defines some spaces of functions defined on time variable domains. Let q ∈ C 2 ([0, ∞), P(Ω, S 0 )) and let

Ψ : [0, ∞) × R 3 → R 3 be such that for every t ∈ [0, T ], the map x → Ψ(t, x) is a C ∞ diffeomorphism from F 0 to F q(t) .
Let v(t, •), t ≥ 0 be a family of functions with v(t, •) : F q(t) → R 3 for every t 0. Denote v Ψ (t, y) = v(t, Ψ(t, y)), for all t 0 and for all y ∈ F 0 . With the above notation we introduce the following function spaces:

C k ([0, T ], H m (F q )) = {v | v Ψ ∈ C k ([0, T ], H m (F 0 ))}, C k ([0, T ], N m (F q )) = {v | v Ψ ∈ C k ([0, T ], N m (F 0 )},
where k ∈ {0, 1}, m 0 is an integer and H m are the usual Sobolev spaces. It is not difficult to check that the above definitions are independent of the choice of the diffeomorphism Ψ.

We can now state the main result in this paper.

Theorem 1.3. Let m 3 be an integer. Let S 0 ⊂ Ω be as in the notational preamble of this work and let

h 0 = 1 vol (S 0 ) S 0 x dx. Let k 0 ∈ R 3 , ω 0 ∈ R 3 and u 0 ∈ H m (F 0 , R 3 ) satisfy: div u 0 (x) = 0 x ∈ F 0 , u 0 (x) • n(x) = (k 0 + ω 0 ∧ (x -h 0 )) • n(x) x ∈ ∂S 0 , u 0 (x) • n(x) = 0 x ∈ ∂Ω.
Then there exists T 0 > 0 such that (1.1) admits a unique solution (q, u, p) with

q ∈ C 2 [0, T 0 ), P (Ω, S 0 ) , (1.10) u ∈ C [0, T 0 ), H m (F q ) ∩ C 1 [0, T 0 ), H m-1 (F q ) , (1.11) p ∈ C [0, T 0 ), N m+1 (F q ) .
(1.12)

2 Idea of the proof of Theorem 1.3

As already mentioned, the basic idea of the proof, borrowed from Bourguignon and Brezis [START_REF] Bourguignon | Remarks on the Euler equation[END_REF], consists in reducing (1.1) to an ODE on an infinite dimensional manifold. In this section we briefly describe this reduction process and we give the main steps of the proof of Theorem 1.3. Let q = h R and u be functions satisfying (1.10) and (1.11) for some T 0 > 0, with div u = 0. We introduce the flow η associated to u, which is defined as the solution of ∂η ∂t (t, y) = u(t, η(t, y)), η(0, y) = y for all y ∈ F 0 .

( Notice that u can be expressed in terms of η and ξ by

u(t, x) = ξ(t, η -1 (t, x)) (t ∈ [0, T 0 ), y ∈ F 0 ). (2.3)
In order to express (1.1) as a first-order ordinary differential equation we note that from the formula

∂ξ ∂t (t, y) = ∂u ∂t (t, η(t, y)) + (u • ∇)u(t, η(t, y)), it follows that u satisfies (1.1a) iff ∂ξ ∂t (t, y) = -∇p(t, η(t, y)), (y ∈ F 0 , t ∈ [0, T 0 )). (2.4)
Consider the function k ∈ C 1 ([0, T 0 ), R 3 ) defined by

ḣ(t) = k(t) (t ∈ [0, T 0 )). (2.5) Define ω ∈ C 1 ([0, T 0 ), R 3 ) by Ṙ(t) = A(ω(t))R(t) (t ∈ [0, T 0 )). (2.6)
As it will be shown in Sections 3 and 4, by solving appropriate Neuman problems, the pressure p can be expressed, for each t ∈ [0, T 0 ) as a function of z = (η, q, ξ, k, ω) * , so that, using (2.4), (2.5) and (2.6), the system (1.1) can be written in the equivalent form

ż(t) = L(z(t)), z(0) = z 0 ,
where L :

F m → E m , with E m = H m (F 0 , R 3 ) × R 3 × M 3 (R) × H m (F 0 , R 3 ) × R 3 × R 3 , (2.7)
and F m is a closed subset of E m . The norm of z ∈ E m is defined by

z 2 E m = η 2 H m (F 0 ,R 3 ) + h 2 + R 2 + ξ 2 H m (F 0 ,R 3 ) + k 2 + ω 2 ,
where • stand for the Euclidean norm on R n . Endowed with this norm E m is a Hilbert space. For (q, u, p) satisfying (1.10)-(1.12) we define To define F m we introduce, for every q ∈ P (Ω, S 0 ), the sets

z(t) =       η(t, •) q(t) ξ(t, •) k(t) ω(t)       , ( 2
Diff m (F 0 , F q ) = η : F 0 → F q η invertible, η ∈ H m (F 0 , R 3 ), η -1 ∈ H m (F q , R 3 ) and det[D y (η(y))] = 1 , (2.9) Σ m (Ω, S 0 ) = σ = η q q ∈ P (Ω, S 0 ) and η ∈ Diff m (F 0 , F q ) , (2.10)
where P (Ω, S 0 ) has been defined in Definition 1.1 and F q is given in (1.8). The set Σ m (Ω, S 0 ), simply denoted by Σ m in the sequel, is formed by the admissible positions of the system. The set of admissible velocities from a position σ = q η describes the tangent space to Σ m at the point σ, which is given by

T σ Σ m = (ξ, k, ω) * ∈ H m (F 0 , R 3 ) × R 3 × R 3 u = ξ • (η -1 ) ∈ H m (F q , R 3 ), div u = 0 in F q , u • n = 0 on ∂Ω, u • n = [k + ω ∧ (x -h)] • n on ∂S q } . (2.11)
The subset F m of E m is defined by

F m = {z ∈ E m | σ = (η, q) * ∈ Σ m , (ξ, k, ω) * ∈ T σ Σ m }.
(2.12)

It is not difficult to check that F m is a locally closed subset of E m , in the sense that for every z ∈ F m there exists a closed ball B of F m centered at z such that F m ∩ B is a closed subset of E m .
Moreover, as it will be shown in Section 6, Σ m is an infinite-dimensional manifold and F m is its tangent bundle.

The precise definition of L requires some preparation, so it is postponed to Sections 3 and 4. In order to prove the main result we show in Section 5 that L satisfies the assumptions of the following version of the Cauchy-Lipschitz theorem, which is a particular case of Theorem 2 from Martin [START_REF] Martin | Differential equations on closed subsets of a Banach space[END_REF].

Proposition 2.1. Let F be a locally closed subset of a Hilbert space E and let

L : [0, T ) × F → E be such that a) L is a locally Lipschitz in z and continuous in t; b) lim s→0 + 1 s dist z + sL t z ; F = 0 t z ∈ [0, T ) × F .
Then for every z 0 ∈ F there exists T 0 > 0 such that the equation

ż(t) = L(t, z(t)), z(0) = z 0 admits a unique solution z ∈ C 1 ([0, T 0 ), F ).
3

Study of the pressure

The study of the pressure p is the key point in order to reduce (1.1) to a system of ordinary differential equations. In this section we write the pressure as the sum of two terms, each of them satisfying a Neumann problem for the Laplacian. We first introduce some function spaces and we recall classical results related to Neumann problems. Let O be a bounded domain in R 3 with smooth boundary ∂O, m ∈ N and n ∈ N. Recall the definition of N m (O) from (1.9) and let V m (O) the space defined by

V m (O) = (f, g) * ∈ H m (O) × H m+1/2 (∂O) O f (x)dx + ∂O g(x)dσ x = 0 . (3.1)
The following classical result on the wellposedness of the Neumann problem for the Laplace operator can be found in the book of Lions and Magenes [START_REF] Lions | Problèmes aux limites non homogènes[END_REF]Chapter 5].

Theorem 3.1. Let m ∈ N. Then, for every f g ∈ V m (O), the boundary value problem -∆ϕ(x) = f (x) x ∈ O, ∂ϕ ∂n (x) = g(x) x ∈ ∂O, admits a unique solution ϕ ∈ N m+2 (O). Moreover, ϕ satisfies O ∇ϕ • ∇ψ dx = O f ψ dx + ∂O gψdσ x (ψ ∈ H m+2 (O)), (3.2)
and there exists a constant C (depending only on O and m) such that

||∇ϕ|| H m+1 (O) ≤ C(||f || H m (O + ||g|| H m+1/2 (∂O) ) f g ∈ V m (O .
In order to prove that the boundary value problem for the pressure is well posed, we need several technical results. Let q ∈ P (Ω, S 0 ). We first note that, thanks to the smoothness of ∂F q , the map x → n(x), defined on ∂F q , can be extended to F q by a function in H m (F q ). This extension is not unique so that the partial derivatives of n on ∂F q are not uniquely determined. However, it can be easily checked that for every vector field τ which is tangent to ∂F q , the quantity 3 j=1 τ j ∂n i ∂x j , with i ∈ {1, 2, 3} does not depend on the choice of the extension. Proposition 3.2. Let m ≥ 3 be an integer, let q ∈ P (Ω, S 0 ) and assume that

w ∈ H m (F q , R 3 ), w • n = 0 on ∂F q . Then the function x → i,j ∂w j ∂x i ∂w i ∂x j , is in H m-1 (F q ) whereas the function x → i,j w i w j ∂n i ∂x j , is in H m-1/2 (Γ), where Γ is either ∂Ω or ∂S q .
Proof. The first property follows from the fact that, under our assumptions, H m-1 (F q ) is an algebra.

To prove the second property we notice that it suffices to use the fact that H m (F q ) is an algebra, the smoothness of the map x → ∂n i ∂x j (x) defined on F q and the trace theorem.

The above lemma allows us to introduce, for every q ∈ P (Ω, S 0 ), the operators:

F q : H m (F q , R 3 ) → H m-1 (F q ), F q (u) = i,j ∂u j ∂x i ∂u i ∂x j , ( 3.3)
G q : H m (F q , R 3 ) → H m-1/2 (Γ), G q (u) = i,j u i u j ∂n i ∂x j | Γ , (3.4)
where Γ is either ∂Ω or ∂S q . An important ingredient allowing to write (1.1) as an ordinary differential equation is the following result:

Proposition 3.3. Let T 0 > 0, let m ≥ 3 be an integer, let h R = q ∈ C 2 [0, T 0 ], P(Ω, S 0 ) , u ∈ C [0, T 0 ), H m (F q ) ∩ C 1 [0, T 0 ), H m-1 (F q ) , p ∈ C [0, T 0 ), N m+1 (F q ) .
Assume that u satisfies

(div u)(t, x) = 0 x ∈ F q(t) , t ∈ [0, T 0 ), (u • n)(t, x) = 0 x ∈ ∂Ω, t ∈ [0, T 0 ), (u • n)(t, x) = v(t, x) • n(t, x) x ∈ ∂S q(t) , t ∈ [0, T 0 ), where v(t, x) = ḣ(t) + ω(t) ∧ (x -h(t)), for all x ∈ F q(t) , t ∈ [0, +∞). (3.5)
Moreover, assume that u, p and q satisfy (1.1a). Then, for very t ∈ [0, T 0 ), we have

-∆p(t, x) = ρ F F q(t) (u)(t, x) (x ∈ F q(t) ), (3.6
)

∂p ∂n (t, x) = ρ F G q(t) (u)(t, x) (x ∈ ∂Ω), (3.7) ∂p ∂n (t, x) = ρ F G q(t) (u -v)(t, x) + 2ρ F (u -v) • (ω(t) ∧ n(t, x)) -ρ F ḧ(t) + ω(t) ∧ (x -h(t)) + ω(t) ∧ (ω(t) ∧ (x -h(t))) •n(t, x) (x ∈ ∂S q(t) ), (3.8)
where F q and G q are defined by (3.3) and (3.4) and v stands for the velocity of the solid defined by (3.5).

Proof. Assume that u, p, q satisfy (1.1a). By applying the div operator to (1.1a) we get that p satisfies, for every t ∈ [0, T 0 ),

-∆p(t, x) = ρ F div ∂u ∂t (t, x) + (u(t, x) • ∇)u(t, x) (x ∈ F q(t) ).
(3.9)

By using the fact that div u ≡ 0, the right-hand side of the above relation can be expressed as

div ∂u ∂t (t, x) + (u • ∇)u(t, x) = i,j ∂u j ∂x i ∂u i ∂x j (t, x) + u.∇div (u)(t, x) = i,j ∂u j ∂x i ∂u i ∂x j (t, x).
The above formula and (3.9) imply (3.6).

On the other hand, by taking normal traces of all the terms in (1.1a) we obtain

∂p ∂n (t, x) = ρ F - ∂u ∂t (t, x) -(u • ∇)u(t, x) • n(t, x) (x ∈ Γ), (3.10)
where Γ = ∂Ω or Γ = S q(t) . The above boundary conditions can be expressed in terms of the velocity and of the position of the solid. First note that

n(t, x) = n(0, x) (x ∈ ∂Ω, t ∈ [0, T 0 )).
Additionally, note that, for every t ∈ [0, T 0 ) and y ∈ ∂S 0 , we have

n(t, Ψ(t, y)) = R(t)n(0, y),
and u(t, Ψ(t, y)) • (R(t)n(0, y)) = v(t, Ψ(t, y)) • (R(t)n(0, y)), (3.11)
where

x = Ψ(t, y) = h(t) + R(t)(y -h 0 ) (y ∈ ∂S 0 ),
and v is the solid velocity given in (3.5). By taking the derivative with respect to t of the two sides of (3.11), we obtain that for every t ∈ [0, T 0 ) and x ∈ ∂S q(t) we have:

∂u ∂t (t, x) + [(v • ∇)u(t, x)] • n(t, x) + u(t, x) • (ω(t) ∧ n(t, x)) = ∂v ∂t (t, x) + [(v • ∇)v(t, x)] • n(t, x) + v(t, x) • (ω(t) ∧ n(t, x)). (3.12)
Using in the above formula the fact (easy to check) that

[(v(t, x) • ∇)v(t, x)] • n(t, x) = -v(t, x) • (ω(t) ∧ n(t, x)),
we obtain that for every t ∈ [0, T 0 ) and x ∈ ∂S q(t) we have:

∂u ∂t + (u • ∇)u • n (t, x) = [(u -v) • ∇](u -v) • n(t, x) -2(u -v) • (ω ∧ n(t, x)) + (v • ∇)v + ∂v ∂t • n(t, x) (3.13)
Using again the relation (u -v) • n = 0 on S q(t) , we have To obtain (3.7) it suffices to apply (3.10) with v = 0 (so that ω = 0).

[(u -v) • ∇](u -v) • n(t, x) = -G q(t) (u -v)(t, x). ( 3
From Proposition 3.3 (more precisely from (3.8)) we note that the pressure depends on ḧ and on ω. In order to make this dependence more precise we introduce, for every q ∈ P (Ω, S 0 ), the potential functions Φ i for i = 1, . . . , 6 which are solutions of the Neumann problems:

-∆Φ i (q; x) = 0 x ∈ F q , (3.15a) ∂Φ i ∂n (q; x) = 0 x ∈ ∂Ω, (3.15b) ∂Φ i ∂n (q; x) = K i (q; x) x ∈ ∂S q , ( 3.15c)
where

K i (q; x) = n i (x) for i = 1, 2, 3, K i (q; x) = [(x -h) ∧ n(x)] i-3 for i = 4, 5, 6.
(3.16)

Denote Φ = (Φ 1 , • • • , Φ 6 ) * .
These functions have been introduced in the book of Lamb [START_REF] Lamb | Hydrodynamics[END_REF] and they were used, in particular, in the work of Houot and Munnier [START_REF] Houot | On the motion and collisions of rigid bodies in an ideal fluid[END_REF] to describe the motion of rigid bodies in a perfect fluid undergoing a potential flow. From Theorem 3.1 on the Neumann problem, it is easy to check that Φ ∈ C ∞ (F q ; R 6 ). Moreover, the following properties are proved in [START_REF] Houot | On the motion and collisions of rigid bodies in an ideal fluid[END_REF].

Proposition 3.4. For every q 0 ∈ P (Ω, S 0 ), there exists a neighborhood O of q 0 of such that

• the mapping q → Φ(q; •) from O to C ∞ (F q ; R 6 ) is of class C 2 ;
• for all i, j ∈ {1, . . . , 6} the mappings

q → I i,j (q) = Fq ∇Φ i (q; x) • ∇Φ j (q; x)dx, are of class C 2 on O.
We also need a potential µ, defined as follows. For z = (η, q, ξ, k, ω) * ∈ F m , where F m has been defined in (2.12), we set u(z; x) = ξ(η -1 (x)), with x ∈ F q . The potential µ is defined as the solution of the boundary value problem

-∆µ(z; x) = F q (u(z; x)) (x ∈ F q ), (3.17a) ∂µ ∂n (z; x) = G q (u(z; x)) (x ∈ ∂Ω), (3.17b) ∂µ ∂n (z; x) = G q (u -v)(x) + 2(u(z; x) -v(z; x)) • (ω ∧ n(x)) -[ω ∧ (ω ∧ (x -h))] • n(x) (x ∈ ∂S q ), (3.17c)
where

v(z; x) = k + ω ∧ (x -h),
and F q , G q are defined in (3.3), (3.4).

Remark 3.5. With the above notation for Φ and µ, if (u, p, q) satisfy (3.6)-(3.8) and z(t) is defined by (2.8), then the pressure can be written

p(z(t); x) = ρ F µ(z(t); x) -ρ F Φ(q(t); x) • ( ḧ(t), ω(t)) * , (3.18)
where • stands for the inner product in R 6 .

An equivalent form of the governing equations

Throughout this section we assume that m 3 and

q ∈ C 2 ([0, T ), P(Ω, S 0 )), u ∈ C([0, T ), H m (F q(•) ) ∩ C 1 ([0, T ), H m-1 (F q(•) ), p ∈ C([0, T ), N m+1 (F q(t) )).
At this point we need the virtual mass of the solid (see, for instance, [START_REF] Houot | On the motion and collisions of rigid bodies in an ideal fluid[END_REF]) which is the six by six matrix K(q) defined, for every for every q ∈ P (Ω, S 0 ), by

K(q) = K S (q) + K F (q), K S (q) = m s Id 3 0 0 J , K F (q) = ρ F F q ∇Φ i (q; x) • ∇Φ j (q; x)dx 1≤i,j≤6
, (

where J = J(q) is the inertia matrix of the solid (1.3). It is easy to check that K S (q) is strictly positive and K F (q) is positive so that K(q) is invertible.

The result below shows that equations (1.1e) and (1.1f) can be rewritten as equations giving ḧ(t) and ω(t), in terms of z(t) defined in (2.8). Proposition 4.1. Assume that u, p, q satisfy (1.1). Then

ḧ(t) ω(t) = [K(q(t))] -1 0 (J(t)ω(t)) ∧ ω(t) + ρ F F q(t)
∇µ(z(t); x) • ∇Φ(q(t); x)dx .

(4.2)

In the above formula, the notation ∇µ • ∇Φ stands for the six dimensional vector of components (∇µ

• ∇Φ i ) 1 i 6
, where µ is the solution of (3.17) and (Φ i ) 1 i 6 are defined by (3.15).

Proof. The decomposition of the pressure (3.18), the formulae (1.1e) and (1.1f) imply that, for every t ∈ [0, T 0 ) we have

m s ḧj (t) = ρ F ∂S q(t) µ(z(t); x)K j (q(t); x)dσ x -ρ F 3 i=1 ḧi (t) ∂S q(t) Φ i (q(t); x)K j (q(t); x)dσ x -ρ F 3 i=1 ωi (t) ∂S q(t)
Φ i+3 (q(t); x)K j (q(t); x)dσ x , (4.3)

3 i=1 J i,j (t) ωi (t) = (J(t)ω(t) ∧ ω(t)) j + ρ F ∂S q(t) µ(z(t); x)K j+3 (q(t); x)dσ x -ρ F 3 i=1 ḧi (t) ∂S q(t) Φ i (q(t); x)K j+3 (q(t); x)dσ x -ρ F 3 i=1 ωi (t) ∂S q(t)
Φ i+3 (q(t); x)K j+3 (q(t); x)dσ x , (4.4)

where K j have been defined in (3.16).

On the other hand, using (3.15), (3.17) and Green's formula we get

∂S q(t) Φ i (q(t); x)K j (q(t); x)dσ x = ∂S q(t) Φ i (q(t); x) ∂φ j ∂n (x)dσ x , = F q(t) ∇Φ i (q(t); x) • ∇Φ j (t, x)dx, ∂S q(t) µ(z(t); x)K j (q(t); x)dσ x = ∂S q(t) µ(z(t); x) ∂φ j ∂n (x)dσ x , = F q(t)
∇µ(z(t); x) • ∇Φ j (q(t); x)dx,

Using the last two formulas in (4.3) and (4.4) we obtain the conclusion (4.2).

Recall the definition of E m and F m from (2.7) and (2.12), respectively, and let L F : F m → H m (Ω, R 3 ) be defined by

L F (z)(y) = ρ F ∇Φ(q; η) • L S (z)(y) -ρ F ∇µ(z; η(y)) (y ∈ F 0 ), (4.5)
for every z = (η, q, ξ, k, ω) * ∈ F m , where Φ is the solution of the Neumann problem (3.15), µ is solution of (3.17) and

L S (z) = [K(q)] -1 0 (J(ω) ∧ ω + ρ F Fq ∇µ(z; x) • ∇Φ(q; x)dx . (4.6)
Let L : F m → E m be defined by

L(z) =       ξ k A(ω)R L F (z) L S (z)       . (4.7)
In the last part of this section we show that the system (1.1) is equivalent to the ordinary differential equation

dz dt (t) = L(z(t)), z(0) = (Id F q 0 , h 0 , Id 3 , u 0 , k 0 , ω 0 ) * . (4.8)
In the following Proposition we prove that every solution of (1.1) generates a solution of (4.8).

Proposition 4.2. Let m 3 an integer, assume that (h 0 , Id 3 ) * ∈ P (Ω, S 0 ) and (u 0 , k 0 , ω 0 ) * ∈ T σ 0 Σ where σ 0 = (Id F q 0 , h 0 , Id 3 ). Moreover, assume that

q ∈ C 2 [0, T ), P (Ω, S 0 ) , u ∈ C [0, T ), H m (F q ) ∩ C 1 [0, T ), H m-1 (F q ) , p ∈ C [0, T ), N m+1 (F q )
satisfy the system (1.1). Then z defined by (2.8) satisfies (4.8).

Proof. The equations for η, h, R in (4.8) are nothing else but the definitions of ξ, k and ω from (2.2), (2.5) and (2.6), respectively. The fact that the equations for ξ, k and ω hold follows from (2.4), Proposition 4.1 and from (3.18).

We still have to show that every solution of (4.8) generates a strong solution of (1.1). Proposition 4.3. Let m 3 an integer, assume that (h 0 , Id 3 ) * ∈ P (Ω, S 0 ) and (u 0 , k 0 , ω 0 ) * ∈ T σ 0 Σ where σ 0 = (Id F q 0 , h 0 , Id 3 ) * . Moreover, assume that

z = (η, h, R, ξ, k, ω) * ∈ C([0, T 0 ); F m ) ∩ C 1 ([0, T 0 ); E m ),
is a solution of (4.8). Let q, u, p be defined by q = (h, R) * ,

u(t, x) = ξ(t, η -1 (t, x)), t ∈ [0, T 0 ), x ∈ F q(t) ,
and let the pressure p be defined by (3.18). Then q, u, p satisfy the smoothness conditions (1.10)-(1.12) and the system (1.1).

Proof. First remark that, since z ∈ C([0, T ); F m ) ∩ C 1 ([0, T ); E m ), we have (div u)(t, x) = 0 (t ∈ [0, T 0 ), x ∈ F q(t) ), u(t, x) • n(t, x) = 0 (t ∈ [0, T 0 ), x ∈ ∂Ω), u(t, x) • n(t, x) = v(t, x) • n(t, x) (t ∈ [0, T 0 ), x ∈ ∂S q(t) ),
so that equations (1.1b), (1.1c), (1.1d) are satisfied. From the definition (4.8) of L we obtain that Ṙ = A(ω)R and

ḧ(t) ω(t) = [K(q(t)] -1 0 (J(t)ω(t)) ∧ ω(t) + ρ F F q (t) ∇µ(t, x) • ∇Φ(t, x)dx , ξ(t) = ρ F ∇Φ(t, η(t, y)) • L S (t, z(t)) -ρ F ∇µ(t, η(t, y)),
where K(q) is given by (4.1). The Newton's laws (1.1e) and (1.1f) come from the definition of the pressure (3.18) in the same way as in the proof of Proposition 4.1. Finally, using the relation ξ = u • η and (3.18), we obtain that (1.1a) also holds.

Locally Lipschitz property of L

In this section we tackle a key point of our approach, which consists in proving that the map L is locally Lipschitz. We frequently use below results and methods from [START_REF] Bourguignon | Remarks on the Euler equation[END_REF] combined with techniques specific to our problem, which require to compare functions defined on two different open sets.

Recall that the manifold F m is defined by

F m =    z = σ ν ∈ E m σ = η q ∈ Σ m ,   ξ k ω   ∈ T σ Σ m    ,
where Σ m is defined by (2.10) and T σ Σ m by (2.11). For an element z ∈ F m , the first three components σ = (η, h, R) * ∈ Σ m define the "position" of the system whereas ν = (ξ, k, ω) * ∈ T σ Σ m defines the velocity. The key point in this section is the study of the application µ from (3.17). Recall the notation q = h R .

The main new issue we need to tackle is the study the dependence of the solution µ of (3.17) with respect to the geometric parameter q. The dependence of µ with respect to ξ and η is studied using the ideas in [START_REF] Bourguignon | Remarks on the Euler equation[END_REF].

We first introduce several functions which are useful for the remaining part of this section. Let α, β 0 , β and τ be the mappings on F m defined by

α(z; y) = F q (ξ • η -1 )(η(y)) (y ∈ F 0 ), (5.1)
β 0 (z; y) = G q (ξ • η -1 )(η(y)) (y ∈ ∂Ω), (5.2) β(z; y) = G q (ξ • η -1 -v)(η(y)) (y ∈ ∂S 0 ), (5.3)
τ (z; y) = 2(ξ(y) -v(η(y))) • ω ∧ RN ( y) -[(v • ∇)v(η(y)) -ω ∧ k] • RN ( y) (y ∈ ∂S 0 ), (5.4)
where y = h 0 + R * (η(y) -h), F q and G q have been defined in (3.3) and (3.4), whereas N is a smooth extension of the unit normal vector of ∂F 0 to F 0 .

Remark 5.1. According to a result from Takahashi [START_REF] Takahashi | Analysis of strong solutions for the equations modeling the motion of a rigid-fluid system in a bounded domain[END_REF] and Cumsille and Tucsnak [START_REF] Cumsille | Wellposedness for the Navier-Stokes flow in the exterior of a rotating obstacle[END_REF], for every q = (h, R) * ∈ P (Ω, S 0 ) and ε > 0 small enough there exists a C ∞ diffeomorphism Ψ q : F 0 → F q such that det[DΨ q (y)] = 1 for all y ∈ F 0 and

Ψ q (y) = y if d(y, ∂Ω) ε, Ψ q (y) = h+R(y -h 0 ) if d(y, ∂S q ) ε. (5.5)
Using Ψ q the unit normal vector field on ∂F q can be extended to F q such that Dn(x) = DN (x)

x ∈ ∂Ω, (5.6)

Dn(x) = RDN (h 0 + R * (x -h))R *
x ∈ ∂S q .

(5.7)

Moreover, the construction of Ψ in [START_REF] Takahashi | Analysis of strong solutions for the equations modeling the motion of a rigid-fluid system in a bounded domain[END_REF] shows that Ψ is C ∞ with respect to q.

Proposition 5.2. Let m 3 be an integer. Then the mappings α, β 0 , β and τ are locally Lipschitz (with respect to z) from

F m to H m-1 (F 0 ), H m-1/2 (∂Ω), H m-1/2 (∂S 0 ) and H m-1/2 (∂S 0 ), respectively. Proof. Let z 0 = (σ 0 , ν 0 ) * ∈ F m with σ 0 =   Id F 0 h 0 Id 3   ∈ Σ m , ν 0 =   u 0 k 0 ω 0   ∈ T σ 0 Σ m . For r > 0 we define B m (r) ⊂ F m by B m (r) = {z ∈ F m | z -z 0 E m r} .
We first note that, by the chain rule, we have

α(z; y) = tr ([Dξ(y)][Dη(y)] -1 ) 2 (y ∈ F 0 ).
Since H m-1 (F 0 ) is a Banach algebra, to show that α is Lipschitz on B m (r) it suffices to check that the maps

z → Dξ, z → [Dη] -1 , are Lipschitz from B m (r) to [H m-1 (F 0 )] 9 .
The first map above is obviously Lipschitz whereas for the second one it suffices to use the fact that for every 3 × 3 matrices A, B of determinant equal to 1 we have

A -1 -B -1 = A -1 (B -A)B -1 , A -1 = cof(A) t ,
where cof(A) is the signed cofactors matrix of A.

For β 0 we remark that, using again the chain rule combined with (5.6), we have

β 0 (z; y) = [DN (η(y)))]ξ(y) • ξ(y) (y ∈ ∂Ω).
By applying Lemma A.3 from [START_REF] Bourguignon | Remarks on the Euler equation[END_REF] it follows that the mapping z

→ DN • η is Lipschitz B m (r) to [H m (F 0 )] 9 .
Using again the fact that H m (F 0 ) is a Banach algebra it follows that β 0 is Lipschitz from B m (r) to H m (F 0 ). By the trace theorem it follows that β 0 is Lipschitz from B m (r) to H m-1/2 (∂Ω).

For β we note that

β(z; y) = R[DN (η(y))]R * (ξ(y) -v(η(y))) • (ξ(y) -v(η(y))) (y ∈ F 0 ).
In the the remaining part of this proof, z 1 and z 2 are generic points in B m (r) and we denote by K(r) any Lipschitz constant obtained in Proposition 5.2. With the notation from this section, it is not difficult to check that the Neumann problem (3.17) can be rewritten as:

-∆µ(z; x) = α(z; η -1 (x)) (x ∈ F q ), (5.10a) ∂µ ∂n (z; x) = β 0 (z; η -1 (x)) (x ∈ ∂Ω), (5.10b) ∂µ ∂n (z; x) = β(z; η -1 (x)) + τ (z; η -1 (x)) (x ∈ ∂S q ), (5.10c)
where α, β 0 , β and τ have been defined in (5.1), (5.2), ( 5.3) and (5.4). The main difficulty consists in the fact that the functions µ 1 = µ(z 1 ; •) and µ 2 = µ(z 2 ; •) are not defined on the same domain. Using Lemma A.4 from [START_REF] Bourguignon | Remarks on the Euler equation[END_REF] we have

χ(z 1 ) -χ(z 2 ) H m (F 0 ,R 3 ) = ∇µ 1 • η 1 -∇µ 2 • η 2 H m (F 0 ,R 3 ) K(r) ∇µ 1 • η -∇µ 2 H m (Fq 2 ,R 3 ) , where η = η 1 • η -1 2 .
By applying Proposition 5.3 we obtain

∇µ 1 • η -∇µ 2 H m (F q 2 ,R 3 ) ≤ K(r)(I 1 + I 2 + I 3 + I 4 ),
where I i , with i ∈ {1, 2, 3, 4}, are given by

I 1 = div (∇µ 1 • η -∇µ 2 ) H m-1 (F q 2 ) , I 2 = curl (∇µ 1 • η -∇µ 2 ) H m-1 (Fq 2 ,M 3 (R)) , I 3 = (∇µ 1 • η -∇µ 2 ) • n H m-1/2 (∂Fq 2 ) , I 4 = |||∇µ 1 • η -∇µ 2 |||,
where

||| • ||| is the norm on H m-1 (F q 2 , R 3 ) defined by |||u||| = sup F q 2 u(x) • γ(x)dx | γ ∈ C m (F q 2 , R 3 ), γ(x) = 0 for all x ∈ ∂F q 2 , γ C m (F q 2 ,R 3 ) 1 , (5.11)
which clearly satisfies (5.8). Using Lemma A.4 from [START_REF] Bourguignon | Remarks on the Euler equation[END_REF] and (5.10) we have

I 1 ≤ K(r) div (∇µ 1 • η) -∆µ 1 • η H m-1 (Fq 2 ,R 3 ) + α(z 1 ) -α(z 2 ) H m-1 (Fq 0 ,R 3 ) .
Using Lemmas A.4 and 4 from [START_REF] Bourguignon | Remarks on the Euler equation[END_REF] we obtain

div (∇µ 1 • η) -∆µ 1 • η H m-1 (F q 2 ) K(r) η -Id Fq 2 H m (F q 2 ) ∇µ 1 H m (F q 1 ) , K(r) η 1 -η 2 H m (F 0 ) ∇µ 1 H m (F q 1 ) .
On the other hand, using again Lemma A.4 and Theorem 3.1 we have

∇µ 1 H m (Fq 1 ) ≤ K(r)( α(z 1 ) H m-1 (F 0 ) + β 0 (z 1 ) H m-1/2 (∂Ω) + β(z 1 ) + τ (z 1 ) H m-1/2 (∂S q 0 ) K(r) (z 1 ∈ B m (r)).
(5.12)

The last two estimates and the locally Lipschitz property of α proved in Proposition 5.2 imply that

I 1 ≤ K(r) z 1 -z 2 E m (z 1 , z 2 ∈ B m (r)). (5.13)
Using the fact that curl(∇f ) = 0 together with arguments completely similar to those used for I 1 , we obtain a constant K(r) such that

I 2 ≤ K(r) z 1 -z 2 E m (z 1 , z 2 ∈ B m (r)). (5.14)
To tackle I 3 , let n i the unit normal vector to ∂F q i , with i ∈ {1, 2}. We have

I 3 ≤ (∇µ 1 • η) • n 2 -n 1 • η) H m-1/2 (∂F q 2 ) + ∂µ 1 ∂n 1 • η - ∂µ 2 ∂n 2 H m-1/2 (∂F q 2 )
.

Using trace inequalities, estimate (5.12) and Lemma A.4 from [START_REF] Bourguignon | Remarks on the Euler equation[END_REF] we obtain that

I 3 ≤ K(r) n 2 (η 2 ) -n 1 (η 1 ) H m-1/2 (∂F 0 ) + ∂µ 1 ∂n 1 • η 1 - ∂µ 2 ∂n 2 • η 2 H m-1/2 (∂F 0 ) = K(r) n 2 (η 2 ) -n 1 (η 1 ) H m-1/2 (∂F 0 ) + β 0 (z 1 ) -β 0 (z 2 ) H m-1/2 (∂Ω) + β(z 1 ) + τ (z 1 ) -β(z 2 ) -τ (z 2 ) H m-1/2 (∂S 0 ) .
Applying Lemma A.3 from [START_REF] Bourguignon | Remarks on the Euler equation[END_REF] to the extensions of n i to F q i (these extensions have been defined in Remark 5.1, it follows that

n 2 (η 2 ) -n 1 (η 1 ) H m-1/2 (∂F 0 ,R 3 ) ≤ K(r) z 1 -z 2 E m (z 1 , z 2 ∈ B m (r)).
The last two estimates and the Lipschitz properties of α, β 0 , β and τ imply that

I 3 ≤ K(r) z 1 -z 2 E m (z 1 , z 2 ∈ B m (r)). (5.15)
To study I 4 we first note that, for every γ ∈ C m (F q 2 , R 3 ) with γ = 0 of ∂F q 2 we have

F q 2 [∇µ 1 • η(x) -∇µ 2 (x)]•γ(x)dx = F q 1 ∇µ 1 •γ•η -1 (x)dx- F q 2 ∇µ 2 (x)•γ(x)dx.
(5.16)

Consider the functions ψ k : F q k → R, defined as the solutions of the Neumann problems:

-∆ψ 1 = -div (γ • η -1 ), in F q 1 , (5.17a) ∂ψ 1 ∂n = 0 on ∂F q 1 .
(5.17b)

-∆ψ 2 = -div γ, in F q 2 , (5.18a) ∂ψ 2 ∂n = 0 on ∂F q 2 .
(5.18b)

Taking the inner product in L 2 (F q 1 ) (respectively in L 2 (F q 2 )) of the first equation in (5.17) (respectively in (5.18)) by µ 1 (respectively by µ 2 ) and the subtracting side by side, we obtain that

F q 1 ∇ψ 1 •∇µ 1 dx- F q 2 ∇ψ 2 •∇µ 2 dx = F q 1 ∇µ 1 •γ•η -1 (x)dx- F q 2 ∇µ 2 (x)•γ(x)dx.
The above formula and (5.16) yield that

Fq 2 [∇µ 1 • η(x) -∇µ 2 (x)] • γ(x)dx = Fq 1 ∇ψ 1 • ∇µ 1 dx - Fq 2 ∇ψ 2 • ∇µ 2 dx.
Using the variational formulation of the Neumann problem (5.10) we obtain that, for i ∈ {1, 2}, we have

Fq i ∇µ i • ∇ψ i dx = Fq i α(z i ; η -1 i (x))ψ i (x)dx + ∂Ω β 0 (z i ; η -1 i (x))ψ i (x) dσ x + ∂Sq i β(z i ; η -1 i (x)) + τ (z i ; η -1 i (x)) ψ i (x)dσ x .
The last two formulas imply that

Fq 2 [∇µ 1 • η(x) -∇µ 2 (x)] • γ(x)dx = Fq 1 α(z 1 ; η -1 1 (x))ψ 1 (x)dx - Fq 2 α(z 2 ; η -1 2 (x))ψ 2 (x)dx + ∂Ω β 0 (z 1 ; η -1 1 (x))ψ 1 (x) dσ x - ∂Ω β 0 (z 2 ; η -1 2 (x))ψ 2 (x) dσ x + ∂S q 1 β(z 1 ; η -1 1 (x)) dσ x - ∂S q 2 β(z 2 ; η -1 2 (x)) dσ x + ∂S q 1 τ (z 1 ; η -1 1 (x))ψ 1 (x)dσ x - ∂S q 2 τ (z 2 ; η -1 2 (x))ψ 2 (x)dσ x . (5.19)
To estimate the difference of the first two terms in the right-hand side of the above formula, we note that

F q 1 α(z 1 ; η -1 1 (x))ψ 1 (x)dx - F q 2 α(z 2 ; η -1 2 (x))ψ 2 (x)dx = F q 0 [α(z 1 ; y)(ψ 1 • η 1 )(y)dy -α(z 2 ; y)(ψ 2 • η 2 )(y)] dy α(z 1 ; •) -α(z 2 ; •) ψ 1 • η 1 + α(z 2 ; •) ψ 1 • η 1 -ψ 2 • η 2 , (5.20)
where all the norms above are in L 2 (F q 0 ). The first term in the right-hand side of the above relation is readily estimated by using Proposition 5.2 to get

α(z 1 ; •) -α(z 2 ; •) ψ 1 • η 1 K(r) z 1 -z 2 E m γ C m (F q 2 ;R 3 ) , (5.21)
for every z 1 , z 2 ∈ B m (r). To estimate the second term in the right-hand side of (5.20) we remark that, using the variational formulations of (5.17) and (5.18) and a simple change of variables we have, for k ∈ {1, 2},

Fq 0 (∇ψ k • η k ) • (∇ϕ k • η k ) dy = Fq 0 (γ • η 2 ) • (∇ϕ k • η k ) dy (ϕ k ∈ H m (F q k )).
Denoting

ψ k = ψ k • η k the last formula becomes Fq 0 (Dη -1 k )(Dη -1 k ) * ∇ ψ k •∇ϕ dy = Fq 0 (Dη -1 k )(γ •η 2 )•∇ϕ dy, (ϕ ∈ H m (F q 0 )).
Subtracting side by side the formulas corresponding to k = 1 and k = 2 it is not difficult to see that, for every

γ ∈ C m (F q 2 , R 3 ), γ = 0 on ∂F q 2 ,
we have

Fq 0 |∇ ψ 1 -∇ ψ 2 | 2 dy K(r) z 1 -z 2 2 E m (z 1 , z 2 ∈ B m (r)).
The above estimate, combined to (5.20) and (5.21) imply that

F q 1 α(z 1 ; η -1 1 (x))ψ 1 (x)dx - F q 2 α(z 2 ; η -1 2 (x))ψ 2 (x)dx K(r) z 1 -z 2 2 E m (z 1 , z 2 ∈ B m (r)).
The other terms in the right-hand side of (5.19) can be estimated in a similar way. In order to keep this paper of reasonable length we skip the proof of the corresponding estimates.

We are now in position to prove that L is locally Lipschitz. Proof. We begin by showing that L S is locally Lipschitz in B m (r) for a given r. From Proposition 3.4 the mapping q → K(q) is C 2 from P (Ω, S) to M 6 (R) (recall that K(q) is the virtual mass matrix defined in (4.1)). Using Proposition 3.4 together with Lemmas A.2 and A.3 from [START_REF] Bourguignon | Remarks on the Euler equation[END_REF], it follows that the mapping

z → ∇Φ • η is Lipschitz from B m (r) to H m (F 0 , M 3×6 (R))
, where (Φ k ) k∈{1,...,6} satisfy (3.15). Moreover, the mapping z → (0 3 , (Jω)∧ω) * is Lipschitz from B m (r) to R 6 . Using the notation in Proposition 4.1, the last term in the right-hand side of (4.2) writes

F q ∇µ(x) • ∇Φ(x) dx = F 0 ∇µ(η(y)) • ∇Φ(η(y)) dy,
so that, using Proposition 5.4 and Proposition 3.4, we obtain that that this term defines a Lipschitz function from B m (r) to R 6 . Using next the smoothness of the map q → K(q), it follows that L S is Lipschitz from B m (r) to R 6 . Finally, the fact that L F is locally Lipschitz readily follows from (4.5) and the corresponding properties of Φ, L S and µ.

L is tangent to F m

In this section we show that the vector field defined by the operator L from (4.7) is tangent to the closed set F m which has been defined in (2.12). More precisely, the main result of this section is Proposition 6.1. Let m 3 an integer and let z 0 ∈ F m . Then

lim r→0 1 r dist(z 0 + rL(z 0 ); F m ) = 0.
In order to prove the above proposition we need some notation and several auxiliary results. Throughout this section e 0 denotes the identity map on F 0 and q 0 = h 0 Id 3 . Moreover z 0 = σ 0 ν 0 denotes a generic element of F m , where

σ 0 =   η 0 h 0 R 0   ∈ Σ m and ν 0 =   ξ 0 k 0 ω 0   ∈ T σ 0 Σ m . Let σ =   η h R   ∈ H m (F 0 , R 3 ) × P (Ω, S 0 ),
and recall the properties (1.6) and (1.7) of ∂Ω. We define the map

ϑ(σ) = ϑ 1 ϑ 2 (σ) (σ ∈ H m (F 0 , R 3 ) × P (Ω, S 0 ) , ( 6.1)
where

ϑ 1 (σ)(y) = det(Dη)(y) - 1 |F 0 | F 0 det(Dη(x)) dx - 1 |F 0 | ∂Ω δ 0 (η(x)) dσ x - 1 |F 0 | ∂S 0 δ(h 0 +R * (η(y)-h)) dσ x σ ∈ H m (F 0 , R 3 ) × P (Ω, S 0 ), y ∈ F 0 , (6.2) ϑ 2 (σ)(y) = δ 0 (η(y)) (σ ∈ H m (F 0 , R 3 ) × P (Ω, S 0 ), y ∈ ∂Ω) , δ(h 0 + R * (η(x) -h)) (σ ∈ H m (F 0 , R 3 ) × P (Ω, S 0 ), y ∈ ∂S 0 ) . (6.
3) Since we obviously have

F 0 ϑ 1 dy + ∂F 0 ϑ 2 dσ y = 0,
it follows that ϑ maps H m (F 0 , R 3 ) × P (Ω, S 0 ) into V m-1 (F 0 ) (see (3.1) for the definition of this space). Lemma 6.2. The function ϑ defined above is of class C k for all integer k ≥ 1 and we have

∂ η ϑ e 0 q 0 (V ) = div V -V • n | ∂F 0 V ∈ H m (F 0 , R 3 ) . (6.4) Moreover, ∂ η ϑ e 0 q 0 maps H m (F 0 , R 3 ) onto V m-1 (F 0 ). Proof. Since m ≥ 3, it follows that H m-1 is an algebra so that the map η → det(Dη) is of class C k from H m (F 0 , R 3 ) to H m-1 (F 0 , R) for every k ≥ 1.
It is easy to check that the other terms in the definition of ϑ 2 and are smooth functions so that ϑ is of class C k for every k ≥ 1.

Using (1.7) it follows that, for every σ ∈ H m (F 0 , R 3 )×P (Ω, S 0 ), V ∈ H m (F 0 , R 3 ), we have

(∂ η ϑ 2 )(σ)(V )(y) = -V (y) • n(η(y)) (y ∈ ∂Ω) , -R * V (y) • n(h 0 + R * (η(y) -h)) (y ∈ ∂S 0 ) . (6.5)
On the other hand, using the fact that the differential of A → det(A) is the linear map H → tr(cof(A)H), where cof(A) is the signed cofactors matrix of A, we obtain that, for every σ ∈ H m (F 0 , R 3 ) × P (Ω, S 0 ), V ∈ H m (F 0 , R 3 ), we have

∂ η ϑ 1 (σ)(V )(y) = tr(cof(Dη)DV )(y) - 1 |F 0 | F 0 tr(cof(Dη)DV ))(y) dy + 1 |F 0 | ∂Ω V (y) • n(η(y))dσ y + 1 |F 0 | ∂S 0 R * V (y) • n(h 0 + R * (η(y) -h)) dσ y .
Taking h = h 0 , R = Id 3 and η = e 0 in the above formula and by using (6.5) we obtain (6.4). Finally, the fact that the right-hand side of (6.4) defines a map from The above lemma can be used, in particular, to show that Σ m is an infinitedimensional manifold over H m (F 0 , R 3 ) × R 3 × M 3 (R) and to compute its tangent space at σ 0 . Proposition 6.3. We have

H m (F 0 , R 3 ) onto V m-1 (F 0 ) is classical,
Σ m = η q ∈ H m (F 0 , R 3 ) × P (Ω, S 0 ) ϑ q η = 0 . (6.6)
Moreover, the tangent space to Σ m at every σ

= η q ∈ O ∩ Σ m (Ω, S 0 ) is the space T σ Σ m defined in (2.11).
Proof. The fact that the set in the left-hand side of (6.6) is a subset of the set int the right-hand side is obvious. To prove the converse inclusion, we first note from ϑ 1 (q, η) = 0 it follows that det (Dη) is constant in F 0 . On the other hand, from ϑ 2 (q, η) = 0 it follows that η(q)(∂F 0 ) ⊂ ∂F q . These assertions, combined to the fact that F 0 and F q have the same volume, imply that det (Dη) = 1 in F 0 . Moreover, the above properties enable us to apply the global inverse mapping theorem of Caccioppoli (see, for instance, Zeidler [18, Theorem 4.G. page 174]) to obtain that η ∈ Diff m (F 0 , F q ). This concludes the proof of (6.6).

In order to prove the second assertion in the proposition, we first note that for every σ

= (η, h, R) * ∈ Σ m and every (ξ, k ω) * ∈ H m (F 0 , R 3 ) × R 3 × R 3 , we have Dϑ 1 (σ)   ξ k ω   (y) = tr Dη -1 Dξ (y ∈ F 0 ), Dϑ 2 (σ)   ξ k ω   (y) = -ξ • n(η(y)) (y ∈ ∂Ω) , -R * [ξ(y) -k -A(ω)(η(y) -h)] • n(h 0 + R * (η(y) -h)) (y ∈ ∂S 0 ) ,
where A(ω) has been defined in (1.2). Since η is a diffeomorphism from F 0 to F q , denoting u = ξ • η -1 and making the change of variable x = η(y), we obtain

Dϑ 1 (σ)   ξ k ω   (η -1 (x)) = (div u)(x) (x ∈ F q ), Dϑ 2 (σ)   ξ k ω   (η -1 (x)) = -u • n(x) (x ∈ ∂Ω) , -[u(x) -k -ω ∧ (x -h)] • n(x) (x ∈ ∂S q ) .
From the above formulas it follows that the kernel of Dφ(σ) is T σ Σ m so that we obtain the second assertion in the proposition.

Proposition 6.4. Let σ 0 =   e 0 h 0 Id 3   ∈ Σ m , ν 0 =   u 0 k 0 ω 0   ∈ T σ 0 Σ m and let γ 0 =   Γ L M   ∈ H m (F (0), R 3 ) × R 3 × R 3 such that div(Γ)(y) = F(u 0 (y)) (y ∈ F 0 ), (6.7a) Γ • n(y) = -G(u 0 (y)) (y ∈ ∂Ω), (6.7b) Γ • n(y) = -G(u 0 (y) -v 0 (y)) -2(u 0 -v 0 ) • (ω 0 ∧ n) + [L + M ∧ (y -h 0 ) + ω 0 ∧ (w 0 ∧ (y -h 0 ))] • n (y ∈ ∂S 0 ), (6.7c)
where v 0 (y) = k 0 + ω 0 ∧ (y -h 0 ). Then there exists ε > 0 and a curve

σ =   η h R   ∈ C 2 ([0, ε]; Σ m (Ω, S 0 )) satisfying σ(0) = σ 0 , σ(0) =   u 0 k 0 A(ω 0 )   σ(0) =   Γ L A(M ) + [A(ω 0 )] 2   . (6.8) Proof. According to Proposition 6.3 the curve σ is contained in Σ m iff ϑ(σ(t)) = 0 for every t ∈ [0, ε].
We begin by constructing the "rigid displacement part" h R of the curve σ.

To do that, we define h, ω : R → R 3

h(t) = h 0 + tk 0 + t 2 2 L, ω(t) = ω 0 + tM t ∈ R, (6.9)
and R : R → SO 3 (R) is defined as the solution of the initial value problem

R(0) = Id 3 , Ṙ(t) = A(ω(t))R(t), ( 6.10)
where A(ω) is the skew-adjoint matrix defined in (1.2). Note that

Ṙ(0) = A(ω 0 ), R(0) = A(M ) + [A(ω 0 )] 2 . (6.11)
The above functions being continuous, it follows that there exists ε > 0 such that

q(t) = h(t) R(t) ∈ P (Ω, S 0 ) t ∈ [0, ε ].
In order to construct the "fluid part" η of the curve σ we first note that, since e 0 q(0) ∈ Σ m , we have ϑ e 0 q(0) = 0. Therefore, by combining Lemma 6.2 with a version of the implicit function theorem (see, for instance, Zeidler [18, Theorem 4.H. page 171]) it follows that there exists ε ∈ (0, ε ] and a function η : [0, ε] → H m (F 0 , R 3 ) such that η(0) = e 0 and ϑ η(t) q(t) = 0, P η(t) -tu 0 -t 2 2 Γ = 0 (t ∈ [0, ε]), (6.12) where P is the orthogonal projector from H m (F 0 , R 3 ) onto Ker ∂ n ϑ e 0 q(0)

. Note that, according to Lemma 6.2, we have Ker ∂ n ϑ e 0 q(0)

= u ∈ H m (F 0 , R 3 ) | div u = 0 u • n = 0 on ∂F 0 (6.13)

In the remaining part of the proof we show that, with the above choice of η, h and R, the curve

σ(t) =   η(t) h(t) R(t)   (t ∈ [0, ε])
satisfies (6.8). We first note that from (6.9) it follows that h(0) = h 0 , ḣ(0) = k 0 , ḧ(0) = L.

From the above formula combined with (6.10) and (6.11) we see that, in order to prove (6.8), we have only to check that η(0) = u 0 , η(0) = Γ. (6.14)

Taking the derivative with respect to t of the formula det (Dη(t)) = 1 we obtain that tr((Dη(t)) -1 D η(t)) = 0 (t ∈ [0, ε]). (6.15) Using next the fact that δ 0 (η(t)) = 0 on ∂Ω it follows that η(t) • n(η(t)) = 0 (on ∂Ω). (6.16)

Moreover, since δ(h 0 + R * (t)(η(t) -h(t))) = 0 on ∂S 0 , we have ( η(t)-ḣ(t)-ω(t)∧(η(t)-h(t)))•R(t)n(h 0 +R * (t)(η(t)-h(t))) = 0 (on ∂S 0 ). (6.17) On the other hand, taking the derivative of the second formula in (6.12) with respect to t we obtain P( η(t) -u 0 -tΓ) = 0 (t ∈ [0, ε]). (6.18)

Taking t = 0 in (6.15)-(6.18) we obtain div ( η(0)) = 0, (in F 0 ), η(0) • n = 0 (on ∂Ω), η(0) • n = (k 0 + ω 0 ∧ (y -h 0 )) (y ∈ ∂S 0 ), P( η(0)) = P(u 0 ).

The above relations clearly imply that the first equality in (6.14) holds.

In order to prove the second equality in (6.14) we take the derivative of (6.15)-(6.18) and then we make t = 0. In this way we obtain div η(0) = F(u 0 ) (in F 0 ), η(0) • n = -G(u 0 ) (on ∂Ω), η(0) • n = -G(u 0 -v 0 ) -2(u 0 -v 0 ) • (ω 0 ∧ n) + [L + M ∧ (y -h 0 ) + ω 0 ∧ (ω 0 ∧ (y -h 0 ))] • n (y ∈ ∂S 0 ), P(η(0)) = P(Γ), where v 0 has been defined in the statement of this proposition. Using (6.7) it follows that the second equality in (6.14) also holds.

We are now in a position to prove the main result of this section.

Proof of Proposition 6.1. Recall the notation for z 0 from the beginning of this section. We first note (by using an appropriate change of variables) that it suffices to prove the result for η 0 = e 0 and R 0 = Id 3 . This will be done by constructing a curve Z(•) in F m such that The main tool of the proof is Proposition 6.4, with an appropriate choice of Γ, L and M . More precisely, u 0 , k 0 and ω 0 are chosen to be those in (1.1) and we take

Γ = L F (z 0 ), L M = L S (z 0 ), (6.20)
where L F and L S have been defined in (4.5) and (4.6), respectively. The fact that Γ, L and M chosen above satisfy the assumptions in Proposition 6.4 follows from (3.15) and (3.17). Define

Z(t) =     σ(t) η(t) ḣ(t) ω 0 + tM     ,
where σ(t) = (η(t), h(t), R(t)) * is the curve constructed in Proposition 6.4. By combining (6.8) and (6.20) it follows that Z(0) = z 0 , Ż(0) = L(z 0 ), which imply (6.19).

The proof of our main result in Theorem 1.3 can be now written as follows.

Proof of Theorem 1.3. The assumptions in Theorem 1.3 imply that z 0 = (e 0 , h 0 , Id 3 , u 0 , k 0 , ω 0 ) * ∈ F m .

Therefore, we can combine Propositions 5.5, 6.1 and 2.1 to obtain that the initial value problem ż = L(z), z(0) = z 0 , (6.21) admits an unique solution z = (η, h, R, ξ, k, ω) * ∈ C 0 [0, T 0 ); F m ) ∩ C 1 [0, T 0 ); E m ) .

According to Proposition 4.3, q, u defined by q = (h, R) * , u(t, x) = ξ(t, η -1 (t, x)) (t ∈ [0, T 0 ), x ∈ F q(t) ), and the pressure p defined by (3.18) define a strong solution of (1.1). We have thus shown the announced existence result.

To prove the uniqueness, it suffices to note that, according to Proposition 4.2, any strong solution of (1.1) defines a solution of (6.21) and to apply Proposition 2.1.

. 14 )

 14 By combining (3.13) and (3.14) and (3.10) we obtain (3.8).

Proposition 5 . 5 .

 55 The mappings L S , L F and L are locally Lipschitz on F m .

  lim

r→0 1 r

 1 dist(z 0 + rL(z 0 ); Z(r)) = 0. (6.19)

  see for instance, Lemma 2.4.1 in Sohr [16, page 79].

The fact that β is Lipschitz from B m (r) to H m-1/2 (∂Ω) can now be proved in the same way as for β 0 .

Finally for τ we notice that for every x ∈ F q we have

Inserting the above formula in (5.4) and applying again Lemma A.3 from [START_REF] Bourguignon | Remarks on the Euler equation[END_REF], the claimed Lipschitz property of τ easily follows.

We also need the following classical result (see, for instance, [1, Lemma 5]):

Proposition 5.3. Let m 3 be an integer and let Ω be bounded domain of R 3 with smooth boundary. Then for every u ∈ H m (Ω, R 3 ) there exists a constant K, which depends on m and on Ω, such that

where

for some C > 0 then there exists a constant K > 0 such that

We are now in a position to give the main ingredient needed to prove that L is locally Lipschitz. This result concerns the potential µ introduced in (3.17).

Proposition 5.4. For every integer m 3, the function χ defined on F m by

)

We use again the notation from the proof of Proposition 5.2, i.e., for r > 0 we set


Nicolas Quaegebeur
email: nicolas.quaegebeur@gmail.com

Antoine Chaigne
email: antoine.chaigne@ensta.fr

Guy Lemarquand
email: guy.lemarquand@univ-lemans.fr

Transient modal radiation of axisymmetric sources: application to loudspeakers

Keywords: Nicolas Quaegebeur, Antoine Chaigne, Guy loudspeaker, nonlinear vibrations, transient, radiation, Rayleigh PACS: 43.20.Px, 43.20.Tb, 43.20.Rz, 43.40.-r, 43.40.Rj

Transient modal radiation of axisymmetric sources:

application to loudspeakers.

The aim of an electrodynamic loudspeaker is to transform an electrical signal into sound. Such a transduction is expected to be linear. However, for high levels of vibrations, nonlinear phenomena appear and are responsible for audible distortions. The various sources of nonlinearities can be separated into two parts [START_REF] Klippel | Nonlinear large signal behavior of electrodynamic loudspeakers at low frequencies[END_REF][START_REF] Quaegebeur | Mechanical resonances and geometrical nonlinearities in electrodynamic loudspeakers[END_REF]: "electrical" nonlinearities due to the large displacement of the coil in the permanent magnet, and "mechanical" nonlinearities due to large displacements of the moving parts of the system (geometrical nonlinearities). In order to describe the complex dynamics of the system subjected to large amplitude motion, the displacement field must be expressed in the timedomain. In that case, a common way to compute the transient displacement field is to expand the transverse displacement of the structure onto the linear modes. This method has received much attention in the literature in the case of circular sources [START_REF] Sridhar | Nonlinear resonances in the forced responses of plates. Part 1 : symetric responses of circular plates[END_REF], axisymmetric sources [START_REF] Touzé | Effect of Imperfections and Damping on the Type of Nonlinearity of Circular Plates and Shallow Spherical Shells[END_REF] and recently in the case of loudspeaker like-structure [START_REF] Quaegebeur | Nonlinear vibrations of loudspeaker-like structures[END_REF] but none of the cited studies include the acoustic radiation of the structure. In order to compute the acoustic field radiated by thin structures subjected to large amplitude motion, the present article proposes a transient model of acoustic radiation that takes advantage of the modal expansion.

Transient acoustic radiation of planar sources vibrating in an infinite planar baffle has received much attention in the literature. The Rayleigh surface integral is in that case an exact expression for predicting the acoustic pressure produced by such sources at any observation point. An extensive review of the various approaches which have been used to evaluate this radiation integral has been given by Harris [START_REF] Harris | Review of transient field theory for a baffled planar piston[END_REF]. The most common way to compute the transient radiation from such sources is to evaluate the so-called Spatial Impulse Response (S.I.R.) as described by Stepanishen [7]. This formulation has been developed initially in the case of an uniform velocity for a circular source, and has then been extended to handle the case of nonuniform velocity patterns [START_REF] Jensen | A new calculation procedure for Spatial Impulse Responses in ultrasound[END_REF][START_REF] Harris | Transient field of a baffled planar piston having an arbitrary vibration amplitude[END_REF] using a modal expansion of the displacement field.

For nonplanar sources, the Rayleigh surface integral is no longer valid, since diffraction effects appear. An analytical formulation of the transient radiation of a body of arbitrary shape has been derived by Hu and Wu [START_REF] Hu | An explicit integral formulation for transient acoustic radiation[END_REF][START_REF] Wu | An alternative formulation for predicting sound radiation from a vibrating object[END_REF], who expressed the acoustic pressure as a sum of integrations of simple and doublet source distributions and their couplings induced by the normal and tangential components of the particle velocity. Under that consideration, the Rayleigh integral represents the first term of this decomposition. This approximation is commonly used to compute transient radiation of spherical radiators [START_REF] O'neil | Theory of focusing radiators[END_REF][START_REF] Arditi | Transient fields of concave annular arrays[END_REF][START_REF] Kates | Radiation from a Dome[END_REF][START_REF] Verhoef | The impulse response of a focused source with an arbitrary axisymmetric surface velocity distribution[END_REF][START_REF] Ketterling | Acoustic field of a wedge-shaped section of a spherical cap transducer[END_REF], but the range of validity (in terms of geometrical configuration and wavenumber) is rarely specified.

Only the case of spherical sources has been treated and only few studies concern the acoustic radiation of axisymmetric sources. Farn [START_REF] Farn | Transient acoustic fields generated by a body of arbitrary shape[END_REF] proposed a numerical method based on source-density method (approximation of the radiating surface by triangular surface elements) and Guyomar [START_REF] Guyomar | Transient radiation from axially symmetric sources[END_REF] used series expansion of the source velocity spatial distribution. More recently, Stepanishen [START_REF] Stepanishen | Acoustic axisymmetric radiation and scattering from bodies of revolution using the internal source density and Fourier methods[END_REF] adapted a mean-square error method to predict harmonic radiation of bodies of arbitrary shape but the radiation of transients remains a problem. Other numerical approaches using retarded potentials combined to boundary elements in the time domain have can also be found in the literature [START_REF] Chappell | A stable boundary element method for modeling transient acoustic radiation[END_REF][START_REF] Bluck | Analysis of three-dimensional transient acoustic wave propagation unsing the boundary integral equation method[END_REF] and give appropriate results but calculation costs remain a problem.

In the present study, a time-domain formulation is adopted to predict the modal sound radiation from axisymmetric sources with a nonuniform vibratory distribution. The pressure field is approximated by the Rayleigh integral corresponding to a monopole source distribution over the nonplanar vibrating surface. The diffraction effects are neglected in the present approach and let for future work. The displacement field is expanded onto the linear modes of the structure and a change of variables in the Rayleigh integral is proposed in the case of a monotonic profile function to compute the Spatial Impulse Response associated to each mode of vibration efficiently. The results are compared to the exact formulation obtained in the case of planar sources [START_REF] Harris | Review of transient field theory for a baffled planar piston[END_REF] and to numerical results obtained by Suzuki and Tichy [START_REF] Suzuki | Sound radiation from convex and concave domes in an infinite baffle[END_REF][START_REF] Suzuki | Radiation and Diffraction Effects by Convex and Concave Dome[END_REF] with the spherical harmonics method in the case of spherical radiators. The method of calculation is then derived in the case of a typical loudspeaker profile (association of a truncated cone with a spherical cap). Finally, the present approach is used to estimate the nonlinear radiation pattern of a prototype loudspeaker and predictions are compared to measurements in anechoic room. In the present study, attention is paid to the acoustic radiation of a baffled axisymmetric body S as represented in Fig. 1. Let a denote the external radius, h 0 (r) the axisymmetric profile function of the considered structure and h max its depth (defined as the maximum of the profile function). The profile function only depends on the radial coordinate r since S is assumed axisymmetric. For the same reason, the observation point x depends on two coordinates (y, z) in the Cartesian space only. The profile function h 0 (r) is called monotonic if its derivative ḣ0 (r) has a constant sign and called non-monotonic otherwise. A concave structure refers to a negative profile function h 0 (r) < 0 while a convex structure relates to a positive profile function. In the following equations, x S denotes a point of the structure with cylindrical coordinates (r, θ, h 0 (r)).

Nonlinear vibration of loudspeaker-like structures

Previous studies [START_REF] Quaegebeur | Mechanical resonances and geometrical nonlinearities in electrodynamic loudspeakers[END_REF][START_REF] Quaegebeur | Nonlinear vibrations of loudspeaker-like structures[END_REF] mention that the nonlinear vibrations of loudspeaker-like structure can be modeled using a nonlinear lumped parameters model coupled with a nonlinear model of the diaphragm vibrations. This formulation allows to include typical electrical nonlinearities that appear in the low frequencies domain but also geometrical nonlinearities that appear around mechanical resonances of the structure. The approach is based on an axisymmetric modal decomposition of the displacement field of the structure:

w(r, t) = ∞ p=0 Φ p (r)q p (t). ( 1
)
where the functions Φ p (r) in Eq. ( 1) represent the modal shapes (that can be either measured or calculated) and the functions q p (t) represent the time function associated to mode p. Using a state-space formulation of the nonlinear electromechanical problem [START_REF] Quaegebeur | Mechanical resonances and geometrical nonlinearities in electrodynamic loudspeakers[END_REF], one obtains:

Ẏ = AY + B(Y )u(t) + N L e (Y ) + N L g (Y ) ( 2
)
where Y is the state vector (time functions q p (t) and their first derivative qp (t)), A denotes the linear dynamics of the problem, B the excitation vector, u(t) the input signal, N L e (Y ) and N L g (Y ) contain respectively all the electrical, mechanical and geometrical nonlinear terms of the electromechanical transduction. Eq. ( 2) is solved numerically using the fixed point method and Runge-Kutta algorithms, so that the global displacement field is obtained from the calculation of the state vector Y using modal reconstruction Eq. ( 1). This global nonlinear formulation describing the complex dynamics of a typical loudspeakers needs to be directly solved in the time-domain. In order to compute the associated radiated sound pressure field for high amplitudes of vibration, a time-domain formulation of sound radiation from axisymmetric sources with a nonuniform vibratory distribution has to be adopted, as presented below.

Integral formulation

The sound pressure field P (x, t) is solution of the well known Kirchhoff-Helmholtz integral [START_REF] Junger | Sound, structures and their interaction[END_REF], that can be reduced to the Rayleigh's integral when considering weakly curved structures, i.e. for h 0 (r) ≪ a:

P (x, t) = ρ t 0 S G(x, x S , t, t ′ ) ∂V n (x S , t ′ ) ∂t ′ dx S dt ′ (3)
Using the axisymmetric modal decomposition Eq. ( 1), the Rayleigh's integral is directly computed in the time-domain by separating the space and the time integration using the Spatial Impulse Response principle [START_REF] Stepanishen | Acoustic axisymmetric radiation and scattering from bodies of revolution using the internal source density and Fourier methods[END_REF]:

P (x, t) = ρ t 0 ∞ p=0 H p (x, t, t ′ ) qp (t ′ ) dt ′ , (4)
where the Spatial Impulse Response (S.I.R.) H p associated to mode p is introduced as follow:

H p (x, t, t ′ ) = S G(x, x S , t, t ′ ) Φ p (x S ) -→ u z . --→ dx S . (5)
where -→ u z .

--→ dx S denotes the dot product between these 2 vectors. The analytical formulation of those functions can be obtained through an appropriate change of coordinates in the case of planar [START_REF] Stepanishen | Acoustic axisymmetric radiation and scattering from bodies of revolution using the internal source density and Fourier methods[END_REF] or spherical structures [START_REF] Arditi | Transient fields of concave annular arrays[END_REF]. However, in the general case, the direct calculation in the time-domain is not obvious due to the presence of the Dirac delta function. The following section proposes an appropriate change of variable in order to compute directly the S.I.R. functions for an axisymmetric source.

CALCULATION OF THE RAYLEIGH INTEGRAL IN THE TIME-DOMAIN

Geometrical considerations

The aim of the following section is to simplify the calculation of the impulse response defined in Eq. ( 5). Geometrical considerations are firstly made in order to separate 3 areas of calculations where the properties of the sound pressure field differ. Then, a new local basis is proposed using those geometrical considerations and the S.I.R. are simply expressed using the new set of coordinates. This method, developped in the case of planar and spherical sources is extended in the present study to the case of axisymmetric sources and offers the advantage of suppressing the Dirac delta function in Eq. (5) (using the convolution theorem), so that the final formulation is reduced to a simple integration.

The Spatial Impulse Response H p in Eq. ( 5) has a compact support [t i ; t o ] where t i denotes the propagation time between the receiver x and the nearest point of the structure x i and t o represents the propagation time between the receiver and the farthest point of the structure x o . The points x i = (0, r i ) and

x o = (0, r o ) are solutions of: |x -x i | = min (r, θ) (|x -x S |) , (6a) |x -x o | = max (r, θ) (|x -x S |) . (6b)
Let x h and x b denote the points located respectively at the top and at the bottom of the structure, as shown in Fig. 2. Depending on the structure and on the location of the receiver x, one can define 3 different zones defined in Tab. 1, where the radiation properties differ.

zone Conditions on x i and x o zone I

x i = {x h , x b } and x o = {x h , x b } zone II x i = {x h , x b } and x o = {x h , x b } zone III x i = {x h , x b } and x o = {x h , x b } Table 1
Definition of the 3 calculation zones.

Change of coordinates

If the present section, the aim is to transform the integration variables of Eq. ( 5) by 2 variables |xx S | and an angle φ that correspond to the geometry of the problem.

Let θ i denote the angle between axis u z and the line (xx i ). A new local orthonormal basis B i = ( t ix , t iy , n i ) centered on point x i is introduced, as presented in Fig. 3. It is defined in the Cartesian basis B c = ( u x , u y , u z ) as follows:

             -→ n i = -sin θ i u y -cos θ i u z , -→ t iy = cos θ i u y -sin θ i u z , -→ t ix = u x . (7)
x ( 0 , y , z ) x ( 0 , y , z ) . . In this basis, a set of cylindrical coordinates is introduced for any point located on the surface x S = (r 2 , φ) and the coordinates in the new local basis are defined in appendix A. Introducing the tangent vectors to the surface ----→ T |x-x S | and -→ T φ with respect to the variables |xx S | and φ (detailed in Appendix A), the integration vector of Eq. ( 5) becomes in the new local basis [START_REF] Kaplan | Advanced Calculus[END_REF]:

--→ dx S = -→ T φ × ----→ T |x-x S | d|x -x S | dφ. (8)
where (.) × (.) denotes the cross product between 2 vectors. The expression of tangent vectors and integration vector can be found in Appendix A. The principle is to replace the integration vector in Eq.( 5) by the new one defined in Eq. ( 8) and to integrate over |xx S | in order to eliminate the Dirac delta function contained in the Green's function (using the convolution theorem [START_REF] Bracewell | Convolution Theorem : The Fourier Transform and Its Applications[END_REF]). Introducing the time delay τ = t -t ′ , one obtains the final formulation of the Rayleigh integral:

P (x, t) = - ρc 2π t 0 p H p (x, τ )q p (t -τ )dτ (9)
where:

H p (x, τ ) = φ Φ p (r) F A + B cos φ -C sin 2 φ |x-x S |=cτ dφ (10)
and:

                             F = 1 + dh 0 dr h 0 (r) -z r , A = 1 -(sin θ i sin φ) 2 cos θ i , B = F r i -y r 2 , C = dh 0 dr y r sin θ i . (11)
In this expression, the Spatial Impulse Response explicitly depends on the profile function h 0 (r), the modal shape Φ p (r) (also dependent on the shape of the source), the position of the receiver x and the integration time τ . The notation |xx S | = cτ means that for each time τ , a distance |xx S |/c is associated for computing the integrand of Eq. [START_REF] Hu | An explicit integral formulation for transient acoustic radiation[END_REF]. In that case, for each time step τ , the intersection between the structure and the sphere centered on observation point x and of radius cτ is sought. This path (determined by the boundary values of integration path over φ) is unique in the case of a monotonic profile function (when the sign of ḣ0 (r) is constant) and can be closed

if φ = [0 : 2π] or open if φ = [φ min : φ max ]
depending on the location of the receiver, the shape of the structure and the integration time τ . The calculation of those boundaries and the application in the case of a spherical cap are expressed in Appendix B for each zone of calculation.

The formulation of Eq. ( 10) allows to reduce the order of integration compared to Eq. ( 5): only one spatial coordinate on variable φ is required and since spatial and time-domain integration are separated, the spatial discretization is independent on the time step (and thus the maximal simulated frequency). In Fig. 4, the Spatial Impulse Response is computed for a hemispherical shell and a receiver located at (y, z) = (a, a) i.e. in the near-field region and for different discretization steps on variable φ in Eq. [START_REF] Hu | An explicit integral formulation for transient acoustic radiation[END_REF]. The solutions in each case is compared to the solution of the Rayleigh integral [START_REF] Arditi | Transient fields of concave annular arrays[END_REF] represented by points. It appears that less than 20 points for φ are required to attain a convergence error below 1% of the integral defined by Eq. ( 10). Fig. 4. Spatial Impulse Responses in the case of a hemispherical transducer and for a receiver located at (y, z) = (a, a). The exact solution of the Rayleigh's integral in the time-domaine derived from [START_REF] Arditi | Transient fields of concave annular arrays[END_REF] is represented by points (• • •) and the results obtained using Eq. ( 10) in the case of different discretizations of the integral over φ are represented by lines : 3 points (dashed light grey line), 6 points (dashed dark grey line), 9 points (dashed black line) and 12 points (solid black line). The convergence is ensured in any case with at less 20 points.

VALIDATION

Planar sources

In the case of a planar structure, the profile function is equal to zero: h 0 (r) = 0.

In that case, Eq. ( 10) becomes [START_REF] Stepanishen | Transient radiation from pistons in an infinite planar baffle[END_REF]:

H p (x, τ ) = φ Φ p (r)| |x-x S |=cτ dφ (12)
with the boundaries derived from Eq. (C.1) in Appendix B in accordance with the results of Harris [START_REF] Harris | Review of transient field theory for a baffled planar piston[END_REF] and Jensen [START_REF] Jensen | A new calculation procedure for Spatial Impulse Responses in ultrasound[END_REF]. The classical example [START_REF] Harris | Review of transient field theory for a baffled planar piston[END_REF] of Spatial Impulse Responses obtained for a planar piston with a uniform velocity distribution is shown in Fig. 5. In this figure, S.I.R. functions are calculated using an integral discretization of 20 points. The in-axis impulse response is equal to 1 in a compact support which means that the effect of radiation filter is equivalent to a low pass filter whose cut-off frequency is determined by the size of the support of the Spatial Impulse Response. The off-axis response decreases with respect to the off-axis position y while its compact support increases. In that case, the radiation filter is also equivalent to a low pass filter whose cut-off frequency is below the the in-axis cut-off frequency. This result is general for axisymmetric structures : the radiation filter is equivalent to a low-pass filters whose cut-off frequency decreases with respect to the off-axis position y.

Spherical sources

In the present section, the approximated formulation of Eq. ( 9) is compared to the numerical results obtained by Suzuki [START_REF] Suzuki | Sound radiation from convex and concave domes in an infinite baffle[END_REF] using 40 spherical harmonics in the wavenumbers domain, taking into account the diffraction effects. In the case of spherical sources, the profile function h 0 (r) equals:

     h 0 (r) = √ R 2 -r 2 - √ R 2 -a 2 for r < a, h 0 (r) = 0 for r > a, (13)
where R denotes the radius of curvature of the spherical shell and a the external radius of the structure. In order to compare equivalent results in the frequency domain, one needs to take the Fourier transform of the Spatial Impulse Response computed using Eq. [START_REF] Hu | An explicit integral formulation for transient acoustic radiation[END_REF] in the case of a rigid body motion (one mode of vibration: p = 1 and uniform displacement of the structure: Φ p (r) = 1 for r < 1). In the work of Suzuki, the calculations are limited by the order of the spherical harmonics (40 for his study), so that his model is restricted to wavenumbers up to ka = 10. In the present formulation, only the time discretization of Spatial Impulse Responses influences the upper ka . On-axis acoustic pressure (up: magnitude in dB, down: phase in degrees) in far-field versus nondimensionalized wavenumber ka computed by Suzuki [START_REF] Suzuki | Sound radiation from convex and concave domes in an infinite baffle[END_REF] (points) and using Eq. ( 9) (lines). The results are presented for 3 different concave structures: planar piston (black), medium spherical cap (dark grey) and hemispherical cap (light gray). The wavenumbers domain results are obtained by taking the Fourier transform of the Spatial Impulse Response computed in the time domain.

limit, so that no limitation of frequency range is observed. In practical cases, the Spatial Impulse Responses are computed using 20 nondimensionnalized time steps τ = tc/a and then completing with zeros, the precision in the wavenumbers domain is adaptable.

The results are presented for 3 different structures: a planar piston, a spherical cap defined by a radius of curvature R = 1.5a and an hemispherical cap . On-axis acoustic pressure (up: magnitude in dB, down: phase in degrees) in far-field versus wavenumber ka computed by Suzuki [START_REF] Suzuki | Sound radiation from convex and concave domes in an infinite baffle[END_REF] (points) and using Eq. ( 9) (lines). The results are presented for 3 different convex structures: planar piston (black), medium spherical cap (dark grey) and hemispherical cap (light gray).

defined by a radius of curvature R = a. The Fourier Transform of the Spatial Impulse Responses computed for an on-axis receiver located in the far-field (y, z) = (0, 100a) are presented in Figs. 6 and7 for concave and convex caps and a good agreement between the present approach and the results of Suzuki is obtained for low wavenumbers (ka < 1) and high wavenumbers (ka > 3).

For low wavenumbers (ka < 1), the monopole approximation is valid and the equivalent radiation filter has a flat response. For high wavenumbers (ka > 3),

δ Fig. 8. Explanation of the interference phenomenon that appears when the wavelength k is above the path difference between the nearest and the farther point seen from the observer point. In the frequency domain, this effect is responsible for a slope of the on-axis sound pressure level of -20 dB/decade and an extension of the directivity compared to planar radiators.

the local curvature induces a decrease of the on-axis sound pressure level and the radiation filter is thus equivalent to a low-pass filter with a slope of -20 dB/decade above a cut-off wavelength depending on the curvature R of the profile and its radius a. In Figs. 6 and 7 the cut-off frequency of the equivalent low-pass filter decreases when the curvature radius R increases. This is explained by interferences phenomena that appears when the wavelength λ = 1/k is below the path difference between the nearest and the farther point seen from the observer point as presented in Fig. 8. Thus the more the structure radius R is important, the more the difference path increases and the more the cut-off frequency of the equivalent low-pass filter decreases.

Around ka = 1, the diffraction taken into account in the calculations of Suzuki [START_REF] Suzuki | Sound radiation from convex and concave domes in an infinite baffle[END_REF] induces an increase of on-axis sound pressure level up to +6 dB in the case of concave caps (a decrease in the case of convex structures) that is not taken into account in the present model. The effect of diffraction appears in the non-dimensionnalized wavenumbers domain 1 < ka < 3 i.e. when the wavelength approaches the maximal depth of the structure. The influence of the diffraction is hardly dependent on the shape and the validity of the present approach is thus no longer viable in the wavelength range 1 < ka < 3.

Figs. 9 and 10 represent the directivity patterns in linear scale computed for the considered spherical caps in the far-field (z = 100a) for wavenumbers such as ka = 1, ka = 3 and ka = 10. Again, the results obtained by taking the Fourier Transform of Eq. ( 10) are compared to results obtained by Suzuki [START_REF] Suzuki | Sound radiation from convex and concave domes in an infinite baffle[END_REF]. For low wavenumbers (ka = 1) the radiation patterns are almost omnidirectional and become more directive for high wavenumbers (ka > 3). The curvature extends the directivity in both convex and concave cases and in the concave case, the maximum of sound pressure is obtained for a precise obser-vation angle depending on the curvature of the considered shell.

In a general manner, the curvature (not necessarily in the case of a spherical shell) induces a decrease of on-axis sound pressure field due to interference for high wavenumbers. The effects of diffraction, not taken into account in the present model, are responsible for an increase of on-axis sound pressure level around 1 < ka < 3 in the case of concave sources (a decrease in the case of convex structures). The on-axis high wavenumbers losses due to interferences are compensated by an extent of directivity related to the local curvature of the source. This result is general for monotonic profile functions and is extended in the following section in the case of non-monotonic profile function and in the case of a typical loudspeaker shape.

Complex sources

The previous formulation is only valid for a monotonic profile function h 0 (r). In the case of a nonmonotonic profile, i.e. when the sign of ḣ0 (r) is not constant, the integration path (intersection between the sphere centered on observation point x and of radius cτ used in the calculation of Eq. ( 10)) is made of multiple paths corresponding to each monotonic part of the global profile.

Neglecting the multiple reflections over the emissive surface, the global radiation can be expanded as the sum of contributions of each monotonic portion of the profile. For example, for a global structure composed by 2 monotonic profiles S = S 1 S 2 , we have:

P (x, t) = - ρc 2π p t 0 H 1 p (x, τ 1 )q p (t -τ 1 )dτ 1 + t 0 H 2 p (x, τ 2 )q p (t -τ 2 )dτ 2 (14)
where H i p (x, τ 1 ) corresponds to the Spatial Impulse Response of mode p associated to the i th monotone part of the profile. This procedure is applied in the case of a typical loudspeaker shape composed as the junction of a truncated cone with a spherical cap. The profile function is then described by:

             h 0 (r) = α(b -a) + √ R 2 -r 2 - √ R 2 -b 2 for r < b, h 0 (r) = α(r -a) for r > b, h 0 (r) = 0 for r > a, (15)
where R denotes the curvature of the inner spherical shell, b the radius of the ka=1 : junction between both sub-structures and α the depth of the truncated conical shell. The following results are obtained for a typical loudspeaker α = 1, b = 0.5a and R = 1.5b. The profile function and the principle of calculation of the on-axis Spatial Impulse Response are presented in Fig. 11 In Fig. 11, the respective contribution of each sub-structure in the time-domain for an on-axis receiver is represented. In the frequency domain, the global ra- diation filter is as previously equivalent to a low-pass filter whose cut-off frequency is not only related to the depth of the truncated conical shell. Indeed, between ka = 1 and ka = 10 in Fig. 12, the comparison between the response of the truncated conical shell only and the global structure shows that the cutoff frequency is increased when the inner spherical cap is present. The effect of the convex spherical cap compensate the interference effect around ka = 1 and an extension of the flat acoustic response is achieved and a global increase of +3dB with respect to the conical source is observed in the interference regime.

The contribution of the inner cap appears in the tail of the Impulse Response and is also responsible for high wavenumber comb-filtering (for ka > 10) that is not present in absence of the inner spherical cap (see Fig. 12). Around ka = 25, a decrease of 10dB with respect to the conical source is observed and can be interpreted as interferences between the conical and the spherical cap because this phenomenum is not present in the case of a spherical shell only neither a conical shell only as presented on Fig. 12. This effect is then amplified above ka = 25 and is responsible for comb filtering at multiples of that nondimensionalized wavenumber ka = 25, ka = 50 and ka = 75. Fig. 12. On-axis acoustic pressure level (magnitude in dB versus wavenumber ka) obtained in the case of a loudspeaker (junction of a truncated cone and an hemispherical cap). The solution (solid black line) is compared to the truncated cone radiation only (dashed grey line). The use of a typical loudspeaker profile increases the low-pass filter cut-off frequency and induces oscillations in the on-axis acoustic pressure level, due to complex interference phenomena.

APPLICATION TO A PROTOTYPE OF LOUDSPEAKER

Presentation

In order to illustrate the influence of the diaphragm shape of a loudspeaker on its acoustical response for high amplitudes of vibration, a prototype of loudspeaker has been developed by G. Lemarquand in the LAUM [START_REF] Lemarquand | New Structure of Loudspeaker. 120th AES Convention[END_REF] as presented on Fig. 13. Its moving part is an aluminium convex spherical shell of 50 mm of diameter. Experiments have been done with a curvature equals to a/R = 0.3 which is defined as the ratio between the external radius and the radius of curvature of the shell. For the prototype, the moving part is not made in one piece, since the emissive surface is glued on the supporting cylinder.

In the studied prototype, attention is paid to the motor and suspensions, in order to minimize the "electrical" nonlinear phenomena. The sources of nonlinearities in the motor are the variations of the fore factor (Bl), and the Eddy currents. As a remedy to these defects, an ironless motor was used, which is made only out of neodymium iron boron permanent magnets. The classical suspensions of loudspeakers are mostly made of rubber, impregnated fabric or molded plastic. They act as a spring, but have a nonlinear behavior. This means that their compliance depends on the movements amplitude and, above all, that the induced damping depends greatly on both amplitude and frequency. These phenomena are the sources of the classical nonlinear behavior of the loudspeakers suspensions at low frequencies. In the studied prototype, the almost perfect compressibility properties of the air were used to create a pneumatic stiffness. Indeed, the air is compressed in the closed box, whose volume tunes the stiffness value. The cabinet is thus a cylindrical pipe which is closed and filled up with an absorbing material. Classical electrical nonlinearities are thus not observed in the present prototype and only geometrical nonlinearities remain and are due to large amplitudes of vibration of the moving part. This kind of nonlinearities has been pointed out mechanically but the influence on the radiated sound pressure field had not been cleared explained.

Electromechanical characterization

The first step of the experimental analysis consist in understanding and modeling the electromechanical behavior of the prototype. Fig. 14 represents the impedance curves in both magnitude. Below 1 kHz, the behavior is mostly resistive and become inductive above 2 kHz. As mentioned above, the inductive part is not modified by eddy currents, due to the ironless conception of the motor, and then a linear dependency of impedance magnitude with respect to frequency is observed (not observable on Fig. 14 because of the logarithmic scale for the frequencies). the prototype. Two tendencies are easily identified on Fig. 14 : • for frequencies below 5 kHz, the resonances are due to stationary waves in the rear cavity : the first resonance is around 70 Hz. Under this assumption, the vibration pattern is then assimilated to a plane piston (rigid body mode).

• for frequencies above 9 kHz, the resonances correspond to mechanical reso-nances of the spherical cap, as described in the following part. The properties of those resonances (eigenfrequencies, modal damping, excitation and shapes) have been analyzed using Laser Doppler Velocity and a linear modal analysis software (IDEAS c ) and are presented on Tab. 2. Axisymmetric modal parameters (eigenfrequency f p , modal damping µ p and excitation T p associated to mode p) measured for the spherical dome. The first mode corresponds to a piston mode (rear cavity resonance) and is very damped compared to the structural modes (above 9 kHz).

Mode

Nonlinear radiation : on-axis response

For high amplitudes of vibration, nonlinearities appear and give rise to harmonic distortion in the velocity and sound pressure signals. Due to the ironless conception of the motor, the electrical nonlinearities classically observed in electrodynamic loudspeakers vanish and only geometrical nonlinearities (due to large deformations of the moving part) remain.

This type of nonlinearities can be modeled using the modal formulation proposed in [START_REF] Quaegebeur | Nonlinear vibrations of loudspeaker-like structures[END_REF]. This approach is absolutely suitable in the present case, since few modes radiate in the audible range. The principle is to model the prototype as a 4 degrees of freedom system (since 4 modes radiate in the audible range, as presented in Tab. 2) and to include geometrical nonlinearities in the modal equations through cubic and quadratic nonlinearities [START_REF] Quaegebeur | Mechanical resonances and geometrical nonlinearities in electrodynamic loudspeakers[END_REF]. In order to illustrate the typical effects of geometrical nonlinearities on radiation of the prototype, only quadratic and cubic nonlinearity terms have been added for the second mode of vibration (at 9650 Hz). The displacement field is calculated for high amplitudes of vibrations using the State-Space formulation described by Eq. ( 2), and the radiation problem is solved using the Spatial Impulse Response (S.I.R.) approach developed in the present study. In this practical case, the measured modal shapes Φ p (r ′ ) and profile function h 0 (r ′ ) are introduced in Eq. ( 10) and convolutions are performed at the sampling frequency of 48 kHz. It is important to notice that the calculations of the displacement field and radiated sound pressure are computed directly in the time-domain and the results are presented in the frequency domain for clarity (taking the Fourier transform of the transient responses). Fig. 15 represent the on-axis sound pressure measurement and prediction at 1 m of the fundamental and of the harmonics 2 and 3 in the case of the convex spherical dome. Measurements have been performed in anechoic room (approximate size 1000m 3 ) using sinusoidal input signals. A good agreement is observed in the audible bandwidth and the geometrical nonlinearities are correctly predicted around discrete frequencies such as :

• f = f a : around a mechanical resonance frequency, the amplitude of vibration increases and harmonics are generated. For example, on figure 15, harmonics 2 and 3 are generated for an excitation frequency around 9.6kHz (first mechanical resonance).

• f = f a /p with p integer : for the submultiple p of a resonance frequency, the frequency of harmonic p corresponds to an eigenfrequency and this harmonic is then amplified. For example, on Fig. 15, an increase of harmonic 3 at 3.2 kHz which corresponds to a third of the first mechanical resonance is observed (the same effect appears with harmonic 2 at 4.8 kHz).

The Spatial Impulse Response appear to model efficiently the complex radiation pattern for high amplitudes of vibration. The influence of the geometrical nonlinearities on the acoustic response of the prototype has been observed experimentally and modeled properly using the present approach. Other sources of nonlinearities could also be taken into account in order to predict the nonlinear acoustic response of classical loudspeakers.

CONCLUSIONS

In order to compute acoustic radiation from axisymmetric structures subjected to large amplitude motion, an explicit integral formulation is derived for predicting the modal time-domain radiation. The radiated acoustic pressure is assumed to be expressed as an integration of a simple source distribution (Rayleigh integral), which is valid except in the nondimensionalized wavenumber range 1 < ka < 3 because the diffraction effects are neglected in the present study and will be detailed in a future work. The displacement is expanded onto the linear modes of the structure and the acoustic pressure is expressed as a sum of modal contributions, by introducing the Spatial Impulse Responses depending on the shape of the source, the position of the receiver and the modal shapes (also dependent of the shape of the source).

The solution is compared to the literature in the cases of a piston and of spherical caps and gives quite good results and low calculation costs compared to frequency-domain formulations. In a general way, the radiation induces a lowpass filtering whose cut-off frequency is related to the radius of the source and the local curvature. The method is extended for predicting the radiation of an axisymmetric source whose profile function is complex and corresponds to a typical loudspeaker shape. Finally, predictions and experimental measurements of on-axis radiation of a prototype loudspeaker subjected to large amplitudes motion are compared and a good agreement is achieved on the audible bandwidth.

A Definition of the new local basis

In the new local basis B i = ( t ix , t iy , n i ), the coordinates of point x S are defined by:

--→ x i x S = h 0 (r i ) -h 0 (r) cos θ i -r 2 cos φ tan θ i -→ n i + r 2 cos φ -→ t iy + r 2 sin φ -→ t ix .

(A.1)

where r 2 is defined by the Cartesian coordinates of integration point x S = (x, y):

B Change of variable

One point located in the integration surface S is parameterized in the Cartesian basis by its 3 coordinates x S = (x, y, h 0 ( √ x 2 + y 2 )). In order to compute the tangent vectors expressed in Eq. (A.3), one needs to express the dependency of |xx S | and φ with respect to those 3 coordinates: where F , A, B and C are defined in Eq. [START_REF] Wu | An alternative formulation for predicting sound radiation from a vibrating object[END_REF].

        

C Boundary values of integration

Depending on the zone of calculation, the boundaries of integration in Eq. ( 10) differ and the integration has to be computed on a closed or an open path (an example of the two cases is represented on

Fig. 1 .

 1 Fig.1. Geometry of the considered structure. The points x and x s represent the receiver and an integration point on the structure respectively. The profile function is denoted by h 0 (r) and corresponds to the projection in the plane z = 0 of the u z coordinate of the point x s .

Fig. 2 .

 2 Fig. 2. Definition of the 3 different areas of calculation. For each position of the receiver, the associated points x i (dotted lines) and x o (dashed lines) are represented . The zones I and III only are defined in front and in back of a convex structure respectively. The intersection point F of the 3 zones is called the sweet spot (or geometrical focal point).

Fig. 3 .

 3 Fig.3. Definition of the new local basis B i = ( t ix , t iy , n i ) and the new coordinates (r 2 , φ). The points x, x i and x S denote the observation point (located in the zone I for the present example), the projection point and an integration point located on the surface S, respectively.

Fig. 5 .

 5 Fig.5. Spatial Impulse Responses obtained in the case of a plane piston with an uniform velocity distribution. The solutions are computed for a receiver located in z = a and an off-axis position y from 0 to 5a and for 30 points of discretization in Eq. (12).

Fig. 6

 6 Fig.6. On-axis acoustic pressure (up: magnitude in dB, down: phase in degrees) in far-field versus nondimensionalized wavenumber ka computed by Suzuki[START_REF] Suzuki | Sound radiation from convex and concave domes in an infinite baffle[END_REF] (points) and using Eq. (9) (lines). The results are presented for 3 different concave structures: planar piston (black), medium spherical cap (dark grey) and hemispherical cap (light gray). The wavenumbers domain results are obtained by taking the Fourier transform of the Spatial Impulse Response computed in the time domain.

Fig. 7

 7 Fig.7. On-axis acoustic pressure (up: magnitude in dB, down: phase in degrees) in far-field versus wavenumber ka computed by Suzuki[START_REF] Suzuki | Sound radiation from convex and concave domes in an infinite baffle[END_REF] (points) and using Eq. (9) (lines). The results are presented for 3 different convex structures: planar piston (black), medium spherical cap (dark grey) and hemispherical cap (light gray).

Fig. 9 .

 9 Fig.9. Directivity functions (using a linear scale) for the 3 different concave caps: plane piston (black), medium spherical cap (dark grey) and hemispherical cap (light gray). The results obtained by the present model (right) are compared to the results obtained by Suzuki[START_REF] Suzuki | Sound radiation from convex and concave domes in an infinite baffle[END_REF] (left) for 3 wavenumbers: ka = 1 (up), ka = 3 (middle) and ka = 10 (bottom).

Fig. 10 .

 10 Fig.10. Directivity functions (using a linear scale) for the 3 different convex caps: plane piston (black), medium spherical cap (dark grey) and hemispherical cap (light gray). The results obtained by the present model (right) are compared to the results obtained by Suzuki[START_REF] Suzuki | Sound radiation from convex and concave domes in an infinite baffle[END_REF] (left) for 3 wavenumbers: ka = 1 (up), ka = 3 (middle) and ka = 10 (bottom).

Fig. 12

 12 Fig. 12 taking the Fourier Transform of computed time function.

Fig. 11 .

 11 Fig. 11. Calculation procedure of the Spatial Impulse Response of a typical loudspeaker (junction of a cone and a spherical cap) for an uniform velocity distribution and for a on-axis receiver located in z = a. The S.I.R. corresponding to the spherical cap H 1 (x, τ ) is represented by a dashed black line, the S.I.R. corresponding to the truncated cone H 2 (x, τ ) is represented by a dashed grey line and the global S.I.R. corresponds to the solid black line.

Fig. 13 .

 13 Fig. 13. Photography of the studied prototype in anechoic room.

However, the impedanceFig. 14 .

 14 Fig. 14. Electrical impedance of the loudspeaker in magnitude (up) and phase (down).

Fig. 15 .

 15 Fig.[START_REF] Verhoef | The impulse response of a focused source with an arbitrary axisymmetric surface velocity distribution[END_REF]. Measurement (points) and prediction (solid lines) of on-axis sound pressure level for a 10 W input. For each frequency, the energy of fundamental (black) and harmonics 2 (dark grey) and 3 (light grey) are measured and predicted.

r 2 2 =

 2 x 2 + cos 2 θ i [y -r i -(h 0 (r) -h 0 (r i )) tan θ i ] 2 . (A.2)Under those considerations, the tangent vectors to the surface ----→ T |x-x S | and -→ T φ with respect to the variables |xx S | and φ are described by:T |x-x S | = ∂ -→ x S . -→ u x ∂|xx S | -→ u x + ∂ -→ x S . -→ u y ∂|xx S | -→ u y + ∂ -→ x S . -→ u z ∂|xx S | -→ u z , -→ T φ = ∂ -→ x S . -→ u x ∂φ -→ u x + ∂ -→ x S . -→ u y ∂φ -→ u y + ∂ -→ x S . -→ u z ∂φ -→ u z , (A.3)

  ∂x ∂|xx S | = |xx S | F x ∂y ∂|xx S | = |xx S | F y -r, change of variables in Eq.(8) leads to:-→ u z . --→ dx S = 2|xx S | F A + B cos φ -C sin 2 φ d|xx S |dφ, (B.4)

Fig. C. 1 . 2 - 1 )Fig. C. 1 .

 1211 Fig. C.1. Example of a closed contour (left) and an open contour (right) of integration in the case of a spherical cap and a receiver x located in the zone I. The integration path corresponds to the intersection of the structure with a sphere of radius cτ centered on x.

  = ( u x , u y , u z ) Cartesian basis B i = ( t ix , t iy , n i ) local basis

	Nomenclature
	a	external radius	x i	impact point
	c	celerity of sound in air	θ i	angle of the impact point
	h 0 (r)	profile function	Φ p	p th natural mode for the
	P (x, t) Acoustic pressure		flexural displacement
	q p	time-function associated	ρ	density of air
		to the p th natural mode	σ	ratio of densities
	r R	radial coordinate curvature radius for the spherical cap	B c (x, y, z)	coordinates in the basis B c
	S t	axisymmetric structure time variable	(r 2 , φ) -→ T α	coordinates in the basis B i tangent vector with respect
	t ′	integration time		to the variable α
	x	observation point	H p	Spatial Impulse Response
	x S	integration point		associated to mode p
	1 INTRODUCTION