Keywords: Locomotion, biomimetic, eel-like robot, hyperredundant robot, Newton-Euler algorithms, Lie groups

This paper presents the dynamic modeling of a continuous three-dimensional swimming eel-like robot. The modeling approach is based on the "geometrically exact beam theory" and on that of Newton-Euler, as it is well known within the robotics community. The proposed algorithm allows us to compute the robot's Galilean movement and the control torques as a function of the expected internal deformation of the eel's body.

I. INTRODUCTION

n this paper, we will present the preliminary results of a multi-disciplinary research project supported by the French CNRS. The goal of the project is to design and control a three dimensional eel-like robot. As many authors working in the biomimetic robotics community have noted, eel-like robots present an interesting perspective for improving the efficiency and maneuverability of underwater vehicles [START_REF] Triantafyllou | Optimal Thrust Development in Oscillating Foils with Application to Fish Propulsion[END_REF][START_REF] Mason | Experiments in Caranguiform Robotic Fish Locomotion[END_REF][START_REF] Morgansen | Trajectory Stabilization for a Planar Caranguiform Robot Fish[END_REF][START_REF] Mcisaac | A Geometric Approach to Anguilliform Locomotion Modelling of an Underwater Eel Robot[END_REF][START_REF] Kato | Guidance and Control of Fish Robot with Apparatus of Pectoral Fin Motion[END_REF]. The prototype we are designing will be a hyper-redundant robot made by connecting many parallel platforms. Moreover in order to guarantee efficient propulsion, it will be covered with a continuous deformable organ, which will mimic the eel's skin. This paper essentially deals with the macroscopic modeling of the future prototype. By "macroscopic" we mean a "high-level" model, which can be used as the basis for a preliminary design of the system and its control strategy. In particular, the macroscopic model does not take into account the detailed technology of the prototype but rather an ideal dynamic behavior useful to fix the guidelines of the project. The new results reported in this paper are multiple. Firstly, and contrary to most previous research on the same topic, the investigated robot is capable of three-dimensional swimming. Secondly, it is based on a continuous model adapted to the macroscopic modeling of the future hyperredundant prototype and to the continuous character of its skin. Some authors, using the concept of backbone curves  Manuscript received July 27, 2004. This work was supported in part by the French CNRS within the research programme "ROBEA". Frédéric Boyer is with Institut de Recherche en Communications et Cybernétique de Nantes (IRCCyN, 1 rue de la Noë, B.P. 92101 44321 Nantes cedex 3, France), école des Mines de Nantes (e-mail : frederic.boyer@emn.fr). Mathieu Porez is with Institut de Recherche en Communications et Cybernétique de Nantes (IRCCyN), école centrale de Nantes (e-mail: mathieu.porez@irccyn.ec-nantes.fr). Wisama Khalil is with Institut de Recherche en Communications et Cybernétique de Nantes (IRCCyN), école centrale de Nantes (e-mail: wisama.khalil@irccyn.ec-nantes.fr). [START_REF] Chirikjan | An Obstacle Avoidance for Hyper-Redundant Manipulators[END_REF][START_REF] Chirikjian | The Kinematics of Hyper-Redundant Robot Locomotion[END_REF][START_REF] Choset | A Follow-The-Leader Approach to Serpentine Robot Motion Planning[END_REF], have previously studied continuous modeling of hyper-redundant manipulators. In order to apply this kind of idea to the dynamics of the 3-dimensional framework, we adopted here the geometrically exact theory of beams in finite deformation, originally by J.C. Simo [START_REF] Simo | A Finite Strain Beam Formulation. The Three-Dimensional Dynamic Problem. Part I: Formulation and Optimal Parametrization[END_REF][START_REF] Simo | A Three-Dimensional Finite-Strain Rod Model. Part II: Computational Aspects[END_REF][START_REF] Simo | On the Dynamics in Space of Rods Undergoing Large Motions -A Geometrically Exact Approach[END_REF]. The idea consists in considering the eel robot as a beam defined by a continuous assembly of rigid cross sections and controlled through distributed laws of internal strains or torque. With this choice, just as in works dealing with articulated locomotion systems [START_REF] Ostrowski | The Geometric Mechanics of Undulatory Robotics Locomotion[END_REF][START_REF] Ostrowski | The Mechanics of Undulatory Locomotion: the Mixed Kinematic and Dynamic Case[END_REF][START_REF] Kelly | Geometric Phases and Robotic Locomotion[END_REF], the head dynamics of the continuous eel are derived on a fiber bundle. But in our case, while the fiber is still the group SE(3) of the head displacements, the shape space is no longer a finite dimensional manifold but rather a functional space of curves in a Lie algebra. In fact, the shape space will be parameterized by the field of the infinitesimal transformations of the cross sections along the eel's backbone. Moreover, in accordance with the works of Simo, the eel's body dynamics will be written on the space of position-orientation of the beam cross sections with respect to the earth frame, i.e. a functional space of curves in a Lie group. However, contrary to the numerical approach proposed by Simo to integrate the dynamics of passive beams, the dynamics problem considered here is not solved with the standard numerical tools of nonlinear structural dynamics but using the "Newton-Euler philosophy" of rigid robotics [START_REF] Luh | On-Line Computational Scheme for Mechanical Manipulator[END_REF][START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF][START_REF] Featherstone | The Calculation of Robot Dynamics Using Articulated-Body Inertias[END_REF]. Finally, the proposed approach turns out to be a generalization of the Newton-Euler based algorithm of Luh and Walker [START_REF] Luh | On-Line Computational Scheme for Mechanical Manipulator[END_REF] applied to the case of a continuous robot with a mobile base (here imitating the eel's head). The algorithm gives the motion of the eel and the control torque evolution as outputs in terms of the deformation time-law of its body as inputs. As is well known from rigid robotics, the recursive nature of the Newton-Euler approach allows us to obtain efficient and fast algorithms, which are very simple to implement. Moreover, it gives us a straightforward link to the modeling of the future poly-articulated prototype. Finally, as far as the interaction of the fluid with the eel is concerned, for control requirements, we need to model the contact in a simple manner regarding the robustness of our future closed loop controllers. Two simple analytical models suited to our purposes exist. Both are based on the Fluid Mechanics' theory of the slender body [START_REF] Katz | Low Speed Aerodynamics[END_REF]. The biomechanics community suggests the first one, the second is offered by the ocean-engineering community. The first is a result of the "Large amplitude elongated body theory of fish locomotion" by Lighthill [START_REF] Lighthill | Mathematical Biofluid-dynamics[END_REF]. This model is based on the basic assumption of the existence of some slices of fluid transversally transported with the cross sections of the eel.

Macro-continuous Computed Torque Algorithm

for a Three-dimensional Eel-like Robot I Frédéric Boyer, Mathieu Porez and Wisama Khalil (IEEE senior member)

Then from kinetic conservation laws, the undulation of the eel's body generates the propulsion by reaction. Nevertheless, this model has been restricted until now to planar swimming. Hence, in order to investigate threedimensional swimming, we use the second model, which is, today, devoted to the dynamics of underwater flexible cables [START_REF] Burgess | Bending Stiffness in a Simulation of Undersea Cable Deployment[END_REF]. In this second model, like in Lighthill's planar solution, fluid forces are introduced through a local analytical model written for each transverse slice of the cable. Moreover this model takes into account, not only some inertial terms (like the Lighthill model does), but also some drag (transversal) and viscous (tangential) forces.

The paper is organized as follows: we will start by briefly presenting the relationships between nonlinear beam theory and hyper-redundant robots designed by the assembly of parallel platforms (section II). Then, some basic definitions and notations about the macroscopic beam model of the prototype are given (section III). Based on this model, the fluid-structure interaction is modeled (section IV). Section V is devoted to the continuous kinematic model of the eel. In section VI, we give the eel's head dynamics based on the principle of virtual works. In the Newton-Euler algorithm, these dynamics give the boundary conditions of the continuous kinematic models. Then the dynamic model of the eel's body is presented (section VII). It is based on the calculus of variations applied to a Cosserat medium [START_REF] Cosserat | Théorie des corps déformables[END_REF], and leads us to a set of "partial differential equations" directly linked to those of Reissner's beam [START_REF] Reissner | On a One-Dimensional Large Displacement Finite-Strain Theory[END_REF]. Finally, all these results are combined in section VIII in order to obtain the computed torque algorithm. The paper ends with numerical examples (section IX) and some concluding remarks (section X).

II. BEAM'S KINEMATICS AND HYPER-REDUNDANT ROBOTS

