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

Bruno Pinçon
email: bruno.pincon@iecn.u-nancy.fr

Locomotion of articulated bodies in an ideal fluid: 2d model with buoyancy, circulation and collisions

Articulated solid bodies (ASB) is a basic model for the study of shapechanging underwater vehicles made of rigid parts linked together by pivoting joints. In this paper we study the locomotion of such swimming mechanisms in an ideal fluid. Our study ranges over a wide class of problems: several ASBs can be involved (without being hydrodynamically decoupled), the fluid-bodies system can be partially or totally confined and fluid circulation, buoyancy force and possible collisions between bodies are taken into account. We derive the Euler-Lagrange equation governing the dynamics of the system, study its well-posedness and describe a numerical scheme implemented in a Matlab toolbox (Biohydrodynamics Toolbox).

Introduction

In the last decade, much work has been done by mathematicians to better understand the dynamics of swimming in a fluid. This interest has grown from the observation that fish and aquatic mammals evolved swimming capabilities far superior to what has been achieved by human technology and consequently provide an attractive model for the design of biomimetic robots. Such swimming mechanisms propelled and steered by shape-changes would be more efficient, stealthier and more maneuverable than if propeller-driven.

So-called autonomous underwater vehicles (AUV) or unmanned undersea vehicles (UUV) have been used extensively to carry out varied missions: for instance, by the oil industry for the maintenance of off-shore installations, by researchers for oceanographic data collection or by the military for mine hunting. The need to improve AUV performance to meet the demand of increasingly more challenging missions has led to intensive research effort in the exploration of biological principles that can be adapted for underwater vehicle engineering applications. The biologically-inspired methods have been envisioned to improve AUVs' low speed maneuvering capabilities including hovering, smallradius turning, sinking and precision station keeping all of which are natural capabilities of aquatic animals. In the area of nano-technologies, the design of nano-robots able to perform basic tasks is a challenge for the forthcoming years.

Significant contributions to the understanding of the biomechanics of swimming have been made by Lighthill [START_REF] Lighthill | Mathematical biofluiddynamics[END_REF], Taylor [START_REF] Taylor | Analysis of the swimming of microscopic organisms[END_REF][START_REF] Taylor | Analysis of the swimming of long and narrow animals[END_REF] Childress [START_REF] Childress | Mechanics of swimming and flying[END_REF] and Wu [START_REF] Wu | Mathematical biofluiddynamics and mechanophysiology of fish locomotion[END_REF]. An interesting survey on the general theme of fish locomotion written by Sparenberg [START_REF] Sparenberg | Survey of the mathematical theory of fish locomotion[END_REF] is worth being mentioned as well.

Experiments have shown that the vortices generated by the tail fins of fish play a crucial role in their locomotion. Indeed, in nature, fish do interact with these vortices. Exploiting the circulation in the flow allows them to reduce the locomotory cost. Therefore, understanding how fish behave in the presence of vortices is essential in studying aquatic propulsion and some models incorporate artificially produced 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][START_REF] Kanso | Swimming due to transverse shape deformations[END_REF]. If we do not neglect the viscosity effects, the relevant model incorporates the non-stationary Navier-Stokes equations for the fluid coupled with Newton's laws for the fish-like swimming object. This perspective is adopted by Carling, Williams and Bowtell 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], Liu and Kawachi in [START_REF] Liu | A numerical study of undulatory swimming[END_REF], Galdi in [START_REF] Galdi | On the steady self-propelled motion of a body in a viscous incompressible fluid[END_REF] or San Martín, Scheid, Takahashi and Tucsnak in [START_REF] San Martin | An initial and boundary problem modeling fish-like swimming[END_REF]. However, and contrary to some common beliefs, forces and momenta acting on the fish body by shed vortices are not solely responsible for the net locomotion and among numerous mathematical articles studying fish locomotion, most of them address the case of a potential flow which is by definition vortex-free: let us mention here the works of Kelly and Murray [START_REF] Kelly | Modelling efficient pisciform swimming for control[END_REF], Kozlov and Onishchenko [START_REF] Kozlov | Motion of a body with undeformable shell and variable mass geometry in an unbounded perfect fluid[END_REF], Kanso, Marsden, Rowley and Melli-Huber [START_REF] Kanso | Locomotion of articulated bodies in a perfect fluid[END_REF], Melli, Rowley and Rufat [START_REF] Melli | Motion planning for an articulated body in a perfect planar fluid[END_REF] and Munnier [START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF][START_REF] Munnier | Locomotion of deformable bodies in an ideal fluid: Newtonian versus lagrangian formalism[END_REF]. This is also the point of view we have chosen in this article. Our model incorporate circulation which is an improvement with respect to some prior vortex free models. However, this circulation is purely due to initial conditions and not continuously generated from the body as in [START_REF] Kanso | Swimming due to transverse shape deformations[END_REF] for instance. To sum up, our model is not valid as representative of real fish swimming but can be seen as a step toward more realistic models and a contribution to the development of mathematical theories involved in fluid-structure interactions.

Thus, in this paper we are interested in studying the dynamics of a set of articulated solid bodies (ASB) immersed in an ideal fluid. The relative positions of the rigid solids composing each ASB are prescribed as functions of time. These functions play the role of controls. They describe thoroughly the shape-changes by which the bodies propel and steer themselves. The study of the locomotion of more general shape-changing swimming bodies is addressed in [START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF].

Since the works of Thomson, Tait and Kirchhoff (summarized in the book of Lamb [START_REF] Lamb | Hydrodynamics. Cambridge Mathematical Library[END_REF]), the standard method to determine the Euler-Lagrange equation driving the motion of rigid solids in an ideal fluid, has consisted in invoking the least action principle of Lagrangian Mechanics for the fluid-bodies system. It has been proved recently in [START_REF] Munnier | Locomotion of deformable bodies in an ideal fluid: Newtonian versus lagrangian formalism[END_REF] that an equivalent method (i.e. leading to the same Euler-Lagrange equation) consists in invoking Newton's laws for the set of solids, coupled with Euler's equations for the motion of the fluid (providing that the flow is irrotational at the initial time). This result is important because it proves first that by considering a potential flow, there is no physical simplification but the assumption that the flow is vortex-free at the time t = 0, and second that the creation of vortices by the swimming bodies is not allowed in a perfect fluid (without incorporating vortex shedding mechanisms). In the present article, we will use the standard method based on Lagrangian Mechanics.

One of the main ingredient in the modeling of fluid-structure problems is the notion of mass matrix arising in the expression of the Lagrangian. Using classical notation we denote by q and q the generalized coordinates and velocities; they parameterize all of the degrees of freedom of the system and allow us to write the Lagrangian of the fluid-bodies system in the short form:

L := 1 2
q M(q) q -P (q), where M(q) is precisely the mass matrix (its elements depend on q and are homogeneous to masses) and P (q) stands for the overall potential energy of the system. The Euler-Lagrange equation driving the free evolution of q with respect to the time is an ODE (ordinary differential equation) resulting from the application to L of the least action principle:

d dt ∂L ∂
q -∂L ∂q = 0, (t ≥ 0).

(1.1) Some degrees of freedom can be prescribed as functions of the others and of time. This idea writes: q = F (t, p),

where p stands for the remaining degrees of freedom of the system and F is a given function. Such constraints on the mechanical system are termed holonomic and well suited to describe the kinematics of swimming robots made of rigid articulated parts. Considering now the Lagrangian as a function of t, p and ṗ, the equation of motion of the ASB whose shape-changes are driven by relation (1.2) is:

d dt ∂L ∂ ṗ - ∂L ∂p = 0, (t ≥ 0). (1.3)
Next, the game consists in expanding the left hand side term of this equality and turn the equation into a convenient form allowing to study the locomotion induced by the shape-changes. The tricky point is the computation of the derivative of the mass matrix with respect to p (or more generally with respect to q) because this task actually requires to differentiate the potential of the fluid with respect to the positions of the solids (parameterized by q). The potential of the fluid being defined as the solution of a Neumann boundary value problem (NBVP) set in the fluid domain (which depends on the positions of the solids), we have to compute so-called shape derivatives.

Nevertheless this difficulty can be avoided when the mass matrix does not depend on p. That means that all of the positions and directions in the fluid are equivalent as seen by an observer attached to the bodies. This isotropy (or symmetry) is ensured providing that: (H1) There is only one body in the fluid.

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

(H3) The body is neutrally buoyant.

(H4) The bodies-fluid system is at rest at the time t = 0. Hypothesis (H1) can equivalently be replaced 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. This approximation makes sense when the immersed bodies are far apart but is not relevant any longer when they are close. Most of the studies on swimming shape-changing bodies in an ideal fluid are carried out under hypotheses (H1-4) or (H1 -4).

One of our aims in this paper is to study the dynamics of a set of ASBs relaxing all of these hypotheses. Thus, in Section 2, we address the case of n rigid solids evolving freely in a fluid with potential flow, the fluid-solids system being either unbounded, either partially or totally confined and we take into account the buoyancy force. We compute the Lagrangian, expand the Euler-Lagrange equation (1.1) and obtain a second order ODE in q that writes: M(q)q + Γ(q), q, q + ∂P ∂q (q) = 0, (t ≥ 0), (1.4) where M(q) is the mass matrix and Γ(q) is a rank-3 tensor, usually called Christoffel symbol in Lagrangian Mechanics. If we denote by M ij (q) the elements of the mass matrix we have:

Γ k ij (q) := 1 2 M ki ∂q j (q) + M kj ∂q i (q) - M ij ∂q k (q) .
We next compute explicitly the elements Γ k ij (q), this task requiring formula of shape sensitivity analysis. The resulting expressions, well suited for numerical simulations, is the first novelty brought in this paper.

In Section 3, we consider a flow with non zero circulation. Circulation is considered as extra degrees of freedom of the fluid flow. We show that the ODE governing the circulation can be decoupled from that driving the motion of the solids and explicitly solved. Therefore, the generalized coordinates relating to the circulation can be canceled out in the Euler-Lagrange equation which finally keeps the same general form as (1.4): M(q)q + Γ(q), q, q + ∂ P ∂q (q) = 0, (t ≥ 0), (1.5) but where the potential energy contains an extra term due to the non zero fluid circulation. The computation of this term requires also the use of shape derivative formula. To our knowledge, such a model involving circulation is new as well and is the second novelty of this paper. We are next concerned with the well posedness of ODE (1.5). This is proved in Section 4, assuming that the boundaries of all of the solids involved are smooth enough. Further, we show that the solution is analytic and can be continued indefinitely unless a collision between two bodies or between a body with a fixed boundary of the fluid domain occurs. Such a result has already be obtained in a more general case in [START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF] but the proof we give here is far more simple.

Unlike what happens in a viscous fluid, collisions are allowed with our model as proved in [START_REF] Houot | On the motion and collisions of rigid bodies in an ideal fluid[END_REF]. How to handle collisions is the aim of Section 5. It is done by adding to the Lagrangian function an electrostatic-like potential energy whose electric charges are located along the boundaries of the solids and generate a repulsive force. The potential energy has to be set for the repulsive force to be neglectful when the solids are far apart but very strong when they get close in order to avoid collisions. This refinement in the model is another novelty of our work.

Then, in Section 6, we address the case of swimming ASBs as described above by introducing holonomic constraints in the motions of the solids. We use the work already done for a set of solids to derive the expanded form of the Euler-Lagrange equation (1.3). The equation we obtain eventually is a generalization of that given in [START_REF] Kanso | Locomotion of articulated bodies in a perfect fluid[END_REF] when hypotheses (H1-4) are relaxed and the buoyancy force taken into account. We study the well posedness of this equation and get results roughly similar to the previous ones, relating to rigid solids: existence of smooth solutions up to an hypothetic collision.

