Alexandre Munnier
email: alexandre.munnier@iecn.u-nancy.fr

Locomotion of deformable bodies in an ideal fluid: Newtonian versus Lagrangian formalisms

Keywords: Biodynamics, ideal fluid, Lagrangian and Newtonian Mechanics, PDE-ODE coupled system, shape sensitivity analysis AMS Subject Classification: 74F10, 70S05, 76B03, 49Q10

This paper is concerned with comparing Newtonian and Lagrangian methods in Mechanics for determining the governing equations of motion (usually called Euler-Lagrange equations) for a collection of deformable bodies immersed in an incompressible, inviscid fluid whose flow is irrotational. The bodies can modify their shapes under the action of inner forces and torques and are endowed with thrusters, what means that they can generate fluid jets by sucking and blowing out fluid through some localized parts of their boundaries. These capabilities may allow them to propel and steer themselves.

Our first contribution is to prove that under smoothness assumptions on the fluid-bodies interface, Newtonian and Lagrangian formalisms yield the same equations of motion. However and rather surprisingly this is no longer true for nonsmooth shaped bodies.

The second novelty brought in in this paper is to treat for the first time a broad spectrum of problems in which several bodies undergoing any kind of deformation can be involved and to display the Euler-Lagrange equations under a form convenient to study locomotion and to perform numerical simulations.

These equations have been used to develop a Matlab toolbox (Biohydrodynamics Toolbox) that allows to study animal locomotion in a fluid or merely the motion of submerged rigid solids. Examples of such simulations are given in this paper.

1 Introduction

History

In the last decade, significant efforts have been done by mathematicians in studying problems relating to the general theme of locomotion in a fluid by shape-changing bodies [START_REF] Lighthill | Mathematical biofluiddynamics[END_REF][START_REF] Childress | Mechanics of swimming and flying[END_REF][START_REF] Galper | Dynamic equations of motion for a rigid or deformable body in an arbitrary non-uniform potential flow field[END_REF][START_REF] Leonard | Stability and drift of underwater vehicle dynamics: mechanical systems with rigid motion symmetry[END_REF][START_REF] Holmes | Dynamics of the Kirchhoff equations. I. Coincident centers of gravity and buoyancy[END_REF][START_REF] Kelly | Modelling efficient pisciform swimming for control[END_REF][START_REF] Wu | Mathematical biofluiddynamics and mechanophysiology of fish locomotion[END_REF][START_REF] Kozlov | Motion of a body with undeformable shell and variable mass geometry in an unbounded perfect fluid[END_REF][START_REF] Patrick | Stability by KAM confinement of certain wild, nongeneric relative equilibria of underwater vehicles with coincident centers of mass and buoyancy[END_REF][START_REF] Kanso | Locomotion of articulated bodies in a perfect fluid[END_REF][START_REF] Melli | Motion planning for an articulated body in a perfect fluid[END_REF][START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF]; see also [START_REF] Sparenberg | Survey of the mathematical theory of fish locomotion[END_REF] for an interesting survey on this topic and further references. The underlying motivation of these works is to better understand the impressive swimming capabilities evolved by fishes and aquatic mammals.

Experiments have shown that the vortices generated by the tail of fishes play a crucial role to understand their locomotion and some models incorporate artificially vortices [START_REF] Mason | Propulsion and control of deformable bodies in an ideal fluid[END_REF][START_REF] Zhu | Threedimensional flow structures and vorticity control in fish-like swimming[END_REF][START_REF] Triantafyllou | Hydrodynamics of fishlike swimming[END_REF]]. If we do not neglect the viscosity effects, the relevant model consists in the non stationary Navier-Stokes equations for the fluid coupled with Newton's laws for the fish-like swimming object. This perspective is the one adopted in [START_REF] Carling | Self-propelled anguilliform swimming: simultaneous solution of the two-dimentional navier-stokes equations and newton's laws of motion[END_REF][START_REF] Liu | A numerical study of undulatory swimming[END_REF][START_REF] Galdi | On the steady self-propelled motion of a body in a viscous incompressible fluid[END_REF][START_REF] San Martin | An initial and boundary problem modeling fish-like swimming[END_REF].

Nevertheless, among numerous mathematical articles studying fish locomotion, most of them address the case of a potential flow which is by definition vortex-free. It is the point of view we have chosen in what follows and on which we focus from now on.

Newtonian versus Lagrangian formalims

Among all of the papers referenced above, one can essentially find two methods to derive the governing equations of motion for the fish-like swimming bodies: either authors invoke Newton's laws (and use Newtonian formalism) or they apply Hamilton's principle (and use Lagrangian formalism). Let us describe succinctly these two approaches.

Classical Mechanics

With Newtonian formalism, the dynamics of the fluid and the dynamics of the bodies are in a first time treated separately and next coupled.

Each body is handled as an independent mechanical system subject to exterior hydrodynamic forces and torques expressed by means of the fluid stress tensor. Newton's laws for linear and angular momenta apply and yield a system of ordinary differential equations (ODE's) whose unknowns are the degrees of freedom of the bodies. The dynamics of the perfect fluid, for its part, is driven by Euler's equations.

The coupling between bodies and fluid is realized not only through the hydrodynamic forces but also through kinematic constraints. Indeed, the region occupied by the fluid as well as its velocity along the bodies' boundaries depend on the locations and velocities of the bodies.

With this approach and in term of complexity, it does not matter whether the fluid is confined or not or how many the submerged bodies are.

Analytic Mechanics

Governing equations of motion can also be obtained by applying Hamilton's principle to the system fluid-bodies in its entirety. This idea goes back to the works of Thomson, Tait and Kirchhoff. In the book [18, chap. VI, page 160], Sir H. Lamb explains that: The cardinal feature of the methods followed by these writers consists in this, that the solids and the fluid are treated as forming together one dynamical system, and thus the troublesome calculation of the effect of the fluid pressures on the surfaces of the solids is avoided.

With this approach, on the contrary to what happens with Newtonian formalism, the more bodies are involved, the more involved is the problem.

Indeed, once the Lagrangian function of the system has been set, the derivation of the equations of motion requires it to be differentiated with respect to the socalled generalized coordinates corresponding to the degrees of freedom of the system. Obviously, more bodies means also more degrees of freedom. Within this process, one usually needs to differentiate the solution of an initial/boundary value problem (namely, equations of Fluid Dynamics) with respect to parameters relating to the geometry of the fluid's domain (i.e. to compute shape derivatives). Although quite involved, this task can be carried out using classical tools of shape optimization [START_REF] Moubachir | Moving shape analysis and control[END_REF][START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF].

Notice however that this difficulty can be overcame when the following restrictive hypotheses hold, what is assumed in most of the papers: (H1) There is only one body in the fluid.

(H2) The overall system body-fluid fills the whole space.

(H3) The body is neutrally buoyant.

(H4) The system body-fluid is at rest at time t = 0.

The method to obtain the Euler-Lagrange equations in this case is sketched in Subsection 1.6. An other way to avoid the shape-derivative difficulty consists in replacing (H1) by: (H1 ) The dynamics of the immersed bodies are hydrodynamically decoupled, which means that each body is handled as if being alone in the fluid. It can be relevant when the immersed bodies are very far one from the others but is not any longer when they are close as shown in Section 7.

Equivalence between both formalisms?

It is classical for monogenic mechanical systems (i.e. systems having forces which derive from a potential) involving both holonomic and nonholonomic constraints, that either Newton's laws or D'Alembert's (or Hamilton's) principle suffice in providing Euler-Lagrange governing equations of motion [3, §13 pages 64-67] (or [11, §2.3 pages 44-45]). More recently, the equivalence has been also shown for multibody systems in [START_REF] Valasek | The complete equivalence of Newton-Euler equations of motion and Lagrange's equations of mixed type[END_REF]. However, to our knowledge, this problem has not be investigated so far for more complex systems like the coupled fluid-structure systems we are dealing with.

Main results

In this article, we consider a set of submerged deformable bodies, not necessarily neutrally buoyant, self-propelled partly by shape-changes and partly by thrusters, inside a potential fluid flow. The system fluid-bodies occupies either a bounded, partially bounded or unbounded domain of R N (N = 2, 3). Shape-changes as well as the running of the thrusters are assigned as functions of time and generically termed controls. We denote by F the domain of the fluid, by B the domain of the n bodies and we denote:

• (SN) the system of coupled equations in which the dynamic of the fluid is governed by Euler's equations and the dynamics of the bodies by Newton's laws.

• (SL) the system of equations obtained by applying Hamilton's principle to the Lagrangian function of the overall system fluid-bodies.

We say that both systems are equivalent if each one can be deduced from the other. We shall prove:

Theorem 1.1 If the boundary of the fluid domain is smooth enough and the thrusters inactive then (SN) and (SL) are equivalent.

The smoothness assumption is fundamental. Indeed, consider the simplest case in which a rigid solid is alone in a fluid of infinite extend. Denoting by T the stress tensor of the fluid and by n the unitary vector normal to ∂B directed toward the interior of the solid, then the net hydrodynamic forces exerted by the fluid on the solid is:

∂B Tn dσ.
We claim that:

Theorem 1.2 There exists a nonsmooth-shaped rigid solid for which:

• The Lagrangian function of the system fluid-solid and next (SL) are well defined and lead to d'Alembert's paradox (the solid moves along a straight line with constant velocity).

• For the motion obtained above |Tn| / ∈ L 1 (∂B). Therefore, Newton's law for the solid's linear momentum fails to apply and (SN) cannot be written.

Going back to the general smooth case, we denote by c the controls variable (i.e. the variable controlling the shapes of the bodies and the running of the thrusters) and by q the generalized variable (which corresponds to the degrees of freedom of the system fluid-bodies, the unknowns of our problem). We will show that for any given smooth function t → c(t), the Euler-Lagrange equation of motion takes the form:

d dt ∂L ∂ q - ∂L ∂q = F t , (1.1)
where L = L(q, q, c, ċ) is the Lagrangian function and F t a generalized force due to the work of the trusters. We will show that, once expanded, equation (1.1) turns into a second order non-linear ODE in q:

Mq = Γ(q, q, c, ċ, c), (1.2)
and we will compute explicitly the mass matrix M and the Christoffel symbol Γ to end up with an expression which involves only data of the problem on ∂B.

Numerical simulations

Equation (1.2) is the core of the Biohydrodynamics Matlab Toolbox (BhT). BhT is free, distributed under licence GNU GPL. A complete description of the toolbox is available on its home page:

http://bht.gforge.inria.fr/

In Section 7 we use it to perform numerical simulations. We show how two rigid moving solids deeply interact in a fluid when they are close what proves that the hydrodynamically decoupled bodies assumption is not always relevant. Further, we illustrate how BhT allows to study in terms of power and torques the forward motion of a pair of fish-like swimming articulated rigid bodies.

Organization of the paper

In the next subsection, we recall the method described for instance in [START_REF] Kanso | Locomotion of articulated bodies in a perfect fluid[END_REF] to derive the Euler-Lagrange equations of motion under hypotheses (H1-4) listed in Subsection 1.2. Section 2 is devoted to the setting of the general problem: topological and geometrical data are made precise, shape-changes and thrusters actuations are properly described as well as the configuration and phase spaces.

In Section 3, we use Lagrangian formalism to obtain formally the Euler-Lagrange equations. For this purpose, we compute the energy of the system and introduce the definitions of allowable controls and Lagrangian solution for the self-propelled motion problem.

In Section 4, we use Newton's laws to get the Euler-Lagrange equations. We end up with a system of coupled equations (some of them being PDE's and other ODE's). We define symmetrically to what has been done in the preceding section what is a Newtonian solution for the self-propelled motion problem.

Newtonian and Lagrangian solutions are next compared in Section 5, distinguishing between models with and without thrusters. The section ends up with a counterexample of a nonsmooth-shaped solid for which Newtonian formalism fails in providing the equations of motion whereas they can be obtained with Lagrangian approach.

In Section 6, we consider back the formal equations obtained at the end of Section 3 and we expand it in the form of a system of second order ODE's whose unknowns are the degrees of freedom of the system fluid-bodies. These ODE's are used in the core of BhT which is used in Section 7 to perform numerical simulations: in a first subsection, we illustrate how distant moving rigid solids interact in a fluid and a in the second subsection we study a pair of fish-like swimming articulated rigid bodies.

The appendix contains technical results.

A particular case