We will consider here the case of a hyper-redundant robot produced by a serial assembly of parallel platforms. The robots investigated like so are of "snake-like" or "elephant trunk" type. In this context the rigid platforms mimic the vertebrae while the kinematics between any two platforms play the role of inter-vertebral kinematics. We develop the dynamics analysis using the geometrically exact theory of nonlinear beams, as developed in the eighties by J.C. Simo [START_REF] Simo | A Finite Strain Beam Formulation. The Three-Dimensional Dynamic Problem. Part I: Formulation and Optimal Parametrization[END_REF][START_REF] Simo | A Three-Dimensional Finite-Strain Rod Model. Part II: Computational Aspects[END_REF][START_REF] Simo | On the Dynamics in Space of Rods Undergoing Large Motions -A Geometrically Exact Approach[END_REF]. In this theory, the beam is considered as a continuous assembly of rigid sections of infinitesimal thickness, i.e. a one-dimensional Cosserat medium [START_REF] Cosserat | Théorie des corps déformables[END_REF]. In the context of robotics, the sections imitate the robot's vertebrae and the beam kinematics, that of the inter-vertebral kinematics. Pushing the analogy forward, the centroidal line of the beam plays the role of the backbone of the snake-like robot. In this framework, one of the first tasks is to relate parallel platforms with the corresponding beam kinematics, where some of them are not related to any standard beam theory. In order to present these correspondences briefly, let us first introduce a few definitions related to the geometrically exact beam theory. Firstly, the material abscissa along the beam axis is denoted by X, which positions a particular cross section in a reference configuration o Σ (see Figure 1). The reference configuration is considered to be straight and represents a Galilean reference to which we fix the earth frame: ( ), ( )) t X t X , where G(X) is the center mass of the X section, and 1 t supports the beam axis 1 E when it is in the configuration o Σ . Let us point out that X X F is not the field of Frenet-Serret frames since each X F is actually attached to the X cross section and not deduced from the geometry of the deformed backbone curve.

The advantage of this choice is that the torsion is no longer geometric but actually related to the torsion strain field of the beam. Now, let us introduce the rigid transformations of SE(3) mapping the cross sections before deformation onto their configurations after deformation (defining the whole deformed configuration of the eel ( ) t Σ (see Figure 1)):

( ) ( ) [0,1] ( ) 0 1 R X d X X gX   ∈ =     (1) ,
where R is the rotation matrix mapping the X mobile basis before deformation onto that after, and d is the vector displacement field of the section centers. With the adopted parametrization (1), the deformation of the beam is defined by the field of twist

[ ] 0,1 ( ) (3) X X s e η ∈ ∈
modeling the infinitesimal changes of the mobile frames situation when sliding along the beam's material line of centroids, i.e., s.t.:

/ . g X gη ∂ ∂ = (2)
, which can be expressed as:

/ / . 0 0 0 1 0 0 R X d X R d K   ∂ ∂ ∂ ∂     Γ =               (3) ,
where we introduced the skew-symmetric tensor:

ˆ( / ) T K R R X =
∂ ∂ associated to the axial vector K. The last two components 2 K and 3 K of K stand for the curvatures of the beam in the two planes 1 3 ( , , )( ) G t t X and 1 2 ( , , )( ) G t t X , while 1 K is the rate (per unit of material length) of rotation of the section around its normal vector, i.e. the torsion strain field. As for the infinitesimal translations between two sections, we introduced the vector ( / )

T R d X Γ = ∂ ∂ , whose
first component 1 Γ is related to the stretching of the beam while the two others are related to its transverse shearing [START_REF] Simo | On the Dynamics in Space of Rods Undergoing Large Motions -A Geometrically Exact Approach[END_REF]. An eel-like robot can then be considered as a beam controlled by a desired time evolution of η :

( , ), d t X t X η ∀ , under the assumption that local controllers are able to impose the desired strains Γ and K instantaneously. We shall see later, in the particular case of our eel-like robot, how to relate the strain law ( , ) K Γ with that of the motor torques. Now, let us relate the general beam kinematics (2)-(3) to the corresponding hyperredundant robots.

To do this, we shall introduce a few possible parallel structures, starting from the most general case to the most specific one, by constraining the time evolution of η more and more. In the most general case, all six components of the vectors K and Γ can be defined by some arbitrary time evolutions. The corresponding parallel platform is of "Gough-Stewart type" [START_REF] Merlet | Parallel Robots[END_REF], and the beam kinematics is that of Timoshenko-Reissner [START_REF] Reissner | On a One-Dimensional Large Displacement Finite-Strain Theory[END_REF][START_REF] Timoshenko | Théorie Elémentaire et Problèmes[END_REF]. Another interesting case consists in imposing 2 3 0 Γ = Γ = , in this case the intervertebral kinematics allows a three d.o.f. rotation and stretching along the backbone. The corresponding continuous model is that of an extensible Kirchoff beam [START_REF] Boyer | Finite Element of Slender Beams in Finite Transformations -a Geometrically Exact Approach[END_REF]. The next interesting kinematics is deduced by imposing 1 1 Γ = . In this case we have 1 / r X t ∂ ∂ = (where r is the position field of the beam axis (see Figure 1)), and the Kirchoff beam is non-stretchable while the corresponding parallel platform is of "spherical joint type". Finally for the purpose of three-dimensional manipulation or locomotion, the minimal kinematics consists in imposing 1 0 K = , in this case the robot cannot twist around its backbone, nevertheless it can roll around it, by combining the two bending curvatures 2 K and 3 K . In this case, the platform is of "universal joint type". If we go further, by imposing 2 0 K = the robot is a planar one as those designed for planar swimming in [START_REF] Mcisaac | A Geometric Approach to Anguilliform Locomotion Modelling of an Underwater Eel Robot[END_REF] Finally, note that the constraints imposed by robot architecture design can be interpreted as specific time evolutions, and that from a dynamics point of view, these constraints will induce internal reaction forces and torque fields playing the role of Lagrange multipliers. In the same manner, time evolution of the internal d.o.f. will be imposed by the corresponding internal control torques and forces. This point will be discussed at the end of the paper taking case of our eel-like robot.

III. BASIC DEFINITIONS AND NOTATIONS

First and foremost, note that throughout the paper, any tensor fields depend on time and the space variable X. Moreover, time dependency can be explicit when the fields are known through their time evolution, or implicit, if their time evolutions require us to integrate the dynamics. In the first case, time appears explicitly as an argument of the field, but not in the second case. Finally the derivative operators " ./ X ∂ ∂ " and " ./ t ∂ ∂ " will be indicated by a "prime" and a "dot" respectively. In this section, we will give all the geometric and inertia characteristics of the macroscopic model of the prototype. Firstly we will introduce the following basic definitions from the continuous approach of the geometrically exact beam theory.

A. Basic definitions

In accordance with Figure 1, we will use the following frames:

• The Galilean frame (fixed to the earth), written as

1 2 3 ( , , , ) g F O E E E = ,
• The mobile frame of the section (vertebra) X of the eel's body, indicated as ( ( ), ( ), ( ), ( )

) X F G X t X t X t X = ,
• The mobile frame of the eel's head, which is written as

1 2 3 ( , , , ) o o o o o F O t t t =
, (it is the mobile frame of the section 0 X = of the body).

With these three sets of frames we define the set of curves, which defines the configuration space of the eel (summation convention on repeated indexes is adopted):

3 [0,1] (3)
( ( ) ( ) , ( ) ( ) ) i i i i SO X r X r X E R X t X E → × = = ⊗ R (4) ,
where " ⊗ " is the tensor product.

In the following we will note more simply:

( 0) o R X R = = , ( 0) o r X r = =
(5) , which describe the orientation and position of the head in relation to the earth.

B. Geometric and inertial characteristics of the eel-like robot

We suppose that a continuous beam of rigid elliptic shape cross sections represents the prototype. The mass distribution is supposed to be homogeneous and of unit density. Using these assumptions, the center of mass of any section coincides with the geometric center of the ellipse. For the same reasons, the inertia principal axes of the inertia tensor density along the beam coincide with the ellipse axes. Thus, taking the mobile vectors basis 2 t and 3 t aligned with the small and great ellipse axis respectively, forces the inertia density to be as follows:

1 1 1 2 2 2 3 3 3 / dI dX J t t J t t J t t ρ ρ ρ = ⊗ + ⊗ + ⊗ (6) ,
where ρ is the mass per unit of volume of the material and , 1, 2,3 i J i = are the geometric second order moments of the section around , 1, 2,3 i t i = . Finally, with these assumptions, the centers of mass and buoyancy are collocated and the whole robot is neutrally buoyant. This means that the simulated prototype is assumed to be able to control its density to allow neutral buoyancy at all depths.

C. Inter-vertebral kinematics

Based on an analysis of three-dimensional swimming of fish, we adopted a spherical inter-vertebral kinematics. With this choice, the corresponding beam theory is that of nonstretchable Kirchoff beams [START_REF] Boyer | Finite Element of Slender Beams in Finite Transformations -a Geometrically Exact Approach[END_REF], while the parallel platforms are of "spherical joint type". Thus, the control inputs of the dynamics problem is the time evolution of the twistcurvature tensor field:

, ( , ) ( / ) T T d X K X t R R X R R′ ∀ = ∂ ∂ = (7) ,
In the following, for any 3 V ∈ R , V will denote the skewsymmetric tensor s.t.:

3 , U V U V U ∀ ∈ = ×

R

, and conversely for a skew-symmetric tensor ξ of so( 3), ξ ∨ will denote the associated pseudo-vector of

3 s.t. V V ξ ξ ∨ = × .
Finally, the "spherical kinematics" imposes the mechanical design constraint:

1 / r X r t ′ ∂ ∂ = = (8) ,
which forces the beam to verify the Kirchoff and nonextensibility assumptions.

IV. FLUID-STRUCTURE INTERACTION MODEL

To model the contact forces between the fluid and the eel's body, we will adopt a standard generalization of the Morison model [START_REF] Morison | The Force Exerted by Surface Waves on Piles[END_REF] applied to the case of under-water flexible cables of circular cross section [START_REF] Burgess | Bending Stiffness in a Simulation of Undersea Cable Deployment[END_REF]. This model takes into account not only the transverse reactive (inertial) forces like Lighthill's model [START_REF] Lighthill | Mathematical Biofluid-dynamics[END_REF], but also the resistive (viscous) tangential ones. However, it excludes the influence of the vortexes shaded in the eel's wake. In order to apply this model to the case of the eel, we added a slight modification to it that takes into account the elliptic profile of the eel's sections (see Appendix). In this model, based on the "slender body approximation" of Fluid Mechanics [START_REF] Katz | Low Speed Aerodynamics[END_REF], the fluid forces are modeled "slice-by-slice" through the field of (left invariant) wrench density along the eel's backbone ( * (3) so is the dual space of the Lie algebra of SO(3), i.e. the space of torques):

