Mathieu Porez
email: mathieu.porez@mines-nantes.fr

Frédéric Boyer
email: frederic.boyer@mines-nantes.fr

Jan Auke

Ijspeert
email: auke.ijspeert@epfl.ch

‡ A J Ijspeert

Improved Lighthill fish swimming model for bio-inspired robots -Modelling, computational aspects and experimental comparisons

Keywords: swimming dynamics, bio-inspired locomotion, Newton-Euler algorithmic, mobile multibody system dynamics

 ranging from hydrodynamic modelling to mobile multibody system dynamics. This article marks a further step towards the Lighthill theory. The LAEBT is applied to one of the best bio-inspired swimming robots yet built: the AmphiBot III, a modular anguilliform swimming robot. To that end, we apply a Newton-Euler modelling approach and focus our attention on the model of hydrodynamic forces. This model is numerically integrated in real time by using an extension of the Newton-Euler recursive forward dynamics algorithm for manipulators to a robot without a fixed base. Simulations and experiments are compared on undulatory gaits and turning manoeuvres for a wide range of parameters. The discrepancies between modelling and reality do not exceed 16% for the swimming speed, while requiring only the one-time calibration of a few hydrodynamic parameters. Since the model can be numerically integrated in real time, it has significantly superior accuracy compared with computational speed ratio, and is, to the best of our knowledge, one of the most accurate models that can be used in real-time. It should provide an interesting tool for the design and control of swimming robots. The approach is presented in a self contained manner, with the concern to help the reader not familiar with fluid dynamics to get insight both into the physics of swimming and the mathematical tools that can help its modelling.

Introduction