In this section only, let us assume that:

(H1) There is only one body in the fluid.

(H2) The overall system body-fluid fills the whole space (identified with R 3 ).

(H3) The body is neutrally buoyant.

(H4) The system body-fluid is at rest at time t = 0.

The shape of the body is described in a moving frame R M attached to the body by a so-called shape variable s.

The location of R M with respect to a Galilean fixed frame is given by g ∈ G := SE(3) (the Euclidian group) and its velocity (expressed in R M ) by a pair (Ω, v), where v is the translational velocity and Ω the angular velocity (we use the same notation as in [23, chap 15] we refer to for details).

Because of (H3), the Lagrangian of the system fluid-body is merely equal to the kinetic energy:

L := 1 2 (Ω, v, ṡ)M(s)   Ω v ṡ   ,
where M(s) is the mass matrix of the system (see e.g. [START_REF] Mason | Propulsion and control of deformable bodies in an ideal fluid[END_REF][START_REF] Kanso | Locomotion of articulated bodies in a perfect fluid[END_REF][START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF]). The fundamental point is that under assumptions (H1-3), the mass matrix depends on s only. It is a consequence of the isotropy of the configuration (as seen by an observer attached to the body, all of the positions and directions in the fluid are equivalent). In [START_REF] Kanso | Locomotion of articulated bodies in a perfect fluid[END_REF], this idea is expressed by saying that the Lagrangian is invariant under superimposed rigid motion.

The mass matrix is symmetric (as being the polarization of a quadratic form) and can be displayed into the following bloc form:

M(s) = M r (s) N(s) N T (s) P(s) ,
so that the Lagrangian can be expanded:

L = 1 2 (Ω, v)M r (s) Ω v + (Ω, v)N(s)ṡ + 1 2 ṡT P(s)ṡ.
Following [START_REF] Lamb | Hydrodynamics. Cambridge Mathematical Library[END_REF], we introduce the angular and translational impulses:

Π P := M r (s) Ω v ,
as well as the impulses related to the deformations:

Λ L := N(s)ṡ.
For any pair (g, ġ) of T G (the tangent bundle of G), we denote (q, q) its coordinates in a chart U ×E and (Ω, v) stand for the related angular and translational velocities. We compute that:

d dt ∂ ∂ q Ω v • q * - ∂ ∂q Ω v • q * = Ω * × Ω Ω * × v -Ω × v * , ∀ q * ∈ E,
where the pair (Ω * , v * ) corresponds to the coordinates (q, q * ). We next easily obtain that:

d dt ∂L ∂ q • q * - ∂L ∂q • q * = d dt Π + Λ P + L • Ω * v * + Π + Λ P + L • Ω * × Ω Ω * × v -Ω × v * , ∀ q * ∈ E.
Since the Euler-Lagrange equations of motion are (see e.g. [START_REF] Marsden | Introduction to mechanics and symmetry[END_REF]Theorem 7.3.3

page 187]): d dt ∂L ∂ q • q * - ∂L ∂q • q * = 0, ∀ q * ∈ E,
we get:

d dt (Π + Λ) = (Π + Λ) × Ω + (P + L) × v, d dt (P + L) = (P + L) × Ω.
Under assumption (H4), this system can be integrated to obtain:

Π + Λ P + L = 0, what leads to: Ω v = -(M r (s)) -1 N(s)ṡ.
This expression is the one obtained in [START_REF] Kanso | Locomotion of articulated bodies in a perfect fluid[END_REF]. It is very convenient to study the motion of the shape-changing body since it gives the velocity with respect to the shape variable. Nevertheless getting the Euler-Lagrange equations under such a simple form rests deeply on assumptions (H1-4) which ensure the isotropy (or symmetry using formalism of [START_REF] Marsden | Introduction to mechanics and symmetry[END_REF]) of the configuration. This isotropy is broken when either:

• the fluid is partially or totally bounded,

• there is several bodies in the fluid,

• or the body is not neutrally buoyant.

In such cases, the mass matrix depends also on q turning out the problem in a much more involved one. All of these cases are addressed in this article.

2 Setting of the problem

Topological and geometrical data

Let (e 1 , . . . , e N ) denote a reference Galilean frame by which we identify the physical space to R N (N = 2 or 3). At any time t > 0 we denote by B := ∪ n i=1 B i the open region of R N occupied by the n deformable bodies (n ∈ N, n ≥ 1) and by F the domain of the surrounding fluid which is assumed to be also open and connected. The set M, independent of time, is next defined as the disjoint union:

M := B 1 ∪ . . . ∪ B n ∪ F.
The boundary ∂F of the fluid is assumed to be a lipschitz continuous surface (N = 3) or curve (N = 2). This regularity allows one to define almost everywhere on ∂F the unitary normal directed towards the exterior of the fluid:

n : ∂F → S N -1 ,
where S N -1 stands for the unitary sphere of R N .

In the sequel, the exterior boundary ∂F \ ∂B is assumed to be bounded (so the fluid domain is either bounded or an exterior domain).

Kinematics of the shape-changing bodies

Systems of coordinates

The coordinates in the frame (e j ) are denoted with lowercase letters x = (x 1 , . . . , x N ) and are commonly called the spatial coordinates (see [23, chap 15]).

The material coordinates of the i-th body, denoted X = (X 1 , . . . , X N ), are given relatively to an orthogonal frame (E i 1 , . . . , E i N ) whose origin is chosen to be at the center of mass of the body at time t = 0. In this frame and at time t = 0, the i-th body occupies the domain B 0 i (usually termed reference configuration). Attached to the i-th body, we define also a moving frame (e i * 1 , . . . , e i * N ). We choose it such that its origin coincides at any time with the center of mass of the body and we denote by x * = (x * 1 , . . . , x * N ) the related so-called body coordinates. In this frame and at any time the body occupies the region B * i . We can assume in addition that at time t = 0, (E i j ) = (e i * j ) and hence that B * i = B 0 i (see figure 2).

Shape space

The deformation of the i-th body is described with respect to the body-fixed frame (e i * j ) by a diffeomorphism φ i :

B 0 i → B * i of class C m (m ∈ N, m ≥ 1)
equal to the identity outside a large ball B i (see figure 1). We denote by D i the set consisting of all such diffeomorphisms. It is thoroughly described in Apppendix, Section A. We have therefore with our notation:

x * = φ i (X).
Let Dφ i stands for the Jacobian matrix of φ i and denote Jφ i := | det Dφ i |. Then the related volume elements satisfy the identity:

dx * = Jφ i (X)dX.
Likewise, the surface elements dσ * and dσ 0 of ∂B * and ∂B 0 respectively are linked by the formula:

dσ * (x * ) = J σ φ i (X)dσ 0 (X), B i φ i B 0 i e i * 3 E i 1 E i 3 E i 2 B i e i * 2 B * i e i * 1 X x * := φ i (X)
Figure 1: Reference configuration B 0 i and domain of the deformed body

B * i := φ i (B 0 i ).
where J σ φ i := |(Dφ i ) -1T n|Jφ i is the tangential Jacobian.

The density ρ * i of the deformed i-th body is deduced from the given density ρ 0 i of the body at time t = 0 by the conservation-of-mass principle:

ρ * i = ρ 0 i (φ -1 i )Jφ -1 i . (2.1)
We denote then dm 0 i := ρ 0 i dX, the element of mass in B 0 i and likewise:

dm * i := ρ * i dx * ,
is the element of mass in B * i . Since we assume that the deformations are due to inner forces and torques only, Newton's laws ensure that the center of mass and the angular momentum of the body with respect to its attached frame (e i * j ) remain unchanged when it undergoes shape-changes. We get:

B * i x * dm * i = B 0 i X dm 0 i = 0,
which can be rewritten, upon a change of variables and taking into account (2.1):

B 0 i φ i (X) dm 0 i = 0. (2.2a)
Differentiating with respect to time, we get:

B 0 i φi dm 0 i = 0, (2.2b)
or equivalently:

B * i φi (φ -1 i ) dm * i = 0.
The same arguments yield for the angular momentum:

r i B 0 i B i e 3 E i 3 e 2 F e i * 3 E i 2 E i 1 e i * 1 e 1 X e i * 2 x := g i φ i (X)
Figure 2: The fixed frame (e j ) and the moving frame (e i * j ) with e i * j = R i E i j . The domain of the body is

B i := R i φ i (B 0 i ) + r i .
• When N = 3:

B 0 i φi × φ i dm 0 i = 0, (2.2c)
or equivalently:

B * i φi (φ -1 i ) × x * dm * i = 0.
• When N = 2:

B 0 i φi • φ ⊥ i dm 0 i = 0, (2.2d)
or equivalently:

B * i φi (φ -1 i ) • x * ⊥ dm * i = 0, where x ⊥ := (-x 2 , x 1 ) for all x = (x 1 , x 2 ) ∈ R 2 .
It will be made clear later on that in the following cases:

(C1) when N = 2, for bounded and unbounded domains F, (C2) for bounded F only when N = 3, we must have in addition:

B * i ∇ • φi dx * = ∂B * i φi (φ -1 i ) • n dσ * = 0, (2.2e)
expressing the conservation of the overall volume of the fluid.

Rigid motion

As already mentioned, G := SE(N ) stands for the rigid body group (or Euclidian group). The elements g ∈ G consists of pairs (R, r) where R ∈ SO(N ) is a rotation matrix and r ∈ R N . They would sometimes be identified with the matrix:

g = R r 0 1 ,
and for any vector x = (x 1 , . . . , x N ) T ∈ R N , we denote

gx 1 := R r 0 1 x 1 . (2.3)
In particular we will made no difference in the notation between the pair (R, r) and the associated isometry of R N defined by (2.3). At any time and for any i = 1, . . . , n, there exists

g i = (R i , r i ) ∈ SE(N ) such that e i * j = R i E i j ,
and r i gives the location of the center of mass of the i-th body (see figure 2). With our notation: x = g i x * , and B i = g i (B * i ). We deduce, since Jg i = 1 and J σ g i = 1, that: dx = dx * , and dσ = dσ * .

With evident notation we can define the density of the i-th body by:

ρ i (x) := ρ * i (x * ),
and the element of mass by: dm i = ρ i dx.

Velocities of the bodies

According to the relation:

x = g i φ i (X),
we deduce that the Eulerian velocity at a point x of B i is:

v i (x) = ġi g -1 i x + R i φi (φ -1 i (g -1 i x)),
where ġi := ( Ṙi , ṙi ), (g i , ġi ) ∈ T G and (φ i , φi ) ∈ S i . Classically, we define the vector ω i (when N = 3) or the scalar ω i (when N = 2) by: Ṙi

R T i x = ω i × x when N = 3, ω i x ⊥ when N = 2,
and we get for all x ∈ B i :

ġi g -1 i x = ω i × (x -r i ) + ṙi when N = 3, ω i (x -r i ) ⊥ + ṙi when N = 2. (2.4a)
We can rewrite the Eulerian velocity:

v i (x) = (ω i × (x -r i ) + ṙi ) + R i φi (φ -1 i (g -1 i x)) when N = 3, (ω i (x -r i ) ⊥ + ṙi ) + R i φi (φ -1 i (g -1 i x)) when N = 2. (2.5)
The body (or convective) velocity reads, for all x * ∈ B * i :

v * i (x * ) = (Ω i × x * + R T i ṙi ) + φi (φ -1 i (x * )) when N = 3, (Ω i x * ⊥ + R T i ṙi ) + φi (φ -1 i (x * )) when N = 2,
where

Ω i := R T i ω i when N = 3, Ω i := ω i when N = 2.
Observe that v i is given in the fixed basis (e j ) whereas v * is given in the bodyattached basis (e i * j ).

Thrusters

The running of the thrusters of the i-th body is described by a function ν i ∈ H -1/2 (∂B 0 i ) (the dual of the Sobolev space H 1/2 (∂B 0 i ) defined for instance in [21, §7.3 pages 38-42]). The function ν i stands for the normal component of the velocity of the fluid blown out by the thrusters. It must satisfies in cases (C1-2) only:

ν i , 1 H -1/2 (∂Bi)×H 1/2 (∂Bi) = 0. (2.6)
When ν i is more regular, namely in L 2 (∂B i ), this condition turns into:

∂Bi ν i dσ = 0. (2.7)
The running of the thrusters can equivalently be defined by a function ν 0 i ∈ H -1/2 (∂B 0 i ) such that, in the cases (C1-2) mentioned above only:

ν 0 i J σ φ i , 1 H -1/2 (∂B 0 i )×H 1/2 (∂B 0 i ) = 0.
Both functions ν i and ν 0 i are linked through the identity:

ν i (x) = ν 0 i (X).

Configuration and phase spaces

The configuration space

The fluid is assumed to be incompressible and inviscid. We denote ρ f > 0 its constant density and we set:

dm f := ρ f dx,
the element of mass. The classical configuration space for such a fluid is the set Diff vol (F 0 , F) of the smooth, volume-preserving diffeormorphisms from F 0 onto F. Notice that this space depends on F and hence on all of the elements φ i and g i . Thus, the fluid cannot be described independently of the shapes and positions of the bodies.

We deduce that the configuration space Q of the system fluid-bodies is the set of all the triplets (φ, g, η) such that:

1. φ := (φ 1 , . . . , φ n ) ∈ D 1 ×. . .×D n , g := (g 1 , . . . , g n ) ∈ G n and η ∈ Diff vol (F 0 , F).
2. For any pair of indices (i, j) such that i = j:

Bi ∩ Bj := g i φ i ( B0 i ) ∩ g j φ j ( B0 j ) = ∅, Bi ∩ (R N \ M) := g i φ i (B 0 i ) ∩ (R N \ M) = ∅.
The latter condition means that the bodies are not allowed to touch or overlap each other or to go out of M.

The phase space

For any (φ, g, η, φ, ġ, η) ∈ T Q, we denote

u := η • η -1 , (2.8a)
the Eulerian velocity of the fluid. The slip boundary condition for inviscid fluid yields:

u • n = v i • n + ν i on ∂B i , (2.8b) u • n = 0 on ∂F \ ∂B, (2.8c)
where v i is defined by (2.5).

Definition 1

The phase space P consists of all the septuplets (φ, g, η, φ, ġ, η, ν) such that:

1. (φ, g, η, φ, ġ, η) ∈ T Q. 2. ν := (ν 1 , . . . , ν n ) ∈ H -1/2 (∂B 1 ) × . . . × H -1/2 (∂B n ).
3. For any index i, conditions (2.2) hold.

4.

For any index i, condition (2.6) holds.

Identities (2.8) hold.

The phase space is a vector bundle whose base space is Q and related projection π Q is defined by:

π Q (φ, g, η, φ, ġ, η, ν) := (φ, g, η).

Potential flow, modified configuration and phase space

We seek the Eulerian velocity of the fluid as the gradient of a potential function ψ:

u := ∇ψ in F. (2.9)
The incompressibility of the fluid entails ∇ • u = 0 and hence: ∆ψ = 0 in F.

(2.10a)

The boundary conditions (2.8) lead to Neumann boundary conditions for ψ, which read:

∂ n ψ = v i • n + ν i on ∂B i , (2.10b)
with the convention that ∂ n ψ = 0 on unspecified boundaries. When (C1-2) hold, we have to add a condition to ensure the uniqueness of ψ, namely:

∂F \∂B ψ dσ = 0. (2.10c)
Boundary value problem (2.10) is well-posed, as detailed in Appendix, Section D. Therein, it is also explained where conditions (C1-2) come from. The fluid has no own degree of freedom left because ψ is completely determined by System (2.10). The element η can be canceled in the definition of the configuration space Q and the definition of the phase space can be simplified as follows:

Definition 2 The phase space P in the case of a potential flow consists of all the quintuplets p := (φ, g, φ, ġ, ν) such that:

1. (φ, g, φ, ġ) ∈ T Q. 2. ν := (ν 1 , . . . , ν n ) ∈ H -1/2 (∂B 1 ) × . . . × H -1/2 (∂B n ).
3. For any index i, conditions (2.2) hold. 4. For any index i, condition (2.6) holds.

The projection π Q must be modified accordingly.

3 Lagrangian description

Mass and inertia tensors

We denote by m i > 0 the mass of the i-th body while its inertia tensor in the moving basis (e i * j ) reads, when N = 3:

I * i := B * i |x * | 2 Id -x * ⊗ x * dm * i ,
or equivalently upon a change of variables:

I * i := B 0 i |φ i (X)| 2 Id -φ i (X) ⊗ φ i (X) dm 0 i .
It can also be expressed with respect to the fixed frame (e j ):

I i = R i I * i R T i .
When N = 2 it is a scalar:

I i := B * i |x * | 2 dm * i ,
or equivalently:

I i := B 0 i |φ i (X)| 2 dm 0 i .

Kinetic and potential energies

Due to relations (2.2), the kinetic energy of the i-th body is, when N = 3:

K b i := 1 2 m i |ṙ i | 2 + 1 2 Ω i • I * i Ω i + 1 2 Bi |R i φi (φ -1 i (g -1 i x))| 2 dm i ,
and when N = 2:

K b i := 1 2 m i |ṙ i | 2 + 1 2 |Ω i | 2 I i + 1 2 Bi |R i φi (φ -1 i (g -1 i x))| 2 dm i ,
the last term being the kinetic energy of deformation. It can be rewritten:

1 2 Bi |R i φi (φ -1 i (g -1 i x))| 2 dm i = 1 2 B 0 i | φi (X)| 2 dm 0 i .
We denote

K b := i K b i ,
the overall kinetic energy of the bodies. Setting the potential gravity function

G(x) := -gx N ,
(where g stands here for the standard gravity), and by virtue of Archimedes' principle, the potential energy associated with the buoyant force is, for the i-th body:

P i := - Bi G(x)(ρ i -ρ f ) dx, (3.1a)
or, upon a change of variables:

P i := - B * i G(g i x * )(ρ * i -ρ f ) dx * . (3.1b)
The overall potential energy related to the buoyant force is:

P := i P i .
At last, the kinetic energy of the fluid is:

K f := 1 2 F |u| 2 dm f = 1 2 F |∇ψ| 2 dm f ,
and the Lagrangian function, defined on the phase space P reads:

L := K f + K b -P. (3.2)

Euler-Lagrange equations of motion of the self-propelled motion problem

As already mentioned, shape-changes and thrusters actuation are referred to as controls.

Definition 3 (Allowable controls) We call allowable controls, any pair (φ(t), ν(t)) of smooth functions that satisfy, for all t: A compatibility condition for the initial positions and velocities of the bodies is also needed for the problem to be well-posed.

1. φ(t) := (φ 1 (t), . . . , φ n (t)) ∈ D 1 × . . . × D n , 2. ν(t) := (ν 1 (t), . . . , ν n (t)) ∈ H -1/2 (∂B 1 ) × . . . × H -1/2 (∂B n ), 3
Definition 4 A pair (g 0 , ġ0 ) ∈ T G n is said satisfying the compatibility condition at a given time t 0 relatively to a pair of allowable controls (φ(t), ν(t)) when (φ(t 0 ), g 0 , dφ dt (t 0 ), ġ0 , ν(t 0 )) ∈ P.

(3.3)
Consider now the projection π T G n : P → T G n , defined by:

π T G n (p) := (g, ġ).
Let U × E be:

• a local chart of T G n in a neighborhood of (g 0 , ġ0 ),

• included in the image by π T G n of a local chart of the phase space P nearby an element p 0 := (φ 0 , g 0 , φ0 , ġ0 , ν 0 ).

We denote (q, q) the variables in U × E and ∂L/∂q and ∂L/∂ q are the related partial derivatives of L.

Definition 5 (Lagrangian solution)

We call Lagrangian solution of the selfpropelled motion problem with given allowable controls (φ(t), ν(t)) and initial data (g 0 , ġ0 ) satisfying the compatibility condition (3.3), any twice continuously differentiable function t → g(t) ∈ G n defined in a neighborhood of t 0 and such that:

1. g(t 0 ) = g 0 and dg dt (t 0 ) = ġ0 .
2. The coordinates of (g, ġ) are (q, q).

3. For all t nearby t 0 , (q, q) solve the EDO:

d dt ∂L ∂ q (p(t)) • q * - ∂L ∂q (p(t)) • q * = F t (p(t)) • q * , ∀ q * ∈ E, (3.4)
where p(t) := (φ(t), g(t), ∂φ ∂t (t), dg dt (t), ν(t)), and F t stands for generalized forces induced by the thrusters actuation.

Since the actuation of the thrusters has been modeled as non-holonomic constraints (namely as constraints on the velocity of the fluid), the Lagangian formulation required the use of generalized forces (see for instance [11, §2.4 pages 45-50]). The expression of these forces will be given later on, in Proposition 5.1.

Newtonian description

We assume that the fluid flow is driven by Euler's equations. In such a fluid, the stress tensor is

T := (-p + ρ f G)Id, (4.1a)
where p stands for the pressure of the fluid. Remind that we have denoted by u the Eulerian velocity of the fluid and by ρ f its constant density. The equations of motion read:

ρ f d dt u = ∇ • T in F, (4.1b)
∇ • u = 0 in F, (4.1c)
where

d dt u := ∂ t u + (u • ∇)u, (4.1d)
is the particular derivative. We add the already mentioned slip boundary condition for inviscid fluid:

u • n = v i • n + ν i on ∂B i , (4.1e) u • n = 0 on ∂F \ ∂B, (4.1f)
where v i is the velocity of the i-th body defined by (2.5). Newton's laws apply to the bodies and yield:

d dt (m i ṙi ) = - ∂Bi Tn dσ -g(m i -m i f )e 3 , (4.2a)
d dt (I i ω i ) = - ∂Bi (x -r i ) × Tn dσ + gm i f (r i f -r i ) × e 3 when N = 3, (4.2b) d dt (I i ω i ) = - ∂Bi (x -r i ) ⊥ • Tn dσ + gm i f (r i f -r i ) ⊥ • e 2 when N = 2, (4.2c)
where

m i f := Bi dm f , (4.2d)
r i f := 1 m i f Bi x dm f . (4.2e)
We add initial conditions for the bodies and the fluid:

u(0) = u 0 , g(0) = g 0 , ġ(0) = ġ0 , (4.3)
with (g 0 , ġ0 ) satisfying the compatibility condition (3.3). The system of equations (4.1), (4.2) and (4.3) allows to described thoroughly the self-propelled motion of the bodies.

Potential flow

Seeking now the Eulerian velocity of the fluid as the gradient of a potential function, u = ∇ψ, System (4.1) must be replaced by (2.10). The pressure (and hence also the stress tensor T) is obtained by means of Bernoulli's formula:

p = p 0 -ρ f ∂ t ψ + |∇ψ| 2 2 , (4.4)
where p 0 is constant in F (but depends on time). As with the Lagrangian approach, the fluid has no own degree of freedom left. In particular, the initial data for the fluid in (4.3) must be canceled.

Definition 6 (Newtonian solution)

We call Newtonian solution of the self-propelled problem with allowable given controls (φ(t), ν(t)) and compatible initial data (g 0 , ġ0 ), any continuously differentiable function t → g(t) = (g 1 (t), . . . , g n (t)) ∈ G n defined in a neighborhood of t 0 and such that:

1. g(t 0 ) = g 0 and dg dt (t 0 ) = ġ0 .

2. All of the pairs (ω i , ṙi ) associated with (g i , ġi ) := (g i , dg i /dt) by relations (2.4) solve the system of coupled equations made up of (2.10), (4.1a), (2.5), (4.2) and (4.4).

Newtonian versus Lagrangian modeling

This section is concerned with comparing Lagrangian solutions of Definition 5 with Newtonian solutions of Definition 6. Remind that the thrusters have been modeled as non-holonomic constraints on the system fluid-bodies (more precisely as constraints on the velocity of the fluid). For such constraints and within Lagrangian formalism, the use of generalized forces is required. Generally, the expression of such forces is not easy to determine. To avoid this difficulty, we shall assume in a first time that there is no thruster and show that Lagrangian solutions and Newtonian solutions match. In a second time, we will use the Newtonian approach to derive the expression of the generalized forces induced by the thrusters in the Lagrangian model.

Equivalence for the model without thruster

In this subsection only, we assume that ν = 0.

Theorem 5.1 Assume that ∂F is of class C 1,1 . Then Newtonian solutions of Definition 6 and Lagrangian solutions of Definition 5 (with

F t = 0) match.
The proof of this theorem rests on a set of lemma. Remind that ψ is the potential of the fluid defined by (2.10) and define, for all p ∈ P the mapping Λ p ∈ H 1 (R N ) (the dual of the Sobolev space H 1 (R N )) by:

Λ p (ϕ) := F ∇ϕ • ∇ψ dm f , ∀ ϕ ∈ H 1 (R N ). (5.1)
Let U × E be the chart of T G n already defined earlier, above Definition 5. Using the same notation as therein we claim:

Lemma 5.1 Under the hypothesis of Theorem 5.1 and for any smooth function

t → p(t) := (φ, g, dφ dt , dg dt , ν) ∈ P,
we have:

d dt ∂Λ p ∂ q • q * - ∂Λ p ∂q • q * = 0 in H 1 (R N ) , ∀ q * ∈ E.
Proof: We give the proof for N = 3, the case N = 2 being similar. We denote (g, ġ) and (g, ġ * ) the elements of T G n of coordinates (q, q) and (q, q * ) respectively. We associate to (g, ġ) and (g, ġ * ) the pairs (ω i , ṙi ) and (ω * i , ṙi * ) like in (2.4). Likewise, to any p := (φ, g, φ, ġ) ∈ P , we associate the velocity fields v i defined on ∂B i by (2.5).

All of the coordinates being fixed but q, the mapping q ∈ E → ψ ∈ H 1 (F) is linear (because q → ∂ n ψ is). We deduce that the function ψ * := (∂ψ/∂ q)• q * is harmonic in F and satisfies the following Neumann boundary conditions (homogeneous on unspecified boundaries):

∂ n ψ * = (ω * i × (x -r i ) + ṙ * i ) • n on ∂B i .
Invoking Reynolds's formula and next Green's formula (see (E.8) in the Appendix), we deduce that:

∂Λ p ∂ q • q * (ϕ) = F ∇ϕ • ∇ψ * dm f (5.2a) = i ∂Bi ϕ(ω * i × (x -r i ) + ṙ * i ) • n dσ f . (5.2b)
On the other hand, we have also, according to Reynolds's formula again:

∂ ∂q F ϕ dm f • q * = i ∂Bi ϕ(ω * i × (x -r i ) + ṙ * i ) • n dσ f ,
which leads to:

∂Λ p ∂ q • q * (ϕ) = ∂ ∂q F ϕ dm f • q * .
Differentiating with respect to t this identity and exchanging both derivatives in the right hand side, it comes:

d dt ∂Λ p ∂ q • q * (ϕ) = ∂ ∂q d dt F ϕ dm f • q * .
Reynolds's formula yields:

d dt F ϕ dm f = i ∂Bi ϕ(v i • n) dσ f = ∂F ϕ∂ n ψ dσ f . (5.3)
It remains to apply Green's formula once more time to the right hand side:

d dt F ϕ dm f = F ∇ϕ • ∇ψ dm f = Λ p (ϕ),
and we get the conclusion of the lemma.

The proof is based on the extensive use of both Green's and Reynolds's formula. These formula required the domain F to be smooth enough (Lipschitz continuous for Green's formula and C 1,1 i.e. continuously differentiable with Lipschitz continuous derivative for Reynolds's formula).

We use this lemma to prove:

Lemma 5.2 Under the hypotheses of Lemma 5.1, we have, for all q * ∈ E:

d dt ∂K f ∂ q • q * - ∂K f ∂q • q * = i ∂Bi ∂ t ψ + |∇ψ| 2 2 (ω * i × (x -r i ) + ṙ * i ) • n dσ f .
Proof: From the identity:

K f = 1 2 Λ p (ψ),
we deduce that, for all q * ∈ E:

∂Λ p (ψ) ∂ q • q * = 1 2 ∂Λ p ∂ q • q * (ψ) + 1 2 Λ p (ψ * ).
According to (5.2a), both terms in the right hand side are equal. We obtain therefore after some algebra:

d dt ∂K f ∂ q • q * - ∂K f ∂q • q * = d dt ∂Λ p ∂ q • q * - ∂Λ p ∂q • q * (ψ) + ∂Λ p ∂ q • q * (∂ t ψ) + 1 2
∂Λ p ∂q • q * (ψ) - 1 2 Λ p (ψ † ),
where we have defined:

ψ † := ∂ψ ∂q • q * .
Lemma 5.1 tells us that the first term in the right hand side vanishes and according to (5.2b) the second one can be rewritten as:

∂Λ p ∂ q • q * (∂ t ψ) = i ∂Bi ∂ t ψ(ω * i × (x -r i ) + ṙ * i ) • n dσ f .
Going back to the definition of Λ p and applying Reynolds's formula, we get:

∂Λ p ∂q • q * (ψ) = i ∂Bi |∇ψ| 2 (ω * i × (x -r i ) + ṙ * i ) • n dσ f + Λ p (ψ † ),
and the proof is completed.

Concerning the buoyant force, we use the same notations as in the proof of Lemma 5.1 and Lemma 5.2 and we deduce easily from the expressions (3.1):

Lemma 5.3 For any i ∈ {1, . . . , n}:

d dt ∂P i ∂ q • q * - ∂P i ∂q • q * = - ∂P i ∂q • q * , ∀ q * ∈ E, and
- ∂P i ∂q • q * = -g(m i -m i f )e 3 • ṙ * i + gm i f (r i f -r i ) × e 3 • ω * i .
The last lemma we need is a straight consequence of Lemma B.2 (in the Appendix):

Lemma 5.4 For all i ∈ {1, . . . , n}:

d dt ∂K b i ∂ q • q * - ∂K b i ∂q • q * = d dt (m i ṙi ) • ṙ * i + d dt (I i ω i ) • ω * i , ∀ q * ∈ E.
We have now all the material for the proof of Theorem 5.1.

Proof of Theorem 5.1: Remind that

L = i K b i + K f - i P i,
and apply successively Lemma 5.2, 5.3 and 5.4 to obtain Newton's laws (4.2).

The regularity required in the theorem for ∂F is crucial as proven in Subsection 5.3.

Model with thrusters

Let us define ψ t as being the potential of the flow related to the thrusters. It solves:

-∆ψ t = 0 in F,
and satisfies Neumann boundary conditions (homogeneous on unspecified boundaries):

∂ n ψ t = ν i on ∂B i .
Consider next the mapping Λ t p ∈ H 1 (R N ) defined like in (5.1) by:

Λ t p (ϕ) := F ∇ϕ • ∇ψ t dm f , ∀ ϕ ∈ H 1 (R N ).
(5.4)

Proposition 5.1 Assume that ∂F is of class C 1,1 . Then, the generalized force arising in Definition 5 is defined by:

F t • q * = - ∂Λ t p ∂q • q * (ψ), ∀ q * ∈ E. (5.5)
The proof of this proposition is similar to the proof of Theorem 5.1 after Lemma 5.1 has been replaced by: Lemma 5.5 Under the hypothesis of Proposition 5.1 and for any smooth function

t → p(t) := (φ, g, dφ dt , dg dt , ν) ∈ P,
we have:

d dt ∂Λ p ∂ q • q * - ∂Λ p ∂q • q * = - ∂Λ t p ∂q • q * in H 1 (R N ) , ∀ q * ∈ E.
Proof: We have no longer

∂ n ψ = v i • n on B i ,
but rather:

∂ n ψ = v i • n + ν i on B i ,
where the velocity fields v i is defined by (2.5). It entails that identity (5.3) turns out to be:

d dt F ϕ dm f = ∂F ϕ(∂ n ψ -∂ n ψ t ) dσ f .
Apart from this modification, the proof is similar to the one of Lemma 5.1.

Non equivalence counterexample

In this section, we show that the conclusion of Theorem 5.1 does not hold any longer when the smoothness assumption on ∂F is missing. We exhibite a two dimensional nonsmooth solid moving vertically in a fluid and for which Newtonian solution of Definition 6 does not exist although Lagrangian solution of Definition 5 be well defined.

B 0 ξ 2 ξ 1 D -3i/2 -2 2 i 1 0 w 1 0 3i/2 w 2 w-complex plane F 0 ξ-complex plane w := f (ξ) ξ := f -1 (w)
Figure 3: The function f is a conformal mapping from C \ D onto F 0 .

Theorem 5.2 Consider the solid pictured on the right of Figure 3 and assume that the system fluid-solid fills the whole space. Then for suitable initial data, Lagrangian formalism applies and leads to d'Alembert's paradox (uniform straight motion). However, with notation of Section 4:

|Tn| / ∈ L 1 (∂B),
and the net hydrodynamic force exerted by the fluid on the solid:

∂B Tn dσ,
is not well defined. Consequently, Newton's laws (4.2) fail to apply.

Proof: We use complex analysis to compute explicitly the potential function, so we identify the space with C (the set of the complex numbers). However we will mix the complex notation ξ = ξ 1 + iξ 2 (i 2 = -1) with the real one ξ = (ξ 1 , ξ 2 ) T and even with the polar coordinates (r, θ), r = |ξ| and θ = Arg (ξ) (i.e. ξ = re iθ ).

The solid has an axis of symmetry and we assume that at the time t = 0 this axis coincides with the imaginary axis and that the center of mass of the rigid body lies at the origin. If we assume in addition that the density of the rigid body is constant, then, for symmetry reasons, its center of mass will remain on the imaginary axis for all t > 0 whenever the initial velocity of the solid is vertical. Therefore, within this configuration, the system has only one degree of freedom.

We denote by ih the position of the center of mass of the solid and i ḣ stands for its velocity. Remind that we denote B 0 the domain of the solid at time t = 0 and F 0 is the domain of the surrounding fluid. The boundary ∂B 0 is the image of the unitary circle by the mapping:

f (ξ) := 3ξ 2 -1 2ξ + 2ξ 3ξ 2 -1 , (5.6)
which is a conformal mapping from C \ D, the exterior of the unitary disk D, onto F 0 . We denote w := f (ξ) for all ξ ∈ C \ D. Following [9, Chap XIII], we introduce, for all w ∈ F, ψ(w) := ψ 1 (w) + iψ 2 (w) the complex holomorphic potential function.

Commonly, ψ 1 is the actual real valued potential function while ψ 2 is the stream function. The potential can be rewritten in the form:

ψ(w) = ḣϕ(w -ih), (5.7)
where ϕ(w) := ϕ 1 (w) + iϕ 2 (w), the functions ϕ 1 and ϕ 2 being both defined and harmonic in the fixed domain F 0 and independent of h. In order to compute ϕ, we introduce the function:

φ := ϕ • f,
and we denote φ 1 := (φ) and φ 2 := (φ). According to [9, formula (5) page 222], φ 2 solves a Dirichlet boundary value problem:

-∆φ 2 = 0 in C \ D, φ 2 = -(f ) on ∂D, |∇φ 2 | → 0 as |ξ| → +∞.
The last condition tells us that the fluid is quiescent at infinity. We can extend f (ξ) and f (1/ξ) as Laurent's series at ξ = 0. We get:

f (ξ) = 3 2 ξ + 1 6 ξ -1 + k≥1 2 3 k+1 ξ -2k-1 ∀ ξ ∈ C, |ξ| > 1/3, f (1/ξ) = 3 2 ξ -1 + 1 6 ξ + k≥1 2 3 k+1 ξ 2k+1 ∀ ξ ∈ C, |ξ| < 3.
Writing then that φ 2 (e iθ ) = -(f (e iθ ) + f (e -iθ ))/2 we obtain the expression of φ 2 (e iθ ) in the form of a Fourier's series:

φ 2 (e iθ ) = - 5 6 cos(θ) - k≥1 3 -k-1 cos((2k + 1)θ), ∀ θ ∈ [-π, π[,
and we deduce that, (using here the polar coordinates (r, θ), due to the geometry of the domain):

φ 2 (r, θ) = - 5 6 cos(θ) r - k≥1 3 -k-1 cos((2k + 1)θ) r 2k+1 , ∀ θ ∈ [-π, π[, ∀ r ≥ 1. (5.8)
Observe that φ 2 can actually be extended as an harmonic function outside the disk of center 0 and radius 1/ √ 3.

Lagrangian formalism

Let us compute now the Lagrangian function L(h, ḣ) of the system fluid-solid. We assume that the density of the solid and the density of the fluid, ρ f > 0 are equal. Hence the solid is neutrally buoyant and the Lagrangian function reduces to the kinematic energy of the system, namely:

L(h, ḣ) = 1 2 m| ḣ| 2 + F ρ F |∇ψ 2 | 2 dw ,
where m > 0 stands for the mass of the solid. Notice that |∇ψ 2 | = |∇ψ 1 | for ψ is an holomorphic function. This expression can be turned into, upon a change of variables and applying Green's formula:

L(h, ḣ) = 1 2 | ḣ| 2 m + ∂D φ 2 ∂ n φ 2 dσ f ,
and we deduce that the Lagrangian function does actually not depend on h. Euler-Lagrange's equation of motion is merely:

d dt ∂L ∂ ḣ ( ḣ) = m + ∂D φ 2 ∂ n φ 2 dσ f ḧ = 0, (5.9)
what, after integrating, leads to so-called D'Alembert's paradox:

ḣ(t) = ḣ0 , ∀ t > 0.
The solid has a straight uniform velocity.

Newtonian formalism

Let us now invoke Newton's laws with the same configuration. Making use of Bernoulli's formula, we get:

m ḧ = - ∂B ∂ t ψ 1 + 1 2 |∇ψ 1 | 2 n 2 dσ f , (5.10)
where n = (n 1 , n 2 ) T is the unitary normal to ∂B directed toward the interior of B.

According to (5.7), we obtain that:

∂ t ψ 1 (w) = ḧϕ 1 (w -ih) -| ḣ| 2 ∂ϕ 1 ∂w 2 (w -ih).
Plugging this expression into (5.10), and upon a change of variables, we get:

m + ∂D φ 2 ∂ n φ 2 dσ f ḧ = -| ḣ| 2 ∂B 1 2 |∇ϕ 2 | 2 - ∂ϕ 1 ∂w 2 n 2 dσ f , (5.11)
where we are use again that |∇φ 1 | = |∇φ 2 | and |∇ϕ 1 | = |∇ϕ 2 | for φ and ϕ are both holomorphic functions. This equation fits with (5.9) if and only if the right hand side term vanishes. We are going to show that it is actually undetermined, in the sense that the amount inside the integral does not belong to L 1 (∂B). We have, on the one hand:

∂B ∂ϕ 1 ∂w 2 n 2 dσ f ≤ ∂B |∇ϕ 1 |dσ f ,
and arguing again that |∇ϕ 1 | = |∇ϕ 2 |, it comes, upon a change of variables:

Γ ∂ϕ 1 ∂w 2 n 2 dσ f ≤ (mes (∂B)) 1/2 ∂B |∇ϕ 2 | 2 dσ f 1/2 ≤ (mes (∂B)) 1/2 ∂D |∇φ 2 | 2 dσ f 1/2 < +∞.
On the other hand, going back to (5.11), we have, applying the same change of variables: In addition, one readily obtains that:

∂B |∇ϕ 2 | 2 n 2 dσ f = ρ f π -π |∇φ 2 (e iθ )| 2 |f (e iθ
|f (e iθ )| -2 (-e iθ f (e iθ )) = - 1 8 θ -1 + O(θ).
(5.14) Summarizing (5.13), (5.14) we get

|∇φ 2 (e iθ )| 2 |f (e iθ )| -2 (-e iθ f (e iθ )) = -γ 2 θ -1 + O(θ),
what proves that (5.12) is undetermined.

The velocity of the fluid is infinite nearby the singular points of ∂B as well as the hydrodynamic forces although the energy of the system be always finite.

This phenomenon is well known in airfoils theory (see [9, Chap. XII, §c. pages 211-212]). The problem can be overcame by introducing circulation in the fluid flow in such a way that the velocity be finite at the singular point. In other words, circulation is used as a new degree of freedom of the fluid flow and its value is set according to Kutta-Joukowsky's law.

However, this method does not apply with our problem since the solid under consideration has two singular points.

Explicit Euler-Lagrange equations

The aim of this section is to repeat in the general case (i.e. after relaxing hypotheses (H1-4)) the method carried out in Subsection 1.6. We start from the Euler-Lagrange equations under their abstract form (3.4)-(5.5) and we expand it into a more explicit form, convenient for numerical simulations.

Kirchhoff 's law

For all (φ, g, φ, ġ, ν) in the phase space P , we define the elementary potentials ψ r , ψ d and ψ t as being harmonic functions in F satisfying Neumann boundary conditions (homogeneous on unspecified boundaries):

∂ n ψ r = (ω i × (x -r i ) + ṙi ) • n on ∂B i when N = 3, (ω i (x -r i ) ⊥ + ṙi ) • n on ∂B i when N = 2, ∂ n ψ d = R i φi (φ -1 i (g -1 i x)) • n on ∂B i , ∂ n ψ t = ν i on ∂B i .
According to Kirchhoff 's law, we get the decomposition for the overall potential:

ψ = ψ r + ψ d + ψ t . (6.1)
When N = 3, Euler's angles, defined in the Appendix, Section B, provide a local charts of SO(3) and TSO(3) in which we denote

θ k = (θ 1 k , θ 2 k , θ 3 k ) and θk = ( θ1 k , θ2 k , θ3
k ) the coordinates of the pair (R k , Ṙk ) (remind that R k is the rotation associated with the rigid motion of the k-th body). It allows us to introduce for i = 1, . . . , 6n the functions ψ r i , also harmonic in F and such that: ∂ n ψ r i = e j • n on ∂B k for i = 6(k -1) + j, (6.2a)

∂ n ψ r i = ω e j × (x -r k ) • n on ∂B k for i = 6(k -1) + 3 + j, (6.2b)
for all k ∈ {1, . . . , n}, j ∈ {1, 2, 3} where the rotation vectors ω e j are also defined in the Appendix, Section B. Denoting as well r k = (r 1 k , r 2 k , r 3 k ) and ṙk = ( ṙ1 k , ṙ2 k , ṙ3 k ) the coordinates of r k and ṙk in the basis (e j ), we get a further decomposition of ψ r :

ψ r = k j ṙj k ψ r 6(k-1)+j + θj k ψ r 6(k-1)+3+j .
When N = 2, the pair (R k , Ṙk ) is merely parameterized by a pair (θ k , ω k ) ∈ R 2 (actually, we have in this case ω k = Ω k = θk , all of them scalar). We get only 3n elementary potentials ψ r i and the boundary conditions (6.3a) turn into:

∂ n ψ r i = e j • n on ∂B k for i = 3(k -1) + j, (6.3a) ∂ n ψ r i = (x -r k ) ⊥ • n on ∂B k for i = 3k, (6.3b)
for all k ∈ {1, . . . , n}, j ∈ {1, 2}. We deduce that in this case we have:

ψ r = k j ṙj k ψ r 3(k-1)+j + ω k ψ r 3k .

Mass matrices

The symmetric added mass matrix associated to the rigid motion of the bodies is defined, when N = 3, by:

M f r :=    F ∇ψ r 1 • ∇ψ r 1 dm f • • • F ∇ψ r 1 • ∇ψ r 6n dm f . . . . . . F ∇ψ r 6n • ∇ψ r 1 dm f • • • F ∇ψ r 6n • ∇ψ r 6n dm f    .
It is a matrix of size 3n × 3n when N = 2. The generalized coordinates read:

• When N = 3: q := (r 1 , θ 1 , . . . , r n , θ n ) and q := (ṙ 1 , θ1 , . . . , ṙn , θn ).

• When N = 2: q := (r 1 , θ 1 , . . . , r n , θ n ) and q := (ṙ 1 , ω 1 , . . . , ṙn , ω n ), and we denote from now on:

p := (q, φ, q, φ, ν), the elements of P . We get the following expression for the part of the kinetic energy of the fluid due to the rigid motion of the bodies:

K f r := 1 2 q • M f r q.
Observe that M f r is a matrix but can also be thought of as a bilinear form in q. To express the overall kinetic energy of the fluid, we introduce:

K f r,d := F ∇ψ r • ∇ψ d dm f , K f r,t := F ∇ψ r • ∇ψ t dm f , K f d := 1 2 F |∇ψ d | 2 dm f , K f d,t := F ∇ψ d • ∇ψ t dm f , K f t := 1 2 F |∇ψ t | 2 dm f ,
what leads to the decomposition:

K f = K f r + K f r,d + K f r,t + K f d + K f d,t + K f t . (6.4)
By analogy with M f r , we introduce as well the following bilinear forms on the phase space P , obtained as the polarization of the different parts of the kinetic energy in (6.4). For all p a := (q, φ, qa , φa , ν a ) and p b := (q, φ, qb , φb , ν b ), we set:

M f r,d ( φb ) qa := i F ∇ψ d ( φb ) • ∇ψ r i dm f qi a , (6.5a) M f r,t (ν b ) qa := i F ∇ψ t (ν b ) • ∇ψ r i dm f qi a , (6.5b) M f d ( φa , φb ) := F ∇ψ d ( φa ) • ∇ψ d ( φb ) dm f , (6.5c) M f d,t ( φa , ν b ) := F ∇ψ d ( φa ) • ∇ψ t (ν b ) dm f , (6.5d) M f t (ν a , ν b ) := F ∇ψ t (ν a ) • ∇ψ t (ν b ) dm f . (6.5e)
In these definitions, we have denoted for instance ψ d ( φ) instead of ψ d to emphasize the dependance of the function with respect to φ. As already mentioned, this dependance is linear (since ∂ n ψ d depends linearly on φ). Rewriting (6.4) by means of the mass matrices we get:

K f = 1 2 q • M f r q + M f r,d ( φ) q + M f r,t (ν) q + 1 2 M f d ( φ, φ) + M f d,t ( φ, ν) + 1 2 M f t (ν, ν). (6.6)
At last, we define the bloc-diagonal matrix:

M b k := diag(m k Id, Ω T k I * k Ω k ), (6.7a)
where Ω k := Ṙk R T k and:

M b := diag(M b 1 , . . . , M b n ). (6.7b)
The kinetic energy of the bodies reads:

K b = 1 2 q • M b q.
(6.7c)

Explicit computation of the Christoffel symbols

The expression of the Lagrangian function is, according to (3.2):

L := K b + K f -P, (6.8)
and the Euler-Lagrange equations of motion are:

d dt ∂L ∂ q - ∂L ∂q = F t , (6.9)
where F t is the generalized force induced by the actuation of the thrusters and computed in Proposition 5.1. Combining (6.8) and (6.9), it comes:

d dt ∂K b ∂ q - ∂K b ∂q + d dt ∂K f ∂ q - ∂K f ∂q = d dt ∂P ∂ q - ∂P ∂q + F t . (6.10)
The game consists now in expanding each term in this equality. It is done in a succession of lemma.

To shorten notation, all of the quantities in what follows are equal to 0 for unspecified indices. Lemma 6.1 (Free rigid motion) The following identity holds:

d dt ∂K b ∂ q - ∂K b ∂q = M b q + Γ b r ( q, q) + Γ b r,d ( φ) q when N = 3, M b q + Γ b r,d ( φ) q when N = 2. (6.11)
• The term Γ b r (called Christoffel symbol in [23, §7.5 pages 196-197]) is a rank-3 tensor defined by:

(Γ b r ) k ij := 1 2 ∂Ω e k ∂θ i - ∂Ω e i ∂θ k • I * l Ω e j + ∂Ω e k ∂θ j - ∂Ω e j ∂θ k • I * l Ω e i + ∂Ω e i ∂θ j + ∂Ω e j ∂θ i • I * l Ω e k ,
if i, j, k ∈ {6(l -1) + 3 + i , 6(l -1) + 3 + j , 6(l -1) + 3 + k } for some l ∈ {1, . . . , n} and i , j , k ∈ {1, 2, 3}.

• The other Christoffel symbol Γ b r,d ( φ) is a rank-2 tensor: -When N = 3:

(Γ b r,d ( φ)) k j := Ω e k • ∂I * l ∂φ l • φl Ω e j ,
if j, k ∈ {6(l -1) + 3 + j , 6(l -1) + 3 + k } for some l ∈ {1, . . . , n} and j , k ∈ {1, 2, 3}.

-When N = 2:

(Γ b r,d ( φ)) k j := ∂I l ∂φ l • φl , if j = k = 3l.
Proof: Starting from the definitions (6.7), straight computations together with the results of Section B yield the expressions of the Christoffel symbols.

We restate Lemma 5.3 under another form: Lemma 6.2 (Potential energy) For all k ∈ {1, . . . , n}:

• When N = 3:

d dt ∂P ∂ q - ∂P ∂q i :=      -g(m k -m k f ) if i = 6(k -1) + 3, g m k f Ω e i • (r k f -r k ) × e 3 if i = 6(k -1) + i , i ∈ {1, 2, 3}. (6.12a) • When N = 2: d dt ∂P ∂ q - ∂P ∂q i := -g(m k -m k f ) if i = 3(k -1) + 2, g m k f (r k f -r k ) ⊥ • e 2 if i = 3k. (6.12b)
Both remaining terms in (6.10) are more involved. We need some additional definitions to pursue. We set on ∂B elementary rigid velocity fields related to the generalized coordinates:

• When N = 3:

w i := e j • n on ∂B k for i = 6(k -1) + j, w i := ω e j × (x -r k ) • n on ∂B k for i = 6(k -1) + 3 + j,
for all k ∈ {1, . . . , n}, j ∈ {1, 2, 3}.

• When N = 2:

w i := e j • n on ∂B k for i = 3(k -1) + j, w i := (x -r k ) ⊥ • n on ∂B k for i = 3k,
for all k ∈ {1, . . . , n}, j ∈ {1, 2}.

We introduce also the velocity field related to the deformations:

w d := R k φk (φ -1 k (g -1 k x)) on ∂B k ,
for all k ∈ {1, . . . , n}. Each w i or w d can be decomposed on ∂B into its normal and tangential components:

w i = w n i n + w τ i and w d = w n d n + w τ d .
We also need the following quantities, homogeneous to accelerations:

• When N = 3:

-For i = 6(k-1)+3+i , j = 6(k-1)+3+j , k ∈ {1, . . . , n}, i , j ∈ {1, 2, 3}:

γ ij := ∂ω e i ∂θ j × (x -r k ) + ω e i × (ω e j × (x -r k )) on ∂B k .
-For i = 6(k -1) + 3 + i , k ∈ {1, . . . , n}, i ∈ {1, 2, 3}:

γ id := ω e i × R k φk (φ -1 k (g -1 k x)) on ∂B k .
• When N = 2:

-For i = j = 3k, k ∈ {1, . . . , n}:

γ ij := -(x -r k ) on ∂B k .
-For i = 3k, k ∈ {1, . . . , n}:

γ id := R k φk (φ -1 k (g -1 k x)) ⊥ on ∂B k .
As usual, the values of these quantities for unspecified indices or/and on unspecified boundaries are zero. Like for the velocities, we denote γ n ij and γ n id the normal components of γ ij and γ id on ∂B.

At last ∇ σ stands for the tangential gradient on ∂B, H for the mean curvature and II for the second fundamental form (see Section C in the Appendix for further details).

To not overload notation, we display the results for the case N = 3 only (the case N = 2 being quite similar). Lemma 6.3 (Euler-Lagrange equation for K f ) Consider a pair of allowable controls (φ(t), ν(t)), assume furthermore that φ(t) is twice continuously differentiable and ν(t) continuously differentiable and denote:

φ := d 2 φ dt 2 and ν := dν dt .
Then:

d dt ∂K f ∂ q - ∂K f ∂q = M f r q + M f r,d ( φ) + M f r,t ( ν) + 1 2 Γ f r ( q, q) + Γ f r,d ( φ) q + Γ f r,t (ν) q + 1 2 Γ f d ( φ, φ) + Γ f d,t ( φ, ν) + 1 2 Γ f t (ν, ν), (6.13)
where M f r,d ( φ) and M f r,t ( ν) are rank-1 tensors defined in (6.5). The Christoffel symbol Γ f r is a rank-3 tensors that reads:

(Γ f r ) k ij := - ∂B (∇ σ ψ r k •∇ σ ψ r i )w n j dσ f - ∂B (∇ σ ψ r k •∇ σ ψ r j )w n i dσ f + ∂B (∇ σ ψ r i •∇ σ ψ r j )w n k dσ f -2 ∂B H∂ n ψ r i ∂ n ψ r j + II(w τ i , w τ j ) -γ n ij + ∇ σ ∂ n ψ r j • w τ i + ∇ σ ∂ n ψ r i • w τ j ψ r k dσ f + ∂B ∂ n ψ r k ∂ n ψ r i ∂ n ψ r j dσ f ∀ i, j, k ∈ {1, . . . , 6n}.
The term Γ f r,d ( φ) is a rank-2 tensor defined by:

(Γ f r,d ( φ)) k i := - ∂B (∇ σ ψ r k •∇ σ ψ r i )w n d dσ f - ∂B (∇ σ ψ r k •∇ σ ψ d )w n i dσ f + ∂B (∇ σ ψ r i •∇ σ ψ d )w n k dσ f -2 ∂B H∂ n ψ r i ∂ n ψ d + II(w τ i , w τ d ) -γ n id + ∇ σ ∂ n ψ d • w τ i + ∇ σ ∂ n ψ r i • w τ d ψ r k dσ f + ∂B ∂ n ψ r k ∂ n ψ r i ∂ n ψ d dσ f ∀ i, k ∈ {1, . . . , 6n}.
Γ f r,t (ν) is also a rank-2 tensor:

(Γ f r,t (ν)) k i := - ∂B (∇ σ ψ r k • ∇ σ ψ t )w n i dσ f + ∂B (∇ σ ψ r i • ∇ σ ψ t )w n k dσ f - ∂B H∂ n ψ r i ∂ n ψ t + ∇ σ ∂ n ψ t • w τ i )ψ r k dσ f + ∂B H∂ n ψ r k ∂ n ψ t + ∇ σ ∂ n ψ t • w τ k )ψ r i dσ f ∀ i, k ∈ {1, . . . , 6n},
Last remaining three terms Γ f d ( φ, φ), Γ f d,t ( φ, ν) and Γ f t (ν, ν) stand for rank-1 tensors:

(Γ f d ( φ, φ)) k := -2 ∂B (∇ σ ψ r k • ∇ σ ψ d )w n d dσ f + ∂B |∇ σ ψ d | 2 w n k dσ f -2 ∂B H|∂ n ψ d | 2 + II(w τ d , w τ d ) + 2∇ σ ∂ n ψ d • w τ d ψ r k dσ f + ∂B ∂ n ψ r k |∂ n ψ d | 2 dσ f ∀ k ∈ {1, . . . , 6n}. (Γ f d,t ( φ, ν)) k := - ∂B (∇ σ ψ r k • ∇ σ ψ t )w n d dσ f + ∂B (∇ σ ψ d • ∇ σ ψ t )w n k dσ f - ∂B H∂ n ψ d ∂ n ψ t + ∇ σ ∂ n ψ t • w τ d )ψ r k dσ f + ∂B H∂ n ψ r k ∂ n ψ t + ∇ σ ∂ n ψ t • w τ k )ψ d dσ f ∀ k ∈ {1, . . . , 6n}, (Γ f t (ν, ν)) k := - ∂B |∇ σ ψ t | 2 w n k dσ f - ∂B H∂ n ψ r k ∂ n ψ t + ∇ σ ∂ n ψ t • w τ k )ψ t dσ f + ∂B ∂ n ψ r k |∂ n ψ t | 2 dσ f ∀ k ∈ {1, . . . , 6n},
Remark 6.1 Observe that it suffices to compute the elementary potentials on the boundaries of the bodies only to evaluate all of the quantities (i.e. mass matrices and Christoffel symbols) of the lemma. It will be useful for numerical simulations.