3 * ]0,1[ (3) so → × R , ( / , / )
ext ext X df dX dc dX , with:

3 3 1 1 1 1 2 2 / ext ld ldi n i i lmi i i i i df dX C V V t C v V t C t = = = - - - ϒ ∑ ∑ (9) ,
where we introduced the velocity and acceleration of each section mass center in the corresponding section frame:

3 1 1 1 i i n i r Vt Vt v = = = + ∑ , 3 1 i i i r t = = ϒ ∑ (10) ,
with n v , the section velocity normal to the eel axis, and ldi C and lmi C , are coefficients depending on the mass per unit of volume of the fluid, the shape and size of the profile (here elliptic) and the Reynolds number of the moving profile in the fluid. At this point, it is worth noting that the local model [START_REF] Simo | A Finite Strain Beam Formulation. The Three-Dimensional Dynamic Problem. Part I: Formulation and Optimal Parametrization[END_REF] generates some drag and lift resultant forces on the global eel dynamics. For instance, let us consider the simplified case where the eel is straight (rigidified) and fixed at rest in a steady flow of velocity v with an angle of attack α (see Figure 2). In these conditions ( 1 , 0 r vE r = -= ), the resultant of local forces (9) turns out to be after computations: where d C and l C are respectively two α dependent drag and lift coefficients, expressed as:

1 2 2 1 3 0 ( / ) ext d l F df dX dX C v E C v E = = - - ∫
1 2 2 1 3 0 ( c c s s ) d l d l d C C C d X = - + ∫ α α α α , 1 1 3 0 ( c s )s c l l d l d C C C dX = - ∫ α α α α .
Note that these expressions of global drag and lift coefficients are approximations within the slender body theory, which do not take into account the influence of the eel's wake. However contrarily to the Lighthill model [START_REF] Lighthill | Mathematical Biofluid-dynamics[END_REF], they take into account tangential viscous resistance as it is shown by imposing 0

α = in d C .
Finally, the acceleration terms of (9) are the added mass forces of [START_REF] Burgess | Bending Stiffness in a Simulation of Undersea Cable Deployment[END_REF] as they are extended to elliptic cross sections in [START_REF] Techet | Fluid Forces on Bodies[END_REF].

As far as the field of fluid torque is concerned, we take into account the drag and added inertia torques generated by the rotation of the planar elliptic cross sections around their normal axis:

1 1 11 1 11 / ext ad am dc dX C t C t = - Ω Ω - Ω (11) ,
where we introduced the angular velocity and acceleration of the section in its mobile frame:

, i i i i t t ω ω = Ω =Ω (12) .
We can summarize all these terms by the following contact law:

/ / 0 . / / 0 drag ext lm drag ext am df dX df dX c r dc dX dc dX c ω        = -                (13) ,
where we introduced the density of drag wrench, which generates some global drag and lift (see the previous simplified example):

1

E 3 E 2 E 1 E 3 E 2 2 E t = 1 t 3 t v α cross section eel 3 1 1 1 1 2 1 1 11 / / ld ldi n i i drag i drag ad C V Vt C v Vt df dX dc dX C t =   - -     =         - Ω Ω   ∑ (14) ,
and also the sectional added mass tensors depending on the geometry of the section profile and the fluid characteristics:

2 2 2 3 3 3 lm lm lm c C t t C t t = ⊗ + ⊗ , 1 1 1 am am c C t t = ⊗ (15)
,

where with an X dependent profile of sections (as in our case), the parameters of (15) will vary along the eel's backbone. It is worth noting that as in the case of lift and drag resultant forces, the added mass contributions of ( 13) will be superimposed slice-by-slice. For instance, if we consider the case where the eel is a rigid ellipsoid (where the lengths of its three axis supported by 1 t , 2 t and 3 t are respectively 2A, 2B, 2C), then the resultant of all the sliceby-slice added mass contributions when calculated at the ellipsoid center mass G, can be broken down as the following 6 6 × matrix in the eel's frame centered on G:

ˆ( ) 0 ˆ0 ( ) T ai a a ai a a G diag m m ms diag I ms I     =         ,
where with standard expressions of . Finally note that these expressions are nothing else than those of the added masses and inertia of a general ellipsoid [START_REF] Susbielles | Vagues et Ouvrages Pétroliers en Mer[END_REF] where / B A and / 0 C A → , i.e. within the "slender body theory". In the case of the eel, this approximation essentially models the reactive transverse fluid forces on which undulatory swimming is based [START_REF] Lighthill | Mathematical Biofluid-dynamics[END_REF].

In order to complete this approximated model let us now consider the forces applied to the two ends of the eel, i.e. its nose and its tail. Because the nose, i.e. section 0 X = is geometrically reduced to a material point, we only take into account some dragging and added mass forces onto the head aligned with its axis, i.e.:

, , 1 1 1
(0) (0) depend on the shape of the head. Finally, because the wake influence is neglected, the wrench applied to the terminal section 1 X = is assumed to be zero.

d head d head o f C V V t = - , , , 1 1 (0)
m head m head o f C t = - ϒ (16)

V. CONTINUOUS KINEMATIC MODEL OF THE EEL

In this section, we develop the continuous version of the kinematic models of open robotic chains. All these models can be interpreted as p.d.e.s in space and time. Nevertheless, when we deal with the dynamic algorithm, these p.d.e.s will be interpreted as spatial o.d.e.s and numerically integrated at each step in a time loop.

A. Continuous geometric model

From ( 7) and ( 8) we obtain the two following models.

Model of orientations:

ˆ , (0) [START_REF] Featherstone | The Calculation of Robot Dynamics Using Articulated-Body Inertias[END_REF] . Model of position:

d o R RK R R ′ = = (
1 , (0) o r t r r ′ = = (18)
.

B. Continuous Kinematic model of velocity

Model of angular velocity:

Time-differentiating [START_REF] Featherstone | The Calculation of Robot Dynamics Using Articulated-Body Inertias[END_REF] and introducing the field of angular velocities:

ˆ:[0,1] (3)
( ) ( ) ( ) T so X R X R X X ω ω → = (19) ,
we find (see [START_REF] Cardona | A Beam Finite Element Non-Linear Theory with Finite Rotations[END_REF] for more details about the demonstration):

ˆˆ[ , ] k k ω ω ′ = + (20) , with ˆˆT d k RK R =
and [.,.] , the Lie bracket of so( 3). Finally, we can rewrite [START_REF] Burgess | Bending Stiffness in a Simulation of Undersea Cable Deployment[END_REF] in terms of axial vectors as:

d k k RK ω ω ′ = + × = (21) ,
with the boundary condition given by the eel's head angular velocity ( 0) ( )

T o o o X RR ω ω ∨ = = =
.

Model of linear velocity:

By time-differentiating (18):

1 r t ω ′ = × (22) ,
with the boundary condition given by: ( 0) o r X r = = .

C. Continuous kinematic model of acceleration:

Time-differentiating ( 21) and [START_REF] Reissner | On a One-Dimensional Large Displacement Finite-Strain Theory[END_REF], gives the two following models.

Model of angular acceleration:

d d k k k RK RK ω ω ω ω ′ = + × + × = + × , ( 0) ( . . ) T T o o o o o X RR RR ω ω ∨ = = = + (23)
.

Model of linear acceleration:

1 1 ( ) r t t ω ω ω ′ = × + × × , ( 0) o r X r = = (24) .
Let us note that ( 23) and ( 24) are two p.d.e.s whose space integration at each step of time will give the acceleration field along the eel's backbone. Thus, they represent two o.d.e.s, which are the continuous counterparts of the forward recurrences of acceleration of a multibody system, where X replaces the body index. Nevertheless, contrary to the case of manipulators, ( 23) and ( 24) are not initialized by some imposed boundary conditions but rather by the acceleration of the eel's head. Hence, we have to derive the head's dynamics in order to compute (at each step in time) the head's acceleration (see section VI).

D. Decomposition of acceleration

As we shall see in the next section, the computation of the head's dynamics requires us to break down acceleration along the eel's body into two components, one depending on the head's acceleration, the other containing the terms relative to velocity and internal strain (curvature and twist) acceleration. Such decomposition is based on the following parametrization of the eel's configuration:

( ) ( ) o e R X R R X = (25) , ( ) ( ) o e r X r r X = + (26) ,
where ( ) e R X and ( ) e r X are the orientation and position of the X section respectively, in relation to the head frame. Finally, the spatial differentiation of ( 25) and [START_REF] Morison | The Force Exerted by Surface Waves on Piles[END_REF] shows that e R and e r are, at any time, the solutions of ( 17) and [START_REF] Katz | Low Speed Aerodynamics[END_REF] , with ( 0) 1 R X = = and ( 0) 0 r X = = .

Decomposition of angular acceleration

If we insert ( 25) into ( 19), we find, after time differentiation: Then, differentiating [START_REF] Cardona | A Beam Finite Element Non-Linear Theory with Finite Rotations[END_REF] with respect to the space variable, and taking (24) into account with the fact that 1 e r r t ′ ′ = = , we find:

( ) ( ) o e X X ω ω ω = + (27)
1 1 ( )) e e t t
γ ω ω ω

′ = × + × × (30) .
Thus, e γ is the solution of ( 24) with the boundary condition ( 0) 0 r X = = , and with ω replaced by e ω .

Finally, the decomposition of acceleration ( 27) and ( 29) can be summarized by the following matrix relation:

γ ω ω ω -        = +                (31) .
Hence, from [START_REF] Meirovitch | Dynamics and Control of Structures[END_REF] we see that e γ and e ω are two acceleration terms containing relative acceleration due to the deformation of the body and all the Coriolis-centrifugal terms.

VI. THE EEL'S HEAD'S DYNAMICS

In this section we adopt the following definition of the eel configuration space: , the shape space of twist-curvatures along the beam. The dynamics of the head are easily derived from the principle of virtual works:

: (3) { : [0,1] (3)
, , δ are the virtual angular and linear displacement fields applied to the cross sections while the time is maintained. The left-hand side of (32) represents the virtual work of acceleration quantities while its right counterpart is the virtual work of external forces here due only to the contact with fluid. Note that the control torque field does not produce any virtual work, since at this point, the internal deformations are directly imposed through the time strain law [START_REF] Chirikjian | The Kinematics of Hyper-Redundant Robot Locomotion[END_REF]. Applying this principle with [START_REF] Meirovitch | Dynamics and Control of Structures[END_REF] and the following form of the virtual displacement field:

dyn ext W W r δ δ δ φ δ = ∀ (32)
1 . 0 1 o e o r r r δ δ δφ δφ -       =            (33) ,
which is compatible with (31) and with the "fixed time constraint" imposed by the application of the principle [START_REF] Goldstein | Classical Mechanics[END_REF]. Finally, we obtain the eel's head's dynamics:

ω         +                   e e a d l e e a d l f f f c c c         + - =                   (34)
, which correspond to the Newton-Euler equations of the whole robot driven by the forced deformations. In [START_REF] Müller | How Body Contributes to the Wake in Undulatory Fish Swimming: Flow Fields of Swimming Eel (Anguilla Anguilla)[END_REF] we introduced the following definitions:

the inertia tensor:

 -   =     -     ∫ , 0 0 0 lm head c   +     (37) ,
the wrench produced by added mass submitted to Coriolis, centrifugal, and strain acceleration: γ ω

     =           ∫ (38) ,
the wrench of drag and lift forces generated by the density of local drag wrench:

1 , 0 / 1 0 ˆ/ 1 0 drag dl d head drag dl e df dX f f dX dc dX c r         = +                ∫ (39) .

Remark: "body deformation dependency"

Note that ( 35)-( 39) take into account the influence of body deformation. For instance, transforming the local drag wrench from the X cross section to the eel's head ( 0 X = ) produces the term under the integral of (39). Hence, (39) first depends on body deformation through the relative position vector e o r r r = -. Moreover, ( 35)-(39) are intrinsic equations, i.e. are not projected in any frame. For instance, projecting (39) in the eel's head frame (this is what is done for numerical purpose), we have:

dX r R R c         = +                   ∫
, where X T denotes the matrix of the components of the tensor T in the mobile base of the X cross section. In particular, o T is the matrix of the components of T in the mobile basis of the section X=0, i.e. in the vector basis of the eel's head. Hence, the resultants (35)-(39) also depend on body deformation through the orientation field e R . Finally, it can be noted that except the basic assumptions about the local contact model, no simplification is made in (35)-(39). These are the characteristics of the geometric exact theory. For instance, plunging the eel into the void and imposing some time varying deformations should conserve six "first integrals", i.e. the components of the linear and angular moments in the earth frame. This has been verified on the simulator with a precision of 6 10 -.

VII. THE DYNAMICS OF THE EEL'S BODY

In this section we adopt the following definition of the eel configuration space:

3 : {( , ):[0,1] (3)} r R SO = × C R , i.
e. a space of curves in a Lie group. We have previously supposed that the internal strains were instantaneously imposed on the vertebrae through their time evolutions [START_REF] Chirikjian | The Kinematics of Hyper-Redundant Robot Locomotion[END_REF]. In order to obtain the time control torque law required to verify this evolution, we are now going to derive the dynamics of internal forces and momentum. This can be done in a straightforward manner from a variational calculus applied to the following augmented Lagrangian:

1 1 0 1 ( , ) (( ) , ( ) ) 2 T T T T d Ar n L r r t K K d X I M ρ ω ρ ω     ′ = + - -         ∫ (40)
, where the first term stands for the kinetic energy of the eel; the second takes into account the constraints imposed to the inter-vertebral kinematics by the design (8) and the control [START_REF] Chirikjian | The Kinematics of Hyper-Redundant Robot Locomotion[END_REF]. In order to force these constraints, we introduced the following field:

3 * [0,1] (3)
( ( ) , ( ) )

T TT so X n X M X → × R ,
where, n is the field of internal forces in the earth frame and M is the field of control torque in the vertebra (cross section) frames. Then posing the "extended Hamilton principle" [START_REF] Meirovitch | Dynamics and Control of Structures[END_REF]:

, , , ,

f f i i t t ext t t L dt W dt r n M δ δ δ δ φ δ δ = ∀ ∫ ∫ (41) ,
where the external contribution is due to the fluid (13), we obtain the partial differential equations of the eel's body which represent the Newton-Euler equations of each section:

]0,1[ X ∀ ∈ :

drag lm drag am df dX A c r n dc dX I c I RM t n ρ ρ ω ω ρ ω ′   -    + =      ′ - + × + ×       1 / ( ) / ( ) ( ) ( )
, which has to be completed with the constraints:

, ( ) ( , ) d X K X K X t ∀ = , 1 0 r t ′ -= (43) ,
and the boundary conditions on the head given by ( 16):

, ,

(0) 0 d head lm head o n f c r M +     =         (44) , (0)
and on the tail:

(1) 0 (1) 0

n M     =         (45) .
Finally, before detailing our algorithm, let us note that the set of equations ( 42)-(45) defines a closed formulation enabling us to solve the direct dynamic problem of the eel, i.e. to compute the motion of the sections based on the knowledge of the internal torque law. This could be achieved for instance through the geometrically exact approach of finite elements as proposed in [START_REF] Simo | On the Dynamics in Space of Rods Undergoing Large Motions -A Geometrically Exact Approach[END_REF]. Instead, this set of equations will be used to solve the inverse problem, consisting in computing the internal control torque law corresponding to some expected deformation. In this case, we shall use (42) as spatial o.d.e.s integrated at each step of a simulation time loop. Let us note that in such a case, only the equations (42) with one of the two boundary conditions (44)-(45) are required since the constraints (43) are implicitly taken into account through continuous kinematics, while the second of the two boundary conditions is taken into account via the eel's head dynamics.

VIII. COMPUTATIONAL ALGORITHM

Based on the previous results, we can derive an algorithm whose inputs are the state variables of the eel's head, and the desired time evolution of its body's curvature and twist (and their time derivatives), i.e.: ( ,,, )

o o o o r r R ω , and: , ( ( , ), ( , ), ( , )
)

d d d X t K X t K X t K X t ∀ (46) .
The outputs of the algorithm are head acceleration and the required control torque for imposing (46), i.e.:

( ( ), ( ), ( , ) )

o o d t r t t M X t X ω ∀ (47) .
The algorithm is based on three integration numerical loops, the first one is in time and allows us to update the state of the head, and the two others are in space and are included in the first. The first space integration loop computes head acceleration. The second space integration loop calculates the control torque distribution applied to the vertebrae.

A. First space integration loop

The first space loop starts with the following spatial o.d.e.'s deduced from the p.d.e.'s of section V, now considered at the fixed current time t of the time loop.

Computation of the configuration

• In orientation: ˆd R RK ′ = , (0) o R R = , • in position: 1 , (0) o r RE r r ′ = = (49) .

Computation of the velocities

• In orientation:

d RK ω′ = , (0) o ω ω
= (50) , • in position: r r ω ′ ′ = × , (0) o r r = (51) .

Computation of the acceleration:

• In orientation:

Computation of the head dynamics:

Let us note that the tensors and wrenches of (34) can be computed at each step of the time loop by integrating the following system of o.d.e.'s with respect to the space variable from the head to the tail: And where (54) and ( 56) are initialized by the following boundary conditions compatible with (44) and ( 16), with

      ′ = +               I 1 ( 1 ) ˆˆ( 1 ) ( 1 )
      ′ = -               F 1 0 ( ) ˆ1 ( ) ( ) lm e e a m e A c r I c I ρ γ ρ ω ω ρ ω +    =    + + ×    (55)
, , 1 1 m head m head o o c C t t = ⊗ : , 0 (0) 0 0 m head c   =     I , , (0)
0

d head dl f   =     F (57) .
While:

0 (0) 0 e   =     F (58) .
For simulation purpose, (54)-( 56) are projected into the head frame. Finally, once the equation (48) replaced by its parameterization in terms of a quaternion field (noted q ), all the equations from (48) to (57) can be solved by integrating the following first order system: ( , ( , ), ( , ), ( , )) ( ,,,,,,( ) ,, )

d d d x f x K X t K X t K X t ′ = , 1 1 ( 0)
T T T T T T T T T T e e e d l
x X r q r vect ω γ ω = I F F , and:

1 ( , , , , 0,0, ( (0)) ,0, (0))

T T T T T T T o o o o o d l x r q r vect ω = I F (

B. Second space integration loop

From the following o.d.e.'s: Next integrating (63) in relation to X at each step time, gives the field of internal forces applied to the vertebrae. At this step, inserting n and ω in the second p.d.e. of (42) and integrating : Finally, once completed by ( 48)-( 51), ( 61)-( 64) are solved by integrating the following first order system: ( , ( , ), ( , ), ( , ))

d d RK R K ω ω ′ = + × , (0)
(( ) ( )) T d a m M K M R I c I ρ ω ω ρ ω ′ + × =- - + × 1 ( / ) T drag R t n dc dX - × - (64)
d d d x f x K X t K X t K X t ′ =
(65) , where: ,,,,,, )

2 ( ) ( ,
T T T T T T T T T
x X r q r r n M ω ω = , and: 2 ( , , , , , , (0),0)

T T T T T T T T o o o o o o o
x r q r r n ω ω = (66) , with (0) n , given by (44). Once the second space loop is completed, the head acceleration is time-integrated twice over in order to update the state of the head. Then the time is shifted by one step and the algorithm resumes…

Remark:

Let us note that (63) and (64) play the role of the backward recursive equations of inter-body wrenches of a serial manipulator [START_REF] Luh | On-Line Computational Scheme for Mechanical Manipulator[END_REF]. Nevertheless, contrary to the case of a manipulator, where the boundary conditions on wrenches are known at the end of the structure and are unknown at the base, in our case we know the boundary wrenches at both ends, so forward or backward integration is equivalent. In the following, we choose to adopt a forward recursive equation on wrenches initialized by (44). It follows that the second boundary condition (45) will be a verification test.

IX. NUMERICAL RESULTS

In this section, we report some numerical results obtained by the algorithm of section VIII. The objective of these examples is to prove that despite the non-linear character of the eel dynamics we can easily realise three-dimensional gaits by combining some elementary curvature laws, while maintaining the twist at zero. Before presenting the numerical results, let us describe the robot characteristics.

A. Geometric description of the eel-like prototype:

In order to take into account the continuous property of the skin of the prototype, we consider the previous macroscopic modeling with the following geometry (see Figure 3):

where the total length of the robot is of 2 L = meters, its total mass is 21,7 kg, its density is 1. The cross section is of elliptic shape (of great axis maximal length equal to 15 cm) on the interval 1 3

[ , ] X X , while the head (

1 [ , ] o X X X ∈
) and the tail (

3 [ , ] X X L ∈
) are respectively half of a sphere (of a diameter of 15 cm), and that of an ellipsoid (of great and small axis length equal to 10 cm and 5 cm, respectively). Finally, denoting the small and great axis' lengths of the elliptic section by b and c, the model ( 9)-( 13), is used with:

1 1 (1/ 2) ( ) / 2 ld C C b c ρ π = + , 2 2 (1 / 2) 2 ld C C c ρ = , 3 3 (1 / 2) 2 ld C C b ρ = , 2 2 2 1 4
(1 / 2) ( )

B. First example: Nominal planar propulsion:

Following the standard uses of biomechanics literature about anguilliform locomotion [START_REF] Gillis | Environmental Effects on Undulatory Locomotion in the American Eel Anguilla Rostrata: Kinematics in Water and on Land[END_REF][START_REF] Müller | How Body Contributes to the Wake in Undulatory Fish Swimming: Flow Fields of Swimming Eel (Anguilla Anguilla)[END_REF], we started our numerical investigations with planar forward propulsion. This planar motion is produced by a curvature law of the following form:

1 2 3 0, ( ) ( ) K K K X ht α = = = , ( 67) ,
with α given in terms of X, by:

2 1 2 3 X X α α α α = + + ,
where the i α 's are arbitrary constant coefficients. The time varying law h is defined by: 1

( , ) ( ) ( , ) o h X t h t h X t = ,
where o h is a starting function that allows the eel to start smoothly from rest to nominal propulsion mode: for:

( ) 1 , for:

o i o i f i o if i i f o f h t
= <  = - -   - - -   ≤ ≤   = >  (68) .
This term multiplies the propulsion term, inspired by the experimental observation of the animal [START_REF] Gillis | Environmental Effects on Undulatory Locomotion in the American Eel Anguilla Rostrata: Kinematics in Water and on Land[END_REF]:

1 ( , ) sin 2 X t h X t T π λ     = -         (69) .
This second term represents the propagation of curvature waves from the head to the tail with a constant time frequency 1/ T and wavelength λ . Simulation is achieved with the following numerical values:

1 1 α = , 2 0.75 α = - , 3 1 α = , 0 i t s = , 1 f t s = , 1 m λ = , 1 T s = .
In Figure 4 we report several configurations of the eel in the plane obtained every 1.25 sec., while Figure 5 represents the time evolution of the eel's head twist in the earth frame.

C. Second example: Plane turning law:

In order to make the eel turn in the plane (Figure 6), we add a constant offset to the previous propulsion term. This is the continuous version of the joint law adopted in [START_REF] Mcisaac | A Geometric Approach to Anguilliform Locomotion Modelling of an Underwater Eel Robot[END_REF]. The constant curvature offset is imposed progressively through the starting law (68). Hence, the curvature law is now:

1 2 3 3 , 0, ( ) ( ) ( ) of o K K K X ht K h t α = = = + .

D. Third example: Submergence

The goal of this example is to achieve submergence from one given altitude to another (see Figure 7). This is accomplished by adding to the propulsive law of the first test, a curvature law 2 ( ) t K t around the second axis of the eel section frames. Hence the resultant curvature law is:

1 2 2 3 0, ( ), ( ) ( ) K K K t K X h t α = = = , with: 2 1 2 1 2 2 1 2 2 2 3 2 2 3 2 2 3
( ) 0 , for:

( ) ( , , 0, ) , for:

( ) ( , , , 0) , for:
( ) 0 , for: with first and second order derivatives equal to zero, while guaranteeing a second order time continuity.

E. Fourth example: Spiral

Lastly, we surmised a three-dimensional plunging following a spiral. This is achieved using the following curvature law, which superimposes the turning curvature law of example 2 with, a 2 K 's time evolution given by: ( )

of o K K h t =
, where o h is the starting time function. The Figure 9 give the time evolutions of the three components of the control torque field M evaluated at / 2 1 X L m = =

.

Figure 9: Time evolution of the internal torque for the section at X=1m in its mobile frame.

Finally, Figure 10 validates the whole dynamic balance of the eel. The computation of internal force n and torque M is computed all over the back-bone (here at 10 t s =

) using a forward space integration and we obtain ( ) ( ) 0 n L M L = = as pointed out in the remark of section VIII.

X. CONCLUSION

This paper presents a dynamic model of a swimming eellike robot. Contrary to the prior results about the biomimetic robotics of the eel, the proposed solution to the dynamics is capable of modeling an ideal three-dimensional continuous prototype. The design of the prototype will be based on the serial assembly of many parallel platforms of "spherical joint type", and the multi-body system will be covered with a continuous flexible organ copying the eel's skin. The continuous approach has the advantage of providing a macroscopic model of the prototype, and in particular giving an overview of swimming dynamics without entering into the details of the modeling of complex internal hybrid (parallel and serial) kinematics and their dynamics. This advantage has already been tested in the framework of our project. Macro-continuous modeling allows us to rapidly investigate the locomotion and the control of the ideal prototype. Moreover it is extensively used today to assist in designing the future prototype. In order to get an exact model (from kinematics to dynamics), we based our development on the "geometrically exact theory" of nonlinear beams, a new promising paradigm of non-linear structural mechanics. Based on the assumption of Cosserat medium, this theory gives an exact model of finite rotations of some rigid "micro-solids" right from the beginning of the analysis. In our case, these micro-solids are the beam cross sections copying the vertebrae of the fish. This choice, allows us to obtain a three dimensional dynamic model which will be used in the future to control three-dimensional swimming, a problem to our knowledge never investigated in robotics. Moreover, based on the literature of fluid mechanics, a simplified, but quite complete model of the fluid-structure contact is adopted. The geometrically exact approach is developed in the Newton-Euler formalism, as it is well known among the robotics dynamics community. Based on these modeling choices, the proposed algorithm allows one to compute the motion of the eel and the control torque distribution from the knowledge of the desired internal deformation imposed to its body. Finally, this algorithm constitutes a generalization of the computed torque approach of articulated manipulators. The generalization here concerns the continuous character of the hyper-redundant robot's model and also the mobility of the base, which plays the role of the eel's head.

Figure 1 :

 1 Figure 1: Frames and parameterisation.

Figure 2 :

 2 Figure 2: Eel in a straight rigidified configuration fixed in a steady flow.

, where oω

  is the angular acceleration of the head and e ω is the solution of the following p.d.e. obtained by spatial differentiation of (27): Decomposition of linear acceleration Now inserting (26) into the definition of r leads us to the introduction of a new acceleration term e γ defined by:

  the added mass tensor induced by the head acceleration:

  where I denotes inertia and added mass tensor, e -F is the wrench of Coriolis-centrifugal and strain forces produced by material and added mass, and dl F is the wrench generated by the local drag forces and torques applied along the eel.

  . At the end, from this first space loop we obtain the head acceleration ( , ) , which initializes the second space loop.

  fields of Galilean acceleration of the eel's body. Then inserting r into the first p.d.e. of (42) gives:

Figure 3 :

 3 Figure 3: Geometric shape qof the eel-like prototype,

Figure 8

 8 Figure 8 shows the results of this test with: 2, . of K =

Figure 7 :

 7 Figure 7: Eel's configuration obtained every 0.5 sec..

Figure 6 :

 6 Figure 6: Eel's configuration every 1.25 sec..

Figure 4 :

 4 Figure 4: Eel's configuration every 1.25 sec.

Figure 8 :

 8 Figure 8: Eel's configuration every 1sec..

Figure 5 :

 5 Figure 5: Head's velocity in the fixed frame.

Figure 10 :

 10 Figure 10: X profile of the second component of the internal force and torque at 10 t s = (in the section frame).

Acknowledgment:

We would like to thank the French CNRS for supporting this work through the project "Robot Anguille" of the Interdisciplinary Research Program ROBEA. We acknowledge the constructive remarks of the reviewers, which have contributed to improve the quality of this paper.

XI. REFERENCE

APPENDIX

The purpose of this Appendix is to show that in the case of an eel's circular cross section, the basic model ( 9) and [START_REF] Simo | On the Dynamics in Space of Rods Undergoing Large Motions -A Geometrically Exact Approach[END_REF], exactly reduces to the model of fluid forces of circular cross sectional cable of [START_REF] Burgess | Bending Stiffness in a Simulation of Undersea Cable Deployment[END_REF]. As a matter of fact, introducing the coefficients of section IX.A into (9) and (11) (with P the cross section perimeter) first gives:

Then, let us consider a circular cross section by making: : added mass coefficient of the circular cross section, gives:

which is quite simply the model of [START_REF] Burgess | Bending Stiffness in a Simulation of Undersea Cable Deployment[END_REF] with 1 t C C = : the "tangential drag coefficient".


Mathieu Porez

Auke Vincent Lebastard

Jan Ijspeert

Frédéric Boyer

Multi-physics model of an electric fish-like robot : numerical aspects and application to obstacle avoidance

The paper deals with the modeling of a fishlike robot equipped with the electric sense, suited to study sensorimotor loops. The proposed multi-physics model merges a swimming dynamic model of a fish-like robot with an electric model of an embedded electrolocation sensor. Based on a TCP-IP and threaded framework, the resulting simulator works in real time. After presenting the modeling aspects of this work, this article focuses on two numerical studies. In the first, the interactions between body deformations and perception variables are studied and a current correction process is proposed. In the second study, an electric exteroceptive feedback loop based on a direct current measurement method is designed and tested for obstacle avoidance.

I. INTRODUCTION

Lissmann in the 1950's [START_REF] Lissmann | The mechanism of object location in gymnarchus niloticus and similar fish[END_REF] was among the first scientists to clearly demonstrate the electric nature of the perception of the weakly electric fish. He assessed that "the electric organ discharges belong to a full sensorial system and are used for scanning the environment and for the interactions with the other electric fishes". After this discovery, scientists begun to study in detail how the environment was electrically interpreted by the electric fish. Brian Rasnow in 1996 [START_REF] Rasnow | The effects of simple objects on the electric field of apteronotus[END_REF] developed the first model of interactions between the environment and the electric currents flowing through the skin of the fish, that we call the electric image. His model derived from simple electromagnetism conditions is dedicated to the study of the effects of a sphere placed in the vicinity of the fish. His simple model showed the relation between the shape of the electric image and the distance and dimensions of spherical objects. Recently Solberg & al [START_REF] Solberg | Active electrolocation for underwater target localization[END_REF] have designed a robotic detection device based on the electric sense. They performed with their device an automatic detection of a sphere. In its current form, it seems that their device is more suited for the design of perception algorithms rather than for an implementation on an existing autonomous robot. More recently, a new project was started in Europe: ANGELS 1(for ANGuilliform robot with ELectric Sense). Its objective is to build an eel-like robot equipped with an electric sensor. The ANGELS' robot would be capable to navigate using the electric sense and to divide itself in several autonomous individual modules for exploration.

M. POREZ and A. J. IJSPEERT are with BioRob, EPFL -STI -IBI -BIOROB, Station 14, 1015 Lausanne, Switzerland. mathieu.porez@epfl.ch and auke.ijspeert@epfl.ch Vincent LEBASTARD and Frédéric BOYER are with IRCCyN, EMN -La Chantrerie -4, rue Alfred Kastler,B.P. 20722, 44307 Nantes Cedex 3 France. vincent.lebastard@mines-nantes.fr and frederic.boyer@mines-nantes.fr

In the ANGELS' context, the article deals with the actionperception loops. In the case of an electric fish-like robot, this problem consists in deriving the laws that rule: 1 • ) the fish and fluid dynamics; 2 • ) electromagnetic phenomena; 3 • ) sensorimotor feedback; involved by the locomotion, perception and their couplings. Thus, in this paper, we propose a model and a simulator in order to study the action-perception problem.

The article is structured as follows. The modeling of an electric fish-like robot is first presented in section II. In section III, the models are implemented in a modular simulation architecture. Then, the resulting simulator is exploited for obstacle avoidance in section IV in the case of the AmphiBot robot of BioRob Lab2 . Lastly, the article ends with section V by some concluding remarks. The fish-like robot is a two-dimensional swimmer, which is composed of N actuated modules, a head and a caudal fin. At any time t, the robot configuration is defined by the joint positions q together with the orientation matrix R 0 and the position vector P 0 of a mobile frame attached to the robot head F 0 = (O, t 1 , t 2 , t 3 ) with respect to F g . Finally, the robot is equipped by discharge electrodes on the tail and measurement electrodes placed on its head in order to electro-sense its environment. In the following, the electrodes are denoted by ǫ i (where i denotes the electrode index).

II. THE ELECTRIC FISH MODELING

A. The problem statement

On the base of this statement, to study the actionperception interactions, we must model: 1 • ) the swimming dynamics of a self-propelled fish-like robot. 2 • ) the electric currents crossing the robot's electrodes in function of the electric activity of its discharge electrodes and the environment; 3 • ) the sensorimotor feedback coupling the current measurements with the parameters of the anguilliform swimming gaits or transient maneuvers. Let us now detail each of them. In the general case, the swimming problem consists in the following sequence of physical causes and effects: 1 • ) in order to move, any swimming robot actuates its internal (shape) degrees of freedom; 2 • ) then, the new body shape disturbs the media which generates contact forces on the robot's skin; 3 • ) finally, these forces trust the external degrees of freedom of the head module. In this string of causalities : the following coupled dynamics appear: the external dynamics for the locomotion, the media dynamics for the contact forces and the internal dynamics for the locomotion control. To model properly a fish-like robot, we must resolve these coupled dynamics. For this purpose, we propose to exploit the three following recent contributions on this topic: 1 • ) the Central Pattern Generator (CPG) by Crespi & al [START_REF] Crespi | Online optimization of swimming and crawling in an amphibious snake robot[END_REF] for the locomotion control; 2 • ) the recursive algorithms based on the Newton-Euler (N-E)'s equations by Khalil & al [START_REF] Khalil | Dynamic modeling and simulation of 3-d serial eel-like robot[END_REF] for the robot dynamic modeling; 3 • ) the analytical hydrodynamic model of a 3-D self-propelled fish swimming by Boyer & al [START_REF] Boyer | Fast dynamics of an eel-like robot -comparisons with navier-stokes simulations[END_REF], [START_REF] Boyer | Poincaré-cosserat equations for the lighthill three-dimensional large amplitude elongated body theory: application to robotics[END_REF] for the hydrodynamic force modeling. As shown on Fig. 2, the proposed swimming model is composed by: 1 • ) a locomotion controller; 2 • ) a direct dynamic model; 3 • ) a hydrodynamic model.