The locomotive performance of fish attracts great interest in submarine robotics. In particular, the quest for underwater vehicles that are both efficient and manoeuvrable has led to the design of a new generation of underwater robots inspired by fish [START_REF] Barrett | Propulsive efficiency of a flexible hull underwater vehicle[END_REF][START_REF] Mcisaac | Motion planning for anguilliform locomotion[END_REF][START_REF] Yu | Development of a biomimetic robotic fish and its control algorithm[END_REF][START_REF] Yamada | Development of amphibious snake-like robot ACM-R5[END_REF][START_REF] Morgansen | Geometric methods for modeling and control of free-swimming fin-actuated underwater vehicles[END_REF][START_REF] Crespi | Online optimization of swimming and crawling in an amphibious snake robot[END_REF][START_REF] Boyer | Fast dynamics of an eel-like robot -comparisons with Navier-Stokes simulations[END_REF][START_REF] Liu | Biological inspiration: from carangiform fish to multi-joint robotic fish[END_REF][START_REF] Salumäe | A flexible fin with bio-inspired stiffness profile and geometry[END_REF][START_REF] Stefanini | A compliant bioinspired swimming robot with neuro-inspired control and autonomous behavior[END_REF]. Whether it be for mechanical design or for control, researchers involved in the field need efficient and accurate models of swimming. In the context of bio-robotics these models have to be able to predict the relations between the internal joint (or shape) motions or forces and the external (or net) motions. Furthermore, this has to be achieved with computational times compatible with their exploitation for control, or fast simulation (as required by optimisation or path planning). Strictly speaking, deriving these relations requires solving the Navier-Stokes (N-S) equations which govern the time evolution of a real (i.e. viscous) fluid at any Reynolds number1 . However, the intrinsic nonlinearities of these equations are such that, to date, their numerical resolution (when it is possible, i.e. for sufficiently simple geometries at low Reynolds numbers [START_REF] Carling | Self-propelled anguilliform swimming: simultaneous solution of the two-dimensional Navier-Stokes equations and Newton's laws of motion[END_REF][START_REF] Liu | A numerical study of undulatory swimming[END_REF][START_REF] Farnell | Numerical model of self-propulsion in a fluid[END_REF][START_REF] Kern | Simulations of optimized anguiliform swimming[END_REF][START_REF] Eldredge | Numerical simulations of undulatory swimming at moderate Reynolds number[END_REF][START_REF] Bergmann | Modelling and simulation of fish like swimming[END_REF], or at high Reynolds numbers, with a simplified model of turbulence [START_REF] Leroyer | Numerical methods for RANSE simulations of a selfpropelled fish-like body[END_REF]) is still far from being achievable in real time. To overcome these limits, several simplifications of N-S equations can be applied. Firstly, due to the relatively high Reynolds numbers involved in fish-like swimming, the fluid around a fish (or a robot) can be considered as inviscid, i.e. it can only exert pressure forces (no shearing stress), and the N-S equations can be replaced by those of Euler [START_REF] Milne-Thomson | Theoretical hydrodynamics -fifth edition[END_REF]. Due to its inviscid character, this fluid model cannot predict the separations of the boundary layer which can occur near the body boundaries, which are responsible for the vortex generation, and in particular of the large vortices observed in the wake of fish [START_REF] Videler | Aquatic vertebrate locomotion: wakes from body waves[END_REF]. To address this lack of viscosity (and vorticity), models based on a perfect fluid (i.e. inviscid and incompressible) can be complemented by an exogenous model of vorticity generation derived from experimental observation and the exploitation of the singularities of potential flow theory2 [START_REF] Kelly | Modelling efficient pisciform swimming for control[END_REF][START_REF] Kanso | Swimming due to transverse shape deformations[END_REF]. These additional models produce vorticity which then obeys the usual conservation laws of a perfect fluid3 . Despite these simplifications, and even if one crudely neglects the wake vorticity, the numerical methods able to solve these swimming models (based on the boundary elements method, see [START_REF] Hill | Large amplitude fish swimming[END_REF][START_REF] Wolfgang | Hydrodynamics of flexible-body swimming motions[END_REF]) are still too computationally heavy to be used in robotics. In opposition to these numerical approaches, a few analytical models were developed in the pre-computer era. Both Wu [START_REF] Wu | Swimming of a waving plate[END_REF] and Lighthill [START_REF] Lighthill | Note on the swimming of slender fish[END_REF] based their models on a perfect fluid, with no vorticity apart from in the fish wake. In order to analytically integrate the equations of the 3-D flow around the body, the flow is replaced by a stratification of planar (i.e. 2-D) potential flows structured by the geometry of an undulating infinite height plate in the case of Wu [START_REF] Wu | Swimming of a waving plate[END_REF], and by a finite axisymmetric slender body, in the case of Lighthill [START_REF] Lighthill | Note on the swimming of slender fish[END_REF]. In both cases, analytic methods in planar harmonic functions enable swimming models in a closed analytical form to be obtained. Even if in some respects, as for instance the modelling of the proto-vortices along the body in the plane of swimming [START_REF] Videler | Aquatic vertebrate locomotion: wakes from body waves[END_REF]], Wu's approach has a few advantages, Lighthill's more realistic body model makes it the most relevant for the roboticist. Thus, we will consider the Lighthill model in the rest of the article, starting with a resumé of a few of the milestones that have marked its elaboration and explain its persistence.

The Lighthill theory was developed over several years [START_REF] Lighthill | Note on the swimming of slender fish[END_REF][START_REF] Lighthill | Aquatic animal propulsion of high hydromechanical efficiency[END_REF][START_REF] Lighthill | Large-amplitude elongated body theory of fish locomotion[END_REF]. In its essence, it consists in explains swimming by the timed variations of body shape that accelerate the (initially quiescent) fluid from front to back thus generating thrust. As has been conjectured by J. Gray4 , the swimming of an elongated fish is based on a kinetic amplification mechanism, where the fluid is gradually accelerated along the fish body by curvature waves of increasing amplitude. These travel from head to tail pushing the fluid laterally until it is shed into the wake, thus generating thrust by virtue of the action/reaction principle [START_REF] Breder | The locomotion of fishes[END_REF][START_REF] Gray | Studies in animal locomotion. I. The movement of fish with special reference to the eel[END_REF]. Starting from this principle, Lighthill's approach is based on two adaptations of the perfect fluid theory. Firstly, in order to avoid the modelling of the complex rotational flow in the fish wake, while taking it into account in the thrust, the modelling is restricted to the fluid contained in a half sphere of infinite radius whose the equatorial plane is at any time perpendicular to the trailing edge of the fish caudal fin. Moreover, based on experimental observations, the flow is assumed to be potential in that hemisphere and produced by the body undulations. Secondly, due to the high aspect ratio (length/thickness) of the body fish, the flow surrounding the fish and the forces which are produced by that flow can be modelled by the slender body theory which approximates the 3-D flow around the body by a stratification of potential 2-D flows perpendicular to the longitudinal axis of the body [START_REF] Munk | The aerodynamic forces on airship hulls[END_REF]. Based on these two ingredients, Lighthill made several crucial contributions to modelling fish swimming. Firstly, he considered an elongated body submitted to low amplitude undulations. Through a perturbation series expansion of Laplace equations with respect to the small parameters of the problem (aspect ratio, amplitude of deformations), he derived -via the direct integration of pressure forces -a closed form of the hydrodynamic thrust [START_REF] Lighthill | Note on the swimming of slender fish[END_REF]. Secondly, he reproduced the results from this first approach by using a more global model based on the balance of kinetic momenta applied to the flow inside the control volume defined by the hemisphere containing the body [START_REF] Lighthill | Note on the swimming of slender fish[END_REF].

Much easier to understand, this second more global point of view allowed Lighthill to derive the mechanism of kinetic amplification first conjectured by Gray, which can be translated in to mathematical terms as follows. When swimming, the body cross sections sweep past a fluid slice while transmitting more and more kinetic energy to it. Once the last cross section has swept past the fluid slice, the slice has a maximum kinetic energy and is shed into the wake so producing the reactive thrust. After his initial paper, Lighthill returned on several occasions to his elongated body theory (EBT) and proved that the thrust (deduced from the kinetic balance applied to the half space including the body) can equivalently be interpreted as generated by the vorticity produced at the trailing edge of the caudal fin [START_REF] Lighthill | Aquatic animal propulsion of high hydromechanical efficiency[END_REF]. Finally, he generalised the EBT to the case of finite body deformations so building the foundations of a large amplitude elongated body theory (LAEBT) [START_REF] Lighthill | Large-amplitude elongated body theory of fish locomotion[END_REF]. Since then, biologists have used the LAEBT to estimate the mean thrust and the swimming power of live fish [START_REF] Tytell | The hydrodynamics of eel swimming. I -Wake structure[END_REF][START_REF] Tytell | The hydrodynamics of eel swimming. II -Effect of swimming speed[END_REF]. More recently, Tytell and co-workers have successfully compared the Lighthill's reactive model to planar Navier-Stokes simulations obtained with a full model of lamprey swimming including body stiffness and muscular activation [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF]. As it is capable of modeling high amplitude shape variations such as those of current fish-like robots, the LAEBT is well adapted to roboticists' needs. Furthermore, it is well adapted to any fish morphology (e.g. caranguiform) that uses the body caudal fin swimming mode [START_REF] Breder | The locomotion of fishes[END_REF]. However, Lighthill's starting point in 1971 was the kinetic balance [START_REF] Lighthill | Large-amplitude elongated body theory of fish locomotion[END_REF] rather than the rigorous computation of the pressure field on the body as that proposed in his EBT [START_REF] Lighthill | Note on the swimming of slender fish[END_REF]. As a result, premising his calculation on a kinetic model of the flow (and not a kinematic model), the LAEBT is biased ab initio by a heuristic character which is not present in the EBT. Recently, we managed to overcome this limit, and the LAEBT has been derived by a direct calculation based on the integral of pressure forces along the body boundaries [START_REF] Candelier | Three-dimensional extension of the Lighthill's large amplitude elongated-body theory of fish locomotion[END_REF]. Furthermore, modelling the fish body as a nonlinear Cosserat beam, we extended Lighthill's model in [START_REF] Boyer | Fast dynamics of an eel-like robot -comparisons with Navier-Stokes simulations[END_REF] (numerically) and [START_REF] Boyer | Poincaré-Cosserat equations for the Lighthill three-dimensional large amplitude elongated body theory: application to robotics[END_REF] (theoretically) to: 1) the 3-D swimming; 2) the self-propelled swimming (i.e. the net motions are calculated rather than being imposed); 3) the computation of the internal stress field of the Cosserat beam (which models the forces of the fish muscles or robot motors); 4) the modelling of resistive forces through a Taylor-like resistive model [START_REF] Taylor | Analysis of the swimming of long narrow animals[END_REF]. This analytical model has been validated by comparisons with N-S simulations [START_REF] Boyer | Fast dynamics of an eel-like robot -comparisons with Navier-Stokes simulations[END_REF]. In the case of planar swimming gaits and manoeuvres, the trajectory and speed discrepancies between the analytical model and the N-S simulations do not exceed 10%, drastically reducing the requirements for the CPU time. In fact, the analytical model can run in real time, while N-S simulations take several hours (or days) for only a few seconds of swimming.

In this article, we will take this model one step further by applying it to one of the most advanced fish-like robots yet built: the AmphiBot III5 of the EPFL BioRob6 . We extended the Newton-Euler (N-E) based algorithm of forward manipulator dynamics [START_REF] Featherstone | The calculation of robot dynamics using articulated-body inertias[END_REF] to the case of a locomotion robot without a fixed base. This approach has the advantage of offering simple access to insight into the physics of swimming [START_REF] Khalil | Dynamic modeling and simulation of 3-D serial eel-like robot[END_REF], while minimising the computational complexity and the time execution. Starting from the description of the AmphiBot III in section 2, we will then derive its dynamic model in the N-E formalism in section 3. In the same section, special attention will be paid to the modelling of the external forces exerted by the fluid onto the robot. The model we describe in this article is self-contained and original in itself. In contrast to our earlier model [START_REF] Boyer | Poincaré-Cosserat equations for the Lighthill three-dimensional large amplitude elongated body theory: application to robotics[END_REF] which was based on abstract variational calculus, the model in this article uses the balance of the fluid kinetic momenta contained in a control volume enclosing each of the segments. Once derived, this model is introduced into the N-E model of the AmphiBot III. Based on this model, the extended Featherstone forward dynamics algorithm is described in section 4 and then used in section 5.3 to compute the motion of the AmphiBot III. These results are then compared to those given by a set of experiments carried out with the real robot. The comparisons concern undulatory gaits and turning manoeuvres. With a very small set of hydrodynamic (dimensionless) parameters each tuned only once, the discrepancies between the model and reality did not exceed 16% in terms of averaged forward velocity, confirming the validity of the LAEBT. This model is, to the best of our knowledge, the most accurate that can be run in real-time on a standard computer. The article concludes with a discussion of the results in section 6.

The snake-like robot AmphiBot III

The snake-like robot AmphiBot III (see figure 1) is a 2-D swimming robot composed of 8 segments, the first (from left to right) being the head segment, while the last segment, or "tail segment", supports a caudal fin. This fin is rigidly clamped to the last segment and is sufficiently stiff7 to be considered as rigid in the subsequent modelling. The AmphiBot III is designed to be modular. In particular, its segments are constructed from several identical polyurethane parts. The external casing of each segment consists of three parts: a top, a bottom and a central box, which are fixed together with magnets. Waterproofing of the assembly is achieved by custom O-rings placed between each part. The segments are connected (both mechanically and electrically) by an elastic connector (see figure 1). A segment has an external length of 9.7 cm, and a cross-section of 4.0 cm × 5.7 cm. The total length of the robot (denoted by L r ) is 88.0 cm. The caudal fin has a length of 10.3 cm and a cross-section of 0.3 cm × 5.7 cm. This relatively a j ,a j+1 a j+2 ,a j+3 a j-2 ,a j-1 j -1 j +1

Figure 2: The series elastic actuators of the AmphiBot III.

low thickness allows one to consider its trailing edge sufficiently sharp to justify the use of the LAEBT. The robot's density is slightly less than 1.0 kg/m 3 , so that the robot swims just below the water surface. The center of mass of each segment is below its geometric center, which ensures the vertical stability of the robot during swimming. Each actuated segment contains a power board, a proportional derivative (PD) position controller, a small water detector, a direct current (DC) gear-motor with an integrated incremental encoder and a rechargeable lithium-ion battery. The DC gear-motor used here has a maximum torque of 0.57 N.m and drives the elastic connector forming "a serial elastic actuator" [START_REF] Spong | Modelling and control of elastic joint robots[END_REF] (see figure 2). The head segment, which resembles a body segment without a motor, contains the locomotion control board which sends out the desired positions to body segments in real-time via a serial bus. A radio device embedded in the head segment permits control commands and swimming parameters to be sent from a controlling computer to the robot when it is just below the water. Its relatively low weight/power ratio, compliant joints, and caudal fin, make the AmphiBot III an excellent swimmer. According to Lighthill's theory, the sharp trailing edge of its tail is largely responsible for its good performance, being responsible for the non zero average thrust produced by body undulations. We show in section 5.3.1 a three fold increase in the AmphiBot III's swimming speed compared to the results obtained with no caudal fin [START_REF] Crespi | Online optimization of swimming and crawling in an amphibious snake robot[END_REF], whatever Figure 3: Schematic views of the fish-like robot and its internal kinematics. In the bottom view, we have voluntarily off-centered the axis of the joint j + 1 with respect to the axis of the joint j in order to make appear on the view the parameters of the rigid shaft j.

O j d j O j-1 O 0 F 0 F j

Segments

the swimming parameters used.

3 Dynamic model of the AmphiBot III

Parametrisation and notation

The AmphiBot III's kinematic schema is presented in figure 3. Let us consider it, and more generally, any serial robot, swimming in a fixed horizontal plane in an unbounded volume of initially quiescent water. We attach a fixed spatial frame denoted by F e = (O e , s e , n e , a e ) to the ambient geometric space, where the unit vector a e is normal to the robot swimming plane. The robot is a mobile multibody system (MMS) composed of a sequence of m + 1 rigid segments interconnected through m actuated joints. Because the joints have a compact geometry, the serial discrete assembly of segments is approximated by a continuous boundary in contact with water. Between two segments, is a "serial elastic actuator" joint [START_REF] Spong | Modelling and control of elastic joint robots[END_REF]: a gear-motor serially connected to a rigid shaft and a revolute spring (see figure 2). The m shafts and the m+1 segments are modeled by 2m + 1 rigid bodies, while each "serial elastic actuator" introduces 2 degrees of freedom (one for the motor axis and one for the rotational spring). As a result, the serial kinematic chain of the AmphiBot III has n = 2m revolute joints parameterised by the m angles of gear-motors and the m strain angles of the revolute springs. Overall, the MMS is composed of n + 1 rigid bodies denoted B 0 , B 1 , ..., B n , with B 0 and B n representing the head and tail segments respectively. According to this convention, an even index denotes a segment while an odd index denotes a shaft. We attach to each body B j a mobile frame F j = (O j , s j , n j , a j ), whose center O j coincides with the center of the joint j. The unit vector s j is supported by the line O j O j+1 , and a j supports the joint axis. At any time t, the robot configuration is defined by the vector of joint positions r(t) = (r 1 , ..., r n ) T (t) defining the relative angles between the segments and shafts, together with the orientation matrix e R 0 and the position vector e P 0 of the mobile frame attached to the head segment F 0 = (O 0 , s 0 , n 0 , a 0 ) with respect to F e . Thus B 0 is the reference body and the time evolution of the pose ( e R 0 , e P 0 ) defines the rigid net motion of the robot. At time t = 0 s this reference pose, and its velocity are known and must be updated by time integration of the head segment accelerations using a computational algorithm which is presented in section 4. Throughout this article we will use the following notation convention. For any physical variable modelled as a tensor, the right lower index will represent the body index (to which it is related) while the left upper exponent will indicate the index of the projection frame (e.g. e R 0 , e P 0 ). When the tensor related to a body is expressed in the mobile frame of this body, the upper index is omitted. Moreover, the temporal derivative ∂./∂t will sometimes be denoted by a 'dot'.

Mobile multibody system model

As the AmphiBot III is composed of a relatively high number of serially connected rigid bodies and revolute joints, we chose to model it using a N-E framework described elsewhere [START_REF] Khalil | Dynamic modeling and simulation of 3-D serial eel-like robot[END_REF][START_REF] Boyer | Recursive inverse dynamic modeling of mobile multi-body systems with joints and wheels[END_REF]. This work is generally devoted to modelling serial MMS, with a mobile base (in this case the head segment). Let us start by introducing the geometric model of the robot which relates the posture of the frame F j with that of the frame F j-1 , both expressed in the earth frame and represented by the two (4 × 4) matrices e g j-1 and e g j of SE(3). This model can be expressed as:

e g j = e g j-1 j-1 g j = e g j-1 j-1 R j j-1 P j 0 1 , (1)
where j-1 R j and j-1 P j are the orientation matrix and the position vector of F j with respect to F j-1 . With the AmphiBot III geometry, we have:

j-1 R j =   cos r j -sin r j 0 sin r j cos r j 0 0 0 1   , and j-1 P j =   d j 0 0   ,
where d j is the distance between the frame centers O j-1 and O j (see Figure 3).

Regarding the velocity of the body j, it is a (6 × 1) vector of se(3) denoted η j and related to the velocity of body j -1 through the recursive relation:

η j = V j Ω j = Ad j g j-1 η j-1 + ṙj A j , (2)
where V j and Ω j are respectively the linear and angular Galilean velocities of the body under consideration (both expressed in its mobile frame), A j = (0 T 3 , a T j ) T is the (6 × 1) unit vector supporting the joint axis j, and Ad j g j-1 is the adjoint map operator which permits a change in velocity from F j-1 to F j :

Ad j g j-1 = j R j-1 j R j-1 j-1 P T j 0 j R j-1 . ( 3
)
Let us remark that in (3), we introduced the notation 'hat' which changes a (3 × 1) vector into its associated (3 × 3) skew-symmetric tensor. Thus, for any vectors A and

B in R 3 , Â is defined such that ÂB = A × B.
Once the Galilean velocities are defined, by temporal derivation of (2), the acceleration of the body j is given by the relation:

ηj = Ad j g j-1 ηj-1 + ζ j + rj A j , (4)
where ζ j represents the component of acceleration in (4) which depends on velocities through the detailed expression:

ζ j = ( j V j-1 + j P j-1 × j Ω j-1 ) × ṙj a j ṙj j Ω j-1 × a j . ( 5
)
Finally, by applying, Newton's law and Euler's theorem to the j th body, one obtains the dynamic equations of body j in the Newton-Euler form:

F j = M j ηj + β j + Ad T j+1 g j F j+1 -F ext,j . (6)
Here, we have introduced the following notation:

• for any j, F j is the (6 × 1) force vector (element of se(3) * ) exerted by the body j -1 onto the body j; and

• M j is the (6 × 6) inertia tensor of B j (element of se(3) * ⊗ se(3)), which can be defined as:

M j = M j -M S j M S j I j = ρ b V B j 1 3 -Ô j Q Ô j Q -Ô j Q Ô j Q dV B j ,
where Q is a point of B j , 1 3 is the 3 × 3 unit matrix, ρ b is the body mass ratio while M j , M S j and I j are the tensor of body mass (spherical in the rigid body case), the tensor of first inertia moments (skew-symmetric in the rigid body case) and the tensor of angular inertia of the body j; and

• the (6 × 1) vector of Coriolis and centrifugal forces, where:

β j = -Ω j × (M S j Ω j ) + Ω j × (M j V j ) Ω j × (I j Ω j ) + M S j (Ω j × V j ) , (7)
(we shall see that the fluid generates other velocity-dependent inertia forces that we will also denote β by extension); and

• the (6 × 1) vector of external hydrodynamic forces denoted by F ext,j which will be defined in the next section (3.3).

O j S j (0)

F j f - pres,j X = l j X S j (X) C j (X) s j n j a j f reac,j c reac,j S j (l j ) -f + pres,j

Model of hydrodynamic forces

Following the terminology in [START_REF] Lighthill | Large-amplitude elongated body theory of fish locomotion[END_REF] we distinguish two types of hydrodynamic forces, produced either by the reaction forces exerted by the acceleration of the fluid around the robot body or by the viscous stresses applied in its boundary layer. In the first case, the forces are said to be "reactive" while in the second, they are termed "resistive". According to [START_REF] Boyer | Fast dynamics of an eel-like robot -comparisons with Navier-Stokes simulations[END_REF], the (6 × 1) vectors of body external (hydrodynamic) forces of ( 6) are defined by the following superimposition of forces:

F ext,j = F reac,j + F res,j if j is even, 0 otherwise, ( 8
)
where F reac,j = (f T reac,j , c T reac,j ) T and F res,j = (f T res,j , c T res,j ) T are the (6 × 1) vectors of reactive and resistive forces respectively. We will now to detail both of them. In the first case, we use the model of the perfect fluid, while in the second, we account for the viscous effects.

The model of reactive forces

As in Lighthill's original theory, the flow around the fish (or robot) is assumed to be potential8 except in its wake where vorticity is shed. The influence of boundary discontinuities, introduced by the joints, onto the hydrodynamic forces exerted on the body, are neglected. These approximations, and the high aspect ratio (length/thickness) of the AmphiBot III, allow the 3-D fluid potential flow laterally surrounding the robot to be approximated using a stratification of 2-D potential flows perpendicular to the robot's backbone. In this simplification of the fluid kinematics9 , the 3-D fluid volume is replaced by a 1-D medium of fluid slices which sweep past the robot cross-sections. As illustrated in figure 4, if we denote by X the abscissa of the crosssection C j (X) of the j th segment, a fluid slice will be defined as the part of the fluid contained at each instant in the geometric section S j (X) which extends C j (X) in the fluid. As the flow is potential, the 2-D velocity field of the fluid in S j (X) is governed by the 2-D Laplace equation submitted to the boundary conditions imposed by the body transverse motion of the robot on C j (X) and by those of a quiescent fluid at infinity 10 . Since the fluid is perfect (i.e. non viscous), the fluid slices of the 1-D stratified fluid can only transmit compression forces along the longitudinal segment axis (aligned with s j ). These cross-sectional pressure forces define the field f pres,j : X ∈ [0, l j ] → f pres,j (X) along each of the segments 11 . In the case of 2-D swimming 12 , these cross-sectional pressure forces can be derived by integrating the unsteady Bernouilli's pressure law 13 in each of the fluid slices [START_REF] Lighthill | Large-amplitude elongated body theory of fish locomotion[END_REF] while in the 3-D case, their computation requires further investigation [START_REF] Candelier | Three-dimensional extension of the Lighthill's large amplitude elongated-body theory of fish locomotion[END_REF]. In both cases, the cross-sectional pressure forces take the form:

f pres,j (X) = Tf,j (X)s j , (9)
where Tf,j (X) is the kinetic energy of the fluid contained in S j (X).

Based on these preliminary findings, we are going to state the balance of kinetic momenta of the fluid laterally surrounding the robot in order to derive the reactive force F reac,j (of ( 8)) exerted onto each of the robot segments. To that end, let us recall that for a given geometric volume (named the "control volume" in fluid dynamics) containing material, N-E equations state the balance of the external forces and moments exerted across the boundaries of the geometric volume onto that material with the time derivative of its kinetic linear and angular momenta. In our case, we consider the fluid which laterally surrounds the segment j, i.e. the fluid contained between the two geometric slices S j (0) and S j (l j ). As we have just seen, the external forces applied on that fluid domain are the reactive force -F reac,j that we seek and the two sectional pressure forces applied onto the isolated fluid across S j (0) and S j (l j ), i.e. f pres,j (0) and -f pres,j (l j ) respectively, which can both be calculated using (9). Thus, applying Newton laws and Euler theorem to that fluid domain (as detailed in the appendix in section 8), finally gives:

F reac,j = -M reac,j ηj -β reac,j -β - reac,j -β + reac,j , (10)
where we have introduced the following definitions:

• M reac,j is the (6 × 6) tensor of added inertia of the fluid accelerated by the segment j:

M reac,j = M f,j -M S f,j M S f,j I f,j =   l j Mf,j l 2 j 2 Mf,j ŝT j l 2 j 2 Mf,j ŝj l 3 j 3 ŝT j Mf,j ŝj + l j Īf,j   , (11)
which depends only on the segment's cross sectional added inertia tensors Mf,j and Īf,j along with its geometry;

• β reac,j is the (6 × 1) vector of the reactive forces produced by the fluid added mass accelerated by the Coriolis and centrifugal accelerations of the segment j:

β reac,j = -Ω j × (M S f,j Ω j ) + Ω j × (M f,j V j ) Ω j × (I f,j Ω j ) + M S f,j (Ω j × V j ) + 0 V j × (M f,j V j ) ; ( 12
)
10 In the rest of the article, the effects of the walls of the pool in which the robot is experimented are neglected.

11 By analogy with solid beam theory, fpres,j(X) represents the compression force applied across Sj(X) by the fluid on its left side onto the fluid on its right side.

12 By "2-D" we here mean that the robot moves in a plane while generating a 3-D flow around its body.

13 A potential flow is entirely parameterised by a unique point-function named the velocity potential from which the velocity field derives. The velocity potential is governed by Laplace equations and Neumann boundary conditions on solid boundaries and defines the flow kinematics. As regards internal forces, they are modelled by a pressure field governed in each point by the unsteady Bernoulli equations.

• β - reac,j and β + reac,j are the (6 × 1) boundary vectors (the two indicesand + indicate the left and right planes S j (0) and S j (l j )) given by:

β - reac,j = - V X Pf,j V X Σf,j (0) + f - pres,j 0 , (13)
and

β + reac,j = Ad T l j g j V X Pf,j V X Σf,j (l j ) - f + pres,j 0 , ( 14
)
where V X , Pf,j and Σf,j are detailed in the appendix in section 8 and respectively denote the linear velocity of the X-cross section C j (X) projected onto s j , its cross-sectional linear momentum and its cross-sectional angular momentum. Physically, β - reac,j and β + reac,j model two kinds of forces. The first ones (first terms of ( 13)-( 14)) account for the net flux of fluid kinetic momenta across the control volume boundaries14 S j (0) and S j (l j ), both moving with segment B j while the fluid slices do not move along the segment axis15 . The second terms of ( 13)-( 14) are produced by the pressure exerted by the fluid outside of the isolated control volume onto its two boundaries S j (0) and S j (l j ).

As previously mentioned (see section 3.1), if the serial assembly of segments is approximated as a continuous boundary in contact with the water, the conservation of the mass, and the action/reaction principle requires the convective and pressure forces transmitted across the slices to cancel each other out, except those applied to the first (from left to right) section of the head segment and the last section of the tail segment. In other terms, among all the forces described in ( 13)-( 14), only β - reac,0 and β + reac,n will generate net forces on the robot. Furthermore, the cross-section of the robot nose being a single point, its cross-sectional added mass tensor Mf,0 (0) and Īf,0 (0) are equal to zero and hence β - reac,0 = 0. On the other hand, β + reac,n , is not zero and models: 1) the effects of the flux of kinetic momentum across S n (l n ) (the terms function of V X (l n ) in ( 14); 2) the pressure force -f + pres,n exerted across the plane S n (l n ) on the fluid on its left side by the fluid on its right side. From (9), f + pres,n is directed along s n with a strength equal to the kinetic energy of the fluid contained in