It is classical in the literature to use as controls directly the shape-changes of the bodies although a more realistic and interesting control would be the internal forces and torques causing these deformations. About this issue, we first explain in Section 7 how to compute a posteriori the internal forces, the shape-changes being given. Then we prove that there is a one-to-one correspondence between the shape-changes and the internal forces they result from. This new result entails the equivalence between controlling directly with the shape-changes and with the internal forces. This is of crucial importance for future works on the controllability of our model.

The two last sections of the article deal with numerics. In Sections 8 we describe a numerical scheme implemented in a Matlab toolbox (Biohydrodynamics Toolbox). It consists of an ODE solver coupled with an integral equations solver for the computation of the fluid potential. Integral equations are numerically solved by applying the Nyström method, which is fast and accurate. An important feature is that the elements of the mass matrices and the shape derivatives arising in the Euler-Lagrange equation can all of them be computed only in terms of the values taken by the fluid potential along the boundaries of the solids. It entails that there is no need to compute the fluid potential everywhere in the fluid domain but only on its boundary. Numerically, there is no need to mesh the 2d fluid domain but only its 1d boundary.

In the last section (Section 9) we present some features of the Matlab toolbox. As an example, we show that a result established by Sparenberg in [30, page 63] ensuring that: a body of finite extent, moving periodically through an inviscid and incompressible fluid without shedding vorticity, cannot exert a force with non-zero mean value, is no longer true when the fluid contains fixed obstacles. This result proves that a fish is not only able to propel itself without shedding vorticity but can also generate a non-zero mean thrust.

2 Dynamics of a set of submerged rigid solids

Notation

Let (e 1 , e 2 ) stand for a reference Galilean frame by which we identify the physical space to R 2 . In this frame, the coordinates of the vectors are denoted x := (x 1 , x 2 ) . At any time t > 0, we denote by S i (i ∈ {1, . . . , n}) the open domain occupied by the i-th solid. We set S := ∪ i S i and F stands for the open region of the fluid. With this notation, the domain M := F ∪ S, as being the domain of the overall fluid-solids system, is independent of time. All of the domains are denoted at the initial time t = 0 with subscript naught: S i 0 , S 0 and F 0 . We assume that ∂F 0 is Lipschitz continuous and since F undergoes only rigid deformations, ∂F has the same regularity for all time. We assume also that the solids touch neither each other nor the fixed boundaries of the fluid at the initial time, which means that Si 0 ∩ Sj 0 = ∅ if i = j and Si 0 ⊂ M for all i ∈ {1, . . . , n}. We denote by n = (n 1 , n 2 ) the unitary normal vector to ∂F directed toward the exterior of the fluid. The regularity of ∂F ensures that n is well defined almost everywhere. Likewise, we introduce the unitary tangent vector τ = (τ 1 , τ 2 ) to ∂F oriented such that {n, τ } be a direct basis (see fig 1). When ∂F is more regular, namely at least C 1,1 (continuously differentiable with first derivative Lipschitz continuous), we can also consider the curvature κ defined almost everywhere on ∂F by: ∂ τ n = -κτ .

Attached to each body, we define a moving frame (E i 1 , E i 2 ) whose origin coincides at any time with the center of mass of the solid and we denote with capital letters X i := (X i 1 , X i 2 ) the relating body coordinates.

Rigid motion

Rigid displacements are described by elements g of G := SE(2) (the rigid body or Euclidian group) consisting of pairs (R, r) where R ∈ SO(2) is a rotation matrix and r := (r 1 , r 2 ) ∈ R 2 is a vector. They would sometimes be identified with the matrix:

g = R r 0 1 ,
and for any vector x = (x 1 , x 2 ) ∈ R 2 , we denote gx the vector such that:

gx 1 := R r 0 1 x 1 . (2.1)
In the sequel, we will make no difference in the notation between the pair (R, r) and the associated isometry of R 2 defined by (2.1). For all (R, Ṙ) ∈ TSO(2) (the tangent bundle to SO(2)), the matrix ṘR is skew-symmetric and hence, as it is classical in Solid Mechanics, we introduce ω ∈ R such that:

ṘR x = ωx ⊥ , ∀ x ∈ R 2 , (2.2)
where x ⊥ := (-x 2 , x 1 ) stands the vector x := (x 1 , x 2 ) positively quarter turned. The manifold SO(2) is parameterized by means of matrices R(θ) (θ ∈ R/2π) defined by: R(θ) := cos(θ) -sin(θ) sin(θ) cos(θ) .

Thereby TSO( 2) is parameterized by pairs (θ, θ) ∈ R/2π × R and ω = θ in (2.2). We deduce that G is parameterized by (θ, r) ∈ R/2π × R 2 and TG by (θ, r, ω, ṙ)

∈ (R/2π × R 2 ) × (R × R 2 ).
Going back to our problem, and since the solids undergo rigid motions, at any time there exists g i := (R i , r i ) ∈ SE [START_REF] Amrouche | Dirichlet and Neumann exterior problems for the n-dimensional Laplace operator: an approach in weighted Sobolev spaces[END_REF] such that

E i j = R i e j ,
and r i gives the coordinates of the center of mass of the solid (see figure 1).

With our notation we have furthermore:

x = g i X i , (2.3)
and S i = g i (S i 0 ). We deduce also the relation between the area elements: dx = dX i .

Differentiating with respect to the time identity (2.3), we deduce that the Eulerian velocity of the i-th body is:

v i (x) : = ġi g i-1 x = ω i (x -r i ) ⊥ + ṙi , (x ∈ S i ).
We next use the partition ∪ i S i of S to define in S the overall Eulerian rigid velocity v of the solids by setting, for all x ∈ S:

v(x) = v i (x), if x ∈ S i .

Fluid dynamics

The flow being assumed to be irrotational, at any time t the Eulerian velocity of the fluid is equal to the gradient ∇ψ of a harmonic potential function ψ defined in F. The classical non-penetrating or slip condition for inviscid fluids yields Neumann boundary conditions for the function ψ, namely:

∂ n ψ = v • n on ∂S, (2.4)
∂ n ψ = 0 on ∂F \ ∂S.

(2.5)

0 S i F E j 2 S j ṙj E i 1 E i 2 e 1 e 2 θ i θ j ṙi E j 1 r j ∂F \ ∂S r i τ n τ n Figure 1:
The fluid domain is denoted F, the i-th solid occupies the domain S i , S := ∪ i S i and M := F ∪ S. The coordinates of the center of mass of the i-th solid in the reference frame (e k ) are r i := (r i 1 , r i 2 ) . The frame (E i k ) is attached to the solid and θ i gives its orientation with respect to (e k ). At the time t = 0, the domains and coordinates are denoted with subscripts naught: F 0 , S i 0 , r i 0 and so on...

For this NBVP to be well-posed, we have to add conditions ensuring the uniqueness of the solution. We choose for instance:

∂S ψ dσ = 0.
Results on the well-posedness of such problems are given in the Appendix, Section A.

For any solid i, we define three elementary potential functions: ψ i 1 , ψ i 2 and ψ i 3 , all of them being harmonic in F and satisfying the following Neumann boundary conditions on ∂F (using the convention that ∂ n ψ i j = 0 on every unspecified boundary):

∂ n ψ i 1 = e 1 • n on ∂S i , (2.6a)
∂ n ψ i 2 = e 2 • n on ∂S i , (2.6b
)

∂ n ψ i 3 = (x -r i ) ⊥ • n on ∂S i . (2.6c)
We invoke next Kirchhoff's law to obtain the following decomposition of the overall potential:

ψ := i ṙi 1 ψ i 1 + ṙi 2 ψ i 2 + ω i ψ i 3 .
(2.7)

Lagrangian function

We denote by ρ i > 0 the density of the i-th solid while ρ f > 0 stands for the constant density of the fluid. The mass and the momentum of inertia of the i-th body are respectively:

m i := S i ρ i dx, I i := S i ρ i |x -r i | 2 dx.
The kinematic energy of the i-th solid is:

K i := 1 2 m i |r i | 2 + 1 2 I i |ω i | 2 ,
and its potential energy relating to the buoyancy reads:

P i := g(m i r i -m i f r i f ) • e 2 , (2.8)
where g stands here for the standard gravity, m i f > 0 is the mass of the fluid occupying the same volume as the solid and r i f is its center of buoyancy. They are defined respectively by:

m i f := ρ f S i dx, r i f := ρ f m i f S i xdx.
When ρ i is constant in S i then r i f = r i . The kinematic energy of the fluid is given by:

K f := 1 2 ρ f F |∇ψ| 2 dx,
and we deduce the expression of the Lagrangian function of the fluid-bodies system:

L := K f + i (K i -P i ).

Generalized coordinates and mass matrices

Let us introduce the generalized coordinates:

q i := (r i 1 , r i 2 , θ i ) and q := (q 1 , . . . , q n ) ,
and the generalized velocities qi and q. At the time t = 0, they turn out to be:

q i 0 := (r i 0 , θ i 0 )
and q 0 := (q 1 0 , . . . , q n 0 ) .

The components of q will be subsequently denoted either q i with i ∈ {1, . . . , 3n} or q i j with i ∈ {1, . . . , n} and j ∈ {1, 2, 3} if we want to emphasize the solids' indices. We denote by Q the open subset of (R 2 × R/2π) n consisting of all the allowable q i.e. all the q for which the solids touch or overlap neither each others nor the boundary of the fluid domain. The set Q being an open set of a Banach space, it can also be seen as an analytic manifold and we denote T Q its tangent bundle consisting of the pairs (q, q).

The generalized coordinates will allow us to rewrite the Lagrangian function in a convenient short form. Beforehand, we need to introduce the mass matrices of the system. The mass matrix of the i-th body is defined by:

M i r := diag(m i , m i , I i ),
while the added mass matrix resulting from the remote interaction between the i-th and the j-th solids in the fluid reads:

M f ij := ρ f    F ∇ψ i 1 • ∇ψ j 1 dx . . . F ∇ψ i 1 • ∇ψ j 3 dx . . . . . . F ∇ψ i 3 • ∇ψ j 1 dx . . . F ∇ψ i 3 • ∇ψ j 3 dx    .
(2.9)

The overall virtual mass matrix of the fluid-solids system is next defined as being the bloc matrix:

M :=       M 1 r + M f 11 M f 12 . . . M f 1n M f 21 M 2 r + M f 22 . . . M f 2n . . . . . . . . . . . . M f n1 M f n2 . . . M n r + M f nn       .
The overall kinetic energy of the fluid-solids system can now be rewritten in a very short form as a matrix-vectors product:

K := 1 2 q M q, ( 2.10a)
and if we set furthermore:

P := i P i , (2.10b)
the expression of the Lagrangian function turns out to be:

L := K -P. (2.10c)
Observe that in the definition (2.9), the elements of the matrix M f can be rewritten, upon an integration by parts, as:

F ∇ψ i k • ∇ψ j l dx = 1 2 ∂S j ψ i k ∂ n ψ j l dσ + 1 2 ∂S i ψ j l ∂ n ψ i k dσ, (2.11)
and this expression involves only boundary integrals.

Euler-Lagrange equation

The Lagrangian is a function of q and q only. We invoke now Hamilton's principle which tells us that the Euler-Lagrange equation driving the dynamics of our system is:

d dt ∂L ∂ q - ∂L ∂q = 0, (t ≥ 0).
Expanding this equality and taking into account the expressions (2.10), we get a second order non-linear ODE:

Mq + Γ(q), q, q + ∂P ∂q = 0, (t ≥ 0), (2.12)
in which Γ(q) is a rank-3 tensor usually called Christoffel symbol in Lagrangian Mechanics. Denoting by M ij the elements of the virtual mass matrix M (i, j ∈ {1, . . . , 3n}) and by q i the elements of q (i ∈ {1, . . . , 3n}), it is defined by:

Γ k ij := 1 2 ∂M ki ∂q j + ∂M kj ∂q i - ∂M ij ∂q k , (i, j, k ∈ {1, . . . , 3n}). ( 2

.13)