Proof: Starting from the expression (6.6) of K f , straight computations of the left hand side term in (6.13) lead, after some algebra, to the following formal expressions for the Christoffel symbols:

Γ f r ( q, q) q * := 2 ∂M r ∂q • q ( q, q * ) -∂M r ∂q • q * ( q, q), (

Γ f r,d ( φ)( q, q * ) := ∂M f r ∂φ • φ ( q, q * ) + ∂M f r,d ( φ) ∂q • q q * - ∂M f r,d ( φ) ∂q • q * q, (6.14b) Γ f r,t (ν)( q, q * ) := ∂M f r,t (ν) ∂q • q q * - ∂M f r,t (ν) ∂q • q * q, ( 6.14a)
Γ f d ( φ, φ) q * := 2 ∂M f r,d ( φ) ∂φ • φ q * - ∂M f d ( φ, φ) ∂q • q * , 6.14c)
Γ f d,t ( φ, ν) q * := ∂M f r,t (ν) ∂φ • φ q * - ∂M f d,t ( φ, ν) ∂q • q * , (6.14d)
Γ f t (ν, ν) q * := - ∂M t (ν, ν) ∂q • q * , ∀ q * ∈ R 6 . (6.14e)
We need then to differentiate the mass matrices with respect to q or φ i.e. with respect to perturbations of the domain F. The elements of the mass matrices defined by (6.5) have generically the form:

F ∇ψ a • ∇ψ b dσ f ,
where ψ a and ψ b solve Neumann boundary value problems in F. This problem of shape differentiation is addressed in Section E in the Appendix. We use formula (E.9) together with (6.14) and we get the expressions of the lemma.

It remains only one term to compute in (6.10).

Lemma 6.4 (Generalized force) The generalized force in (6.10) reads:

F t = Γ t r (ν) q + Γ t r,d ( φ, ν) + Γ t r,t (ν, ν), (6.15)
where, for all i, k ∈ {1, . . . , 6n} we have:

(Γ t r (ν)) k i := ∂B (H∂ n ψ r k ∂ n ψ t + ∇ σ ∂ n ψ t • w n k )ψ r i dσ f - ∂B ∂ n ψ t ∂ n ψ r k ∂ n ψ r i dσ f , (Γ t r,d ( φ, ν)) k := ∂B (H∂ n ψ r k ∂ n ψ t + ∇ σ ∂ n ψ t • w n k )ψ d dσ f - ∂B ∂ n ψ t ∂ n ψ r k ∂ n ψ d dσ f , (Γ t r,t (ν, ν)) k := ∂B (H∂ n ψ r k ∂ n ψ t + ∇ σ ∂ n ψ t • w n k )ψ t dσ f - ∂B |∂ n ψ t | 2 ∂ n ψ r k dσ f .
Proof: Use Proposition 5.5 and proceed like in the preceding lemma.

Plugging (6.11), (6.12), (6.13) and (6.15) into (6.10), we obtain a second order ODE in q. Thus, we have realized our aim: the ODE yields the net rigid motion of the bodies with respect to the controls φ and ν (and their derivatives φ, φ and ν).

Numerical simulations

All the simulations displayed in this section have been performed with the Biohydrodynamics Matlab Toolbox (BhT) available at

• http://bht.gforge.inria.fr/ (hosted by INRIA Gforge)

• http://www.mathworks.com/matlabcentral/fileexchange/21872

BhT is free (distributed under licence GNU GPL) and comes with complete html documentation, including examples and tutorials. It allows to realize numerical experiments involving a set of rigid bodies or articulated rigid solids in a perfect fluid. The fluid domain can be partially or totally bounded and the buoyant force and collisions between bodies are supported. The equations of motion used in BhT are those obtained in Section 6.3 for the 2d case.

As already mentioned, to solve these equations it is required only to compute the elementary potential functions on the boundaries of the bodies. For being solutions of Neumann boundary value problems, these potential functions can be computed by means of boundary integral formula, that is, only in terms of the boundary data. As a consequence, from a numerical point of view, the problem can be solved by meshing the fluid's boundary only.

Distant interaction between solids in a fluid

In this subsection we wish to illustrate how submerged bodies interact in a fluid. In particular, we will show that when the bodies are close, the common hydrodynamically decoupled dynamics assumption is not relevant.

Kissing ellipses

Consider two ellipses neutrally buoyant and set at time t = 0 as shown in figure 4(a). We choose horizontal translational initial velocities directed to the left for the upper ellipse and to the right for the other one and null initial rotational velocity for both of them. We let then the system fluid-ellipses evolve freely. Under the hydrodynamically decoupled dynamics assumption, according to d'Alembert's paradox, both ellipses should have a straight uniform motion. Actually, they kiss and half-rotate following the trajectories pictured in figure 5(a).

Aspiration phenomenon

This example illustrate an aspiration phenomenon. Two ellipses are set at time t = 0 as in figure 6(a). The upper one is neutrally buoyant while the density of the other is slightly greater than the fluid's. The system fluid-ellipses is at rest at time t = 0. The lower ellipse falls (see figures 7.1) and drags the other down. One more time, this phenomenon would not occur under the hydrodynamically decoupled dynamics assumption.