S n (l n ), i.e. f + pres,n = f pres,n (l n ) = T f,n (l n )s n . As a result, if V (X)
and Ω(X) denote the linear and angular velocity of the X-cross section C j (X) in the segment frame, we can write (see (26) in appendix):

β + reac,n = Ad T ln gn ( Mf,n V )V X ( Īf,n Ω)V X (l n ) - 1 2 (V T Mf,n V + Ω T Īf,n Ω)(l n )s n 0 . ( 15
)
Finally, in the planar case (where Ω X = 0, V Z = 0), we reproduce the conditions in which the original LAEBT was formatted [START_REF] Lighthill | Aquatic animal propulsion of high hydromechanical efficiency[END_REF]. In particular, the second term of ( 15) is the key explanation of thrust generation. Indeed, it is responsible for the non-zero average thrust and models the reaction force exerted on the caudal fin and generated by the cross-sectional kinetic energy amplified along the robot body by the curvature wave and shed in the wake, or equivalently by the vorticity generated along the sharp trailing edge. This is demonstrated in [START_REF] Lighthill | Aquatic animal propulsion of high hydromechanical efficiency[END_REF] for the 2-D swimming in small body deformation and discussed in [START_REF] Boyer | Poincaré-Cosserat equations for the Lighthill three-dimensional large amplitude elongated body theory: application to robotics[END_REF][START_REF] Candelier | Three-dimensional extension of the Lighthill's large amplitude elongated-body theory of fish locomotion[END_REF] in the case of 3-D swimming with finite deformations.