We wish now to compute explicitly these quantities. Considering the expression of M, we realize that the main difficulty in carrying out this task is to differentiate the elements of the added mass matrix:

M f ij := ρ f F ∇ψ i i • ∇ψ j j dx, (2.14)
(i = 3(i -1) + i , j = 3(j -1) + j , i , j ∈ {1, . . . , n}, i , j ∈ {1, 2, 3}) with respect to the generalized coordinates. Indeed, the quantities (2.14) depend on q not only because the domain of integration F does but also because the elementary potentials are defined as solutions of NBVPs in F. Although quite involved, formulae to compute such derivatives are available. To display them we need to introduce some additional notation. We use in the sequel the results and the notation of the Appendix, Section B. We denote by w i k (i ∈ {1, . . . , n}, k ∈ {1, 2, 3}) the velocity fields defined on ∂S i by:

w i 1 := e 1 , w i 2 := e 2 , w i 3 := (x -r i ) ⊥ , so that the boundary conditions (2.6) can be rewritten as ∂ n ψ i k = w i k
• n for all i and k. Let us introduce also for all i ∈ {1, . . . , n} the functions ξ i , ζ i : ∂S i → R 3 , defined by:

ξ i :=   -n 2 n 1 (x -r i ) • n   , ζ i :=   n 1 n 2 (x -r i ) ⊥ • n   ,
and the matrices, defined on ∂S i as well, by:

S i := κ(ξ i ⊗ ξ i -ζ i ⊗ ζ i ) + ξ i ⊗ f 3 + f 3 ⊗ ξ i -[(x -r i ) • n](f 3 ⊗ f 3 ),
where f 3 := (0, 0, 1) . At last, we recall the definition of the Kronecker symbol:

δ i j := 1 if i = j, 0 otherwise.
We can state now the:

Proposition 2.1
The quantities (2.14) are analytic with respect to q. If ∂F is of class C 1,1 , the Christoffel symbols (2.13) can be computed in terms of the elementary potentials and we get, for all k = 3(k -1) + k , i = 3(i -1) + i and j = 3(j -1) + j with k , i , j ∈ {1, . . . , n} and k , i , j ∈ {1, 2, 3} (which means that k , i and j are the indices of the solids):

Γ k ij = - ρ f 2 ∂S i ∂ τ ψ k k ∂ τ ψ j j (w i i • n)dσ - ρ f 2 ∂S j ∂ τ ψ k k ∂ τ ψ i i (w j j • n)dσ + ρ f 2 ∂S k ∂ τ ψ i i ∂ τ ψ j j (w k k • n)dσ + δ j i ρ f ∂S i ψ k k S i i j dσ + δ j i δ k i ρ f 2 ∂S k (w k k • n)(w i i • n)(w j j • n)dσ. (2.15)
Notice that, like for the elements of the mass matrix, all of the terms in (2.15) are boundary integrals.

Proof: We study the sensitivity of the elements (2.14) of the matrix M f at the time t = 0, the method being similar at any time. We want to apply results of Section B, so we need first to define the mapping φ s arising over there.

Since the solids do no touch each other at the time t = 0, we can consider open domains S i ε and O i such that Si 0 ⊂ S i ε , Si ε ⊂ O i , Ōi ⊂ M for all i ∈ {1, . . . , n} and Ōi ∩ Ōj = ∅ for all i = j. According to Urisohn's Lemma, there exists a function χ of class C ∞ , compactly supported and such that 0

≤ χ ≤ 1, χ = 1 in S i ε and χ = 0 in R 2 \ ∪ i Ōi .
With the notation of Section B, we specialize s = q and we set:

φ s (x) := (1 -χ(x))x + χ(x)(R(θ i )(x -r i 0 ) + r i ) if x ∈ Ōi , x if x ∈ R 2 \ ∪ i Ōi .
For all q close enough to q 0 , the mapping φ s is an infinitely differentiable diffeomorphism from R 2 onto itself and

φ s (x) = R(θ i )(x -r i 0 ) + r i in S i ε .
We next apply Proposition B.1 to get the analyticity of the quantities (2.14).

We wish then to apply Proposition B.3. To this purpose, we need to compute the quantities involved in formula (B.3). Straightforward computations yield, on ∂S i and for all i ∈ {1, . . . , n}:

w i 1 • τ = -n 2 , ∂ τ (w i 1 • n) • τ = κn 2 , w i 2 • τ = n 1 , ∂ τ (w i 2 • n) • τ = -κn 1 , w i 3 • τ = (x -r i ) • n, ∂ τ (w i 3 • n) • τ = (-1 -κ)(x -r i ) • n.
At last, we get: In (2.12), it remains only to compute the derivatives of the potential energy (2.8) with respect to the generalized coordinates and it is straightforward to obtain:

γ ij (x) = -(x -r i ) • n if i = j = 3k, k ∈ {1, . . . ,
∂P i ∂q k j =      g(m i -m i f ) if i = k and j = 2, g(m i r i -m i f r i f ) ⊥ • e 2 if i = k and j = 3, 0 otherwise.
(2.16)

Model with circulation

Adding circulation means adding a certain number of degrees of freedom to the fluid flow. This number depends on p, the number of connected components of ∂F. We denote ∂F i (1 ≤ i ≤ p) these connected components. Obviously, since there are n solids in the fluid, we have p ≥ n and we set the indices such that ∂F i := ∂S i for 1 ≤ i ≤ n. Let us now introduce p -1 elementary stream functions: ϕ i , i ∈ {1, . . . , p -1}, harmonic in F and satisfying Dirichlet boundary conditions:

ϕ i = 1 on ∂F i , ϕ i = 0 on ∂F \ ∂F i .
Observe that with this definition, there is no need to consider the function ϕ i for i = p because it is equal to 1 -p-1 i ϕ i . At any time, the overall stream function relating to the circulation can be written as a linear combination of the functions ϕ i :

ϕ = i c i ϕ i , (3.1)
where the coefficients c i are real for all i ∈ {1, . . . , p -1} and depend on time.

We call c := (c 1 , . . . , c p-1 ) the generalized circulation variable, homogeneous to a generalized velocity. The unknowns of our problem are now q, q and c. The Eulerian velocity of the fluid becomes, when circulation is taken into account:

u := ∇ψ + ∇ϕ ⊥ in F.
We defined an additional symmetric (p -1) × (p -1) mass matrix:

C := ρ f    F ∇ϕ 1 • ∇ϕ 1 dx . . . F ∇ϕ 1 • ∇ϕ p-1 dx . . . . . . F ∇ϕ p-1 • ∇ϕ 1 dx . . . F ∇ϕ p-1 • ∇ϕ p-1 dx    ,
and we point out that its elements can be rewritten, using Green's formula, in the form:

F ∇ϕ i • ∇ϕ j dx = 1 2 ∂F i ∂ n ϕ j dσ + 1 2 ∂F j ∂ n ϕ i dσ, (3.2)
i.e. once more, as boundary integrals. The virtual mass matrix of the fluid-solids system is now the bloc diagonal matrix:

M 0 0 C ,
because all of the cross terms involving both an elementary potential and an elementary stream function are null. Indeed, for any i ∈ {1, . . . , n}, k ∈ {1, . . . , p -1} and j ∈ {1, 2, 3}, we have:

F ∇ψ i j • (∇ϕ k ) ⊥ dx = - F (∇ψ i j ) ⊥ • ∇ϕ k dx = - ∂F k (∇ψ i j ) ⊥ • ndσ = ∂F k ∂ τ ψ i j dσ = 0.
Observe that the elementary stream functions depend on the position of the solids and hence on q. We draw the same conclusion for the mass matrix C. The Lagrangian function given in (2.10) rewrites now:

L := 1 2 q M q + 1 2 c Cc -P.
Once more, Hamilton's principle tells us that (remember that c is homogeneous to a generalized velocity):

d dt ∂L ∂c = 0, (t ≥ 0), which simplifies into: d dt (Cc) = 0, (t ≥ 0).
It means that the impulse of the fluid (see [START_REF] Lamb | Hydrodynamics. Cambridge Mathematical Library[END_REF] for a definition) corresponding to the circulation is constant. We recover in particular a well known result: if there is no circulation at the time t = 0, the circulation remains null for all time. We denote Π f c ∈ R n the initial impulse related to the circulation and we get:

c = C -1 Π f c , (t ≥ 0). (3.3) The matrix C is positive because for any c ∈ R p-1 , c Cc = ρ f F |∇ϕ| 2 dx ≥ 0
where the function ϕ is given by (3.1). Furthermore c Cc = 0 if and only if ϕ = 0 which means that c = 0. We deduce that C is positive definite and hence always invertible. We next introduce the modified potential energy:

P := P + P c , with P c := - 1 2 (Π f c ) C -1 Π f c .
The Lagrangian function turns out to be:

L := 1 2 q M q -P ,
and the Euler-Lagrange equation (2.12) must be modified as follows:

Mq + Γ(q), q, q +

∂ P ∂q = 0, (t ≥ 0). (3.4)
For all k ∈ {1, . . . , n} and k ∈ {1, 2, 3}, we have:

∂P c ∂q k k := - 1 2 (Π f c ) (C -1 ) ∂C ∂q k k C -1 Π f c . (3.5)
The following Lemma allows to go further in the computation of this term:

Lemma 3.1 If ∂F is of class C 1,1
, the elements of the matrix ∂C/∂q k k read:

∂ ∂q k k F ∇ϕ i • ∇ϕ j dx = - ∂S k ∂ n ϕ i ∂ n ϕ j (w k k • n)dσ, (3.6)
for all i, j ∈ {1, . . . , p -1}.

Proof: We just give the outline of the proof, most of the arguments being similar to those of Section B. Using the same notation as in the proof of Proposition 2.1 we first compute ∂ si ϕ j for all i ∈ {1, . . . , 3n} and j ∈ {1, . . . , p -1}. Straightforward computations yield:

-∆∂ si ϕ j = 0 in F,

∂ si ϕ j = 0 on ∂F \ ∂S i ,
where i = 3(i -1) + i , i ∈ {1, . . . , n}, i ∈ {1, 2, 3}. Differentiating with respect to s i the remaining boundary condition:

ϕ j (φ s ) = δ i j on ∂S i 0 ,
we get: ∂ si ϕ j + ∇ϕ j • w i i = 0 on ∂S i . But since ϕ j is constant along ∂S i , we have ∇ϕ j = (∂ n ϕ j )n and the equality above turns into:

∂ si ϕ j = -∂ n ϕ j (w i i • n) on ∂S i .
Then, Reynolds's formula (see [11, Corollary 5.2.2 page 172]) yields for all i ∈ {1, . . . , 3n} and j, k ∈ {1, . . . , p -1}:

∂ ∂s i F ∇ϕ j • ∇ϕ k dx = ∂S i (∇ϕ j • ∇ϕ k )(w i i • n) dσ + F ∇(∂ si ϕ j ) • ∇ϕ k dx + F ∇ϕ j • ∇(∂ si ϕ k ) dx. (3.7)
Once more, we observe that ∇ϕ j = ∂ n ϕ j on ∂F and we apply Green's formula with the two last terms to obtain formula (3.6).

From equations (3.3) and (3.4) we deduce that the motion of the solids in the fluid cannot generate circulation providing that the circulation is zero at the initial time, but the converse is false: starting from rest, solids immersed in a fluid with non zero circulation are certainly set into motion by the fluid.

Well-posedness of the Euler-Lagrange equations

We study now the well-posedness of equations (2.12) and (3.4). The main result is:

Theorem 4.1 For any initial data (q 0 , q0 ) ∈ T Q, there exists a unique solution t ∈ [0, T ) → q ∈ Q to equation (3.4) such that (q(0), dq dt (0)) = (q 0 , q0 ). Further, the solution is analytic and either T = +∞ or T corresponds to the time of a collision between two solids or between a solid with any boundary of the fluid domain.

Proof: Since F is at least Lipschitz-continuous, we deduce that for all i, j ∈ {1, . . . , n} and k, l ∈ {1, 2, 3}, the quantities:

F ∇ψ i k • ∇ψ j l dx and F ∇ϕ i • ∇ϕ j dx,
are analytic with respect to q. It is a non-trivial result of shape sensitivity analysis proved in [START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF]. It entails that all of the terms in (3.4) are analytic with respect to q and also with respect to q because they are linear or quadratic with respect to this variable. The next argument is that the matrix M is definite-positive and hence invertible for all q ∈ Q. The Cauchy-Lipschitz Theorem applies to EDO (3.4) and ensures that for any initial data, there exists a unique analytic solution local in time.

It remains to prove that the solution (q, q) cannot blow up in finite time unless a collision occurs. To this purpose, we introduce the overall energy of the fluid-solids system:

E := 1 2 q M q + P .
Differentiating with respect to time the identity:

E = ∂L ∂ q • q -L,
and applying the chain rule, we obtain the relation:

d dt E = d dt ∂L ∂ q - ∂L ∂q • q,
which means that E is conserved. The solution t → (q, q) is therefore bounded in any compact subset of T Q. Classical results for ODEs tell that the solution can be continued up to a maximum time T and that either T = +∞ or:

lim t→T (q, q) ∈ ∂(T Q),
which means, since the tangent spaces are boundary-less vector spaces, that:

lim t→T q ∈ ∂Q.
This condition characterizes a collision (or a contact) and the proof is completed.

Collisions

Collisions between immersed rigid solids evolving freely in a potential flow can occur as proved in [START_REF] Houot | On the motion and collisions of rigid bodies in an ideal fluid[END_REF]. In the numerical simulations, collisions and overlapping of solids will cause the computation to fail and hence must be avoided. To this purpose, we add to our model an electrostatic-like potential energy defined by:

V := ε n i=0 n j=i+1 ∂S j ∂S i |x -x | -α dσ x dσ x ,
where we have denoted ∂S 0 := ∂F \ ∂S (i.e. ∂S 0 is the boundary of the fluid not shared with the solids) and ε > 0, α > 0 are two constants. We set them such that the repulsive force induced by the potential V be neglectful when the solids are not very close. The potential V must be subtracted to the right hand side term in the definition of the Lagrangian function and the Euler-Lagrange equation (3.4) becomes:

Mq + Γ(q), q, q + ∂ P ∂q

+ ∂V ∂q = 0, (t ≥ 0). (5.1)
We easily obtain that, for k = 3(k -1) + k , (k ∈ {1, . . . , n}, k ∈ {1, 2, 3}):

∂V ∂q k k = 2ε j =k ∂S j ∂S k |x-x | -α (x -x ) |x -x | 2 • n + κ (w k k •n)dσ x dσ x . (5.2)
Once again, observe that this expression involves only boundary integrals.

Dynamics of ASBs

To simplify, let us assume in this section that there is no circulation. The dynamics of the ASBs can easily be deduced from the dynamics of a set a solids. Indeed, ASBs are defined as a collection of rigid solids linked together by holonomic constraints which are easily handled within Lagangian formalism. We write that:

q = F (p, u),
where p := (p 1 , . . . , p m ) belongs to P, an analytic manifold of dimension m (1 ≤ m < 3n). The variable p and stands for the new degrees of freedom of the system and u := (u 1 , . . . , u d ) ∈ U (an analytic manifold of dimension d ≥ 1 such that m + d = 3n) is a given control variable governing the relative position of the solids composing the ASBs. To simplify we assume that F is also analytic from P × U into Q.

Remember that the Lagrangian L = L(q, q) is a function of q and q only. We define then:

L(p, u, ṗ, u) := L F (p, u), ∂F ∂p ṗ + ∂F ∂u u , (6.1)
and the equation governing the motion of the ASBs is:

d dt ∂L ∂ ṗ - ∂L ∂p = 0, (t ≥ 0). (6.2)
Introducing the matrices:

N p := ∂F ∂p , N u := ∂F ∂u ,
and the rank-3 tensor (the second order derivative of F ):

D 2 F, ṗ, u := ∂ 2 F ∂p 2 , ṗ, ṗ + ∂ 2 F ∂u 2 , u, u + 2 ∂ 2 F ∂p∂u , ṗ, u ,
the chain rule allows us to expand (6.2) as follows:

N p MN p p + N p M D 2 F, ṗ, u + N p Γ(F (p, u)), N p ṗ + N u u, N p ṗ + N u u + N p MN u ü + N p ∂P ∂q + N p ∂V ∂q = 0, (t ≥ 0). (6.
3)

The function F being given, the matrices N p and N u and the tensor D 2 F are explicitly known. We know how to evaluate all of the terms in (5.1) and therefore we draw the same conclusion for (6.3). Some care has to be taken when selecting the function F , for the ODE (6.3) to be well-posed. Let us consider a given control u : R + → U.

Theorem 6.1 Assume that t ∈ R + → u ∈ U is of classe C 2,1 (respectively C 2+k with k ≥ 1 or analytic) and assume that N p is full rank for all (p, u) ∈ P × U. Then, for any initial data (p 0 , ṗ0 ) ∈ T P there exists a unique solution t ∈ [0, T ) → p ∈ P to ODE (6.3) such that (p(0), dp dt (0)) = (p 0 , ṗ0 ). Further, the solution if of class C 2 (respectively C 2+k or analytic) and either T = +∞ or T corresponds to the time of a collision between two ASBs or between a body with any boundary of the fluid domain.

The hypothesis on the rank of N p is not surprising. It means that for any value of the control variable u, the function p → F (p, u) is locally an immersion from P into Q.

Proof: We use the same arguments as in the proof of Theorem 4.1.

The hypothesis on N p ensures that the matrix N p MN p is positive definite for all (p, u) and hence always invertible.

From the definition (6.1), we deduce that:

L := 1 2 ṗ N p MN p ṗ + ṗ N p MN u u + 1 2 u N u MN u u -P (F (p, u)).
We set then:

P * := - 1 2 u N u MN u u + P (F (p, u)),
and we introduce the energy like quantity:

E * := 1 2 ṗ N p MN p ṗ + P * .
Proceeding like in the proof of Theorem 4.1, we prove that E * remains constant for all time and hence that the solution cannot blow up unless a collision occurs.

Computation of internal forces and torques

To govern the shape-changes, in place of the variable u, a more realistic control variable would be the internal forces and torques causing these shape-changes. Actually, under an additional assumption on F , we will show that both choices are equivalent by proving that there is a one-to-one relation between u and the internal forces. The generalized forces, denoted by ν in the sequel, are defined in Lagrangian Mechanics by:

d dt ∂L ∂ u - ∂L ∂u = ν, (t ≥ 0). (7.1)
We can easily expand this identity by copy-pasting (6.3) and exchanging u and p for they play symmetrical roles:

N u MN u ü + N u M D 2 F, ṗ, u + N u Γ(F (p, u)), N p ṗ + N u u, N p ṗ + N u u + N u MN p p + N u ∂P ∂q + N u ∂V ∂q = ν, (t ≥ 0). (7.2)
Therefore, for any given smooth control t ∈ R + → u ∈ U, we can compute first p by solving (6.3) and next use it to compute ν with (7.2). The converse is true if the system of equations ( 6.3)-(7.2), seen as a system of ODEs in (p, u) can be solved for any given right hand side ν. So, let now the function t ∈ R + → ν ∈ R d be given.

Theorem 7.1 Assume that both matrices N p and N u are full-rank for all (p, u) ∈ P × U and that the function t → ν(t) is Lipschitz continuous (respectively of class C k with k ≥ 1 or analytic). Then, for any (p 0 , ṗ0 , u 0 , u0 ) ∈ T P × T U, the system of ODEs (6.3)-( 7.2) admits a unique solution t ∈ [0, T ) → (p, u) ∈ P ×U such that (p(0), dp dt (0), u(0), du dt (0)) = (p 0 , ṗ0 , u 0 , u0 ). Further, the solution is of class C 2 (respectively C 2+k or analytic) and either T = +∞ or T corresponds to the time of a collision between two ASBs or between a body with any boundary of the fluid domain.

The hypotheses on N p and N u entail that F is a local diffeomorphism from P × U onto Q.

Proof: Once more, we proceed like in the proof of Theorem 6.1 and we obtain the existence of a solution t → (p, u) ∈ P × U local in time. We denote by T the maximal time of existence of the solution, we assume that T < +∞ and we define the natural overall energy by:

E := 1 2 ṗ N p MN p ṗ + ṗ N p MN u u + 1 2 u N u MN u u + P (F (p, u)).
We next easily check that:

∂L ∂ ṗ • ṗ + ∂L ∂ u • u -L = E.
Differentiating this equality with respect to time, we get:

d dt ∂L ∂ u - ∂L ∂u • u = dE dt , (7.3)
which means that the variation of energy is equal to the power-like amount ν • u expended by the internal forces. Assume now that F (p, u) remains in a compact subset K of Q for all t ∈ [0, T [. Then, there exist α > 0 and P min ∈ R (the minimum of the potential energy P in K) such that, for all (p, u) ∈ F -1 (K):

α| u| 2 ≤ E -P min .
According to (7.3), we deduce that for all t ∈ [0, T [:

α| u| 2 ≤ λ + t 0 |ν|| u|ds, (7.4)
where λ := |E(0) -P min |. Setting then:

φ(t) := t 0 |ν|| u|ds,
we obtain after some basic algebra, that for all t ∈ [0, T [:

φ (t) φ(t) + λ ≤ |ν| √ α .
Integrating this inequality with respect to time, we get the estimate:

φ(t) ≤ λ + √ λ + 1 2 √ α t 0 |ν|ds 2 .
After plugging this result into (7.4), we just have proved the Gronwall-type inequality:

| u| 2 ≤ 2 α λ + 1 α √ λ + 1 2 √ α t 0 |ν|ds 2 ,
which means that u is bounded for all t ∈ [0, T [. Considering (7.3), this entails that E is also bounded and hence that ṗ is bounded as well. Classical behavior results for solutions of ODEs tell us that this result contradicts the hypothesis T < +∞. We deduce that either T = +∞ or lim t→T F (p, u) ∈ ∂Q, which means a collision in our system at the time T .

One of the interest in computing the internal forces and torques is that it allows us to quantify the efficiency of a locomotion strategy. Indeed, a relevant cost functional associated with a displacement over a time interval [0, T ] could be:

T 0 |ν| 2 ds. (7.5)
This approach was for instance adopted in [START_REF] Kanso | Optimal motion of an articulated body in a perfect fluid[END_REF] to seek optimal strokes.

8 Numerical scheme

Introduction

The equations governing the motion of a collection of solids or ASBs are ODEs. Therefore our numerical scheme will involve an ODE solver. Consider the ODE (5.1) without the circulation term (circulation has not been implemented yet) and let us turn it into the normal form:

d dt q q =   -M -1 Γ(q), q, q + ∂P ∂q + ∂V ∂q q   . (8.1)
Each evaluation of the right hand side term requires the computation of the elementary potentials. Therefore, the numerical scheme will also involve a NBVPs solver.

The amounts we need to compute are more precisely the elements of the mass matrices:

F ∇ψ i k • ∇ψ j l dx,
for all i, j ∈ {1, . . . , n} and k, l ∈ {1, 2, 3} as well as their derivatives with respect to q. As already mentioned before, all of these quantities can be expressed as boundary integrals (see formulae (2.15) and (2.11)). This is also true for the repulsive force ∂V /∂q (see (5.2)).

Similarly, all of the elementary potentials, as being solutions of NBVPs, are completely determined by their (Neumann) data on the boundaries of the solids (see Appendix, Section A). Based on these considerations, the numerical scheme will involve only computations on ∂F. It is a crucial point because one of the main difficulty in dealing with fluid-structure problems is that the fluid domain is not fixed but depends on the unknown positions of the immersed structures. Therefore, if computations have to be performed in this domain, it must be remeshed at each time step. Integral formulations allow this problem to be avoided. Another advantage of this approach is that the discretization of an integral equation set on the 1d fluid boundary leads to a linear system having many less unknowns than the discretization of the same problem set in the 2d fluid domain. At last, when the fluid domain is not bounded, dealing with integral equations allows to avoid the tricky problem of truncating the domain and determining artificial boundary conditions.

How to turn boundary value problems (2.6) into integral equations is a classical process (see [START_REF] Mikhlin | Mathematical physics, an advanced course[END_REF]). For our problem, we have two choices for the integral formulation: the direct or the indirect type. The direct formulation meets the problem requirements because it gives directly the Dirichlet boundary data in terms of the Neumann boundary data, that is, it corresponds to the Neumannto-Dirichlet operator. It reads (see [START_REF] Atkinson | The numerical solution of integral equations of the second kind[END_REF]) for all i ∈ {1, . . . , n}, k ∈ {1, 2, 3} and x ∈ ∂F:

ψ i k (x) - 1 π ∂F (x -x) |x -x| 2 • n(x ) ψ i k (x )dσ x = - 1 π ∂S i log |x -x|b i k (x )dσ x , (8.2)
where we denote by b i k := ∂ n ψ i k = w i k • n the Neumann boundary data of ψ i k and by n, the unitary normal to ∂F directed toward the exterior of the fluid. An additionnal equation is needed to ensure the uniqueness of the solution when F is bounded:

∂S ψ i k (x )dσ x = 0. (8.
3)

The next step consists in selecting a numerical method to solve (8.2) (coupled with (8.3) in the bounded case). In most of the articles, the authors use the so-called panel method [START_REF] Kanso | Optimal motion of an articulated body in a perfect fluid[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 planar fluid[END_REF]. It is a collocation method where the approximation function space consists of piecewise constant functions built from a discretization of the boundary by small curves or straight segments.

However, when the boundaries are regular enough, several authors [START_REF] Atkinson | The numerical solution of integral equations of the second kind[END_REF][START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF] advocate the Nyström method as being the best possible choice because it is almost as simple to implement as the panel method but has far better convergence properties. In the following Subsections 8.2 and 8.3, we describe this method applied to our NVBP. In Subsection 8.4 we give some convergence results and Subsection 8.5 is devoted to the comparison between the Nyström method and the panel method.

The Nyström method

The accuracy of the Nyström method depends on the quadrature rule choosen to approximate the integrals. For any i ∈ {1, . . . , n}, we introduce a C ,1 ( ≥ 2) parameterization of the boundary ∂S i :

γ i : s ∈ [0, 2π[ → γ i (s) ∈ R 2 . (8.4)
It is obtained merely by composing a given parameterization of the fixed curve ∂S i 0 with the rigid displacements g i := (R i , r i ). Similarly, we denote by γ i (i ∈ {n + 1, . . . , N }) the parameterizations of the remaining N -n fixed boundaries (not shared with moving solids). Via these parameterizations, all of the functions defined on the fluid boundaries can be seen as 2π-periodic functions for which the rectangle rule is very well adapted as proved in [START_REF] Kress | Numerical analysis[END_REF]Corollary 9.27 page 210] (see also Appendix, Lemma C.2).

We next define for all j, l ∈ {1, . . . , N } the kernels K l j and R l j by:

K l j (t, s) := (γ j (s) -γ l (t)) |γ j (s) -γ l (t)| 2 • n(γ j (s))|γ j (s)|, R l j (t, s) := log |γ j (s) -γ l (t)||γ j (s)|, (s, t ∈ [0, 2π[),
and we rewrite (8.2) for all i ∈ {1, . . . , n} and k ∈ {1, 2, 3} as:

ψ i k (γ l (t)) - 1 π N j=1 2π 0 K l j (t, s)ψ i k (γ j (s))ds = - 1 π 2π 0 R l i (t, s)b i k (γ i (s))ds, (l = 1, . . . , N ). (8.5)
For any pair of indices (i, k) ∈ {1, . . . , n} × {1, 2, 3}, (8.5) is a linear system of N integral equations with N unknowns. The unknowns are the pieces of the solution ψ i k living on each one of the N boundaries. The kernels K l l (t, s) have no singularity at t = s because:

lim s→t K l l (t, s) = 1 2 κ(γ l (t))|γ l (t)|,
where we recall that κ is the curvature of the boundary.

Notation:

The process we are going to describe now is similar for any pair of indices (i, k), so we will drop them in the notation. For instance, we will merely denote ψ in place of ψ i k .

The Nyström method consists of two steps:

step 1
The domain [0, 2π[ of each function γ j (j ∈ {1, . . . , N }) is uniformly discretized with m j := 2n j +1 points (n j ≥ 1) denoted s j p := 2πp/m j (p = 0, . . . , m j -1). The integrals in the left hand side of (8.5) are approximated for all j, l ∈ {1, . . . , N } by means of the rectangular quadrature rule:

2π 0 K l j (t, s)ψ(γ j (s))ds 2π m j mj -1 p=0 K l j (t, s j p )ψ(γ j (s j p )), (t ∈ [0, 2π[).
This approximation is quite accurate because, as already mentioned, s → K l j (t, s)ψ(γ j (s)) is a smooth 2π-periodic function. The quadrature rule is used also for the right hand side of (8.5) except when l = i because the kernel R i i (t, s) has an (integrable) singularity at t = s. This problem is overcome by using the following decomposition (see [3, page 329]), available for any (s, t) ∈ [0, 2π[×[0, 2π[:

log |γ i (t) -γ i (s)| = log 2e -1/2 sin t -s 2 + B i (t, s), (8.6)
where:

B i (t, s) :=      log |e 1/2 (γ i (t) -γ i (s))| |2 sin((t -s)/2)| , if s = t, log |e 1/2 γ i (t)| if s = t.
The function B i (t, s) is a smooth periodic kernel without singularity at s = t and the corresponding integral is approximated using again the rectangular rule with m i = 2n i +1 points. To integrate the singular part in (8.6), we use the following formula, specializing ϕ(s)

:= |γ i (s)|b i k (γ i (s))|:
-

1 π 2π 0 ϕ(s) log 2e -1/2 sin t -s 2 ds = 1 √ 2π   ϕ 0 + |p|>0 ϕ p |p| e ipt   . (8.7a)
The equality is true for any 2π-periodic smooth function ϕ defined on [0, 2π[. The approximation of the right hand side of (8.7a) is done by

• truncating the series at p = ±n i ;

• replacing the Fourier coefficients ϕ p (|p| ≤ n i ) by their discrete Fourier transforms Φ p (based on the same rectangular quadrature rule using the m i = 1 + 2n i points s i j of [0, 2π[). These two points correspond to the following matrix vector product:

Φ = √ 2π m i F mi Φ, (i = 1, . . . , n), (8.7b)
where:

Φ := [ϕ(s i 0 ), ϕ(s i 1 ), . . . , ϕ(s i mi-1 )] , (8.7c)
Φ := [ Φ 0 , . . . , Φ ni , Φ -ni , . . . , Φ -1 ] , (8.7d)
and F m is defined for any m ≥ 1 as the usual discrete Fourier matrix:

F m := e -2iπkj/m 0≤k,j≤m-1 . (8.7e)
Denoting by ψ l (t) and y l (t) (l ∈ {1, . . . , N }) the approximations obtained by a similar process for respectively ψ(γ l (t)) and for the right hand side of (8.5), the N discretized equations (8.5) read:

ψ l (t) - N j=1 2 m j mj -1 p=0 K l j (t, s j p ) ψ j (s j p ) = y l (t), (l = 1, . . . N ). (8.8)
Observe that this equation allows to define the numerical solution ψ = [ ψ 1 , . . . , ψ N ] everywhere on the boundary, once the quantities ψ j (s j p ) are known. It is the so-called Nyström interpolation formula. step 2 Finally the N equations (8.8) are turned into an equivalent linear square system whose unknowns are ψ j (s j p ), by imposing that the identities (8.8) be exact for t = s l p (p = 0, . . . , m l -1). We get then:

N j=1 A l,j Ψ j = Y l , (l = 1, . . . , N ), (8.9)
where each block matrix A l,j is of size m l × m j . More precisely, we have:

A l,j l =j q,p = - 2 m j K l j (s l q , s j p ), 0 ≤ q ≤ m l -1, 0 ≤ p ≤ m j -1, A l,l q,p = δ q,p - 2 m l K l j (s l q , s l p ), 0 ≤ q, p ≤ m l -1, Ψ j = [ ψ j (s j 0 ), . . . , ψ j (s j mj -1 )] , Y l = [ y l (s l 0 ), . . . , y l (s l m l -1 )] , l = 1, . . . , N.
As already mentioned, Y l for l = i is built using the rectangular rule:

[Y l ] q = - 2 m i mi-1 p=0 R l i (s l q , s i p )b i k (γ i (s i p )), 0 ≤ q ≤ m l -1.
For l = i, the special treatment described throughout the equations (8.7) leads to:

Y i = DΦ + F -1 mi diag(w)F mi Φ, [D] q,p = - 2 m i B i (s i q , s i p ), 0 ≤ q, p ≤ m i -1, w = [1, 1, 1/2, . . . , 1/n i , 1/n i , . . . , 1/2, 1] ,
where the term DΦ corresponds to the rectangular rule applied to the smooth kernel B i . The other part is computed fastly in O(log(m i )m i ) operations using the fast Fourier transform algorithm followed by m i divisions and finally by the inverse fast Fourier transform (this reads simply ifft(w.*fft(Phi)) in Matlab language).

The overall matrix A corresponding to the linear system (8.9) is defined by blocks:

A :=    A 1,1 . . . A 1,n . . . . . . A n,1 . . . A n,n    .
Unlike the right hand side term in system (8.9) which depends on the potential ψ i k , the matrix A is the same for all of the elementary potentials. When F is bounded, the matrix A is not invertible any longer (or very badly conditioned). It is not surprising since the solution of the continuous problem is defined only up to a constant and must be supplemented by (8.3) to recover uniqueness. This problem can be overcome by replacing one of the equations by the linear equation resulting from the discretization of (8.3):

n l=1 2π m l m l -1 p=0 ψ i k (s l p )|γ i (s l p )| = 0. ( 8.10)
For symmetry reasons, the removed equation is first substracted to all of the others before being replaced by (8.10). Although we did not carry out a rigorous analysis, this choice has resulted better, after several numerical tests, than the same method but without substraction and than the classical method consisting in imposing a value to the potential at a given point.

Additional details on the scheme

The linear system (8.9) must be solved for all of the elementary potentials ψ i k (i ∈ {1, . . . , n}, k ∈ {1, 2, 3}), that is to say for 3n different right hand side terms built from the 3n boundary data b i k . Since the matrix A remains unchanged, the more efficient way to carry out this task consists in using a LU factorization rather than an iterative method especially when n is large.

When the computations of the elementary potentials is done, a lot of contour integrals must be approximated in (2.11) and (2.15). We use again the rectangular formula with the same boundary discretization (s j p = 2πp/m j , p = 0, . . . , m j -1). Some terms in (2.11) and (2.15) require the knowledge at the integration points s j p of the tangential derivatives of the potentials. This is done by trigonometric interpolation using the discrete Fourier transform and its inverse.

Eventually, we tried out several Matlab built-in ODEs solvers and selected ode113 which is a variable order Adams-Bashforth-Moulton PECE solver (we refer to the Matlab documentation for details).

Accuracy of the scheme

The purpose of this subsection is to study how accurate is the approximation of the right hand side of ODE (8.1) obtained when the potentials are computed with the Nyström method. More precisely, the quantities to be evaluated are the entries of the mass matrix M f (see (2.11)) having the form:

∂S j ψ i k ∂ n ψ j l dσ = ∂S j ψ i k (w j l • n) dσ, (8.11)
and the elements of the rank-3 tensor Γ k ij which involves terms of the form (see (2.15)):

∂S i ∂ τ ψ k k ∂ τ ψ j j ∂ n ψ i i dσ = ∂S i ∂ τ ψ k k ∂ τ ψ j j (w i i • n)dσ. (8.12)
According to [3, page 333], the accuracy of the Nyström method depends on the smoothness of both the Neumann boundary data and the boundaries parameterization. Denoting by h := max i 2π/m i the size of the discretization, we obtain in our case:

Lemma 8.1 Assume that the boundary ∂F is of class C ( ≥ 2), then we have, for all i ∈ {1, . . . , n} and k ∈ {1, 2, 3}:

ψ k i -ψ k i L ∞ (∂F ) = O(| log(h)|h -1
(h ) + O(h -1 | log h|) = O(h -1 | log h|).
Integrals (8.12) also involve 2π-periodic functions of class C -1,1 but we must estimate first the error of the approximated tangential derivatives. Applying Proposition C.1 with a perturbation of size e L ∞ (∂F ) = Ch -1 | log h|, we obtain an error bounded by:

C| log h|h + C | log h| h Ch -1 | log h| = O((log h) 2 h -2 ),
where C and C are two positive constants. Then, applying Lemma C.2 with this perturbation, we obtain an overall error for (8.12) of order O(h

)+O((log h) 2 h -2 ) = O((log h) 2 h -2 ).
Summarizing the results of Lemma ).

An experimental study of the accuracy of the scheme involving the ODE solver will be presented in Subsection 9.2.

Nyström method versus panel method

As already mentioned before, the panel method is a collocation method of order 1: the error between the exact solution ψ and its approximation ψ h on ∂F satisfies ψ h -ψ ∞ = O(h) with possibly superconvergence property at the collocation points t j : |ψ h (t j ) -ψ(t j )| = O(h 2 ) (see [10, remark 4.4.8 page 85]). However, this superconvergence property requires an exact approximation of the boundary geometry unlike what is done in most of the cases [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 planar fluid[END_REF][START_REF] Kanso | Swimming due to transverse shape deformations[END_REF] where the boundary is only approximated with small segments. In this case, the approximation in O(h) can not be improved, whatever is the regularity of the boundary. On the contrary, the accuracy of the solution computed with the Nyström method increases along with the regularity of ∂F. Indeed, according to Lemma 8.1 we have the estimate:

ψ h -ψ ∞ = O(| log(h)|h -1
), the convergence becoming even exponential with C ∞ boundaries. The elements (8.11) of the mass matrix M f , when the potentials are computed with the panel method and assuming superconvergence property, will next be obtained with an error of order O(| log h|h) (most likely worth than that if segments are used to approximate the boundary) versus O(| log h|h -1 ) with the Nyström method. When the boundary is of class C with > 3, the Nyström method is definitively more accurate (or many less points are needed to obtain the same accuracy).

It does not really make sense to compare both methods for the computation of the elements (8.12) of the rank-3 tensors arising in the rhs of the main ODE (8.1). Indeed, this part of the ODE is due to the fact that hypotheses (H1-4) have been relaxed and is an original point of this paper. 9 Biohydrodynamics Matlab Toolbox

Introduction

The Biohydrodynamics Toolbox (BhT) is a Matlab Toolbox designed to perform numerical simulations of rigid solids immersed in a potential flow. The system fluid-solids can be either confined, partially bounded or unbounded. The buoyant force and collisions between solids are supported (but circulation is not yet). The solids can be free of constraints or linked together in order to constitute ASBs. In the latter case, the constraints are prescribed by the user as functions of time and the induced motion is computed. Post-processing functions are provided to compute the internal torques, the power and the energy expended by the swimming bodies.

BhT is free (distributed under license GNU GPL) and comes with complete html documentation including plenty examples and tutorials. It is available at:

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

• http://www.mathworks.com/matlabcentral/fileexchange/21872 (hosted by Matlab Central).

Going through all of the features of BhT would be too long and is out of the scope of this article. So after giving some results on how the numerical scheme performs, we will rather use BhT to discuss and illustrate the following Sparenberg's affirmation (see [30, page 63]): A body of finite extent, moving periodically through an inviscid and incompressible fluid without shedding vorticity, cannot exert a force with non-zero mean value. In particular, we will show that this statement seems no longer true when the fluid domain contains fixed obtacles.

Convergence of the scheme

The purpose of this subsection is not to present a rigorous numerical analysis of BhT's solver but only to provide some benchmarks and to illustrate how the numerical scheme performs. Keep in mind that BhT's solver is an ODE solver (namely the buit-in MATLAB solver ode113) coupled with a solver of integral equations (based on Nyström's method). We set the relative tolerance of the ODE solver ode113 to 1e-14 with the option odeset('RelTol',1e-14).

We consider two identical neutrally buoyant ellipse-shaped solids of semimajor axis 2 and semiminor axis 0.9. The coordinates of their centers of mass are respectively (2, -1.1) and (-2, 1.1) at the initial time and their initial linear velocities (-0.5, 0) and (0.5, 0). Both have zero initial angular velocity. We study how the number of points used in the discretization of the ellipses boundaries, impacts the accuracy of the trajectory computation. We compute first the trajectories and velocities of the ellipses with 1890 points and use this solution as the reference one. Next, we decrease the number of points and perform again the same computation. The results are displayed in Figures 3456below. Observe that all of these graphs have log scale y-axis. In figure 7 is ploted the CPU time required for the computation of the trajectories with respect to the number of points (with log scale x and y-axis). The left hand side of the graph is not relevant since the time used by MATLAB's interpreter is larger than that used for computations. From 100 points of discretization, the time grows as the number of points to the power 2.3. However, if we assume that for a larger numbers of points most of the CPU time will be used to factorize at each time step a n × n linear system (n being the number of points), it is expected to asymptotically grow as the number of points to the power 3.

Simulation settings

We consider a fish-like swimming body made of four ellipse-shaped rigid solids as pictured in figure 8. The three angles of the joints are driven by the functions drawn in fig 9 over the time interval [0, 6π]. Notice that all of the functions are compactly supported. The fish starts from rest, undergoes shape-changes yielding a net forward motion and stop dead in its initial shape. We are interesting in studying its final velocity. Indeed, denoting by m > 0 the total mass of the fish and by ṙ the velocity of its center of mass, first Newton's law tells us that:

1 6π mṙ(6π) = 1 6π 6π 0 F dt, (9.1)
where F stands for the hydrodynamical force generated by the fish to propel itself. Thereby, the right hand side in (9.1) is precisely the mean value of this force over the time interval [0, 6π]. We are going to perform twice the same simulation. In a first time, the fish will be alone in the fluid and the system fluid-fish will fill the whole space. In a second time, the fluid will contain also two fixed rectangular obstacles and the fish will swim between them in its course (see figure 10).

The movies are available at:

http://bht.gforge.inria.fr/Examples/demos.html.

Results

With BhT we can compute the velocity of the center of mass of the fish at any time. It is pictured in figures 11 and 12.

We clearly observe in figure 11 that the final velocity of the fish is zero, which confirms Sparenberg's statement this case. On the other hand, with the same sequence of strokes, figure 12 tells us that the fish generates non zero mean thrust when the fluid contains obstacles.

An other interesting point that can be observed using BhT is that the torques applied at the joints to produce the shape-changes are not equal in both cases although the deformations are. In figure 13 we have drawn the values of the torques for the first case (fluid without obtacles) and in figure 14 the difference between the values of the torques in the first and second cases.

Computing next the cost of each displacement as quantified by the cost functional (7.5), we obtain that the second strategy is 3.22% cheaper than the first one. These surprising results can be summarized by saying that not only generates the fish non-zero mean thrust by swimming between the obstacles but also it swims more effortlessly. Assume now that ∂F is C 1,1 (continuously differentiable with first derivative Lipschitz continuous) and g ∈ H 1/2 (∂F), then the solution ψ of (A.3) belongs to:

X 2 (F) := {ψ ∈ H 2 (F) : ∆(∂ xi ψ) ∈ L 2 (F), ∀ i = 1, 2}.
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.

A.2 Unbounded domain

When F is not bounded, we define the weight function:

ρ(x) := 1 + |x| 2 log(2 + |x| 2 ), (x ∈ R 2 ),
and the Lebesgue space:

L 2 N (F) := {ψ ∈ D (F) : ρ -1 ψ ∈ L 2 (F)},

B Shape sensitivity analysis

This section is self-contained and independent, including for the notation. We recall results of shape optimization theory, most of them being proved in [START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF]. Let F 0 ⊂ R 2 and define for all s = (s 1 , . . . , s n ) ∈ R n (n ≥ 1) diffeomorphisms φ s of class C m (m ≥ 1) from R 2 onto itself, equal to the identity outside a large fixed compact ball B. We denote by D m the Banach space of such diffeomorphisms, endowed with the norm of the Sobolev space W m,∞ (B, B) (uniform convergence in B of the function and the derivatives up to order m).

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 any i ∈ {1, . . . , n}, let then ψ i be the solution of the NBVP:

-∆ψ i = 0 in F, ∂ n ψ i = w n i on ∂F,
and define for any i, j ∈ {1, . . . , n} the functional:

J ij (s) := F ∇ψ i • ∇ψ j dx.
Providing that F 0 is Lipschitz continuous, we have: The proof is given in [START_REF] Munnier | On the self-displacement of deformable bodies in a potential fluid flow[END_REF] in a more general framework. Assuming additional regularity on ∂F 0 , we can compute the partial derivatives of J ij : Proposition B.2 Assume that F 0 is of class C 1,1 and that the mapping s → Then, for all i, j, k ∈ {1, . . . , n}, we have:

∂ ∂s i F ∇ψ j • ∇ψ k dx = - ∂F (∂ τ ψ j ∂ τ ψ k )w n i dσ - ∂F κw i • w j -γ n ij + (∂ τ w n j • τ )w τ i + (∂ τ w n i • τ )w τ j ψ k dσ - ∂F κw i • w k -γ n ik + (∂ τ w n k • τ )w τ i + (∂ τ w n i • τ )w τ k ψ j dσ + ∂F ∂ n ψ i ∂ n ψ j ∂ n ψ k dσ.
Straightforward but tedious computations yield:

Proposition B.

C Error estimates for the numerical approximations

In this section we recall or slightly adapt some results about Fourier series and interpolating trigonometric polynomial series which are needed in Subsection 8.4. In all the sequel, f will be a 2π-periodic function of class C ,α with ≥ 0 and α ∈ (0, 1]. For any k ∈ Z, we denote by fk its Fourier coefficients: fk := 1/2π 2π 0 f (t)e -ikt dt and for any positive integer n, F n (f ) stands for its truncated Fourier series: F n (f )(t) := |k|≤n fk e ikt . When this series is convergent, its sum is denoted F(f ). We also consider a uniform partition of [0, 2π] built with m := 2n + 1 points t j := jh (j = 0, . . . , m -1, h := 2π/m). The interpolating trigonometric polynomial of f of degree n is denoted I n (f ). It is defined by: |D n (t -t j )|.

I n (f )(t) =
From the definition of the points t j = 2πj/m and the fact that |D n (t)| is an even 2π-periodic function we deduce that:

  Figure 2 displays some screenshots of what the solids motion looks like and the animation can be watched on the web page http://bht.gforge.inria.fr/ Examples/demos.html. The ellipses never touch and their deviations with respect to straight lines are only due to hydrodynamic forces (if hydrodynamically decoupled, the trajectories would be straight lines).

Figure 2 :

 2 Figure 2: Screenshots of the ellipses at different times. The dashed blue lines show the trajectories of the centers of mass.

Figure 3 :

 3 Figure 3: Error for the position of the centers of mass with respect to the number of points

Figure 4 :

 4 Figure 4: Error for the orientation of the ellipses with respect to the number of points

Figure 5 :Figure 6 :

 56 Figure 5: Error for the linear velocities with respect to the number of points

Figure 9 :

 9 Figure 9: Angles of the joints with respect to time (given functions). Joint 1 corresponds to the fish's head and joint 3 to the tail.

Figure 10 :

 10 Figure 10: Fish swimming between two obstacles

Figure 11 :

 11 Figure 11: Velocity of the center of mass of the fish with respect to time (fluid without obstacle).

Figure 12 :

 12 Figure 12: Velocity of the center of mass of the fish with respect to time (fluid with obstacles).

Figure 13 :Figure 14 :

 1314 Figure 13: Values of the torques at the joints with respect to time in the first case.

|k|≤n c k e ikt , c k = 1 m m- 1 j=0fm m- 1 j=0D

 11 (t j )e -2iπkj/m . Another useful expression of I n (f ) given in[START_REF] Atkinson | Theorical Numerical Analysis[END_REF] involves the Dirichlet kernel:I n (f )(t) = 2 n (t -t j )f (t j ), D n (t) := sin((n + 1/2)t) 2 sin(t/2) . (C.1)Subsequently, we will use the notation C and C for positive constants depending on f but not on the discretization parameters (m, n or h) and C will stand for any pure constant. The first result concerns the convergence of the Fourier series and the interpolating trigonometric polynomial series. The proof can be found in[4, pages 128-129].Lemma C.1 For any integer n ≥ 2, we have:f -F n (f ) ∞ ≤ C log n n +α and f -I n (f ) ∞ ≤ C log n n +α .We consider now the quantity I(f ) := 2π 0 f (t)dt and its approximation based on the rectangular rule:I R m (f ) := h m-1 j=0 f (t i ).Lemma C.2 The error between I(f ) and I R m (f ) can be estimated as follows:|I(f ) -I R m (f )| ≤ C m +α = C h +α . (C.2a)If we replace f (t i ) by f (t i ) = f (t i ) + e i , an approximation of f (t i ) with an error e i , then we get:|I(f ) -I R m ( f )| ≤ C h +α + 2π e ∞ . (C.2b)Proof: The estimate (C.2a) is trivial when = 0. When ≥ 1, F n (f ) converges absolutely to F(f ) and F(f ) = f . We compute first the error for the functions e ikt and we deduce that:I(f ) -I R m (f ) = -2π p∈Z * fpm .But the Fourier coefficients of a function of class C ,α satisfy the estimate: | fk | ≤ C|k| --α (see[START_REF] Zygmund | Trigonometric Series[END_REF],p. 45). The conclusion follows. The estimate (C.2b) being obvious, the proof is then completed.We define b(t) := (2/m) m-1 j=0 |D n (t -t j )| for all t ∈ R and we state: Lemma C.[START_REF] Atkinson | The numerical solution of integral equations of the second kind[END_REF] The following estimate holds:max t b(t) ≤ C n log n.Proof: Proceeding as in[28, page 15], we easily obtain that b is an even 2π/mperiodic function. So we can seek its maximum on the interval [0, π/m]. We use the crude bound:

  t∈[jπ/m,(j+1)π/m] |D n (t)|. Using the expression D n (t) = 1/2 + n k=1 cos(kt) of the Dirichlet kernel we get D n (t) = -n k=1 k sin(kt). We deduce that |D n | is increasing on the interval [0, π/m] (because for any t ∈ [0, π/m], nt ≤ nπ/m = nπ/(2n + 1) < π/2) and next that M 0 ≤ M 1 . To estimate M j (j = 1, . . . , m -1) we differentiate the expression (C.1) of the Dirichlet kernel to obtain that:D n (t) = 1 2 sin(t/2) (n + 1/2) cos((n + 1/2)t) -12sin((n + 1/2)t) tan(t/2) .

  ). If the boundary is of class C , then the Neumann boundary data b i k = w i k •n are of class C -1 . We next apply the result given in [3, page 333]. Assume the parameterization of the boundary ∂F if of class C ,1 ( ≥ 3). Then we can approximate the integrals (8.11) with an error of order O(h -1 | log h|) and the integrals (8.12) with an error of order O((log h) 2 h -2 ). The parameterization being of class C ,1 , the Neumann data b i k := w i k •n are C -1,1 . Classical results on the regularity of solutions of NBVPs ensure that the solutions ψ i k are C ,1 on ∂F. So the functions ψ i k (w j l • n) arising in (8.11) are C -1,1 on ∂F. These functions are first approximated by ψ k i (w j l • n), where ψ k i is computed by means of the Nyström method. So, integrals (8.11) involve a 2π-periodic function of class C -1,1 which is known up to a perturbation term bounded by Ch -1 | log h| (C a positive constant). Lemma C.2 gives for these integrals an error of size O

	Proof: Based on the results given in Appendix, Section C, we next show:
	Lemma 8.2 Proof:

  8.1 and Lemma 8.2, we obtain: Assume that the parameterization of the boundary ∂F if of class C ,1 ( ≥ 3), then the rhs of ODE (8.1) is computed with an error of order O((log h) 2 h -2

	Proposition 8.1

3

  Under the assumptions of Proposition B.2 and for all i, j, k ∈ {1, . . . , n}, the following identity holds:ψ k • ∂ τ ψ i )w n j dσ -∂F (∂ τ ψ k • ∂ τ ψ j )w n i dσ + ∂F (∂ τ ψ i • ∂ τ ψ j )w n k dσ -2 ∂F κ w i • w j -γ n ij + (∂ τ w n j • τ )w τ i + (∂ τ w n i • τ )w τ j ψ k dσ + ∂F ∂ n ψ k ∂ n ψ i ∂ n ψ j dσ. (B.1)

	∂ ∂s j	∂F	ψ i ∂ n ψ k dσ +	∂ ∂s i	∂F	ψ j ∂ n ψ k dσ -	∂ ∂s k	∂F	ψ i ∂ n ψ j dσ =
	-

∂F

(∂ τ

A Neumann boundary value problem

This section is dedicated to recalling some results stated in [START_REF] Brezis | Analyse fonctionnelle[END_REF][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], on the well-posedness of NBVPs. We write the problem in the generic form:

-∆ψ = 0 in F, (A.1a)

∂ n ψ = g on ∂F, (A.1b) and we assume that 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]). For problem (A.1) to admit solutions, g has to satisfy the compatibility condition:

A.1 Bounded domain

In this case, we introduce the function space

For any ψ ∈ X 1 (F), its normal trace ∂ n ψ exists in H -1/2 (∂F). Lax-Milgram Theorem ensures that there exists a unique weak solution ψ ∈ X 1 (F) satisfying both identities:

Moreover, we have an estimate: where D (F) is the space of distributions. Classical Sobolev spaces must be replaced by weighted Sobolev spaces:

When g ∈ H -1/2 (∂F), the solution of (A.3a-A.3b) exists in H 1 N (F) and is unique. Further, it can be proved that the solution belongs to X 1 N (F). As in the bounded case, if ∂F is C 1,1 and g ∈ H 1/2 (∂F), then the solution is in X 2 N (F).

One easily shows that, for all t ∈ R:

which implies that

and so M j ≤ m/(2| sin(πj/2m)|) for j = 1, . . . , m -1. Now using the inequality sin x ≥ 2x/π, (x ∈ [0, π/2]), we get 1/sin(πj/2m) ≤ m/j whence we deduce that:

and the proof is completed.

In Subsection 8.4, we need to compute boundary integrals involving the tangential derivatives of some potentials. Tangential derivatives are approximated by the derivatives of interpolating trigonometric polynomials.

Proposition C.1 Assume that ≥ 1. Then we have the error estimate:

If we replace f (t i ) by an approximation f (t i ) = f (t i ) + e i then we get:

Proof: Since we have, for all integer n, F n (f ) = F n (f ) , we can write that:

According to Lemma C.1, the first term is bounded by C log n/n -1+α . For the second term we have, since

Using the expression (C.1), we get:

and therefore:

|D n (t -t j )|.

Using Lemma C.1 and Lemma C.3, we obtain estimate (C.3a). To prove (C.3b), we write that:

D n (t -t j )e j , and the estimate for the last term is obtained by applying again Lemma C.3. The proof is then completed.


Dominique Bergogne
email: dominique.bergogne@insa-lyon.fr

Asif Hammoud

Dominique Tournier

Cyril Buttay

Youness Amieh

Pascal Bevilacqua

Abderrahime Zaoui

Hervé Morel

Bruno Allard

Abderahime Zaoui

Electro-thermal behaviour of a SiC JFET stressed by lightning-induced overvoltages

Keywords: <<JFET>>, <<Silicon Carbide>>, <<Voltage Source Converter>>

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

SiC JFET inverters have been presented at various conferences over the last few years. High temperature and high power densities were the main subjects. A paper from SiCED [START_REF] Friedrichs | Behavior of high voltage sic vjfets under avalanche conditions[END_REF], states the ruggeness of the devices when submitted to repetitive avalanche. Short-circuit and over load behaviour as also been studied in [START_REF] Boughrara | Robustness of sic jfet in short-circuit modes[END_REF]. In some applications the inverter output conductors might be struck by lightning induced over-voltage. The modelling of lightning induced over-voltages is developped in ref [START_REF] Nucci | Lightning-induced voltages on overhead lines. Electromagnetic Compatibility[END_REF]. This work focuses on SiC JFETs under lightning-induced over-voltages according to standards : CEI 801-5 Electromagnetic compatibility of electrical equipment, Part 5, class 3 (latest reference is CEI 61000-4-5). The JFETs are studied in the electrical environment of a voltage source inverter. We have developped a surge generator and a method to estimate the internal temperature of the JFET by using thermo-sensitive electric parameters. The novelty of the paper is the experimentation with JFET carying much more current than the channel saturation current. In such conditions, standard models are out of range and thermal behaviour must be taken into account. This paper describes an experimental setup to provide data for electro-thermal modelling of the JFET under overload.

Lightning surge generator

The european standard for electromagnetic compatibility of electric and electronic equipments [START_REF] Hirschi | Protection d'équipements électriques et électroniques contre les surtensions transitoires. montena emc sa[END_REF] states the references for the normative lightning-induced overvoltages. A so-called 1.2/50µsec pulse, is used in our experiment with a peak voltage of 2kV (open load) and a maximum current of 120A (short-circuit). These values have been verified by simulation and measured experimentally on the prototype generator, see figure1. For information purposes here are the values used in the generator : C=5uF, L=50uH, R1=16.5Ω, R2=11Ω and R3=66Ω. K1 is a high voltage IGBT driven through a fiber-optic link, L1 is an air-coil inductor to get saturation-free operation. The generator circuit and the pre-charge voltage are set to comply with CEI 61000-4-5 class 3 voltage surge immunity tests regulation. This class applies for moderately exposed equipment, with voltages up to 2kV. Higher levels of surge energy depends on the global qualification of the whole electrical environment of the inverter, this is outside the scope of this paper. Inverter, JFET and lightning surge

Assumptions

On a practical voltage fed inverter, the DC bus impedance is low because the stray inductance has been reduced to improve efficiency and bus capacitance is high to sustain the reactive power of the load. In figure 2 this translates as the following assumption : DC voltage source Zr is low enough and DC bus capacitor C is large enough to maintain a constant Ubus voltage despite the lightning stroke. In that case, the voltage source can be replaced by a short-circuit and experimental measurements can be carried out directely on a SiC JFET stuck by a lighning surge. Another consequence is that applying lightning surges on the DC bus is not relevant.

Simplification

In the inverter of figure 2 each JFET, as any power switch, can be driven-On or Off by the gate driver. In the experiment described here, the simplification implies studying one JFET at a time, which is possible if one assumes that the 'other' JFET is always Off. Turning the opposite JFET On during a surge will cause the two JFET to produce a short circuit, this case is not within the scope of this paper. In ref [START_REF] Boughrara | Robustness of sic jfet in short-circuit modes[END_REF], short-circuit behaviour of JFETs is studied in details. In the practical case of an inverter, the JFET gate driver will turn the opposite JFET Off -short circuit protection-so the simplification proposed is valid.

Three cases

The lightning surge may have a positive or negative polarity when applied to the inverter's terminal. Inside, two JFETs forms a leg, each JFET can be controlled On or Off at the obvious exeption of both JFETs On. Table I summarises the 6 possible cases induced by the preceeding statements. But one more information must be included here : a JFET is able to conduct current despite an 'Off' control applied by the Gate driver. A reverse voltage will force a current through the internal diode of the JFET, an Off JFET doesn't mean there is no current. Table I is used to draw the 6 corresponding circuits of figure 3. The inverter's DC link is represented, and, as assumed earlier, the impedance of the DC-bus is low compared to the JFET's ability to conduct current. From that, three cases remain to be studied :

The experiment

The principle is to aply the voltage surge produced by the lighning stroke generator to a JFET, for various biasing conditions -the three cases-and to measure a thermo-sensitive parameter to compute the temperature response. All measurements are transients, single shot, so no averaging or synchonous detection techniques apply.

Thermal sensing Thermo-sensitive parameters

To get the internal temperature of the device during a surge, for the three cases listed before, two separate thermo-sensitive parameter (TSP) are used.

The choice of two thermo-sensitive parameters, according to the JFET Gate biasing, is to simplify the experiment and reduce possible errors. The experiment uses the thermo-sensitive parameter that is directely available on the JFET, without modifing the Gate biasing. Using a single parameter would mean changing Gate bias during the stroke, thus creating one more source of difficult-to-control transients and adding complexity to the circuit to be taken into account for modelling. First, when the Gate of the JFET is negative in such a way that no current flows in the channel, the device is equivalent to a diode. Hence, the choice of the forward voltage drop in a diode at fixed constant current. It is a good linear temperature indicator despite a poor gain (2mV/K). Second, when the JFET is turned On, the Gate to Source voltage is nulled, the internal diode in shunted by the channel conduction. Here, the temperature estimator is the channel resistance at reduced voltage, to make sure no conduction occurs in the diode. Here the gain corresponds to the RDSON of the JFET, which is, again, quite low (0.3 Ohm at room temperature).

Thermal calibration

Thermo-sensitive parameters are calibrated prior to the experiment by applying a fixed temperature for a long enough time on the JFET and measuring the electrical response. A fixed current source is used for both parameters: diode voltage drop and On-resistance. To avoid self heating during the experiment and possible interaction between the stroke and the sensing circuit, the measuring current is applied during a short period through a controlled switch. The full experiment circuit diagram is described in the next paragraph, in figure 5. During calibration a constant current pulse is applied on the device under test. The electrical circuit is not altered between calibration and measurement. The same control unit is used for sequencing the switches (see figure 5), the difference during calibration is that the stroke carries no power (Vc=0) and the temperature is set externally by a thermo-regulated heated slab. The TSP is an image of temperature, fitting the data sets permits to write the equations of temperature as a function of V F for each TSP, see table II.

V F [V] 1,8 2 2,2 2,4 2,6 T
Table II: Equations used to sense the temperature with thermo-sensitive parameters.

Diode @ 200mA T diode = -5, 17.10 2 V F + 1, 37.10 3 Channel @ 300mA T On-resistance = -1, 17.10 2 + 1, 31.10 3 V F -1, 42.10 3 V 2 F + 6, 99.10 2 V 3

F

Thermal sensing during a stroke

Thermal sensing through electrical measurements during a stroke is performed by interrupting the stroke at increasing time steps, isolating the device under test, and then, applying a fixed current to get a voltage response. Each step (during each stroke) produces one temperature value associated with a time delay.

The set of values can then be supperimposed on a single plot to form a 'real time' thermal response. Four switches are used and are controlled by a common programable logic unit, triggered by the push of a button. Now lets give a look at figure 5. K1 is the main switch of the stroke generator, it connects the precharged capacitor C to the passive network that produces the normalized stroke waveform. K2 permits to interrupt the stroke by diverting the current to ground, K2 is operated after a variable duration ∆t1 and remains closed until the total energy of capacitor C is dissipated in the network. K3 isolates the stroke generator from the JFET under test to allow the voltage measurement of Vf, after a deadtime ∆t3 of 1µs. A second deadtime of 1µs is placed before the application of the measuring current Im by K4. This switch is turned On for a fixed duration, ∆t4 of 10µs. This is when the TSP voltage Vf is acquired. It reflects the temperature response of either the diode or the channel. t1 is increased gradually to the entire stroke's duration (up to 200µs) to measure the temperature step by step.

Experimental results

Due to severe Electro-Magnetic-Interference (EMI) it is not possible to get the measured TSP voltage during a period of several µs after the interruption of the stroke. For that reason, linear extrapolation is used to 'go back in time' from acquired data. The acquisitions are visible on the picture of figure 6, blue dots are calculated from each separate acquisition. The error due to the electrical measurement is estimated to 10°C. The error due to extrapolation is more complex to estimate and a filtering effect can be observed.

A first section describes the experiment for cases 3 and 5, as in Table I, a stroke is applied to provoque a direct conduction from Drain to Source, in the channel only. A second section treats both cases 1 and 4 and cases 2 and 6, as the JFET is stressed in reverse conduction.

For all experiments, a specific Gate biasing circuit is used. This circuit is described at the end of this section. To get a conduction in the direct conduction mode, one has to bias the Gate-Source voltage to zero volt. A short circuit is an obvious answer. Nevertheless, care must be taken to control the Gate current in order to prevent destruction. We have set a biasing circuit (see figure 8) that permits a high common mode impedance for EMI reduction. More, the stroke generator produces a strong EMI across the whole circuit. Carefull decoupling of common mode and single ground connexion is compulsory. One more difficulty is the dynamics of the signals: within 1us, measured current and voltages change by a ratio of 1000. The estimated temperature is plotted in figure 6. The temperature rise is not important, this is due to the limitation we have put on the input voltage : 1kV instead of 2kV as stated in the CEI 61000-4-5 class 3 standard. This limit is directed by several facts : increasing the energy on the circuit increases EMI, thus degrading the precision of the measurements. More, we have observed an unexpected voltage variation between Gate and Source as shown on figure 6. The voltage is supposed to be set at zero by the Gate biasing circuit. During the stroke, there is a voltage appearing between Gate and Source, reaching -15V. The breakdown voltage of the Gate is -19V for the device we have used. For safety reasons, we have not increased the energy for this experiment. For the sake of comparison, the energy level, although low, was kept identical for all tests. The difference between the plots of figure 7 is the Gate bias value. To get the internal diode of the JFET to conduct on its own, a negative Gate bias is applied. Both channel and diode may conduct if the Gate is zero biased (in reference to the Source).

Direct conduction

Gate biasing circuit

A specific circuit, (see figure 8) is designed to obtain safe Gate biasing and a good control of the channel conduction. A high common mode impedance is obtained by using high values for Ra and Rb (typically 68 kΩ). Cg is a high value film capacitor for get a low differential impedance in order to maintain a constant voltage, 10 µF. Rg emulates the output impedance of a Gate driver, here a 33 Ω resistor. The key role of Ra-Rb and Cg is to drive safely the JFET Gate-Source into punch-trough in order to be sure that no current flows through the channel. Ra-Rb limit the static Gate current to acceptable values, while Cg maintains the voltage during transients. The dynamic Gate current is limited by Rg.

+ -

Ra

Verification of thermal sensing by simulation

To verify the consistency of the estimated temperature we have developped a thermal model in witch the experimentally measured power is used as an input. The simulation is purely thermal as the building of an electro-thermal model of the JFET under over load is very complex. For simplicity, the thermal model is a classical R-C equivalent thermal circuit as shown on figure 9. The SiC chip is 'cut' into 30 regions getting wider as the distance from the active area increases. We also made simplifications due to the fact that only short durations are observed : a single thermal path from the top of the JFET to the backside. The backside is kept at room temperature. It is not important to increase the precision of the thermal model knowing that the temperature estimation has an error in the order of 10°C. By using experimental measurements of current and voltage accross the JFET, instantaneous power is computed and used as a power source in the thermal model. This technique permits to separate the thermo-electrical behaviour of the JFET from the purely thermal cooling phenomenom in order to verify the temperature estimation. In the next sections, the simulated temperature is plotted on the same graphs as the experimentally estimated temperature for each of the three experiments (the three cases). In all cases, the matching is not close but the curves shows the same behaviour. It is in direct conduction that the maximum temperature is reached. This observation can be linked to the observation of the Gate-Spource voltage. For direct conduction, the Gate is biased at zero volt by the circuit of figure 8, thus letting the channel to be normally-On. During the stroke, the overload current produces a Drain-Source voltage that can be as high as 20V, as a consequence, the Gate potential is moved, and an external Gate-Source voltage of -15 can be observed. We suppose that in such conditions the channel is no more fully On, further investigation is needed at this point.

Direct and inverse conduction

T [°C]

Discussion

On the plot in figure 11, as in the one in direct conduction, the simulated response (solid line) shows a peak that is delayed from the measured temperature using TSPs. As the same triggering signal is used in all experiments and as the reference signal for the thermal simulation is the current and voltage acquired with that trigger, trivial error is not possible. Another information from the plots is that the estimated temperature (measured) response looks like a low-pass filtering has occured. This can be a thermal phenomenom inside the JFET as power dissipation is not located exactely where the TSP region is situated. More, the TSP voltage is extrapolated over a segment of 10µs. For each of the three cases the current path differs hence different regions of the device are stressed differently. In this paper, we do not pretend to address the physical aspect of the problem, nevertheless, important variations on experimental thermal responses asks for an attempt of explanation. As expressed in the previous section, errors on the use of TSP take part in the difference observed between experiment and simulation. Another source of incertainty is the thermal model that as probably too simple for very short periods of time. And more generally, the difference of geographic situation of the heated region versus the TSP region inside the JFET. In figure 12, for example, the Drain-Source diode is situated at the top of the chip while power dissipation occurs in both the channel and the diode in reverse conduction with zero volt Gate bias. In that case the TSP used is the Diode forward voltage.

Further study of the Gate behaviour and its hability to withstand current surges will enable high power strokes to be tested and perhaps the thermal limits of SiC clearly showed?

Figure 12: Simplified cross section of JFET

Figure 1 :

 1 Figure 1: Normative lightning stroke generator : Generator circuit and measured short-circuit current.

Figure 2 :

 2 Figure 2: JFET voltage inverter stressed by lightning-induced surges on it's output terminals

Figure 3 :

 3 Figure 3: Possible current paths of a lighning surge applied at the output of a JFET inverter.

  200 250 300 350 400 VGS=0 V F = Diode forward voltage drop @ 200mA.V F = On-voltage across channel @ 300mA

Figure 4 :

 4 Figure 4: Calibration curves of thermosensitive parameters.

Figure 5 :

 5 Figure 5: Experimental setup circuit and switching sequence.

Figure 6 :

 6 Figure 6: Experimentally estimated thermal response and Gate voltage.

Figure 7 :

 7 Figure 7: Experimentally estimated thermal response

Figure 8 :

 8 Figure 8: Gate biasing circuit.

Figure 9 :

 9 Figure 9: Equivalent electrical circuit for the thermal behaviour of the JFET chip at very short time durations

Figure 10 :

 10 Figure 10: Comparison of simulation versus experimental thermal response.

  Figure 11: Comparison of simulation versus experimental thermal response.

Table I :

 I Possibles cases when applying a lightning surge at the output of an inverter.

	Case number	1	2	3	4	5	6
	Surge Polarity	+	+	+	-	-	-
	High side JFET Off On Off Off On Off
	Lowside JFET	Off Off On Off Off On
	• Conduction of the internal diode only, such as in cases 1 and 4
	• Conduction in the channel only, such as in cases 3 and 5
	• Conduction in both channel and diode, such as in cases 2 and 6
	1		2				3
	4		5				6