Articulated solid bodies

Articulated rigid solids are deformable bodies made up of rigid solids linked together by hinges. It is a particular simple case of deformable bodies. In the following example, we consider two articulated fishes swimming one toward the other.

Remind that the deformations and more specifically for articulated bodies, the relative positions of the solids, are given as functions of time. BhT provides tools that allow to compute the power and the torques at the joints. Torques are actually the relevant quantities we are interested in when studying fish locomotion. We refer to the on-line documentation of BhT for details.

Concluding comments

Lagrangian formalism consists in computing first the energy of the fluid flow: it is obtained by integrating the gradient of the potential (or stream) function over the fluid domain. This is a smoothing process for the resulting kinetic energy is more regular with respect to the generalized coordinates than the potential function. Then, according to the least action principle, Euler-Lagrange equations are obtained partly by differentiating the fluid's kinetic energy. Notice that this process is carried out without difficulties in the example treated in Section 5.3.

Conversely, with Newtonian formalism, forces are computed first, what requires to differentiate the potential (or the stream) function and yields, in the nonsmooth case of Section 5.3 an infinite and nonintegrable quantity. Next these forces are integrated (when it is allowable) along the boundary of the bodies in order to obtain the hydrodynamic resultant force and next to apply Newton's laws.

We have proved that both approaches coincide for smooth bodies. Although Lagrangian formalism is more simple to handle, it does not allow to obtain the expected governing equations of motion when the bodies are endowed with thrusters.

At last, regarding the problem of existence and uniqueness of solutions for the self-propelled motion problem, it is precisely the topic the article [START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF] deals with.

A Shape space

In this section, all the quantities are described in the frame (e i * j ) attached to the i-th body. Let B i ⊂ R N be a large open ball centered at the origin of the frame (in particular, let choose B i large enough to contain B 0 i ; see figure 1) and denote C m 0 (B i , R N ) (1 ≤ m ≤ +∞) the vector space of all the m-times continuously differentiable functions F := (F 1 , . . . , F N ) from R N into R N that vanish outside Bi . This space, endowed with the norm:

F m,∞ := max γ 1 +...+γ N ≤m γ 1 ,...,γ N ∈N max 1≤i≤N ∂ γ1+•••+γ N F i ∂ γ1 x 1 . . . ∂ γ N x N L ∞ (B)
, is a Banach space. Figure 9: Angles of the joints with respect to time (given functions). Joint 1 correspond the fish's head and joint 3 to the tail.

We denote U i the open subset of this space consisting of all the functions such that F 1,∞ < 1. Further, U i := Id + U i is the open subset of the affine space

Id + C m 0 (B i , R N ) := {Id + F : F ∈ C m 0 (B i , R N )},
consisting of the functions Id + F such that F ∈ U i . A fixed point argument ensures that for any such function, (

Id + F )| Bi is a C m diffeomorphism from B i onto B i .
The set U i is a manifold, a chart being (U i , U i , Λ i ) where

Λ i : F ∈ U i → φ := Id + F ∈ U i .
For any φ ∈ U i , the tangent space T φ U i is equal to the Banach space C m 0 (B i , R N ). Let us define now D i , the set of all the C m diffeomorphismes φ of R N such that φ( Bi ) = Bi and φ| R N \Bi = Id. The set D i is a manifold as well since for any φ 0 ∈ D i , a compatible system of local charts is (U i , φ 0 (U i ), φ 0 (Λ i )). We denote θ := (θ 1 , θ 2 , θ 3 ) and ω e i the vector of R 3 such that:

∂R ∂θ i (θ)R(θ) T x = ω e i × x ∀ x ∈ R 3 ,
and Ω e i such that:

R(θ) T ∂R ∂θ i (θ)x = Ω e i × x ∀ x ∈ R 3 .
One easily check that ω e 1 = e 1 , ω e 2 = exp(θ 1 S 1 )e 2 and ω e 3 = exp(θ 1 S 1 ) exp(θ 2 S 2 )e 3 . Since we have R(θ)Ω e i = ω e i and exp(θ i S i )e i = 0, we deduce that Ω e 1 = exp(-θ 3 S 3 ) exp(-θ 2 S 2 )e 1 , Ω e 2 = exp(-θ 3 S 3 )e 2 and Ω e 3 = e 3 . These identities entail:

∂Ω e i ∂θ j - ∂Ω e j ∂θ i = Ω e i × Ω e j . (B.1)
Let θ(t) be a smooth given function and denote θ = ( θ1 , θ1 , θ1 ) the first derivative of θ. Then if Ω is defined by the relation:

R(θ) T d dt R(θ) x = Ω × x ∀ x ∈ R 3 ,
it reads also:

Ω = i θi Ω e i ,
and we have:

Lemma B.1 For all i = 1, 2, 3, the following identity holds:

d dt ∂Ω ∂ θi - ∂Ω ∂θ i = Ω e i × Ω. (B.2) Proof: We have ∂Ω ∂ θi = Ω e i ,
and then:

d dt ∂Ω ∂ θi = j ∂Ω e i ∂θ j θj .
On the other hand:

∂Ω ∂θ i = j ∂Ω e j ∂θ i θj ,
and the conclusion arises with (B.1).

Let now J e be a symmetric, time-dependent matrix and define J := RJ e R T , ω := RΩ and K := Ω • J e Ω/2. Then:

Lemma B.2 For any i = 1, 2, 3, the following identity holds:

d dt (Jω) • ω e i = d dt ∂K ∂ θi - ∂K ∂θ i .
Proof: We have:

d dt ∂K ∂ θi - ∂K ∂θ i = d dt (J e Ω) • Ω e i + (J e Ω) • d dt ∂Ω ∂ θi - ∂Ω ∂θ i ,
and, with (B.2):

d dt ∂K ∂ θi - ∂K ∂θ i = d dt (J e Ω) • Ω e i + (J e Ω) • (Ω e i × Ω).
On the other hand:

d dt (Jω) • ω e i = R d dt (J e Ω) • RΩ e i + ṘJ e Ω • RΩ e i , = d dt (J e Ω) • Ω e i + (J e Ω) • (Ω e i × Ω),
and the proof is completed.

C Second fundamental form

Let us denote S N -1 the unitary sphere of R N . The second fundamental form of ∂F, denoted by II, is defined from the differential Dn of the function

n : ∂F → S N -1 ,
by the relation:

Dn(x) • u, w = -II(x), u, w , ∀ (x, u), (x, v) ∈ T ∂F.
The second fundamental form is bilinear and symmetric, so when N = 3 it admits two real eigenvalues κ 1 , κ 2 called the principal curvatures of ∂F. The associated eigenvectors τ i (i = 1, 2) are the principal directions of ∂F (see [6, page 129] for a more precise definition). Therefore, the following relation holds:

dn • τ i , τ i = -II, τ i , τ i (C.1a) = -κ i , ∀ i = 1, 2. (C.1b)
The signs of κ 1 and κ 2 are chosen such that {n, τ 1 , τ 2 } be a right-handed orthogonal basis of R 3 . We denote by H := κ 1 + κ 2 the mean curvature of ∂F. When N = 2, H coincides with κ, the curvature of F(t) and relation (C.1) turns into:

Dn • τ = -Hτ ,
where τ is a unitary tangent vector to the curve ∂F(t). As for N = 3, the sign of H is chosen such that {n, τ } be a direct orthogonal basis of R 2 .

D Neumann boundary value problem

In this section we recall some results about the well-posedness of Neumann boundary value problems. These results can be found in [START_REF] Amrouche | Weighted Sobolev spaces for Laplace's equation in R n[END_REF], [START_REF] Amrouche | Dirichlet and Neumann exterior problems for the n-dimensional Laplace operator: an approach in weighted Sobolev spaces[END_REF] in the generically form of which is:

-∆ψ = 0 in F, (D.1a) ∂ n ψ = g on ∂F. (D.1b)
We assume that ∂F is Lipschitz continuous and g ∈ H -1/2 (∂F) is given (H -1/2 (∂F) is the dual space of the Sobolev space H 1/2 (∂F); see [21, §7.3 pages 38-42]).

D.1 Bounded domain

In this case, the compatibility condition:

g, 1 H -1/2 (∂F )×H 1/2 (∂F ) = 0, (D.2)
is required for problem (D.1) to be well posed. We introduce the function space:

X 1 (F) := {ψ ∈ H 1 (F) : ∆ψ ∈ L 2 (F)}.
For any ψ ∈ X 1 (F), its normal trace ∂ n ψ exists as an element of H -1/2 (∂F). Lax-Milgram Theorem ensures that there exists a unique weak solution ψ ∈ X 1 (F) satisfying:

F ∇ψ • ∇ϕ dx = ∂F gϕ dσ, ∀ ϕ ∈ H 1 (F), (D.3a) ∂B ψ dσ = 0. (D.3b)
Moreover, we have an estimate:

ψ X 1 (F ) ≤ C g H -1/2 (∂F ) , (D.4)
for some constant C > 0 independent of ψ and g. Assume now that ∂F is C 1,1 (continuously differentiable with first derivative Lipschitz contiuous) and g ∈ H 1/2 (∂F), then the solution ψ of (D.3) is in:

X 2 (F) := {ψ ∈ H 2 (F) : ∆(∂ xi ψ) ∈ L 2 (F), ∀ i = 1, . . . , N }.
For such function, both quantities ∂ 2 n ψ := n•D 2 ψ n and ∆ σ ψ (the Laplace-Beltrami operator) are well defined on ∂F as elements of H -1/2 (∂F). We have again an estimate:

ψ X 2 (F ) ≤ C g H 1/2 (∂F ) ,
for some constant C > 0 independent of ψ and g.

D.2 Unbounded domain

When F is not bounded, the compatibility condition (D.2) is required only when N = 2. We define the weight functions:

ρ(x) := 1 + |x| 2 log(1 + |x| 2 ) if N = 2, (D.5) ρ(x) := 1 + |x| 2 if N = 3, (D.6)
and the Lebesgue space:

L 2 N (F) := {ψ ∈ D (R 3 ) : ρψ ∈ L 2 (R N )},
where D (R 3 ) is the space of distributions. Classical Sobolev spaces must be replaced by Sobolev spaces:

H 1 N (F) : = {ψ ∈ L 2 N (F) : ∂ xi ψ ∈ L 2 (F), ∀ i = 1, . . . , N }, X 1 N (F) : = {ψ ∈ H 1 N (F) : ρ -1 ∆ψ ∈ L 2 (F)}, H 2 N (F) : = {ψ ∈ H 1 N (F) : ρ -1 ∂ 2 xixj ψ ∈ L 2 (F), ∀ i, j = 1, . . . , N }, X 2 N (F) : = {ψ ∈ H 2 N (F) : ρ -2 ∆ψ ∈ L 2 (F)}.
When ∂F is Lipschitz continuous and g ∈ H -1/2 (∂F):

• Then the solution of (D.3a), replacing H 1 (F) by H 1 N (F), exists and is unique in the case N = 3. In particular, condition (D.3b) is not needed in this case to ensure the uniqueness.

• When N = 2, solutions of (D.3a) exist in H 1 N (F) but (D.3b) is required to ensure the uniqueness.

It can be proved that the solutions are in fact in X 1 N (F). As in the case F bounded, if ∂F is C 1,1 and g ∈ H 1/2 (∂F), then the solutions are in X 2 N (F).

E Shape sensitivity analysis

This section is self-contained and independent, including for the notation. We recall results of shape optimization theory. Most of them have been proved in [START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF]. Let F 0 ⊂ R N and define for all s = (s 1 , s 2 , s 3 ) ∈ R 3 diffeomorphisms φ s of class C m (m ≥ 1) in R N and equal to the identity outside a large fixed ball (like in section A). We denote D m the Banach space of all such diffeomorphisms.

Denote

F := φ s (F 0 ), w i := ∂ si φ s • φ -1 s , γ ij := ∂ 2 sisj φ s • φ -1
s and w n i and w τ i the normal and tangential components of w i on ∂F such that:

w i = w n i n + w τ i .
Define likewise γ n ij , the normal component of γ ij .

For i ∈ {1, 2, 3}, let then ψ i be the solution of the Neumann boundary value problem:

-∆ψ i = 0 in F,