The model of resistive forces

As far as the resistive model is concerned, the shear stresses exerted by the fluid in the robot boundary layer are approximated through a Taylor-like resistive model [START_REF] Taylor | Analysis of the swimming of long narrow animals[END_REF]. In this model, the resistive forces take the form of a field of cross-sectional steady drag forces and couples, which resist the lateral motion of the cross sections:

F res,j = 1 2 l j 0 Ad T X g j C 1 (|V X |V X )s j + C 2 (|V Y |V Y )n j + C 3 (|V Z |V Z )a j (C 4 |Ω X |Ω X )s j dX , where V Y (X) = V (X) T n j , V Z (X) = V (X) T a j
and Ω X (X) = Ω(X) T s j are the components of the lateral velocity of the X-cross section C j (X) in the B j mobile frame, while C i for i = 1, 2, 3, 4 are resistive coefficients derived from experimental fluid dynamics which are further described in section 5.2.2.

Newton-Euler model of the swimming robot

Once the model of external forces has been developed, we can model of the robot's dynamics.

According to ( 8) and ( 10), F ext,j = 0 if j is odd (B j is a shaft), and :

F ext,j = -M reac,j ηj -β ext,j , (16)
if j is even (B j is a segment), with:

If j = n: β ext,j = β reac,j + F res,j , while :

β ext,n = β reac,n + β + reac,n + F res,n . (17)
It is worth noting here that all the terms of ( 17) are known functions of the state variable η j .

Finally, inserting this model in ( 6) gives the N-E model of the segments of our swimming robot.

For j = 0, 1, ..., n -1:

F j = Ad T j+1 g j F j+1 + (M j + M reac,j ) ηj + β j + β reac,j + F res,j if j is even, Ad T j+1 g j F j+1 + M j ηj + β j otherwise, (18)
and for j = n:

F n = (M n + M reac,n ) ηn + β n + β reac,n + F res,n + β + reac,n . (19)
Once completed with recursive kinematics (1), ( 2) and ( 4), these equations can be exploited to solve the forward dynamics of a swimming robot as detailed in the next section.

The computational algorithm

To compute the dynamics of our robot, we have used Khalil et al 's recursive forward dynamic algorithm based on the N-E equations [START_REF] Khalil | Dynamic modeling and simulation of 3-D serial eel-like robot[END_REF]. This algorithm is used in a global timeintegration loop. It is a generalisation of Featherstone's N-E algorithm of serial manipulators [START_REF] Featherstone | The calculation of robot dynamics using articulated-body inertias[END_REF] to the case of a system without a fixed base. According to the flow chart in figure 5, it gives the head (net) and joint (shape) accelerations, i.e. η0 and r, as a function of joint torques Γ and the current state (g 0 , η 0 , r, ṙ). At each time step of a global time integration loop, the direct algorithm executes the three following recursions with respect to the body index.

Swimming gait

Joint laws Direct dynamic model

η0 1 s 2 1 s 2 1st forward recursion loop 2nd forward recursion loop 1st backward recursion loop η0 η 0 , e g 0 Direct dynamic model r, η0 η 0 , e g 0 , ṙ, r, Γ ṙ, r r r d Γ e g j , η j , ζ j , M j , β j , M reac,j , β ext,j H j , K j , α j , M * j , β * j η0 = -(M * 0 ) -1 β * 0

The first forward recursion on the kinematics

As the current robot's state ( e g 0 , η 0 , r, ṙ) is known, the algorithm starts by the following forward (from the head to the tail) recursion:

For j = 0, ..., n, compute:

• the j-1 R j , j-1 P j and the body transformations e g j from recursion (1),

• the body velocities η j from recursion (2),

• the terms ζ j of (4) from ( 5),

• the body inertia matrices M j and M n of ( 18)-( 19) from ( 7),

• the body Coriolis and centrifugal forces β j and β n of ( 18)-( 19) from ( 7),

• the added inertia matrices M reac,j from (11),

• the fluid Coriolis and centrifugal forces β reac,j from (12),

• the velocity dependent hydrodynamic forces β + reac,n and F res,j of ( 18)-( 19) from ( 15)-( 16).

End for.

The backward recursion on the "external" dynamics

Once all the state-dependent variables known, the next step of the computational algorithm consists in calculating the rigid net acceleration of the robot, i.e. the head acceleration η0 , as follows:

η0 = -(M * 0 ) -1 β * 0 . ( 20
)
This represents the external dynamics, i.e. the dynamics of the net-motions, with M * 0 (r) and β * 0 (Γ, r, ṙ, η 0 ) respectively denoting the (6 × 6) inertia matrix and the (6 × 1) vector of Coriolis, centrifugal and external (control and hydrodynamic) forces, both expressed in the head (reference) frame with the joints locked in their current configuration. They are given by the following backward recursion (from the tail to the head) initialised by M * n = M n and β * n = β n + β ext,n :