B. The swimming modeling

0 V0 1 s 2 q d q Γ q, q 0 V 0 P 0 R 0 C.P.G P. D. controller Direct dynamic model A ν k α 1 s 2

Hydro. model

1) The locomotion controller: It is composed of a CPG and a proportional derivative PD controller computing, re-spectively, the joint set points q d and the torques Γ applied by the motors on joints. A simple manner to generate the swimming rhythmic motions of fishes (see [START_REF] Breder | The locomotion of fishes[END_REF]) is to use a CPG, i.e. a system of coupled nonlinear oscillators inspired from the locomotor circuits found in the spinal cord of vertebrates (e.g. the lamprey). This bio-inspired controller has several explicit parameters, which can be continuously modified, controlling the body undulation shape for forward and backward swimming (through the wave amplitude A, the wave frequency ν and the number of waves along the robot backbone k) and the average curvature for turning maneuvers (through the backbone average curvature α). In addition to its reduced parameter set, the CPG can adapt quickly and smoothly the joint trajectories after any abrupt parameter change. In this paper, this last property will be extensively used to smoothly modulate the swimming gaits depending on the electro-sense variables.

2) The direct dynamic model: In this work, the computing of the swimming dynamics of the fish-like robot is achieved by the direct recursive N-E algorithm developed in [START_REF] Khalil | Dynamic modeling and simulation of 3-d serial eel-like robot[END_REF]. The algorithm based on Newton's law and Euler's theorem is dedicated to the dynamic modeling of mobile serial robots. It allows to compute, as a function of Γ, the robot motion, i.e. the robot head accelerations 0 V0 with respect to F g together with the joint accelerations q. At each step of a global time integration loop; the direct algorithm solves one after one the three following recursive loops dedicated to : 1 • ) the robot kinematics; 2 • ) the external dynamics taking into account the hydrodynamic forces F h ; 3 • ) the internal dynamics.