∂ n ψ i = w n i on ∂F.
Introduce also for all s ∈ R 3 , ν 1 , ν 2 ∈ H -1/2 (∂F) such that ν i , 1 H -1/2 (∂F )×H 1/2 (∂F ) = 0, when N = 2 or when N = 3 and F bounded. Let then ψ 4 and ψ 5 be defined by:

-∆ψ i = 0 in F,

∂ n ψ i = ν i-3 on ∂F,
for i ∈ {4, 5}. At last, we introduce the functionals, for all s := (s 1 , s 2 , s 3 ) and all i, j ∈ {1, 2, 3, 4, 5}:

J ij (s) := F ∇ψ i • ∇ψ j dx.
The proof of the following lemma can be found in [START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF]. Then J ij is of class C k-1 .

We can compute the derivation of J ij with respect to s, assuming additional regularity on F 0 , φ s and ν l :

Lemma E.2 Assume that:

• F 0 is of class C 1,1 ,
• the mapping s → φ s ∈ D m is of class C 2 with m ≥ 2,

• the mappings s → ν l • φ -1 s ∈ H 1/2 (∂F 0 ) are of class C 1 for l = 1, 2.

Then, for all i = 1, 2, 3 and j = 1, . . . , 5, ∂ si ψ j is well defined in H 1 (F). We prove formula (E.1a). Let us fix s and j, specify i = 1 and denote, for all h ∈ R, φ h := φ (s1+h,s2,s3) (φ -1 s ), F h := φ h (F) = φ (s1+h,s2,s3) (F 0 ) and ψ h the solution of the Neumann boundary value problem:

-∆ψ h = 0 on F h , (E.2a)

∂ n ψ h = b h on ∂F h , (E.2b)
where b h := w n h with w h := ∂ sj φ h (φ -1 h ). When h = 0 we have ψ h = ψ j . We introduce B h := b h (φ h ) and N h := n(φ h ), both being defined on ∂F h . We denote merely ψ , B and N the derivatives of ψ h , B h and N h with respect to h at h = 0. On the one hand, we obtain merely:

-∆ψ = 0 in F.

The derivation of the boundary conditions are more involved. Differentiating with respect to h the identity:

∇ψ h (φ h ) • N h = B h on ∂F,
when h = 0, we obtain:

∂ n ψ + n • D 2 ψ j w i + ∇ψ j • N = B on ∂F,
where D 2 ψ j is the Hessian matrix of ψ j . Denoting ∂ 2 n ψ j := n • D 2 ψ j n and since w i = w n i n + w τ i , we obtain that:

∂ n ψ + ∂ 2 n ψ j w n i + n • D 2 ψ j w τ i + ∇ψ j • N = B on ∂F. (E.3)
We use on the one hand, the identity of [12, Proposition 5.4.12 page 196], valid for any function ψ regular enough in a neighborhood of ∂F:

∆ψ = ∆ σ ψ -H∂ n ψ + ∂ 2 n ψ on ∂F,
where ∆ σ is the Laplace-Beltrami operator and H the mean curvature. On the other hand, computing the tangential derivative in the direction w τ i of the boundary condition of System (E.2), we get:

n • D 2 ψ j w τ i -II(∇ σ ψ j , w τ i ) = ∇ σ v n i • w τ i on ∂F. (E.4)
Proposition 4 in Section 5 of [START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF] yields:

∇ψ j • N = -∇ σ ψ j • ∇ σ w n i -II(∇ σ ψ j , w τ i ) on ∂F. (E.5)
We get hence, plugging (E.4) and (E.5) into (E.3):

∂ n ψ = ∆ σ ψ j w n i -Hw n i w n j + ∇ σ ψ j • ∇ σ w n i + B -∇ σ w n j • w τ i .

It remains to compute the derivative of B h with respect to h. Taking into account its definition, we have:

B = ∂ sj φ h (φ -1 h ) • N h | h=0 . (E.6)
Invoking again formula (E.5), we obtain:

B = γ ij -II(w τ i , w τ j ) -∇ σ w n i • w τ j , (E.7)
which leads to formula (E.1a). We proceed likewise for formula (E.1b).

Reynolds's formula reads: When i, j ∈ {1, 2, 3} and k ∈ {4, 5}, this equality turns into:

∂
∂ ∂s i F ∇ψ j • ∇ψ k dx = - ∂F (∇ψ j • ∇ψ k )w n i dσ - ∂F (H∂ n ψ i ∂ n ψ j + II(w τ i , w τ j ) -γ n ij + ∇ σ ∂ n ψ i • w τ j + ∇ σ ∂ n ψ j • w τ i )ψ k dσ - ∂F (H∂ n ψ i ∂ n ψ k + ∇ σ ∂ n ψ k • w τ i )ψ j dσ + 2 ∂F ∂ n ψ i ∂ n ψ j ∂ n ψ k dσ.
At last, for any i ∈ {1, 2, 3} and j, k ∈ {4, 5}, we have:

∂ ∂s i F ∇ψ j • ∇ψ k dx = - ∂F (∇ψ j • ∇ψ k )w n i dσ - ∂F (H∂ n ψ i ∂ n ψ j + ∇ σ ∂ n ψ j • w τ i )ψ k dσ - ∂F (H∂ n ψ i ∂ n ψ k + ∇ σ ∂ n ψ k • w τ i )ψ j dσ + 2 ∂F ∂ n ψ i ∂ n ψ j ∂ n ψ k dσ.
Combining results of this lemma and after some tedious algebra, we get:

1 . 2 1 . 3

 1213 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Newtonian versus Lagrangian formalims . . . . . . . . . . . . . . . . 3 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.5 Organization of the paper . . . . . . . . . . . . . . . . . . . . . . . . 5 1.6 A particular case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

  )| -2 (-e iθ f (e iθ ))dθ.(5.12) Let us then compute the behavior of |∇φ 2 (e iθ )| 2 |f (e iθ )| -2 (-e iθ f (e iθ )) nearby θ = 0. By virtue of (5.8), we get:

Figure 4 :

 4 Figure 4: Screenshots of the motion of the ellipses All of the movies and data (densities of the solids, of the fluid, initial positions and velocities...) in this section are available at http://bht.gforge.inria.fr/Examples/demos.html.

  Trajectories of the centers of mass of the ellipses.

  Orientations (radians) of the ellipses with respect to time.

Figure 5 :

 5 Figure 5: Kissing ellipses: trajectories of the center of masses and orientations.

Figure 6 :

 6 Figure 6: The lower ellipse falls and drags the other down.

  y-coordinates of the ellipses with respect to time (dashed line corresponds to the lower ellipse).

  Vertical velocities of the ellipses with respect to time (dashed line corresponds to the lower ellipse).

Figure 7 :

 7 Figure 7: The aspiration phenomenon becomes clear with the graphs of the trajectories of the center of masses and the one of the vertical velocities.

(a) t = 12 (

 12 b) t = 35 (c) t = 64

Figure 8 :

 8 Figure 8: Screenshots of the swimming fishes

Figure 10 :Figure 11 :

 1011 Figure 10: Power expanded at the joints

Lemma E. 1

 1 Assume that• F 0 is Lipschitz continuous, • the mapping s → φ s ∈ D m is of class C k with k ≥ 2 and m ≥ 1, • the mappings s → ν l • φ -1 s ∈ H -1/2 (∂F 0 ) are of class C k-1 for l = 1, 2.

When j = 1 , 2 , 3

 123 it solves the Neumann boundary value problem -∆(∂ si ψ j ) = 0 in F, and on ∂B:∂ n (∂ si ψ j ) = ∆ σ ψ j w n i -H∂ n ψ i ∂ n ψ j -II(w τ i , w τ j ) + ∇ σ ψ j • ∇ σ w n i + γ n ij -∇ σ ∂ n ψ i • w τ j -∇ σ ∂ n ψ j • w τ i . (E.1a)When j = 4, 5, the boundary condition on ∂B turns out to be:∂ n (∂ si ψ j ) = ∆ σ ψ j w n i -H∂ n ψ i ∂ n ψ j + ∇ σ ψ j • ∇ σ w n i -∇ σ ∂ n ψ j • w τ i . (E.1b)Proof: The problem of differentiating with respect to the shape of the domain the solution of a Neumann boundary value problem is treated in [12, Subsection 5.5 pages 202-206], in[START_REF] Houot | On the motion and collisions of rigid bodies in an ideal fluid[END_REF] Section 4] and in [27,Section 4]. We give here only the formal computations leading to formula (E.1).

Proposition E. 1 3 ∂F ψ 2 ∂ n ψ 1 dσ + ∂ ∂s 2 ∂F ψ 3 ∂ n ψ 1 dσ - ∂ ∂s 1 ∂F ψ 2 ∂(∇ σ ψ 2 • ∇ σ ψ 3 )w n 1 dσ - 2 ∂FH∂ n ψ 2 ∂ n ψ 3 + II(w τ 2 , w τ 3 ) -γ n 23 + 2 ∂F ψ 4 ∂ 1 ∂F ψ 2 ∂ 1 ∂F ψ 4 ∂

 132231222232323241214 The following identities holds:∂ ∂s n ψ 3 dσ = -∂F (∇ σ ψ 1 • ∇ σ ψ 2 )w n 3 dσ -∂F (∇ σ ψ 1 • ∇ σ ψ 3 )w n 2 dσ + ∂F ∇ σ ∂ n ψ 3 • w τ 2 + ∇ σ ∂ n ψ 2 • w τ 3 ψ 1 dσ + ∂F ∂ n ψ 1 ∂ n ψ 2 ∂ n ψ 3 dσ. (E.9a) ∂ ∂s n ψ 1 dσ -∂ ∂s n ψ 4 dσ = -∂F (∇ σ ψ 1 • ∇ σ ψ 4 )w n 2 dσ + ∂F (∇ σ ψ 2 • ∇ σ ψ 4 )w n 1 dσ -∂F H∂ n ψ 2 ∂ n ψ 4 + ∇ σ ∂ n ψ 4 • w τ 2 )ψ 1 dσ + ∂F H∂ n ψ 1 ∂ n ψ 4 + ∇ σ ∂ n ψ 4 • w τ 1 )ψ 2 dσ, (E.9b) ∂ ∂s n ψ 5 dσ = -∂F (∇ σ ψ 4 • ∇ σ ψ 5 )w n 1 dσ -∂F H∂ n ψ 1 ∂ n ψ 4 + ∇ σ ∂ n ψ 4 • w τ 1 )ψ 5 dσ -∂F H∂ n ψ 1 ∂ n ψ 5 + ∇ σ ∂ n ψ 5 •w τ 1 )ψ 4 dσ + ∂F ∂ n ψ 1 ∂ n ψ 4 ∂ n ψ 5 dσ. (E.9c)

  si ψ j ) • ∇ψ k dx + F ∇ψ j • ∇(∂ si ψ k ) dx. (E.8)It is classical for smooth functions (see for instance[10, formula (85) page 46]) and can be found, in the context of Sobolev spaces in [12, Theorem 5.2.2 page 172]. More precisely, we apply here [12, Corollary 5.2.2 page 172]. Using, for any i ∈ {1, 2, 3} and j, k ∈ {1, . . . , 5}, the identities:ψ j • ∇ σ w n i )ψ k + (∇ σ ψ j • ∇ σ )ψ k w n i dσ,and∇ψ i • ∇ψ j = ∇ σ ψ i • ∇ σ ψ j + ∂ n ψ i ∂ n ψ j on ∂F,straight but tedious computations yield:Lemma E.3 For any i, j, k ∈ {1, 2, 3}, we have:

	∂s i	F	∇ψ j • ∇ψ k dx =	∂F (∇ψ j • ∇ψ k )w n i dσ +
	∂
	∂s i

F ∇(∂ ∂F ∆ σ ψ j w n i ψ k dσ = -∂F (∇ σ F ∇ψ j • ∇ψ k dx = -∂F (∇ψ j • ∇ψ k )w n i dσ -∂F (H∂ n ψ i ∂ n ψ j + II(w τ i , w τ j ) -γ n ij + ∇ σ ∂ n ψ i • w τ j + ∇ σ ∂ n ψ j • w τ i )ψ k dσ -∂F (H∂ n ψ i ∂ n ψ k + II(w τ i , w τ k ) -γ n ik + ∇ σ ∂ n ψ i • w τ k + ∇ σ ∂ n ψ k • w τ i )ψ j dσ + 2 ∂F ∂ n ψ i ∂ n ψ j ∂ n ψ k dσ.

Author partially supported by ECOS C07-E05