For j = n, n -1, ..., 1:

H j = A T j M * j A j ; K j = M * j -M * j (A j H -1 j A T j )M * j , α j = K j ζ j + M * j A j H -1 j Γ j -A T j β * j + β * j , M * j-1 = M j-1 + M reac,j-1 + Ad T j g j-1 K j Ad j g j-1 , β * j-1 = β j-1 + β ext,j-1 + Ad T j g j-1 α j ,
End for.

The second forward recursion on the "internal" dynamics

Finally, the algorithm ends with a third forward recursion initialised by the current state as well as η0 , computed from (20):

For j = 1, 2, ..., n:

rj = H j -A T j M * j Ad j g j-1 ηj-1 + ζ j + Γ j -A T j β * j , ηj = Ad j g j-1 ηj-1 + ζ j + rj A j .
End for.

Finally, the time t, the head state ( e g 0 , η 0 ) and the joint position and speed vectors ( ṙ, r) are updated, and the algorithm resumes according to the flow chart in figure 5. Moreover, let us note that, technically, η0 is numerically integrated with a numerical integrator based on quaternions.

Comparisons between experiments and simulations

In order to assess the dynamic model of the robot (and especially the model of hydrodynamic forces) in section 3, we compared experimental data from the AmphiBot III with simulation results obtained by applying the algorithm of section 4 with the model in section 3.

Experimental setup

The trajectory and speed of the AmphiBot III robot were measured in a swimming pool with a video tracking system, as illustrated in figure 6. The pool measured 6.0 m × 1.5 m × 0.3 m and was filled with water to a depth of approximately 0.2 m. The video tracking system was composed of 3 mm light-emitting diodes (LEDs) fixed to each robot segment and two black and white video BASLER A622F cameras (each with a TAMRON M12VM412 wide-angle lens). Each camera has a resolution of 1280×960 pixels in 8-bit grey scale and a frame rate of 15.0 frames per second. A time trigger synchronises the time of the frame capture. The frames processed using "in-house" software based on the OpenCV library. For each new frame capture, the software assembles the images from the video cameras, corrects the frame distortions due to the wide-angle lenses, extracts from the corrected frames, by a thresholding process, the positions of the segment LEDs with respect to a frame attached to the pool, and saves the computed positions in a file. Then, from (1) and the positions recorded in the dark16 , the robot's head trajectory and the body shape are computed and, by derivation in time, the forward swimming speed is estimated. After calibration, the tracking system measured position to an accuracy that was ±5 mm.

Parameters of the dynamic model

Before presenting the comparison of the results, we have to specify the parameters used in our model in order to fit the dynamics of the model of section 3 with those of the AmphiBot III. This fitting concerns the MMS model parameters (including those of the controllers) along with the hydrodynamic coefficients. This calibration is achieved once for all by starting from physical considerations and, using as desired shape motions, the law defined (for j odd) as:

r d,j (t) = A(cos(2π(νt + jk/m)) + α). ( 21
)
This law is well known from the zoological literature on elongated (anguilliform) fish practicing body caudal fin swimming [START_REF] Breder | The locomotion of fishes[END_REF][START_REF] Gray | Studies in animal locomotion. I. The movement of fish with special reference to the eel[END_REF]. It corresponds to a wave of body curvature propagating along the fish backbone from the head to the tail superimposed on a constant curvature offset. In (21), ν is the wave frequency, A is the wave amplitude, k is the number of waves along the robot backbone and α governs the turning angle.

Parameters of the mobile multibody system model

The AmphiBot III introduced in Section 2 is composed of 8 segments attached together by m = 7 serial elastic actuators. (shafts and segments) connected together. Using the numbering convention introduced in section 3.1, the mechanical parameters of the model are those of table 1. Note that all the segments and shafts have the same characteristics except the last segment (j = 14) which supports the rigid caudal fin (see the left of figure 1). The torques applied to the robot joints are related to the joint motions (see figure 5) through the following physical law when the joints are passive (i.e. when j is even):

Γ j = -K e r j ,
where K e is the stiffness of the passive joints (assumed to be identical for all the passive joints) whose value, (given in table 1), has been measured by static tests. On the other hand (j is odd), the active joints are ruled by the following PD control law (identical for all the active joints):

Γ j = K p (r d,j (t) -r j ) -K d ṙj , ( 22
)
where K p and K d are the proportional and derivative gains of the robot controllers, given in table 1, while r d,j (t) is the desired joint position known by its time evolution (21).

Parameters of the hydrodynamic model

As mentioned in section 2, the AmphiBot III robot is statically balanced in order to swim in an horizontal plane just below the water surface. As the pitching and rolling net motions are very limited, we will ignore their dynamics in the simulations and the reactive model ( 11), ( 12) and ( 15) can be computed with (for j = 0, 1..., n):

Mf,j = C m ρ f πh 2 4 (n j n T j )
, and Īf,j = 0 , where ρ f is the fluid volume mass, h = 5.7 cm is the segment immersed height and C m is a dimensionless shape coefficient given in Table 1.

The resistive model is derived from experimental fluid dynamics [START_REF] Hoerner | Fluid Dynamics Drag[END_REF]. The coefficients C i of ( 16) are defined by:

C 1 = C f ρ f P , C 2 = C d ρ f h , C 3 = 0 ,
and C 4 = 0 , where P = 19.4 cm is the segment cross-section perimeter, C f and C d are the dimensionless friction and drag coefficients stated in table 1, while the remaining drag coefficients are forced to zero because of the planar nature of the motions. It is worth noting here that the hydrodynamic model follows the modular nature of the AmphiBot III. Thus, since all the segments are identical, the entire hydrodynamic model is parameterised by only 3 coefficients (C m , C f and C d ) which can be fixed once in the following comparisons.

Comparisons and discussions

For validation purposes, the algorithm presented in section 4 has been coded in C. The developed simulator uses a fourth-order Runge-Kutta integrator for the time integrations (with a time step of 5 × 10 -4 s) and a numerical integration by Gauss quadrature (with 6 points) for the spatial integration of the resistive forces along the robot's segments ( 16). These choices allow us to compute the swimming motion of our virtual AmphiBot III robot in real time on a laptop (with a 2.0 GHz dual-core CPU and 2048 Mb of RAM). In this section, we present comparison results for two types of shape motions well known to zoological literature [START_REF] Breder | The locomotion of fishes[END_REF][START_REF] Gray | Studies in animal locomotion. I. The movement of fish with special reference to the eel[END_REF]: the forward swimming gait and the turning manoeuvre.

Steady-state forward swimming gaits

In this first study, we systematically tested how the three swimming parameters ν, A and k of ( 21) affect the forward swimming speed of the robot. To that end, 100 simulations were performed, with 5 experiments per simulation, leading to 500 experiments, for 4 values of ν ranking from 0.4 Hz to 1.0 Hz (with a step ∆ν = 0.2 Hz), for 5 values of A ranking from 25 • to 45 • (with a step ∆A = 5 • ) and for 5 values of k ranking from 0.5 to 1.5 (∆k = 0.25). Let us note that for this study, the turning angle α is fixed to 0. For each gait, the (simulated or real) robot, starts from resting and is stopped after 10 s of swimming or when it has travelled 5 m (i.e. when it reaches the opposite side of the pool). Then, the following cruising forward swimming speed:

Vf = || e V0 || , with e V0 = 1 T t+T t e V 0 dτ ,
is computed with T = 1/ν and for a t taken beyond the transitory regime of the gait. The left and middle columns of 7 (a)-(d) illustrate the results obtained numerically and experimentally. The right column is related to the following relative error between simulations and experiments:

ǫ = 100 ( Vf,exp -Vf,sim ) 2 V 2 f,max 1/2 (%).
Undulations with an amplitude greater than 40 • and a wave number less than 0.75 are impossible experimentally because the control torques exceed the intrinsic limits of the servomotors and may damage them. The squares with a white cross bottom right in the panels of figure 7 represent this condition. The variability of the experimental measurements is measured by the averaged standard deviation which is about 0.018 m/s. Depending on the experimental parameters, the cruising speed of the robot varied between 0.05 and 0.59 m/s, i.e. in the range of real fish [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF]. For any swimming frequency ν, the maximum speed was reached when A = 30 • and k = 0.5. This agrees well with the simulation for which the swimming speed varied between 0.07 and 0.55 m/s while for any swimming frequency ν, the maximum speed was reached when A = 25 • and k = 0.5. This good agreement between the simulation and the experiment is also seen in the panels of figure 7 which show a maximal relative error ǫ ≃ 16% and a mean error ǭ ≃ 4%. The discrepancy between the experiments and the simulations can be seen to increase with the frequency ν. This is probably due to the fact that the water in the pool becomes increasingly agitated by the robot, which becomes subjected to surface waves reflecting from the pool's walls, a phenomenon which is not accounted for by the model. These effects increase (upper part of panels) with the number of waves (k) per (robot) length (L r ) and with their amplitude (A). The energy transmitted to the unmodeled free surface waves increases with A and k so amplifying the phenomenon. Roll oscillations which are neglected in our model of the planar swimming are another perturbation effect which could explain the (small) discrepancies in the range of small wave numbers (for example, for k = 0.5 and ν ≥ 0.8 Hz). In this case a small amount of the energy supplied by the motors goes into the rolling dynamics. Finally, as observed on live fish [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF], for sufficient wave lengths (k = 0.5, k = 0.75), the velocity increases monotonically with the frequency while after increasing with the amplitude, it attains a maximum value beyond which it decreases. This first set of tests is related to the dynamics of the external d.o.f. and does not give any indication about the internal d.o.f. (shape) of the robot. However, recent studies in biology tend to show that kinematics and swimming performance are strongly affected by internal parameters, particularly body stiffness [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF]. In particular, in our case the hydrodynamic modelling errors could be hidden by the calibration of the compliant and controlled joint dynamics. As a result, we had to assess the accuracy of the internal dynamic model. To that end, we reconstructed the shape of the robot using the video tracking system and compared these observed shapes to those calculated using the simulation. For the purpose of illustration, figures 8 (a)-(d), show a series of snapshots for 4 arbitrary chosen forward swimming gaits, sampled throughout one swimming cycle (between 3T and 4T , i.e. beyond the transient phase whose duration is about 3T ) comparing the simulator and the robot. Once again, good matching is observed, both in terms of external (head) and internal (shape) d.o.f. To quantify the errors due to internal (shape) dynamics modelling, we should compare the angle of the joint behind each series elastic actuator on the simulator and the robot. Unfortunately, the angular encoders of the AmphiBot are located between the actuators and the gears. To bypass this difficulty, we decided to compute the lateral displacement maxima for each segment in each of the gaits of figure 8. The results are displayed in figure 9. Such a measure accounts not only for the errors brought by the internal (shape) dynamics but also, since the head is free, a small proportion of those due to the external (lateral) dynamics. As a result, the errors indicated in figure 9, though they do not exceed 20%, are greater than those we could deduce from the afore mentioned direct measurement of the joint angles after the gears. Beyond the unavoidable errors introduced by the LED tracking device, the causes of these errors are the imprecisions introduced by the model of hydrodynamic forces, those due to the identified model of compliant joints (mass-spring parameters) and joint servo systems, as well as others due to unmodeled phenomena such as the actuators' saturation or the friction and backlash of the gears. Finally, for the four forward swimming tests of figure 8, the hydrodynamic thrust T h along the unit vector supporting the averaged forward velocity (denoted u) has been calculated from the relation:

T h = u T .   n j=0 Ad T j ge F ext,j   , with u = e V0 || e V0 || .
The results are displayed in figure 10. As expected from the Lighthill model, these plots show that among all the reactive forces, only the caudal term produces a non zero averaged thrust over one cycle. As discussed in section 3.3, this term models the forces exerted by the vorticity shed into the wake at the caudal fin's trailing edge [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF][START_REF] Candelier | Three-dimensional extension of the Lighthill's large amplitude elongated-body theory of fish locomotion[END_REF]. It is dominant over the others and produces the thrust which counterbalances the drag generated by the resistive forces of the Taylor-like model.

Steady-state turning experiments

The second set of tests concerns steady-state turning. We started by studying the external dynamics by assessing how the turning parameter α of (21) influences the curvature of the robot trajectory when it is stabilized in its steady-state motion. To that end, 16 simulations were performed, with 5 experiments per simulation, leading to 80 experiments, for 4 shape motions (21) with (ν, A, k) as displayed in table 2, and 4 values of α ranking from 0.1 to 0.4 (∆α = 0.1). The swimming frequency ν is fixed to 0.6 Hz (i.e. T = 1.67 s) in all cases. For each set of parameters, both the simulated or real robots, start from rest, and are stopped after 10 s (i.e. 6T ) of swimming or when its head collides one of the pool's walls. The curvature radius of the robot trajectory denoted by R a is then computed from the video and simulation data by applying a least squares fitting. The discrepancy between the curvature radius observed in the simulation and the experiment reflects the errors introduced by the model of the lateral hydrodynamic forces (perpendicular to the robot axis) exerted on the robot when its joint dynamics are in their steady-state phase. In contrast to the non-turning experiments in which symmetries average out some of the discrepancies, when turning, any modelling error on the lateral forces can generate a net lateral drift which affects the turning radius. Figures 11 (a)-(d) show the radius of the robot trajectory as a function of the turning parameter α, for the 4 sets of parameters in table 2. The overall tendencies observed in the simulations and experiments are similar. However, in this case the robot steers into its own wake where big vortices are being advected by the flow. As a result, the fluid laterally bounding the robot is far from quiescent and even further from being a potential flow (as postulated by the LAEBT), which probably explains why the drift in turning experiment after one cycle is more significant compared to the forward swimming experiment.

To complete the discussion, a sequence of 10 snapshots of the simulated and real robot sampled along three of the previous turning trajectories over one period (between 3T and 4T ) is displayed in figure 12 in the same way as for the non-turning gait illustrated in figure 8. It can be seen that the simulated and real shapes match as well as in non-turning motions, at least when α is sufficiently small. When α increases (and the turning radius decreases), some discrepancies appear about half way through the cycle. Since this moment at which, the actuators' torque attains its maximal value, these discrepancies are probably due to the nonlinearities introduced by the actuators (saturations) and gears (nonlinear stiffness). However, a good matching of the shape is obtained at the beginning and the end of each cycle. This good match is used in figure 12 to translate and rotate the head of the simulated robot in order to align its configuration with that of the AmphiBot at the instant t = 3T and then let the simulation evolve as displayed in figure 12. Although the external dynamic tendencies are qualitatively the same for the simulation and the robot, the net position of the model progressively drifts with respect to reality. The discrepancies we are interested in measuring are introduced by a dynamic model, i.e. they first affect the computation of accelerations. As a result, because accelerations are integrated twice to reconstruct the position and shape of the robot in space, the discrepancies shown in figure 12 are mostly due to this integration which inevitably accumulates the dynamic errors and amplifies the discrepancies with time. Finally, note that this drift is amplified by increasing α (i.e. decreasing the turning radius). This is probably due to the fact that in this case, because the robot swims into its own wake, Lighthill's basic assumptions on the flow are less and less justified. The shape dynamics introduce other errors when the turning radius is small.

Transient experiments

The behavior of the model during transients is also important if it is to be used for control, since a control law will constantly push the robot in such transient phases in order to react to unmodeled perturbations or to track the time varying control inputs. For the purpose of illustration, figure 13 shows snapshots of the robot and the simulator for one of the forward-swimming gaits and one of the turning maneuvers previously presented. The AmphiBot III is controlled by a connection of central pattern generators (CPGs), i.e. a set of coupled nonlinear oscillators organized in pairs serially connected along the body length, whose outputs define the time-evolution of the desired joint angles r d of the control law (22). According to [START_REF] Crespi | Online optimization of swimming and crawling in an amphibious snake robot[END_REF]), the state of each oscillator is governed by the following dynamic system: with (φ i , θi , θ i ) the state vector of the i th oscillator and x i its output. For each joint, a pair of left-right oscillators generate the desired actuator angle through the relation r d,j = x l j -x r j , l and r denoting the left and the right oscillator of the joint while the parameters w ij , D i , ν i and ψ ij are chosen in order to ensure the asymptotic convergence of the desired actuator angles toward the desired robot shape motion (21). This explains why in all the steady-state tests previously presented the motions have been considered after several cycles (namely 3 cycles), in order to eliminate the transient phenomena from the observations. By contrast, in the transient experiments illustrated in figure 13, motion is recorded from the beginning when the robot is at rest. In these experiments, the dynamic model of the coupled oscillators ( 23) is coupled to the simulator of the robot with the parameters of the AmphiBot III defined as follows. For any connection between oscillators, we have w ij = 10, while between two consecutive pairs of oscillators, we have ψ ij = (2πk/m) or ψ ij = -(2πk/m) depending if the connection is ascending or descending. Between two oscillators of a same pair, we have ψ ij = π, while ν i = ν, D l i = (A -αA)/2, D r i = (A + αA)/2 and a i = 20, for any oscillator. The results illustrated in figure 13 show that even in the context of transients, the model qualitatively (and quantitatively) still matches reality well. In particular, there is no significant discrepancy of net position and orientation occurring in the first few moments of motion although in this case, the velocities being low, the Reynolds number characterizing the flow is smaller than that required by Lighthill's reactive model.

     φi = 2πν i + j θ j w ij sin(φ j -φ i -ψ ij ) , θi = a i a i 4 (D i -θ i ) -θi , x i = θ i (1 + cos φ i ) ,

Conclusion