3) The hydrodynamic model: As far as the hydrodynamic is concerned, we used here the generalization of the large amplitude elongated body theory of Lighthill (see [START_REF] Lighthill | Large-amplitude elongated body theory of fish locomotion[END_REF]) proposed in [START_REF] Boyer | Poincaré-cosserat equations for the lighthill three-dimensional large amplitude elongated body theory: application to robotics[END_REF] and numerically validated in [START_REF] Boyer | Fast dynamics of an eel-like robot -comparisons with navier-stokes simulations[END_REF]. This model appears as a superimposition of a reactive model with a resistive one which depends only on segment motions. As regards the reactive part of the hydrodynamic model, the thrust and the lateral lift are modeled through the effect of the fluid inertia on the undulating fish body. Basically, it is based on the slender body theory of Munk [START_REF] Munk | The aerodynamic forces on airship hulls[END_REF], where the 3-D flow around an elongated body is approximated by a stratification of planar lateral flows, which are then resolved analytically. Concerning the resistive part, the effect of the fluid viscosity on the nose and on the skin of the fish is modeling through a Taylor-like resistive model [START_REF] Taylor | Analysis of the swimming of long narrow animals[END_REF].

C. The electric modeling

From a technical view point, the working of the electricperception sensor embedded in the robot is based on the generation of an electric field E in the water by polarization of measurement electrodes with respect to discharge ones. The robot body having a lower conductivity than the water, hence E is focalized through the measurement electrodes. By this mechanism, the robot can build an electric image of its environment by comparison between the expected currents (measured in a free environment) and those actually measured. From the point of view of electromagnetism, the electric state of the fluid can be considered as quasi steady.