This article presents the dynamic model of a fish-like robot named AmphiBot III. The approach is based on the adaptation of J. Lighthill's large amplitude elongated body theory to a serial mobile multibody system. Based on this model, we used a fast Newton-Euler based algorithm which solves the forward dynamics of the robot. This algorithm, which is a generalisation of the Featherstone algorithm to mobile multi-body systems allows the robot motion to be integrated in real time on a standard computer. The results of this integration are compared to the motion of the real robot for forward swimming gaits and turning manoeuvres. In all cases, the model of internal dynamics, whose parameters correspond to known or measurable mechatronic quantities of the real robot, shows moderate discrepancies (< 20%) with respect to experiments. When turning, in spite of the highly vortical flow, the results show a good match between the simulation and observation in terms of curvature radius. Finally, all these results have been obtained by calibrating only 3 dimensionless hydrodynamic parameters (i.e. C m , C f and C d ) whose values were fixed only once. They accord well with the physics, confirming the LAEBT as a good solution for fish-like robot modelling. Furthermore, beyond its application to the AmphiBot III and any anguilliform fish-like robot, it has a more universal value which should enable the modelling of a wide diversity of robot morphologies including those inspired by caranguiform and thuniform fish or more widely any fish that uses body caudal fin swimming [START_REF] Breder | The locomotion of fishes[END_REF]. It would suffice to adapt the cross-sectional added inertias (and drag coefficients) to each specific case. Similarly, adding roll and vertical cross-sectional added inertias and drag forces would allow the present model to be extended to the three-dimensional swimming [START_REF] Boyer | Poincaré-Cosserat equations for the Lighthill three-dimensional large amplitude elongated body theory: application to robotics[END_REF][START_REF] Candelier | Three-dimensional extension of the Lighthill's large amplitude elongated-body theory of fish locomotion[END_REF]. The numerical efficiency and the modular character of the programming of the model could prove a useful tool for robot design, gait optimization and control approaches that require an online model [START_REF] Porez | Multi-physics model of an electric fish-like robot : numerical aspects and application to obstacle avoidance[END_REF][START_REF] Morel | Nonlinear motion control of cpg-based movement with applications to a class of swimming robots[END_REF]Morel et al., , 2012b,c,a],c,a]. A simplified version of this model (with no caudal term) has already proved useful in controlling the RAAMO eel-like robot (see [START_REF] Boyer | The eel-like robot[END_REF] for more details). More specifically, the model has been used: 1) as a fast simulator to develop 3-D swimming controllers with and without pectoral fins [START_REF] Alamir | Feedback design for 3-D movement of an eel-like robot[END_REF]El-Rafei et al., 2008a]. In particular, it has been used as a real-time interface to control our eel-like robot in a virtual world using a joy-stick [El-Rafei et al., 2008a]; 2) as a simulator to identify a highly-speed black box model [START_REF] El-Rafei | Multi-variable constrained control approach for a three-dimensional eel-like robot[END_REF]. Finally, the modelling framework developed for the RAAMO robot was structured in three hierarchical stages. The first level was a Navier-Stokes resolver [START_REF] Leroyer | Numerical methods for RANSE simulations of a selfpropelled fish-like body[END_REF] used to assess the second stage based on the LAEBT [START_REF] Boyer | Fast dynamics of an eel-like robot -comparisons with Navier-Stokes simulations[END_REF], itself used to identify the third (black box) stage. This last high-speed model was used to integrate a predictive controller in a receding horizon feedback schema [START_REF] El-Rafei | Multi-variable constrained control approach for a three-dimensional eel-like robot[END_REF]. Beyond the RAAMO project, the model has been used in a same manner (fast simulator and identification) to control the AmphiBot III [START_REF] Morel | Nonlinear motion control of cpg-based movement with applications to a class of swimming robots[END_REF][Morel et al., , 2012a]]. Most of these results were restricted to simulation but in future we hope to use this model to control the AmphiBot III.

Acknowledgments

The ANGELS project is funded by the European Commission, Information Society and Media, Future and Emerging Technologies contract number 231845. The authors gratefully acknowledge Alessandro Crespi, Camille Joggi, Steve Berger, Yannick Morel and Douglas Carnall for their help.

8 Appendix: Proof of (10)

This appendix aims at deriving the expression (10) of the reactive force F reac,j applied on the segment j. This derivation is based on the following simple remark. If we apply the kinetic momenta balance to the fluid contained in the control volume between the two boundary planes of the segment j, i.e. S j (0) and S j (l j ) (see figure 4), we find the time rate of kinetic momenta of this piece of fluid on the first term of the balance, and on the second term, the pressure forces applied across S j (0) and S j (l j ) along with those applied by the lateral boundary of the segment j. Therefore, this third force being the opposite (by action/reaction principle) of F reac,j , and since the pressure forces are known functions of the segment state given by ( 9), one just needs to calculate the time rate of kinetic momenta to finally find F reac,j .

From perfect fluid theory, we cannot compute the kinetic momenta of the fluid contained in the volume control by using the basic definitions of mechanics 17 since they lead to divergent integrals contrary to the kinetic energy which is a convergent quantity [START_REF] Batchelor | An introduction to fluid dynamics[END_REF]. As a result, the fluid kinetic momenta are defined as the linear and angular impulses [START_REF] Lamb | Hydrodynamics[END_REF], i.e. the conjugate linear and angular momenta computed from the kinetic energy of all the fluid slices between S j (0) and S j (l j ). Applying the Kirchhoff potentials approach [START_REF] Kirchhoff | Ueber die bewegung eines rotationskörpers in einer flüssigkeit[END_REF] to the 2-D Laplace equations in S j (X), the kinetic energy of the fluid contained in S j (X), can be expressed as a simple quadratic form of the lateral velocities of the segment cross section C j (X):

Tf,j (X) = 1 2 η T Mf,j η (X) , ( 24
)
where η is the field of cross-sectional velocity along the segment j, given for X ∈ [0, l j ] by:

η(X) = (V T , Ω T )(X) = Ad ( X g j ) η j ,
with Ad X g j the adjoint map operator from the cross-sectional frame (O(X), s j , n j , a j ) centered on C j (X) to (O j , s j , n j , a j ) which is defined by:

Ad X g j = 1 3 -X ŝj 0 1 3 . (25)
In ( 25), the rotation part reduces to identity since the segments are rigid and straight. Moreover, in (24), we introduced the cross-sectional (2-D) added mass tensor:

Mf,j (X) = Mf,j 0 0 Īf,j (X) ,

where Mf,j and Īf,j represent the density of linear and angular cross-sectional added inertia per unit of segment length. Now computing the impulses as the conjugate momenta associated to (24), one obtains:

Pf,j Σf,j (X) = ∂ ∂V Tf,j ∂ ∂Ω Tf,j (X) = Mf,j V Īf,j Ω (X) ,
where Pf,j and Σf,j respectively denote the densities of linear and angular cross-sectional fluid impulses along the segment j. Before going any further, let us remark that in the Newton law and Euler theorem, the time-derivatives are applied to the momenta of fixed material particles. As a result, the application of these principles, requires the introduction of the particular derivative [START_REF] Batchelor | An introduction to fluid dynamics[END_REF] which can be written here as 18 :

d dt (f (X, t)) = ∂ ∂t (f (X, t)) -V X (X) ∂ ∂X (f (X, t)) ,
17 i.e. in terms of the fluid velocity field. 18 Let us recall that the particular derivative represents the variation in time as measured by an observer attached to a given material particle. In our case, since the material particles are replaced by slices, (26) stands for a 1dimensional version of the particular derivative. Hence, when applied to any function f (X, t), d./dt (respectively ∂./∂t) is the rate of change of time measured by an observer attached to the fluid slice Sj(X) (respectively, to the segment cross section Cj(X).

where V X (X) represents the axial (along the j th segment length) velocity of the X cross section (C j (X)) in the B j frame, i.e. V X (X) = V (X) T s j = (V j +Ω j ×Xs j ) T s j = V T j s j . As a corollary, for any function f of X and t related to the fluid stratification, we have the following 1-dimensional version of the Reynolds (transport) theorem of fluid mechanics (see [START_REF] Boyer | Poincaré-Cosserat equations for the Lighthill three-dimensional large amplitude elongated body theory: application to robotics[END_REF]):

d dt l j 0 f (X, t)dX = l j 0 ∂f ∂t dX + [V X f ] l j 0 . (26)
In ( 26), the first (integral) term represents the contribution brought by the instantaneous variation of f (X, t) between S j (0)) and S j (l j ), while the second (boundary) term represents the variation in time of f between the two slices brought by the flux of fluid across this control volume attached to the j th segment19 . In order to derive the expression of the (6 × 1) vector F reac,j = (f T reac,j , c T reac,j ) T of reaction forces applied by the fluid onto the segment B j we apply Newton's law and Euler's theorem to the fluid between the two geometric slices S j (0) and S j (l j ). That allows one to obtain:

• linear dynamics (from Newton's law) d dt (P f,j ) = -f reac,j + f pres,j (0) -f pres,j (l j ) , ( 27)

where:

P f,j = l j 0
Pf,j dX , is the linear kinetic momentum of the fluid between the two geometric slices S j (0) and S j (l j ) .

• angular dynamics (from Euler's theorem):

d dt (Σ f,j ) + V j × P f,j = -c reac,j + Xs j × f pres,j (l j ) ,

where:

Σ f,j = l j 0 ( Σf,j + Xs j × Pf,j )dX , is the angular kinetic momentum of the fluid between the two geometric slices S j (0) and S j (l j ).

Then, using ( 26) in ( 27) and ( 28), gives the N-E equations of the isolated fluid between S j (0) and S j (l j ), from which we deduce the expression (10) of (6 × 1) reactive forces.

List of Notations 1 3

The unit 3 × 3 matrix.

A

The joint motion amplitude.

A j

The unit vector supporting the joint axis j.

Ad j gj-1 The adjoint map operator allowing a change in velocity from F j-1 to F j .

C i

The hydrodynamic coefficients of the resistive model.

C d

The drag coefficient.

C f

The friction coefficient.

C m

The added mass coefficient.

F j

The force vector exerted by the body j -1 on to the body j expressed in F j .

F ext,j The vector of external hydrodynamic forces exerted on the body j.

F reac,j The vector of reactive forces exerted on the segment j.

F res,j The vector of resistive forces exerted on the segment j.

H j

The angular inertia of the body j around the axis A j .

I j

The tensor of angular inertia of the body j.

K d

The derivative gains of the joint controllers.

K e

The stiffness of passive joints.

K p

The proportional gains of the joint controllers.

L r

The total length of the robot.

M S j the tensor of first inertia moments of the body j.

M j

The tensor of body mass of the body j.

P

The body cross-section perimeter.

R a

The radius of the robot trajectory.

T

The swimming wave period.

T h

The hydrodynamic thrust.

V j

The linear Galilean velocity of the body j in F j .

V X The velocity field along the segment j.

j-1 P j The position vector of F j with respect to F j-1 .

j-1 R j The orientation matrix of F j with respect to F j-1 .

j-1 g j The transformation matrix of F j with respect to F j-1 .

Ω j

The angular Galilean velocity of the body j in F j .

α

The turning angle ratio.

α j

The reduced vector of Coriolis, centrifugal and external forces exerted onto the body j.

Pf,j

The field of linear kinetic momenta of the fluid along the segment j in F j .

Vf

The cruising forward swimming speed.

Σf,j

The field of angular kinetic momenta of the fluid along the segment j in F j .

Tf,j

The field of kinetic energy of the fluid along the body j.

β j The vector of Coriolis and centrifugal forces exerted onto the body j in F j .

β reac,j The vector of the reactive forces produced by the fluid added mass accelerated by the Coriolis and centrifugal accelerations of the segment j.

ηj

The time derivative of η j .

η j

The velocity of the body j in F j .

B j

The body j.

C j A cross-section of the body j.

F e

The fixed Galilean frame.

F j

The mobile frame attached to the body j.

K j

The reduced inertia tensor of the body j.

M j

The inertia tensor of the body j.

M reac,j The tensor of added inertia of the fluid accelerated by the segment j.

S j

The geometric slice of fluid which prolongs C j .

ν

The joint motion frequency.

φ i
The phase state of the oscillator i.

ρ b

The body mass ratio of the body j.

ρ f

The fluid mass ratio.

θ i

The amplitude state of the oscillator i.

ζ j

The part of the acceleration of the body j which depends on its velocity expressed in F j .

d j

The distance between the frame centers O j-1 and O j along s j-1 .

f pres,j The cross-sectional pressure forces field exerted on the body j.

k

The number of waves along the robot backbone.

m

The number of actuators.

n

The number of joints.

r j

The position of the joint j. r d,j

The set-point of the joint j.

t

The time.

u

The unit vector supporting the forward motion.

Figure 1 :

 1 Figure 1: From the left to the right, a global view of the snake-like robot AmphiBot III and an internal view of two robot segments.

Figure 4 :

 4 Figure 4: Schematic view of the fluid laterally surrounding a robot segment (i.e. when j is even).

Figure 5 :

 5 Figure 5: Schematic view of the direct dynamic model (DDM). The DDM computes the head segment accelerations and joint accelerations as a function of torques Γ, hydrodynamic forces F ext and the head state ( e g 0 , η 0 ).

Figure 6 :

 6 Figure 6: From the left to the right, a global view of the experimental setup and a schematic view of the measurement process.

Figure 7 :

 7 Figure 7: Forward swimming in a straight line: Swimming speed comparisons for all the tested forward swimming gaits as function of ν, A and k. The left, middle and right columns represent, from ν = 0.4 Hz (top row) to ν = 1.0 Hz (bottom row), the numerical results, the experimental results and the relative errors respectively. The white crosses represent the unstable gaits for which the speed could not be estimated.

Figure 8 :

 8 Figure8: Forward swimming in a straight line: The blue lines with circular markers represent the successive midline profiles of the AmphiBot reconstructed from the video tracking device by means of LED lights (blue circles) during a complete swimming cycle (between 3T and 4T ). The magenta lines with stars represent the midline of the simulated robot reconstructed numerically using the dynamic model presented (see section 4). The magenta stars are the positions of the simulated LEDs.

  ν = 1 Hz, A = 40 • , and k = 1.25.

  ν = 0.8 Hz, A = 35 • , and k = 1. ν = 0.6 Hz, A = 30 • , and k = 0.75.

  ν = 0.4 Hz, A = 25 • , and k = 0.5.

Figure 9 :

 9 Figure 9: Forward swimming in a straight line: Maximal lateral displacements (100 × |x|/L r ) of segments along the body (100 × |y|/L r ) measured with respect to the forward direction of swimming (dashed lines: robot; continuous lines: simulator; L r = robot length). The vertical error bars show the standard deviation of experimental measurements.

  ν = 1 Hz, A = 40 • , and k = 1.25. ν = 0.8 Hz, A = 35 • , and k = 1. ν = 0.6 Hz, A = 30 • , and k = 0.75. ν = 0.4 Hz, A = 25 • , and k = 0.5.

Figure 10 :

 10 Figure 10: Forward swimming in a straight line: Magnitude of the resultant of the different components of the hydrodynamic thrust for the motions of figure 8, including the reactive body component (+++), the caudal reactive component (---), the resistive transverse drag component ( ), and the resistive viscous axial component (• • •).

  ν = 0.6 Hz, A = 25 • , and k = 0.5. ν = 0.6 Hz, A = 35 • , and k = 0.5. ν = 0.6 Hz, A = 25 • , and k = 1.0.

  ν = 0.6 Hz, A = 35 • , and k = 1.0.

Figure 11 :

 11 Figure 11: Steady state turning manoeuvre: The radius of the robot trajectory R a is plotted as a function of the turning parameter α for the four desired shape motions defined in table 2. Dashed line: experimental results. Continuous line: simulation results. The vertical error bars indicate the standard deviation of experimental measurements.

  ν = 0.6 Hz, A = 25 • , k = 1.0 and α = 0.1. ν = 0.6 Hz, A = 35 • , k = 0.5 and α = 0.2. ν = 0.6 Hz, A = 25 • , k = 0.5 and α = 0.4.

Figure 12 :

 12 Figure 12: Steady state turning manoeuvre: The blue lines with circular markers represent the successive midline profiles of the AmphiBot reconstructed from digitized video fields by means of LED lights (blue points) during a complete swimming cycle (between 3T and 4T ). The magenta lines with stars represent the midline of the simulated robot reconstructed numerically using the dynamic model presented. The red spots indicate the head position of both the simulated and real robot. Wave form label A (deg) ν (Hz) k 1 25 0.6 0.5 2 25 0.6 1.0 3 35 0.6 0.5 4 35 0.6 1.0Table 2: Parameters of (21) used for the turning manoeuvres.

  ν = 1.0Hz, A = 30 • , and k = 1.0. ν = 0.6Hz, A = 35 • , k = 1.0 and α = 0.2.

Figure 13 :

 13 Figure13: Transient phases of a non-turning (left) and a turning motion (right): The blue lines with circular markers represent the successive midline profiles of the AmphiBot (blue dots) over 3T s starting from rest with a snapshot taken at each period. The magenta lines with stars represent the midline of the simulated robot which now includes a model of the CPGs. The red dots indicate the head position of both the simulated and real robot.

Table 1 :

 1 Thus, it is an open chain of n = 14 angular joints and 15 rigid bodies Parameters of the AmphiBot III dynamic model.

	Parameter Unit	j is odd	j is even	j = 14
	l j	m	-	0.097	0.2
	d j	m	0.097	-	-
	m j	Kg	-	0.218	0.248
	ms j I j	Kg.m Kg.m 2	-2.2 × 10 -3 a j a T j	6.0 × 10 -3 s j 3.37 × 10 -4 a j a T j	9.7 × 10 -3 s j 8.35 × 10 -4 a j a T j
	K e	Nm/rad	-	1.5	1.5
	K p	Nm/rad	50	-	-
	K d ρ b ρ f	Nm.s/rad Kg/m 3 Kg/m 3	0.25 --	-1000.0 1000.0	-1000.0 1000.0
	C m	-	-	0.5	0.5
	C f	-	-	0.006	0.006
	C d	-	-	1.75	1.75

In fluid dynamics, the Reynolds number is a dimensionless number characterising the relative influence of the viscous and inertial forces.

An edge on a body boundary generates a singular potential flow whose discontinuities can model the regions of high shearing strains of the real flow.

For instance, the vorticity shed at a sharp trailing edge can be modeled by a sheet of singularities (of dipoles or vortex) whose strength is deduced at each step of time by a Kutta-like condition[START_REF] Saffman | Vortex dynamics[END_REF], and whose the advection in the wake is governed by Helmoltz laws[START_REF] Milne-Thomson | Theoretical hydrodynamics -fifth edition[END_REF].

whose a meeting with him initiated Lighthill's works on swimming.

This robot is part of a series of amphibious robots that was originally founded by the SwissNational Science Foundation and by the European commission (FP7).

see http://biorob.epfl.ch/

Note that it is possible to attach a compliant fin instead of the rigid one.

By definition, the flow of a perfect fluid (non viscous and incompressible) is said potential if its vorticity field is zero everywhere. In this case the velocity is the gradient of a potential field φ governed by Laplace equations.

Which takes its origin in the slender body theory of Munk[START_REF] Munk | The aerodynamic forces on airship hulls[END_REF].

They come from the boundary terms of (26) in Appendix (see Section 8).

According to the slender body theory, they only move laterally, i.e. in the planes perpendicular to the segment axis.

In order to improve the LED light tracking, all the experiments were conducted in darkness in order to minimise the ambient light reflections on the water's surface.

By the conventions of fluid mechanics, this second term is named "convective" since it is produced by the advection of f by the flow sweeping past the body.