Thus the electric field is irrotational and can be determined by the gradient of the electric potential field φ. Then, the Ohm law and the conservation of electric currents allow one to state the set of partial differential equations ∆φ = 0 named Laplace equations, which rule φ ∈ D (hence E = -∇φ) with respect to the electric conditions imposed on ∂D (i.e. the electrodes, robot body, tank and passive objects). In order to resolve the Laplace equations, we propose to use the 2-D Boundary Integral Equations (BIE) formulation and its discretization using the conventional Boundary Element Method (BEM) [START_REF] Gaul | Boundary element methods for engineers and scientists[END_REF]. The solution to the boundary value problem described by the Laplace equations can be written for any internal point x 0 ∈ D by:

φ(x 0 ) = ∂D G(x 0 , x)E + ⊥ (x)dl - ∂D F (x 0 , x)φ(x)dl , (1)
where E(x) + ⊥ = (∂φ/∂n)(x) (with x ∈ ∂D and n the outward normal to ∂D), F = (∂G(x 0 , x)/∂n)(x) and G(x 0 , x) are the Green's function in 2-D space define by:

G(x 0 , x) = 1 2π ln 1 r and F (x 0 , x) = 1 2πr ∂r ∂n
, where r = ||xx 0 ||. In order to numerically solve (1), the BEM consists in meshing ∂D into M boundary elements L l , i.e. ∂D = M l=1 L l (see Fig. 3). Then, the two layers of singularities φ and E + ⊥ are approximated on each elements through a polynomial interpolation of their nodal values. Hence, the electric potential φ l computing in x l , center of the l th element, can be written as the discrete superimposition of the contributions due to the elements m = 1, ...M , i.e. (2) Then, the M equations ( 2) can be written in the form of an implicit linear system, in applying the boundary conditions to φ and E + ⊥ . After the resolution of (2), the currents I i crossing electrodes ǫ i are computed as follows: I i = γη E + ⊥l L l where η is a 2-D/3-D correlation coefficient and L l is the length of the element L l . Let us note that the last sum operation is realized only on constituent meshes of ǫ i . Finally, we realized an experimental validation of our numerical model. The used set-up is described in [START_REF] Jawad | Sensor model for the navigation of underwater vehicles by the electric sense[END_REF]. The differences between the currents given by our 2-D BEM solver and experiments do not exceed 10% after setting η at 16.5.

III. THE SIMULATOR FRAMEWORK

As mentioned above, the swimming and the electric models are included in a modular simulation framework based on a TCP-IP network. The highlight of this simulator is the real-time computing of the global model. This simulator computes at a frequency of 20 Hz (for M ≤ 1000) the current state of the measurement electrodes with respect to the swimming dynamics and the surrounding environment. In accordance with Fig. 4

IV. APPLICATION TO THE AMPHIBOT

A. Electro-AmphiBot Fig. 5. The three modules AmphiBot robot equipped with the electric sense.

To illustrate the proposed multi-physic model, we address the two following problems: 1 • ) the study of the interactions between the body deformations and the perception variables; 2 • ) the design of an electric exteroceptive feedback loop for obstacle avoidance based on a direct current measurement method. In this context, we used the modular robot family of BioRob Lab as a reference. These robots named AmphiBot are fish-like robots composed by a serial assembling of identical modules which have been recently equipped with the electric sense further to a collaboration with IRCCyN Lab 3 and Subatech Lab 4 in the ANGELS' context (see Fig. 5). For this work as shown on the Fig. 6, the used robot has three modules, i.e. one head (m 1 ), two actuated modules (m 2 and m 3 ) and a caudal fin. The robot has an external length of 390 mm, a cross-section of 40 by 57 mm 2 , and its joint axis is placed at 20 mm from the front edge of segments. Moreover, the robot is neutrally buoyant, i.e. its density is equal to that of the fluid. The electric sensor is composed of five hemispherical electrodes with a diameter of 15 mm. They are placed as follows: three on the head distributed on the left ǫ 1 , right ǫ 3 and front ǫ 2 of the module; two on the caudal fin, one on each lateral face ǫ 4 and ǫ 5 . Moreover, in order to measure on each electrode an electric current, we imposed the electric potentials of -5 V to the head (measurement) electrodes and +5 V of the caudal fin (discharge) ones.

B. Interaction action-perception

This subsection deals with the interactions between the measured currents and the body undulations (set by ν, A and k) together with the average curvature (set by α). Until now, all contributions of the electro-perception in robotic (see [START_REF] Solberg | Active electrolocation for underwater target localization[END_REF], [START_REF] Lebastard | Underwater robot navigation around a sphere using electrolocation sense and kalman filter[END_REF]) have been realized on a rigid body, i.e. the geometric configuration of the electric sensor is constant with respect to the time. Thus, the current variations, which are the inputs of perception algorithms and other controllers, are only due to the environment. Nevertheless, in this paper, the sensor is embedded in a fish-like robot. Thus, under the swimming undulation of the body shape, the sensor is deformed. As Jawad & al have shown in [START_REF] Jawad | Sensor model for the navigation of underwater vehicles by the electric sense[END_REF], the currents crossing the sensor are strongly dependent of the interdistances between electrodes (see ( 8) of [START_REF] Jawad | Sensor model for the navigation of underwater vehicles by the electric sense[END_REF]), which vary in the swimming case. Hence, the geometric configuration of the sensor disturbs the measured currents which in this case is significant and hides the environment effects. In order to cancel these unwanted effects, in the following, we will propose a correction on the measured currents.

To observe the effect of the undulation on the measured currents, we carry out the following numerical test: the robot swims in a straight line towards a wall (we fix ν = 0.8Hz, A = 25 • , k = 0.5 and α = 0 • ) in a fluid with a conductivity of 0.04 S/m . As illustrated in Fig. 7, we place on the robot path two insulated objects (denoted by A and B) which are less conductive than the water. Hence, when the robot swims near to one of them, in accordance with the Ohm law, the local media resistivity increases and the currents through the electrodes decrease. As an illustration of the electrolocation, on Fig. 8, we drew the electric field around the object A. Fig. 9 shows the evolution of I i for i = 1, 2, 3 with respect to the distance d between the robot nose and the wall. We observe that I 2 is disturbed by all the objects. The effect of the wall on I 2 is significant for d ≤ 0.2 m. Near the objects A and B, I 2 decreases on average of 8%. On the other hand, the lateral current I 1 (respectively I 3 ) decreases on average of 15 % (respectively 10 %) under the resistive effect of B (respectively A) but it is not affected by the object A (respectively B) and the wall. In view of these first observations, we can deduce that ǫ 2 detects objects in front of the robot while ǫ 1 and ǫ 3 detect objects on the left or right respectively. As far as the current oscillations are concerned, they are due to body swimming deformations. Their effects are of the order of 7% which can conceal environment effects and reduce the performances of possible perception or control algorithms. Thus, in order to minimize the undulation effects on measurements, a solution will be to work with the mean currents Īi (t) = t t-1/ν I i (t)dt (see dot lines in Fig. 9). As the undulation effect has been highlighted, let us look at to the curvature disturbances. To illustrate them, we carry out a simulation where the robot turns within a corner. As illustrated in Fig. 10, the (preset) motion path of the robot is as follows: 1 • ) the robot goes forward to the point C; 2 • ) between C and E, it turns to the right with a constant curvature α = 20 • ; 3 • ) after E, it goes forward again. In Fig. 11, we have plotted I i and Īi measured on the robot head during the test described above. In view of the results obtained in the straight line test, we could expect that Ī3 is constant during the turning test, but this is not the case. In fact, the disturbance observed on ǫ 3 is due to the body curvature, which moves electrodes closer to each other. To characterize this interaction, we carried out a study in infinite space (i.e. a free environment), where we measured currents denoted by Īi∞ for different values of the curvature α (see Fig. 12). In order to remove the contribution of the curvature from Īi , for any swimming motion, we propose the following correction: I c i = Īi -Īi∞ (α)+ Īi∞ (0), where I c i is the corrected measured current through ǫ i . Thus, in applying this correction to ǫ 3 for the turning test, it is possible to reduce the curvature effect. As shown in Fig. 13, the correction is efficient during the turning phase: I c 3 is approximatively constant between C and E. Let us note that the peaks observed in Fig. 13 are due to the time delay introduced by computing the average. Based on these results, we designed a sensorimotor feedback loop dedicated to the obstacle avoidance.

C. Application to the obstacle avoidance

The proposed electric exteroceptive feedback does not use electric models of the environment but simple mathematic operations on I c i . In studying the signs of these variables, the low-level-perception algorithm can detect an insulating obstacle in the surrounding of the robot, on the left, on the right or in front of it and find a free area to escape. The proposed feedback loop is based on the following variables:

c f = β f (I c 1 -I c
3 ), and c h = β h (I c 2 + Ī2∞ (0)), where, c f is the left-right difference between measured currents on the robot's head, whereas c h indicates the presence or not of an insulated object being in front of the robot, β h and β f are sensibility gains. Then, the idea consists, as function of c f and c h , to drive the average body curvature in order that the robot turns towards a free space, by using α = β α atan2(c f , c h ), where β α is a sensibility gain. In order to illustrate the working of the proposed lowperception algorithm, the first example is dedicated to corner avoidance. In this example, we present the different outputs of the perception and locomotion algorithms. Fig. 14 shows the trajectory of the robot in the tank. This trajectory produced by the obstacle avoidance control law is composed of three parts :1 • ) start of the turn (F); 2 • ) effect of the corner (G); 3 • ) end of the turn (H). Fig. 15 shows the mean currents and the internal variables of the control law. Fig. 15 shows the time evolution of α, q 1 and q 2 . More precisely, on the start of turn, I 3 decreases, thus c f and α decrease while c h increases. Then, the robot under the effect of the feedback control law bends on the left. As regards the effect of the corner, the robot being far from the corner, c h , c f and α are close to zero. The robot goes in strait line to the next wall. The wall being an insulated object and being on the right, the robot bends to the left. At the end of the turn, the robot goes in strait line to the next obstacle. Īi , c h , c f , α, q 1 and q 2 with respect to the time t for the avoidance of a corner. In the second and last examples, two tests have been done. In the first one, the robot moves in an empty tank, while for the second one, the robot avoids an object in the same tank. Fig. 16 shows 5 trajectories of the robot in the tank with and without an insulating object. The test is stopped when the head of the robot goes out of the tank which happens when an obstacle is detected too late.

V. CONCLUSIONS

In this paper, we have presented some results regarding sensorimotor feedback in the case of an electric fish-like robot. For this work, we have designed a multi-physics simulator working in real time. This simulator solves in parallel the swimming locomotion problem and the electrokinetic equations. Using this tool, we studied the interactions between the currents crossing the electrodes and the body shape. On this base, we proposed a current correction process in order to cancel the unwanted effects of the anguilliform swimming on the perception. Thanks to these results, we designed a low-level-perception algorithm to address the obstacle avoidance problem of an electric fish-like robot. With the proposed controller based only on 3 measurements, the robot can swim safely for a quite long duration.

VI. ACKNOWLEDGMENTS

The ANGELS project is funded by the European Commission, Information Society and Media, Future and Emerging Technologies, contract number: 231845. Finally, the authors gratefully acknowledge Pol-Bernard Gossiaux for rewarding exchanges on electromagnetism.

Fig. 1 .

 1 Fig. 1. The problem statement.

F hFig. 2 .

 2 Fig. 2. Schematic view of the swimming model.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Schematic view of the electric model.

Fig. 6 .

 6 Fig. 6. The geometry of the 2-D simulated robot. All distances are in mm.

Fig. 7 .

 7 Fig. 7. The robot trajectory for the straight line swimming towards a wall test. The blue line is the trajectory of the head module.

Fig. 8 .Fig. 9 .

 89 Fig. 8. The electric fields around the insulating object A. The black lines are the tangents to the electric field lines.

Fig. 10 .

 10 Fig.10. The preset robot trajectory for the turning within a corner test. The blue line is the trajectory of the head module.

Fig. 11 .

 11 Fig. 11. The electric currents crossing the electrodes in function of the time t.

Fig. 12 .

 12 Fig. 12. The currents in infinite space with respect to α.

Fig. 13 .

 13 Fig. 13. Ī3 , Ī3∞ and I c 3 with respect to the time t.

Fig. 14 .

 14 Fig.[START_REF] Jawad | Sensor model for the navigation of underwater vehicles by the electric sense[END_REF]. The robot trajectory obtained using the control law for obstacle avoidance. The blue line is the trajectory of the head module.

  Fig.15. Īi , c h , c f , α, q 1 and q 2 with respect to the time t for the avoidance of a corner.

Fig. 16 .

 16 Fig. 16. Several trajectories of the robot with and without object.

http://www.theangelsproject.eu/

http://biorob.epfl.ch/

http://www.irccyn.ec-nantes.fr/

http://www-subatech.in2p3.fr/