Frederic Boyer
email: frederic.boyer@imt-atlantique.fr

Vincent Lebastard
email: vincent.lebastard@imt-atlantique.fr.

Fabien Candelier
email: fabien.candelier@univ-amu.fr.

F Renda

Forward dynamics of continuum and soft robots: a strain parametrization based approach

In this article we propose a new solution to the forward dynamics of Cosserat beams with in perspective, its application to continuum and soft robotics manipulation and locomotion. In contrast to usual approaches, it is based on the non-linear parametrization of the beam shape by its strain fields and their discretization on a functional basis of strain modes. While remaining geometrically exact, the approach provides a minimal set of ordinary differential equations in the usual Lagrange matrix form that can be solved with standard explicit time-integrators. Inspired from rigid robotics, the calculation of the matrices of the Lagrange model is performed with a continuous inverse Newton-Euler algorithm. The approach is tested on several numerical benches of nonlinear structural statics, as well as further examples illustrating its capabilities for dynamics.

I. INTRODUCTION

A nonlinear Cosserat beam is defined by stacking along a single material dimension a continuous set of rigid cross sections of infinitesimal thickness [START_REF] Cosserat | Théorie des corps déformables[END_REF]. In contrast to the socalled Euler-Bernoulli, Rayleigh or Timoshenko models of strength of materials, in this non-linear theory, whose the most general instantiation is due to Reissner [START_REF] Reissner | On a one-dimensional large displacement finite-strain beam theory[END_REF], the cross sections can be subject to finite magnitude transformations [START_REF] Boyer | Poincare's equations for cosserat media: Application to shells[END_REF]. With the recent progresses of continuum and soft robotics [START_REF] Walker | Continuous backbone "continuum" robot manipulators[END_REF], this model is increasingly used today. In medical robotics, it is an efficient modelling tool for designing and analysing new concepts of highly redundant needles and endoscopes consisting of a passive elongated organ internally actuated through different technologies [START_REF] Webster | Nonholonomic modeling of needle steering[END_REF], [START_REF] Mahvash | Stiffness control of surgical continuum manipulators[END_REF], [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF], [START_REF] Chikhaoui | Comparison of modeling approaches for a tendon actuated continuum robot with three extensible segments[END_REF]. In bio-robotics, it has been used for modelling locomotion of continuum elongated animals as snakes, worms and fish with the aim of understanding the secrets of their performance and to reproduce them on bio-inspired artifacts [START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional eel-like robot[END_REF], [START_REF] Candelier | Note on the swimming of an elongated body in a non-uniform flow[END_REF]. More recently, taking inspiration from hydrostats as the octopus, the model of Cosserat beams has attracted attention of the young community of soft robotics for modelling arms made of soft materials as silicon, actuated with cables, pressurized chambers or electro-active materials [START_REF] Renda | Dynamic model of a multibending soft robot arm driven by cables[END_REF], [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF]. There exists a powerful numerical method specifically designed to simulate the Reissner beam theory, namely the Geometrically-Exact Finite Elements Method (GE-FEM) initiated by J.C. Simo in the 80s [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF], and developed by successors as [START_REF] Cardona | A beam finite element non-linear theory with finite rotations[END_REF], [START_REF] Ibrahimbegovic | Finite rotations in dynamics of beams and implicit time-stepping schemes[END_REF], [START_REF] Boyer | Finite element of slender beams in finite transformations: A geometrically exact approach[END_REF], [START_REF] Sonneville | Geometrically exact beam finite element formulated on the special euclidean group se(3)[END_REF]. However, with one exception [START_REF] Grazioso | A geometrically exact model for soft continuum robots: The finite element deformation space formulation[END_REF], it seems that this numerical method does not spread in the robotics community. One of the major reasons of this is probably that the GE-FEM is heavy, and practically difficult to implement for researchers not involved in non-linear structural dynamics. Furthermore, even in this research field, it is not conventional since it requires to cope with the differential geometry of the Lie group SO(3) [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF], its Stieffel fibration S 2 × SO(2) [START_REF] Boyer | Finite element of slender beams in finite transformations: A geometrically exact approach[END_REF], or directly SE(3) [START_REF] Sonneville | Geometrically exact beam finite element formulated on the special euclidean group se(3)[END_REF]. Another obstacle to the penetration of the GE-FEM in robotics, is due to the fact that the beam configuration being parameterized with the absolute pose of its cross-sectional frames, when imposing internal kinematic constraints (as the inextensibility condition), or external ones (e.g. those of a joint between two beams), the dynamics turn out to be ruled by stiff Differential-Algebraic-Equations (DAEs), which are not easy to numerically solve, and even more difficult to exploit for control purposes. Alternatively to the FEM, the oceanic engineering community have addressed in the 90s the forward dynamics of undersea cables using the shearing-free Cosserat model of Kirchhoff rods [START_REF] Burgess | Bending stiffness in a simulation of undersea cable deployment[END_REF], [START_REF] Sun | Modelling and simulation of low-tension oceanic cable/body deployment[END_REF]. In this other context, implicit finite difference schemes in time are used to set the dynamics of the cable in the form of a nonlinear Two Point Boundary Value Problem (TPBVP) in space, that can be changed into a set of algebraic equations with a space finite difference scheme, and numerically addressed with a finite-difference solver in space-time [START_REF] Burgess | Bending stiffness in a simulation of undersea cable deployment[END_REF], or directly solved at each time-step with a shooting method [START_REF] Sun | Modelling and simulation of low-tension oceanic cable/body deployment[END_REF], [START_REF] Gatti-Bono | Physical numerical modelling of the dynamic behavior of fly line[END_REF]. In both cases, the nonlinear TPBVP is linearized to be iteratively solved at each time step of a simulation. Recently, this shooting approach has been improved and applied to the real time simulation of several design of continuum robots [START_REF] Till | Real-time dynamics of soft and continuum robots based on cosserat-rod models[END_REF]. The computational graphics community has also extensively used the Cosserat model, with the aim of designing fast interactive simulations of rods able to be bent and twisted into complex shapes with knots. In this context, authors of [START_REF] Bergou | Discrete elastic rods[END_REF] have proposed an original spatial discretization of a Cosserat (Kirchhoff) rod by integral averaging its rotation field along edges connecting a set of discrete vertices. Coupled with discrete parallel transport of frames, this alternative approach is geometrically exact and has been shown recently to be a physical realistic solution for addressing several forward and inverse dynamics of Cosserat rods [START_REF] Gazzola | Forward and inverse problems in the mechanics of soft filaments[END_REF].

Despite their respective interest, being based on the same absolute parametrization, all these approaches do not provide dynamic models in the form of a minimal set of ODEs, as those governing the joint-parameterized dynamics of our usual manipulators [START_REF] Siciliano | Handbook of Robotics[END_REF]. This is in contrast to the so-called Piecewise-Constant-Curvature (PCC) parametrization, an approach which has rapidly spread in robotics [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF]. In this approach, the beam is first partitioned into a serial connection of pieces along which the curvatures are assumed to be constant and used as configuration parameters for the robot. Originally introduced with a further torsional component in [START_REF] Bertails | Linear time super-helices[END_REF], this kind of parametrization has been recently extended in the framework of Cosserat beams for soft robots through the Piecewise-Constant-Strain (PCS) parametrization [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF]. In this approach, the configuration of a Cosserat beam is no more described with the absolute poses of its cross-sectional frames, but rather with the relative poses between them, or in short, with its strain fields. Being constant, these fields can be used as parameters of a set of finite transformations similar to the joint dependent interbody transformations of a rigid manipulator. Owing to its relative nature, such a parametrization allows setting the dynamic model of soft robots in the usual form of rigid ones, i.e., it provides some minimal sets of ODEs suited to the design of model-based controllers as well as to the use of standard explicit time-integrators for simulation [START_REF] Renda | A geometric and unified approach for modeling soft-rigid multi-body systems with lumped and distributed degrees of freedom[END_REF]. Moreover, the approach allows modelling more easily the effects of actuation and proprioceptive sensing.

In the article here presented, we push forward the strain based parametrization of the PCS approach, by considering that the strains are no more constant along some pieces, but rather discretized on a truncated basis of strain functions according to the Ritz method [START_REF] Gander | From euler, ritz, and galerkin to modern computing[END_REF]. To assess the feasibility of this new parametrization, and in particular its ability to represent the statics and dynamics of non-linear Cosserat beams, we here apply it to the entire domain of a single beam which is kinematically free or clamped at each of its ends, and subject to any set of external and internal forces. While opening promising perspectives for the multi-beam and multi-section case, this first milestone can already be used in any application where a continuum robot can be represented by an actuated single beam. In particular, the internal constitutive law can encompass an activation component able to model the internal actuation of a soft arm [START_REF] Renda | Screw-based modeling of soft manipulators with tendon and fluidic actuation[END_REF], [START_REF] Till | Real-time dynamics of soft and continuum robots based on cosserat-rod models[END_REF], or the muscular activity of an elongated animal, as this is expected in biomechanics of fish swimming or snake creeping [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF]. As far as we know, the approach is also new for non-linear structural dynamics, where the discretisation of passive non-linear beams always concerns the positional (and rotational in the GE case) fields, and not the strain ones. It can easily capture in an explicit form, all the sub-models of the full Reissner model, as the shearing free Kirchhoff model, or that of inextensible Kirchhoff beams. For future application to bioinspired locomotion, we consider the general case of a beam free at its both ends, whose overall rigid motion is parameterized by the pose of its root frame. In this regard, the approach can be seen as a geometrically exact generalisation of the floating frame approach of flexible robots [START_REF] Canavin | Floating reference frames for flexible spacecraft[END_REF], where the root frame stands for the floating frame, and the deformation are no more linearly parameterized on a basis of vibrational modes, but rather non-linearly with a basis of strain modes. Such a ("modal") reduction is of great interest for robotics since it allows shifting from an infinite dimensional model to a finite dimensional one with few modes, to which one can apply control and optimization technics [START_REF] Siciliano | Handbook of Robotics[END_REF]. As a consequence of this new parametrization, two further difficulties arise. First, in contrast to the PCS [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF], the strain fields being not constant along the beam, the kinematics are no more analytically integrable. Second, due to the relative parametrization, all the geometric nonlinearities express through the mass operator which becomes full, highly nonlinear, and difficult to calculate. To address the first issue, the beam kinematics are numerically integrated with a quaternion-based integrator. Regarding the second issue, all the matrices of the inverse Lagrangian dynamic model of the beam are calculated with the Newton-Euler (NE) inverse dynamic algorithm of an equivalent continuum manipulator according to an idea originally proposed in [START_REF] Walker | Efficient dynamic computer simulation of robotic mechanisms[END_REF] for calculating the Lagrangian model of rigid systems. Using this inverse NE algorithm, the approach, which is at that time, not optimized in terms of computational time, is however very easy to implement.

The article is structured as follows. In section 2, the model of Cosserat beam is reminded, with its usual absolute parametrization. In section 3, we show how this model can be adapted to bioinspired, continuum and soft robotics by adding a model of the actuation to the passive constitutive law of the standard model. In section 4, the strain based parametrization of a Cosserat beam and its discretization on a Ritz basis are presented. Based on this relative parametrization, the forward dynamic model of an internally free-free actuated beam is provided in section 5. Its practical numerical computation is performed with a general algorithm presented in section 6. This general algorithm requires using several times the NE inverse dynamic algorithm of an equivalent virtual manipulator here derived in section 7. In section 8, the general algorithm of section 5 is applied to several numerical benches proposed in the past years to validate the GE FEM, as well as to other examples illustrating the approach in dynamics.

II. MODEL OF COSSERAT BEAMS

For the purpose of concision, we use a few of the Lie group notations of [START_REF] Murray | Robotic Manipulation[END_REF], and invite the reader not familiar with them, to refer to Appendix 1. We here consider an hyper-elastic beam subject to finite displacements and small strains. In the Cosserat approach, such a medium is modeled by a continuous set of rigid cross sections stacked along a material line parameterized by a coordinate X ∈ [0, 1] which plays the role of a continuous label for the cross sections. To each X-cross section, a mobile cross-sectional frame F(X) = (O, t 1 , t 2 , t 3 )(X) is attached, where O(X) and t 1 (X) coincide with the center of the cross section and its unit normal vector, respectively. In this context, the configuration of any cross section is defined by the action of an element of g ∈ SE(3) on a unique inertial frame F s . Hence, the beam configuration space is naturally defined as:

C = {g : X ∈ [0, 1] → g(X) ∈ SE(3)}, (1)
which stands for a functional space of curves in SE(3). Throughout the article, partial derivatives ∂./∂X and ∂./∂t are indicated by a "prime" and a "dot", respectively. The field g depending on both X and t, its space-time variations can be entirely described by the two vector fields η and ξ from [0, 1] to se(3) ∼ = R 6 (see Appendix 1 for notations):

η = (g -1 ġ) ∨ , ξ = (g -1 g ) ∨ , (2)
where

η = (Ω 1 , Ω 2 , Ω 3 , V 1 , V 2 , V 3 ) T = (Ω T , V T ) T stands
for the field of the velocity twists of the cross-sections in their mobile frames, while

ξ = (K 1 , K 2 , K 3 , Γ 1 , Γ 2 , Γ 3 ) = (K T , Γ T )
T is the exact geometrical counterpart of η when replacing t by X. By a slight abuse of terminology, this second field of twist will be often named the strain field, since any set of strain tensors has to depend on the current configuration g through this field only [START_REF] Boyer | Poincare's equations for cosserat media: Application to shells[END_REF]. For instance, if g o defines the reference configuration of the beam in which its internal strain energy is zero, a strain field is simply provided by ξ -ξ o with ξ o = (g -1 o g o ) ∨ . The velocity and strain fields η and ξ have their dual counterparts in se( 3

) ∼ = R 6 . Denoted Σ = (Υ 1 , Υ 2 , Υ 3 , P 1 , P 2 , P 3 ) T = (Υ T , P T ) T and Λ = (C 1 , C 2 , C 3 , N 1 , N 2 , N 3 ) T = (C T , N T ) T ,
these two further fields of wrench model the cross-sectional kinetic momenta and stress along the beam respectively. In the small strains assumption, these two pairs of dual fields are related by the (generalized) constitutive relations:

Σ = Mη , Λ = H(ξ -ξ o ), (3)
where M = diag(ρI 1 , ρI 2 , ρI 3 , ρA, ρA, ρA) and H = diag(GI 1 , EI 2 , EI 3 , EA, GA, GA) are two tensors fields in se(3) ⊗ se(3) ∼ = R 6 ⊗ R 6 , which stand for the crosssectional inertia and Hooke tensors along the beam. Applying Hamilton principle in this geometric context, provides the partial differential equation (PDE) of a Cosserat beam subject to a density of state-dependent external wrench F on ]0, 1[ and two tip external wrenches F 0 and F 1 at X = 0 and 1 respectively [START_REF] Boyer | Poincare's equations for cosserat media: Application to shells[END_REF]:

Σ -ad T η Σ = Λ -ad T ξ Λ + F , (4)
with boundary conditions at X = 0:

Λ(0) = -F 0 , or g(0) = g 0 , (5)
and at X = 1:

Λ(1) = F 1 , or g(1) = g 1 , (6)
depending whether, we impose the motion (through some known tip transformations g 0 and g 1 ), or the force at the tips. From (3), introducing Σ = Mη into (4) provides a PDE in acceleration η which governs the time-evolution of η and needs to be supplemented with the reconstruction equation (a kinematic model):

ġ = g η. (7)
Finally, the Cosserat formulation (3-7) stands for the model of Reissner beams [START_REF] Reissner | On a one-dimensional large displacement finite-strain beam theory[END_REF], which is at the basis of the geometrically exact finite element method of [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF]. In this case, (4) is replaced by its equivalent weak-form of virtual works, which once spatially discretized in finite-elements, is solved at each step of a time-loop with respect to a set of nodal velocities. Then the beam configuration is updated by using a discrete version of the reconstruction equation [START_REF] Rucker | Statics and dynamics of continuum robots with general tendon routing and external loading[END_REF], using the exponential map [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF].

Remark 1: In all the subsequent developments, we will consider that the internal beam kinematics which define the motions between cross sections can be constrained by some restrictions. For instance, if we neglect the transverse shearing and the extensibility, one has Γ = (1, 0, 0) T , and the beam kinematics can be reduced from the 6 internal dofs of a Reissner beam to the 3 internal dofs of those of a Kirchhoff beam [START_REF] Antman | Ordinary differential equations of nonlinear elasticity i: Foundations of the theories of non-linearly elastic rods and shell[END_REF]. To generically handle this kind of restriction, we can use the following decomposition of the field ξ:

ξ = Bξ a + Bξ c , (8)
where for any n a ≤ 6, ξ a (respectively ξ c ) is the n a × 1 vector field of the free strains allowed by the beam kinematics (respectively, the (6 -n a ) × 1 vector field of prohibited strains), while B and B are two complementary selection matrices of 1 and 0 such that B T B = 1, B T B = 0 and BT B = 1. Going further into details, ξ a , with the subscript "a" meaning "allowed", is a field of kinematic unknowns governed by the dynamics of the beam, while ξ c , with "c" for "constrained", is defined by some prescribed constants with respect to time. For instance, for an inextensible Kirchhoff beam, we have:

ξ a = (K 1 , K 2 , K 3 ) T , ξ c = (1, 0, 0) T , B = 1 3×3 0 3×3 , B = 0 3×3 1 3×3 . (9)
Finally, the decomposition (8) has its counterpart on the stress:

Λ = BΛ a + BΛ c , (10)
where Λ a is the field of stress components governed by the reduced constitutive law:

Λ a = B T Λ = H a (ξ a -ξ ao ), (11)
with H a = B T HB, the matrix of the reduced Hooke coefficients, and ξ ao the value of ξ a in the reference configuration. Finally, in [START_REF] Candelier | Note on the swimming of an elongated body in a non-uniform flow[END_REF], Λ c stands for a set of Lagrange multipliers fields (internal reaction forces and torques), in charge of imposing the internal constraints BT ξ = ξ c .

III. APPLICATION TO CONTINUUM AND SOFT ROBOTICS

Though being originally motivated by the study of passive beams in finite deformations, the above formulation can be easily adapted to model several systems relevant to bio-inspired and continuum robotics. Indeed, in several fields of robotics, continuum robots can be considered as internally actuated Cosserat beams [START_REF] Renda | Screw-based modeling of soft manipulators with tendon and fluidic actuation[END_REF], [START_REF] Till | Real-time dynamics of soft and continuum robots based on cosserat-rod models[END_REF]. Going further into details, modelling the body of the robot by the above model, and assuming that the inertia of the actuators are negligible with respect to those of this actuated body, the model of such an internally actuated Cosserat beam can be easily deduced from the above passive one, by replacing the reduced constitutive law [START_REF] Renda | Dynamic model of a multibending soft robot arm driven by cables[END_REF], by the active one:

Λ a = Λ ad (t) + D a ξa + H a (ξ a -ξ ao ), (12)
where Λ ad now stands for the autonomous time-dependent internal forces and/or torques produced by the actuation (in all the paper, the subscript "d" denotes some imposed forces defined by prescribed time-laws), while H a , and D a = B T DB model the viscoelastic properties of the actuated passive beam, all being compatible with the allowed internal beam kinematics. We now consider the model of section 2 with the active constitutive law [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF], and illustrate how it can be used for addressing several issues in robotics.

Imposing D a = H a = 0 in [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF], provides the dynamic model of a continuum multibody system constituted of an infinite number of bodies of infinitesimal length, namely the cross sections. Such a model has been used in [START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional eel-like robot[END_REF] for solving the continuum inverse dynamics of hyper-redundant fish robots. In the same context, ( 12) more generally provides a simplified model of the muscles inside an elongated animal, with Λ ad standing for the exogenous rhythmic stress field governed by a central pattern generator, while (D a , H a ) can be used for modelling a certain damping and stiffness of the muscles that are here mounted in parallel with the actuation [START_REF] Ekeberg | A combined neuronal and mechanical model of fish swimming[END_REF]. Note that in this first context, the Cosserat model is used as an idealization of the original continuum animals from which bio-robotics takes its inspiration.

The same law can be used in other fields, such as medical or soft robotics, where roboticists have designed elongated structures which are internally actuated by different types of technologies like pressurized chambers [START_REF] Marchese | Design, kinematics, and control of a soft spatial fluidic elastomer manipulator[END_REF], or routine cables [START_REF] Renda | Screw-based modeling of soft manipulators with tendon and fluidic actuation[END_REF]. In this latter context, the passive organ of the robot (e.g., a flexible rod along which discrete disks are arranged, a piece of silicon through which cables are passed...) can be modelled by a Cosserat beam (or a serial connection of several of them), i.e., with the above model, where in [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF], (D a , H a ) encode the viscoelastic properties of the flexible passive organ, while Λ ad captures the effects of the actuators on the system. In practical applications, the model of Λ ad has to be derived on a case-by-case basis. A possible derivation consists in shifting the actuation forces from the physical space, where they are primarily exerted, to the space of allowed internal strains ξ a , by using the kinematic-invariance of their internal virtual works.

To illustrate such a derivation, we consider the case of a cable driven soft arm as in [START_REF] Renda | Screw-based modeling of soft manipulators with tendon and fluidic actuation[END_REF]. This system is modelled by an inextensible Kirchhoff beam actuated by a set of frictionless cables of negligible inertia and extension indexed i = 1, 2...N , each of them being pulled in a controlled manner. Such cables are attached to the tip cross section X = 1 and routed through any other X ∈ [0, 1[, in such a way that in any beam cross sectional frame F(X), the i th cable crosses the section at a position: D i (X) = (0, D i,Y , D i,Z )(X). In these conditions, each cable has for positional field inside the beam: X → ϕ i (X) = r + RD i , and is able to transmit along its unit tangent ϕ i / ϕ i , a tension t → T i (t) ∈ R + only. The cable length and its variation are defined by:

l i = 1 0 ϕ i dX , δl i = 1 0 δϕ T i ϕ i ϕ i dX. ( 13
)
Introducing the expressions of ϕ i = r + RD i into (13), one can easily show that any variation δl i , can be expressed in terms of the beam shape variations δξ = (δK T , δΓ T ) T , as:

δl i = 1 0 (Γ + K × D i + D i ) T (δΓ + δK × D i ) Γ + K × D i + D i dX. (14)
The above identity defines the Jacobian from the kinematic space of fields δξ a = B T δξ to the discrete set of δl i s (which stands for the starting kinematic space of the cables). Therefore, using duality, i.e., invoking the kinematic-invariance of internal virtual works, we have for N cables:

δW = N i=1 δl i T i = 1 0 δξ T a Λ ad dX. (15)
Introducing ( 14) into [START_REF] Ibrahimbegovic | Finite rotations in dynamics of beams and implicit time-stepping schemes[END_REF] provides by simple identification, the expected vector of internal wrenches exerted on the Kirchhoff beam by a set of N cables:

Λ ad (X, t) = N i=1 B T D i × Γ c,i Γ c,i T i (t) Γ c,i = N i=1 D i (X) × Γ c,i Γ c,i T i (t). ( 16
)
where we used the B-matrix of an inextensible Kirchhoff beam [START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional eel-like robot[END_REF], with notations

Γ c,i = E 1 + K × D i + D i E 1 = (1, 0, 0) T ,
and where the subscript "c" here means "cable". Note that when K is constant over the beam, this field becomes uniform and has been obtained in [START_REF] Renda | Screw-based modeling of soft manipulators with tendon and fluidic actuation[END_REF] by other means. Finally, in all the subsequent developments we consider the model of section 2 with the full constitutive law [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF], and refer to it as that of an internally actuated Cosserat beam, or more concisely as that of a continuum robot.

IV. STRAIN PARAMETRIZATION OF A COSSERAT BEAM

In the original Cosserat model, the configuration space of a beam is defined as (1), i.e. each configuration of the beam is defined as the field of homogeneous transformations g parameterizing the pose of its cross-sectional frames with respect to a unique inertial frame. In [START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional eel-like robot[END_REF], we introduced an alternative definition of the configuration space of a Cosserat beam which is based on the following simple remark. If one knows at each instant the pose g(0)

g 0 of the root cross-sectional frame F(0) F 0 w.r.t. the inertial frame, as well as the field of spatial twist ξ along the beam length, one can then easily reconstruct the field of g along it, by forward integrating (i.e. from X = 0 to 1), the continuous geometric model:

g = g ξ. (17)
with initial condition g(0) = g 0 . Therefore, one can parametrize the Cosserat beam configuration by (g 0 , ξ), where remind that ξ is a field of 6×1 vectors in se(3) ∼ = R 6 . This leads to the second definition of configuration space of a floating beam:

C = SE(3) × S, (18)
where:

S = {ξ : X ∈ [0, 1] → ξ(X) ∈ R 6 }. (19)
Geometrically, such a space defines a principal fiber bundle SE(3) × S, where the fiber SE(3) here stands for the configuration space of F 0 , while the space S is the configuration space of its shapes, or "shape space". The beam being a continuous system, S is a functional space, i.e. an infinite dimensional space of X-parameterized curves in R 6 .

Practically, the reconstruction of the beam configuration in space can be easily performed by numerically integrating [START_REF] Sonneville | Geometrically exact beam finite element formulated on the special euclidean group se(3)[END_REF] as follows. Parameterizing the rotation matrix R of g = (R, r), with a unit quaternion Q ∈ R 4 , one can integrate [START_REF] Sonneville | Geometrically exact beam finite element formulated on the special euclidean group se(3)[END_REF] in the alternative form:

Q r = (2 Q ) -1 A(R(Q)K)Q R(Q)Γ , ( 20
)
starting from Q 0 (R 0 ) and r 0 , both deduced from g 0 = (R 0 , r 0 ). In [START_REF] Sun | Modelling and simulation of low-tension oceanic cable/body deployment[END_REF], we introduced the usual operator :

A(K) =     0 -K 1 -K 2 -K 3 +K 1 0 -K 3 +K 2 +K 2 +K 3 0 -K 1 +K 3 -K 2 +K 1 0     , (21)
as well as R(Q) ∈ SO(3) and Q(R) ∈ R 4 , which stand for the standard quaternion parametrization of rotations and its inverse map, respectively [START_REF] Murray | Robotic Manipulation[END_REF]. Now, referring to remark 1, for any Cosserat beam kinematics, one can adopt [START_REF] Grazioso | A geometrically exact model for soft continuum robots: The finite element deformation space formulation[END_REF] as definition of the configuration space, with now:

S = {ξ a : X ∈ [0, 1] → ξ a (X) ∈ R na }, (22)
which includes the general case of Reissner beams for which n a = 6. Finally, for the purpose of numerical resolution, we now discretize the field ξ a as:

ξ a (X, t) = ξ ao (X) + Φ(X)q(t), (23)
where Φ = (Φ 1 , Φ 2 , ..., Φ n ) defines a n a × n matrix of n basis functions or "shape functions" since they parameterize the shape space S of ( 18) with the n × 1 vector q of strain generalized coordinates. With this further reduction, the configuration space of a Cosserat beam simply reads:

C = SE(3) × R n . ( 24
)
This configuration space, which is still a principal fiber bundle, is well known of rigid robotics locomotion [START_REF] Kelly | Geometric phases and robotic locomotion[END_REF]. In this other context the finite dimensional shape space of ( 24) stands for the coordinates that parameterise the n joints connecting a set of n + 1 rigid bodies B 0 , B 1 ...B n , while the fiber SE(3) is the configuration space of an arbitrarily distinguished reference body B 0 .

Remark 2: In ( 23), remind that ξ ao (X) is the X-rate of variation of g on a reference configuration of the beam. As a result, if the beam is straight in its reference configuration, we have ξ ao = B T ξ o with ξ o = (0, 0, 0, 1, 0, 0) T . In this case the only non-zero component of ξ ao is brought by an allowed stretching field, as this is the case for an extensible Kirchhoff beam. More generally, if the reference configuration is curved, we need to shift from Cosserat beams to arches and (K 2o , K 3o ) will introduce further nonzero components in the time-independent field ξ ao .

V. DYNAMICS OF A CONTINUUM ROBOT PARAMETERIZED BY STRAINS

We can now build on the above geometric picture and seek for a dynamic model governing the configuration (g 0 , q) along time. To that end, one can reconsider the weak form of virtual works in the most general case where Λ a is defined by an active constitutive law of the type [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF]. A straightforward application of this principle in the configuration space [START_REF] Gazzola | Forward and inverse problems in the mechanics of soft filaments[END_REF] provides the dynamic model of a continuum robot in Lagrangian form:

0 Q ad = M 0 M 0 M 0 M η0 q + F v (q, q, η 0 ) Q v (q, q, η 0 ) + F c (q, g 0 ) Q c (q, g 0 ) + 0 K q + D q , (25)
which needs to be supplemented with the reconstruction equation ġ0 = g 0 η0 , i.e., a kinematic model allowing to recover the time evolution of g 0 from that of η 0 governed by the dynamic model [START_REF] Siciliano | Handbook of Robotics[END_REF]. Due to our definition (24) of the configuration space, all these matrices are organized in sub-blocks and the top row of the above matrix equality represents a balance of wrench (i.e. 6 × 1 vectors) w.r.t. to the reference frame F 0 , while the bottom one, stands for a balance of n × 1 generalized forces in the shape space (here parameterized by the strain coordinates of q). Going into details, in the right hand side of (25), we find from left to right, the matrix of generalized inertia of the Cosserat beam, the vector of its accelerations, that of all velocity-dependent forces, which include the Coriolis and centrifugal forces as well as some eventual velocity-dependent external forces as those exerted on a swimming robot. The next vector stands for that of configuration-dependent external forces as the gravitational ones, or a set of imposed time-varying forces appearing first in the two concentrated wrenches F 0 , F 1 and/or the wrench density F . Next, we find the generalized restoring and damping forces, where we have introduced the matrices of generalised stiffness and damping:

K = 1 0 Φ T H a ΦdX , D = 1 0 Φ T D a ΦdX. (26)
Finally, on the left hand side of (25), Q ad stands for the vector of generalised forces of internal actuation:

Q ad (t) = - 1 0 Φ(X) T Λ ad (X, t)dX. ( 27
)
Remark 3: The model ( 25) is enough general to encompass several cases in nonlinear structural dynamics, biomechanics and soft robotics. For instance, when the robot is a soft manipulator, we just have to remove the d.o.f.s of the floating frame F 0 (i.e., the top row of ( 25)). If one wants to consider the ideal case of a continuum rigid robot internally actuated, one has just to remove the restoring and damping forces.

On the other hand, if one wants to consider the case of a passive Cosserat beam, it suffices to remove the activation forces Q ad . Moreover, if in this latter case the stiffness matrix is zero, the same model can be used to describe a three dimensional continuum rigid chain swinging in gravity. More generally, using the full model ( 25) allows describing the locomotion dynamics of an elongated animal as a fish or a snake actuated by the continuous muscular model [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF]. Finally, for specific technological systems as those recently developed for soft robotics and medical needles [START_REF] Renda | Dynamic model of a multibending soft robot arm driven by cables[END_REF], [START_REF] Webster | Nonholonomic modeling of needle steering[END_REF], the right hand side of (25) remains valid, while the modelling effort is confined to the calculation of Λ ad in Q ad . For instance, if one reconsiders the example of the cable driven inextensible Kirchhoff beam of section 3, introducing ( 16) into [START_REF] Bertails | Linear time super-helices[END_REF] allows detailing the left hand side of (25) in the form:

Q ad = - N i=1 1 0 Φ T (D i × Γ c,i )dX T i (t) = LT (28)
where L(q) denotes the n × N actuation matrix, and T = (T 1 , T 2 ....T N ) T is the vector of control variables.

Remark 4: The strains being not subject to boundary conditions of geometric type (these conditions are imposed to the field of transformations g), the shape functions of [START_REF] Bergou | Discrete elastic rods[END_REF] are not required to satisfy any constraints at the boundaries of [0, 1]. Different functional basis can be used to discretize the stain fields as long as the discretization is set in the form [START_REF] Bergou | Discrete elastic rods[END_REF]. For instance, in case where n a = 1 the basis which comes immediately to mind is a polynomial basis of monomials of increasing order and we shall have:

Φ(X) = (1, X, X 2 ..., X n-1 ) P n (X). (29)
Note that the PCS parametrization of [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF] is a special case of [START_REF] Renda | A geometric and unified approach for modeling soft-rigid multi-body systems with lumped and distributed degrees of freedom[END_REF] in which all the non zero powers of X are removed. Continuing with polynomials one could give more structure to them and replace the above basis of monomials by an orthogonal basis of Tchebychev polynomials. In this case which will be numerically tested, the stiffness matrix ( 26) is changed into the identity matrix 1 na×na . Finally, when n a > 1, i.e. if one considers several allowed strains (possibly actuated), one has just to repeat the above scalar discretization for each of the strain components. For instance, if ξ a = (K 1 , K 2 , K 3 ) T , one can use [START_REF] Renda | A geometric and unified approach for modeling soft-rigid multi-body systems with lumped and distributed degrees of freedom[END_REF] for each of the three components and obtain the form [START_REF] Bergou | Discrete elastic rods[END_REF] with the block-diagonal matrix:

Φ(X) = diag(P k1 (X), P k2 (X), P k3 (X)), (30)
where k 1 + k 2 + k 3 = n, i.e. q = (q 1 , q 2 ...q n ), with (q 1 , ...q k1 ), (q k1+1 , ...q k1+k2 ) and (q k1+k2+1 , ...q n ) representing the torsional, and the curvature coordinates (in the planes (t 1 , t 2 ) and (t 1 , t 3 )) respectively.

Remark 5: Thanks to the strain parametrization, the stiffness and damping matrices of ( 26) as well as the actuation matrix of ( 28), are particularly simple and the two first can be calculated once for all outside the time-loop of a simulation. On the other hand, the explicit calculation of the other matrices is much more complex and in fact analytically unfeasible. Going further into details, when shifting from the g-parametrisation of section 2, to the strain parametrisation, the geometric nonlinearities are transferred from the model of the restoring forces to that of the inertial forces. As a result, while our stiffness matrices are easily computable, the mass matrices and velocity dependent forces are highly non-linear and a priori difficult to calculate. In the next section we propose a simple algorithm for calculating these matrices numerically. Inspired from past works in rigid robotics [START_REF] Walker | Efficient dynamic computer simulation of robotic mechanisms[END_REF], the principle of this algorithm consists in using an inverse Newton-Euler dynamic algorithm of the system to numerically reconstruct all the matrices of its Lagrangian model. In a further step, we will see how such an inverse NE algorithm can be derived.

VI. PRACTICAL CALCULATION OF THE FORWARD DYNAMIC MODEL OF CONTINUUM ROBOTS

To practically calculate all the matrices of ( 25), it is convenient, for reasons that will be detailed soon, to consider a fully actuated virtual rigid continuum manipulator associated to our free-free Cosserat beam. By "virtual manipulator", we mean that the reference section X = 0 is no more considered as being free but rather virtually actuated with an external wrench F 0 exerted on it. Moreover, this manipulator is said "rigid" and "fully actuated" since the field ξ a is directly actuated with no internal damping nor stiffness. In these conditions, ( 25) is simply changed into:

F 0 Q a = M 0 M 0 M 0 M η0 q + F v (q, q, η 0 ) Q v (q, q, η 0 ) + F c (q, g 0 ) Q c (q, g 0 ) . (31)
In short, [START_REF] Renda | Screw-based modeling of soft manipulators with tendon and fluidic actuation[END_REF] defines the Lagrangian model of the virtual continuum manipulator associated with the original beam, and except the damping and restoring forces, all the matrices of ( 25) are handled by [START_REF] Renda | Screw-based modeling of soft manipulators with tendon and fluidic actuation[END_REF]. As any Lagrangian model, ( 31) is an inverse dynamic model which is explicit w.r.t. the accelerations. Alternatively, it is well known from rigid robotics [START_REF] Khalil | Dynamic modeling of robots using recursive newton-euler techniques[END_REF] that the same dynamics can be described in an implicit form, thanks to an inverse Newton-Euler (NE) algorithm as that proposed by Luh and Walker [START_REF] Luh | On-line computational scheme for mechanical manipulator[END_REF] or by Featherstone [START_REF] Featherstone | The calculation of robot dynamics using articulatedbody inertias[END_REF] for rigid multibody systems. In the present case, such a Newton-Euler algorithm provides for any state (g 0 , η 0 , q, q) of the manipulator, the outputs F 0 and Q a that accelerate the manipulator with ( η0 , q). Therefore, this inverse NE algorithm defines an input-output map, here named IDM (for Inverse Dynamics Model), of the general form:

F 0 Q a = IDM (q, q, q, g 0 , η 0 , η0 ). (32)
Now, according to [START_REF] Walker | Efficient dynamic computer simulation of robotic mechanisms[END_REF], if one has in hand such an IDM , by feeding it with some specific known inputs, it becomes easy to calculate all the matrices of the Lagrangian model [START_REF] Siciliano | Handbook of Robotics[END_REF] at each instant of a simulation time-loop. Going further into details, let us define δ i ∈ R 6 (respectively δ α ∈ R n ) as a vector with all entries equal to zero except the i th one (respectively, the α th one), which is equal to 1. Then, using the definitions of δ i , δ α , simple identifications of ( 31) and ( 32) provides us with the following algorithm, where for any matrix A, A (k) denotes its k th column (with

k = i = 1, 2...6, or k = α = 1, 2...n):
• Step 1: Calculate:

F c Q c = IDM (q, 0, 0, g 0 , 0, 0) (33)
• Step 2: Calculate:

F v Q v = IDM (q, q, 0, g 0 , η 0 , 0) - F c Q c . (34)
• Step 3: For i = 1, 2...6, calculate:

M (i) 0 M (i) 0 = M 0 δ i M 0 δ i = IDM (q, 0, 0, g 0 , 0, δ i ) - F c Q c . (35)
• Step 4: For α = 1, 2...n, calculate:

M (α) 0 M (α) = M 0 δ α M δ α = IDM (q, 0, δ α , g 0 , 0, 0) - F c Q c . (36)
Finally, when all the matrices of ( 31) are computed, one can simply add [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF] and ( 27) to build the model [START_REF] Siciliano | Handbook of Robotics[END_REF], invert the mass matrix with the Cholesky algorithm, and deduce the forward dynamics of the original actuated free-free Cosserat beam in the explicit form:

d dt     q g 0 η 0 q     =     q g 0 η0 M 0 M 0 M 0 M -1 -F tot -Q tot     , (37)
where we used notations

F tot = F v + F c and Q tot = Q v + Q c + K q + D q -Q ad .
In the above system, the second equation (from top to bottom), is a reconstruction equation which can be integrated with the quaternion-based integration scheme (20) where t replaces X. Finally, time-integrating numerically (37) also provides the time evolution of q from which the configuration of all the cross-sectional frames can be reconstructed at any instant by integrating [START_REF] Sun | Modelling and simulation of low-tension oceanic cable/body deployment[END_REF] with ξ given by [START_REF] Chikhaoui | Comparison of modeling approaches for a tendon actuated continuum robot with three extensible segments[END_REF][START_REF] Bergou | Discrete elastic rods[END_REF].

Remark 6: Note that if we have to deal with a soft manipulator, i.e., an actuated Cosserat beam clamped in X = 0 and free in X = 1, referring to remark 3, its dynamics can be more simply written in the well known Lagrangian form :

Q a = M (q)q + Q v (q, q) + Q c (q), (38)
and the corresponding IDM takes the reduced form Q a = IDM (q, q, q). This reduced IDM which was originally introduced in [START_REF] Walker | Efficient dynamic computer simulation of robotic mechanisms[END_REF] to calculate the matrices of (38) for rigid manipulators is no more than [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF] where g 0 = 1, η 0 = η0 = 0. This confirms on the NE side, that the case of a free-free actuated beam includes that of soft actuated manipulators, and for this reason in all the following we shall continue to consider this more general context.

Remark 7: Finally, the above forward dynamic algorithm is conditioned by the obtention of a NE inverse dynamic algorithm IDM of a continuum manipulator in the form [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF]. Such an algorithm is presented in the next section. It's derivation is based on the generalisation of the Newton-Euler formulation from rigid multibody systems to continuum ones.

VII. INVERSE NE DYNAMICS OF CONTINUUM

MANIPULATORS

In [START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional eel-like robot[END_REF], a continuum NE inverse algorithm has been proposed for solving the inverse dynamics of a free-free Cosserat beam modelling a continuum hyper-redundant locomotor. This algorithm has been extensively used for the study of bio-inspired swimming [START_REF] Boyer | Fast dynamics of an eel-like robot-Comparisons with Navier-Stokes simulations[END_REF]. In this section, we remind it in the case of a continuum manipulator. In a second step, we use this continuous IDM to design the (discrete) IDM required by the forward dynamics algorithm of section 6.

A. Statement of the problem

We address the inverse dynamic problem of a fully actuated Cosserat beam, i.e., a beam actuated by a field of internal wrench Λ a and an external wrench F 0 . This problem can be stated as follows. Imposing any two sets (g 0 , η 0 , η0 ) and (ξ a , ξa , ξa ), find the field Λ a and the wrench F 0 that ensure them. This formally defines the continuous IDM as the input-output map:

F 0 Λ a = IDM (ξ a , ξa , ξa , g 0 , η 0 , η0 ), (39)
where let us remind that in contrast to [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF], the inputs ξ a , ξa , ξa and the output Λ a are fields, i.e., elements of some infinite-dimensional functional spaces.

B. Continuous inverse NE algorithm

For solving [START_REF] Kelly | Geometric phases and robotic locomotion[END_REF], we first remark that introducing Σ = Mη into the PDEs (4) gives:

Λ = ad T ξ Λ + M η -ad T η Mη -F , (40)
which are no more than the NE equations of a rigid continuum robot, where the cross sections stand for the discrete bodies, and the stress Λ, for the interbody wrenches. Furthermore, differentiating twice [START_REF] Sonneville | Geometrically exact beam finite element formulated on the special euclidean group se(3)[END_REF] provides the continuous version of the usual recursions on transformations, velocities and accelerations of the discrete NE model of manipulators [START_REF] Khalil | Dynamic modeling of robots using recursive newton-euler techniques[END_REF]. Exploiting this NE structure of the Cosserat model, one can easily derive a continuum version of the computed torque algorithm of rigid manipulators [START_REF] Siciliano | Handbook of Robotics[END_REF]. To that end, it suffices first to forward X-integrate (from X = 0 to X = 1), the continuum model of transformations, velocities and accelerations along the beam:

  g η η   =   g ξ -ad ξ η + B ξa -ad ξ η -ad B ξa η + B ξa   , (41)
which provides all the kinematic fields from the basis, where they are imposed, to the tip. Once all the kinematics are known over the entire beam, one can introduce them into [START_REF] Renda | Dynamic model of a multibending soft robot arm driven by cables[END_REF], that is X-integrated backward from the tip X = 1, where Λ(1) = F 1 is imposed, to X = 0, where Λ(0) is unknown reaction wrench. This second spatial integration provides the field of stress wrench Λ and finally the expected outputs of (39). This algorithm is summarized as follows and detailed in Appendix 2:

1) Forward X-integrate equations (41) initialized with (g 0 , η 0 , η0 ). 2) Calculate the state-dependent functions F 1 , F , ad T η (Mη).

3) Backward X-integrate (40) with initial condition

Λ(1) = F 1 . 4) Once Λ known along [0, 1], calculate the outputs:

Λ a = B T Λ and F 0 = -Λ(0).

Remark 8:

The above algorithm is nothing more than a continuous version of the inverse algorithm of Luh [START_REF] Luh | On-line computational scheme for mechanical manipulator[END_REF]. As its discrete homologous, it is structured into similar functional steps, where the 2 (1 forward and 1 backward) recursions on body indices of the discrete case, are replaced by 2 ODEs (1 forward and 1 backward) on the continuous label X. Therefore, one can easily hybridise the above continuum inverse algorithm with that of [START_REF] Luh | On-line computational scheme for mechanical manipulator[END_REF], in order to produce a more general IDM holding for systems connecting several soft and/or rigid bodies in a tree-like topology. To that end, the forward kinematic ODEs of the above algorithm (step 1) which holds along soft bodies, need to be composed with the usual discrete recursions on transformations, velocities and accelerations along rigid bodies. In the same manner, the backward ODEs on the internal stress wrenches of step 3 are composed with the discrete recursion on inter-body wrenches of the rigid elements. In these compositions, the final conditions provided by integrating an ODE, or iterating a recursion, need first to be changed of frame with one of the discrete operators g, Ad g or Ad T g , where g here denotes a discrete transformation at the junction between two bodies. Once so transformed, these final conditions feed the initial conditions of the next ODE or recursion. Repeating this process from body to body along the topology, provides the expected IDM which ensures the continuity of all the fields (transformation, velocities, accelerations and stress) throughout the entire structure.

C. Discrete inverse NE algorithm

Now that we have a realisation of the continuous IDM of a continuum manipulator, it is straightforward to build a discrete version of the same algorithm, which is consistent with our strain discretisation. Indeed, it suffices to impose [START_REF] Bergou | Discrete elastic rods[END_REF] along with:

ξa (X, t) = Φ(X) q(t) , ξa (X, t) = Φ(X)q(t), (42)
as inputs of the continuous IDM , and to project its output Λ a on the same functional basis according to:

Q a = - 1 0 Φ(X) T Λ a (X)dX. ( 43
)
This discretisation of the continuous IDM input-output map defines the discrete IDM of a continuum manipulator [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF].

Remark 9: Finally, applying the algorithm of section 6. with this discrete IDM allows reconstructing all the matrices of the Lagrangian model [START_REF] Siciliano | Handbook of Robotics[END_REF], except those of the generalised restoring and damping ones, whose expressions [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF] are indeed enough simple to be added a posteriori (see remark 5). Nevertheless, if one wants to automatically deduce the restoring and damping generalized forces with the reconstruction algorithm, it suffices to replace -Λ a , by -Λ ad +H a (ξ a -ξ ao )+D a ξa in the continuum IDM and to gather the restoring forces (respectively, the damping forces) with the configuration (respectively, the velocity) dependent ones, while Λ ad replaces Λ a a an output of IDM . This output being real if we deal with an internally actuated continuum robot, or fictitious, if we deal with a passive beam.

VIII. NUMERICAL TESTS

We first address a set of benches used in the past to evaluate the GE-FEM of Reissner and Kirchhoff beams [START_REF] Simo | A three-dimensional finite-strain rod model. Part II: Computational aspects[END_REF], [START_REF] Boyer | Finite element of slender beams in finite transformations: A geometrically exact approach[END_REF]. In these examples, the beam is considered as passive and modelled with the general constitutive law [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF] with no activation, and with D a = µH a , where µ is a coefficient tuning the internal damping. We close this section with the case of a cable driven beam free of external forces.

A. Static simulations

In static regime, the BCs of a beam need to impose at least 6 independent geometric constraints at its tips. Moreover, the beam needs to be subject to slow variations of external forces or/and geometric boundary conditions. In this case, the external d.o.f. along with the accelerations (including Coriolis and centrifugal ones) can be removed from the model [START_REF] Siciliano | Handbook of Robotics[END_REF] which turns out to be:

K q = -Q c (q, t), ( 44
)
where Q c here represents some generalized external forces which can be directly imposed or indirectly defined as some reaction forces (Lagrange multipliers) forcing a set of slowly moving geometric constraints. In particular, t in ( 44) stands for a fictitious numerical time that encodes the incremental variations imposed by these external forces or/and geometric constraints. The static balance ( 44) defines a set of nonlinear algebraic equations that can be solved with a Newton algorithm. However, this approach can be advantageously circumvented by adding an internal damping D artificially chosen, its role being to change the set of algebraic equations ( 44) into the set of ODEs:

D q + K q = c (q, t), (45)
which can be integrated explicitly with respect to the numerical time t. Note here that we are only interested in the asymptotic solutions of [START_REF] Simo | A three-dimensional finite-strain rod model. Part II: Computational aspects[END_REF], which are the roots of [START_REF] Boyer | Fast dynamics of an eel-like robot-Comparisons with Navier-Stokes simulations[END_REF] as soon as D is such that the dissipations have enough time to operate between successive incremental excitations.

In short, with a good accordance of the numerical damping matrix and the numerical time step ∆t, the beam will follow a succession of quasi-static equilibriums leading it to the expected static equilibrium configuration. In this section, ( 44) and ( 45) are named the undamped and damped static balances respectively. 1) Test 1: Clamped beam subject to a concentrated tip torque: This first test is very usual since it corresponds to one of the only analytical solutions of Cosserat rods in finite deformations. This exact solution was numerically recovered with the GE-FEM in [START_REF] Simo | A three-dimensional finite-strain rod model. Part II: Computational aspects[END_REF], [START_REF] Boyer | Finite element of slender beams in finite transformations: A geometrically exact approach[END_REF]. The beam is clamped in X = 0 and subject to an external imposed concentrated torque at X = 1. As a result, the forward algorithm has just to calculate Q c with [START_REF] Canavin | Floating reference frames for flexible spacecraft[END_REF], in which F = 0 and F 1 = F 1 (t) = (0, 0, C(t), 0, 0, 0) T . This calculation can be performed numerically with the general algorithm. However, in this particular case, and for the purpose of the illustration, one can entirely integrate by hand the undamped static balance [START_REF] Boyer | Fast dynamics of an eel-like robot-Comparisons with Navier-Stokes simulations[END_REF] as follows. First, using the relation ad T ξ = -Ad T g Ad -T g (see Appendix 1), the ODE ( 40) can be explicitly integrated with F = η = η = 0, which gives:

Q c (t) = - 1 0 Φ T B T Ad T g dXAd -T g(1) F 1 (t). ( 46
)
Second, since the load is planar, one can consider that only the curvature K 3 is excited and take B = (0, 0, 1, 0, 0, 0) T . In these conditions, adopting as functional basis [START_REF] Renda | A geometric and unified approach for modeling soft-rigid multi-body systems with lumped and distributed degrees of freedom[END_REF], the integral ( 46) can be detailed as:

Q c (t) = - 1 0 Φ T dXC(t). (47)
Which, after integration gives:

Q c (t) = -C(t)(1, 1/2, 1/3, 1/4....1/n) T , (48)
which defines the right hand member of the static balance (44) (we assume µ = 0). Regarding the left hand side of this balance, it is straightforward to calculate the stiffness matrix with [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF] for the strain basis [START_REF] Renda | A geometric and unified approach for modeling soft-rigid multi-body systems with lumped and distributed degrees of freedom[END_REF], which gives the static balance:

EI       1 1 2 1 3
..

1 n 1 2 1 3 1 4
..

1 n+1 1 3 1 4 1 5
..

1 n+2 .. .. .. .. .. 1 n 1 n+1 1 n+2
..

1 2n             q 1 q 2 . . q n       = C       1 1 2 1 3
.

1 n      
. Now, since the first column of the above n × n matrix is equal to EI/C times the vector of the right hand side, the above system has for solution q = (C/EI)(1, 0, 0...0) T . Moreover, this solution is unique since K is invertible as soon as all the basis functions are different. Finally, this simple analysis shows that the strain parametrization allows recovering the exact analytical solution of Euler which stipulates that a beam subject to a tip torque enrols around itself with a curvature radius of EI/C. with the general algorithm of section 6, fed with

F 1 = C(t)(e T c , 0 T 1×3 ) T
, where e c = (0, 1/ √ 2, 1/ √ 2) T is the unit vector around which the torque is applied incrementally from 0 to C = 5Nm with a step of C/k, and k = 5 a fixed number of steps. Note that in contrast to the above analytical solution, the tip torque being applied in an arbitrary tilted plane, we used a strain basis defined by Φ = diag(P k2 , P k3 ), with here k 2 = k 3 = 5. The (numerical) damping matrix is such that D = µK with µ = 10 -3 s. The numerical time step between two consecutive increments of load is set to ∆t = 0.02s. Figure 1 shows several equilibrium configurations in three-dimensions (3D). Figure 2 illustrates the typical first-order type convergency of K2 , K3 along loading steps, while figure 3 displays the evolution of the beam's tip position vs the magnitude of the tip torque. As expected, one recovers the solution of Euler in the tilted symmetry plane where the torque is exerted. 2) Test 2: Clamped beam subject to a follower tip force: This test was originally simulated with the GE-FEM applied to Reissner and Kirchhoff beams in [START_REF] Simo | A three-dimensional finite-strain rod model. Part II: Computational aspects[END_REF], and [START_REF] Boyer | Finite element of slender beams in finite transformations: A geometrically exact approach[END_REF] respectively. It consists in imposing to a beam, a follower tip force defined by F 1 = (0, 0, 0, 0, F (t), 0) T . In this case, the integral of ( 46) is no more g-independent and needs to be calculated numerically with the general algorithm of section 6 in the conditions of [START_REF] Canavin | Floating reference frames for flexible spacecraft[END_REF] where the configurationdependent forces are here only due to F 1 . In figures 4 and 5, we reported the results of a simulation of the damped balance [START_REF] Simo | A three-dimensional finite-strain rod model. Part II: Computational aspects[END_REF] obtained with the 5 first monomials of ( 29) as curvature basis. As in the example of [START_REF] Boyer | Finite element of slender beams in finite transformations: A geometrically exact approach[END_REF], we consider a circular cross-sectional beam with diameter 2R b = 0.57m, length l = 100m, and E = 4.015 × 10 8 Pa.

The simulation is performed by incrementing F (t) from 0 to 130kN. The numerical damping matrix is such that D = µK , with µ = 10 -3 s. The numerical time step between two consecutive increments of load is ∆t = 0.02s. Figure 4 displays the equilibrium configurations of the beam for several values of the tip force, while figure 5 shows the evolution of the beam's tip position versus the tip force magnitude. As expected, the results perfectly fit with the GE-FEM of [START_REF] Boyer | Finite element of slender beams in finite transformations: A geometrically exact approach[END_REF]. 3) Test 3: Clamped-clamped beam subject to imposed tip motions: We follow with a three-dimensional static test simulated with the GE-FEM of Kirchhoff rods in [START_REF] Boyer | Finite element of slender beams in finite transformations: A geometrically exact approach[END_REF]. This is a key test, since it allows qualifying the couplings between flexion and torsion which are responsible of most of the geometric nonlinearities in statics. It consists in twisting a rod, and then, to push its two tips toward each other. As this can be observed on the electric cables of our computers, the strain energy initially stored in the torsion state is so converted into flexural energy, while the cable forms a loop. To numerically reproduce this example, one considers an extensible Kirchhoff beam with the strain parametrization of remark 2, i.e.:

ξ = B(ξ ao + Φq) = B((0, 0, 0, 1) T + Φ(X)q), (49)
with Φ(X) = diag(P k1 , P k2 , P k3 , P k4 )(X). The beam is clamped at its two tips, but while in X = 0, F(0) is maintained fixed, in a first phase (Phase 1), one imposes F(1) to rotate around its first vector t 1 (1) of 2πrad, while maintaining its origin fixed. After this first phase, one enters into a second one (Phase 2) along which, while the orientation of F( 1) is maintained fixed, its origin (X = 1) is pushed toward the root X = 0 of a certain displacement. This imposed motion can be modelled by the geometric constraints on transformations:

g(1) = g 1 (t), (50)
where g(1) depends on q only, while g 1 (t) is a prescribed tip transformation consistent with the two phases of the test. Time differentiating (50) provides a model of the constraints in the kinematic form:

J 1 (q) q = η(1) = η 1 (t), (51)
where η 1 (t) = (g -1 1 (t) ġ1 (t)) ∨ is the prescribed time varying twist imposed at the moving tip of the beam, while J 1 is the 6 × n linear map encoding how the frame in X = 1 moves (with twist δψ(1) = (g(1) -1 δg(1)) ∨ ), when one imposes some arbitrary perturbations δq to the strain coordinates. This linear map J 1 can be easily calculated by remarking that the field δψ = (g -1 δg) ∨ produced by any δq, is solution of the ODE on η of (41) with η replaced by δψ, i.e.: δψ = -ad ξ δψ + BΦδq , δψ(0) = 0 6×1 .

(52)

Therefore, using ad ξ = Ad -1 g Ad g (see Appendix 1) allows one to integrate (52) explicitly and to get:

J 1 = Ad -1 g(1) 1 0 Ad g BΦdX, (53)
which, according to the duality between forces and velocities, is no more than the transposal of the matrix appearing in [START_REF] Boyer | Poincare cosserat equations for the lighthill three-dimensional large amplitude elongated body theory: Application to robotics[END_REF]. Now, if g 1 (t) varies enough slowly along time, the dynamic effects are negligible and the evolution of the beam can be described by the damped dynamics [START_REF] Simo | A three-dimensional finite-strain rod model. Part II: Computational aspects[END_REF], in which the (configuration dependent) external forces consist of:

Q c (q) = Q g (q) -J T 1 (q)Λ 1 , (54)
where Q g (q) is the vector of gravity forces and Λ 1 is a wrench of unknown Lagrange multipliers expressed in the cross sectional frame X = 1, whose role is to force the constraints (51). In these quasi-static conditions, the test can be modelled by the closed system of DAEs:

D q + K q = -Q g (q) + J T 1 (q)Λ 1 , J 1 (q) q = η 1 (t). (55)
Instead of using a Newton algorithm, (55) is here addressed with a projective method. To that end, one considers that q is at any time, the sum of a particular and an homogeneous solution of the kinematic constraints, i.e. :

q = J + 1 (q)η 1 (t) + H(q) qr , (56)
where J + 1 is the pseudo-inverse of J 1 , H(q) is a n × m matrix which gathers m (independent) basis vectors of the kernel of J 1 , while qr represents a reduced vector of admissible velocities (compatible with the constraints). Projecting (54) onto H, allows removing Λ 1 and to find:

qr = -(H T D H) -1 H T K q + Q g + D J + 1 η 1 (t)
, which once reintroduced into (56), provides the expected ODEs of the beam:

q = -H(H T D H) -1 H T (K q + Q g (q)) + 1 -H(H T D H) -1 H T D J + 1 (q)η 1 (t).
(57) In these equations, η 1 (t) can be considered as a kinematic feedback control of the form:

η 1 (t) = T SE(3) (Ψ 1 ) Ψ1 , with: Ψ1 = k p Ψ 1 , (58)
and:

Ψ 1 = Ψ 1 (q, t) = log SE(3) (g -1 1 (t)g(1)
), where T SE(3) and log SE(3) denote the tangent and the inverse of the exponential map Ψ → exp SE(3) (Ψ) in SE(3) [START_REF] Sonneville | Geometrically exact beam finite element formulated on the special euclidean group se(3)[END_REF], while k p is a diagonal matrix of proportional gains. When applied to (57) with g 1 (t) incremented step by step along the two phases of the test, the control (58) ensures (50). Note that at each new q, provided by the time-integration of q, one has to X-integrate forward the continuous geometric model g = g(Bξ ao + Φ(X)q) ∧ from g(0) = 1 4×4 , to find the new g(1) required by Ψ 1 in (58). Finally, the algorithm of section 6, is only used to calculate Q g , with [START_REF] Canavin | Floating reference frames for flexible spacecraft[END_REF] in which F = (0 1×3 , (R T (0, 0, γ g ) T ) T ) T and γ g = 9.81ms -2 .

Remark 10: Note that this algorithm only requires to compute log SE(3) and T SE [START_REF] Boyer | Poincare's equations for cosserat media: Application to shells[END_REF] . In contrast, the alternative Newton-algorithm, would require the further operators:

T -1 SE(3)
, and:

∂ ∂Υ T -T SE(3) (Υ)Λ 1 ∈ R 6 × R 6 .
In figure 6, some snapshots of the static equilibria are reported for several values of g 1 (t). The beam considered is the cross-sectional one of [START_REF] Boyer | Finite element of slender beams in finite transformations: A geometrically exact approach[END_REF] with for diameter 0.0176m, length l = 20m, E = 2.1 × 10 11 Pa, G = 4.1 × 10 7 Pa and volume mass ρ = 2.04 × 10 8 kgm -3 . This simulation is carried out with Tchebytchev polynomials, namely 5 for each curvature, 3 for torsion, and 2 for stretching [START_REF] Sun | Modelling and simulation of low-tension oceanic cable/body deployment[END_REF][START_REF] Sun | Modelling and simulation of low-tension oceanic cable/body deployment[END_REF][START_REF] Sun | Modelling and simulation of low-tension oceanic cable/body deployment[END_REF][START_REF] Cosserat | Théorie des corps déformables[END_REF][START_REF] Cosserat | Théorie des corps déformables[END_REF][START_REF] Cosserat | Théorie des corps déformables[END_REF]. With E 1 = (1, 0, 0) T , the prescribed tip motion is defined by g 1 (t) = (exp SO(3) (θ 1 (t) Ê1 ), E 1 ), for Phase 1, and g 1 (t) = (1 3×3 , x 1 (t)E 1 ), for Phase 2. In the two phases, θ 1 (t) (respectively, x 1 (t)) is incremented of 2π/k (respect., of -1.5/k) from 0 to 2π (respect., from l to -0.5l), with k = 18 (respect., k = 15), the number of tip pose increments. The numerical time-step between two consecutive increments of tip transformation is ∆t = 0.1s. As expected the torsional energy is converted into flexural one, while the snapshots fit very well (e.g. discrepancies on the torsional field K 1 do note exceed 0.1%) with those given by the GE-FEM of [START_REF] Boyer | Finite element of slender beams in finite transformations: A geometrically exact approach[END_REF].

(k 1 = 3, k 2 = k 3 = 5, k 4 = 2). Other simulation parameters are µ = 10 -3 s, k p = 10 2 diag

B. Cantilever Kirchhoff beam actuated with cables

In this section, we consider a cable driven Kirchhoff beam as evoked in the example of section 3. The beam is embedded in X = 0, free of any external force over its entire length, and internally actuated with one cable. To simulate the effect of imposing a given tension T along the cable, we need to integrate the damped static balance:

D q + K q = Q ad (q, t), (59)
with Q ad (q, t) = L(q)T (t) the vector of generalized internal actuation forces generated by the cable tension T ∈ R + and L defined by [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF]. Note that simulating (59) does not require using the algorithm of section 6. Indeed, owing to the strain-based parametrization, these matrices are easily computable from ( 26) and [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF]. Below, we reported the numerical solutions of (59) obtained with different routines of the cable defined by X → D(X) = (0, D Y (X), D Z (X)) T , with D Y and D Z taken as indicated in the table I. These routines, noted parallel, convergent, cross, spiral and parallel-truncated, are defined when the beam is in its reference straight configuration. TABLE I: Several routines of a one cable driven beam.

further modes), while due to the discontinuity, the paralleltruncated routine needs ten modes to converge.

C. Dynamic simulations

In all subsequent simulations we use the standard ODE45 integrator of Matlab, which allows one to time-integrate non-stiff ODEs.

1) Test 4: Clamped beam rolled up and released: In this first test, we start from the static equilibrium of the test 1, i.e., we first apply a concentrated fixed torque of magnitude C at the tip of a clamped beam. Then, this load is removed instantaneously at t = 0, and the beam freely moves with no gravity according to the dynamics: M (q)q + Q v (q, q) + Q c (q) + D q + K q = 0, (60) which are deduced from [START_REF] Siciliano | Handbook of Robotics[END_REF] in which the internal activation forces as well as the external dofs are removed (see remark 3 and 7). Except for restoring and damping forces, all the matrices of this model are computed with the full algorithm of section 6. The simulated beam is that of test 1 with volume mass ρ = 10 3 kgm -3 and µ = 10 -3 s. In figures 8(a,b,c), we reported a set of snapshots plotted each ∆t = 0.12s over the first second of simulation for C = 0.8, 2.3, 5.4Nm, while figure 9 shows the time evolution of the coordinates of the beam tip along time. In these figures, the strain basis is a Tchebychev polynomial basis with dimension n = 5. To illustrate the convergence of the Ritz decomposition, we use the following error between n and n + 1 modes:

(n + 1/n) = 1 t f t f 0 1 0 (K (n+1) 2 -K (n) 2 ) 2 dXdt, (61)
where

K (n) 2
stands for K 2 approximated with n strain functions, while t f is the duration of the simulation. Applying (61) to the 3 simulations, shows that (2/1) is always high, while increasing n makes (n + 1/n) rapidly drop up to be lower than 1% of (2/1). This rate of convergence decreases when the deformations increase. In short, convergence needs to increase n with the amplitude of deformations, while n = 5 is enough for three tests. This is confirmed by 2) Test 5: Three-dimensional continuum pendulum swinging in gravity: In this second dynamic test, we consider a string whose root X = 0 is clamped to a moving basis while its opposite tip is equipped with a lumped mass m 1 rigidly fixed to it. The string is modelled by an inextensible and non twistable Kirchhoff beam with EI 2 = EI 3 = 0. According to remark 10, the presence of the tip mass needs to feed the backward ODE (40) of the IDM with Λ(1) = F 1 = -M 1 η1 + F in,1 + F ext,1 , where M 1 and F in,1 + F ext,1 are the 6 × 6 inertia matrix of the tip body and the 6 × 1 wrench of inertia and external forces applied to it, both expressed in F(1). With a lumped mass we have M 1 = -m 1 diag(0, 0, 0, 1, 1, 1), and F in,1 + F ext,1 = m 1 (0 1×3 , (R T (1)(γ g , 0, 0) T ) T -(Ω 1 ×V 1 ) T ) T . Finally, this IDM can be used to calculate all the matrices of the Lagrangian model of the system, which in this case, reads:

M (q)q + Q v (q, q, η 0 (t)) + Q c (q) + D q + K q = -M 0 (q) η0 (t), (62)
where t → (η 0 , η0 )(t) here represent the time-evolution of the velocity and acceleration of the root X = 0 that excite the pendulum. The simulated string has a mass per unit length ρA = 7.85 × 10 -5 kgm -1 and a length l = 1m. The tip mass is m 1 = 0.03kg. In figures 11, 12 and 13, we reported simulation results obtained with g 0 (t) = (1 3×3 , r 0 (t)), and a r 0 (t) defined by a Lissajous curve r T 0 (t) = (0, a sin(ωt), a sin(ωt+ϕ)), with a = 0.15m, ω = 2πrad.s -1 , and ϕ = π/2. In these conditions, η 0 (t) = (0 1×3 , ṙT 0 (t)) T and η0 (t) = (0 1×3 , rT 0 (t)) T . We used 3 Tchebychev polynomials for K 2 and K 3 , which are the only strains allowed by the internal kinematics. Figure 11 displays several snapshots drawn every ∆t = 0.15s. In figure 12, the trajectories of the two tip cross sections are plotted in the y -z plane, while in figure 13 the coordinates of the tip section are plotted along time.

3) Test 6: Internally actuated floating beam in vacuum:

In this third dynamic test, one considers a free-free inextensible and non twistable Kirchhoff beam in vacuum, to which one imposes along K 2 and K 3 , an internal stress field governed by [START_REF] Trivedi | Geometrically exact models for soft robotic manipulators[END_REF] with an activation stress of the form:

X ∈ [0, 1] : Λ ad (X, t) = f r (t)f (2) w (X, t) f r (t)f (3) w (X, t) . (63)
Referring to anterior works on swimming [START_REF] Boyer | Poincare cosserat equations for the lighthill three-dimensional large amplitude elongated body theory: Application to robotics[END_REF], [START_REF] Candelier | Three-dimensional extension of lighthill's large-amplitude elongated-body theory of fish locomotion[END_REF], the law of each internal torque field is built with two spacetime functions. The first, is a sinusoid ramp time-function noted f r of the form:

t ∈ [0, t i [ : f r (t) = 0, t ∈ [t i , t f [ : f r (t) = t -t i t f -t i - 1 2π sin 2π t -t i t f -t i , t ≥ t f : f r (t) = 1, (64)
which is twice continuously differentiable and thereby guarantees smooth time transitions from 0 to 1. The second function, is a space-time wave function f

(i) w , i = 2, 3, defined ∀(X, t) ∈ [0, 1] × [0, t f ] by: f (i) w (X, t) = a i sin 2π X λ (i) - t T (i) , (65)
in which λ (i) is the wave length, T (i) is the period, and a i is the (constant) amplitude envelop for each of the two curvature components i = 2, 3. In these conditions, we need to consider the general dynamics (25) whose all matrices are calculated by the algorithm of section 6, except Q ad which is given by ( 27) with (63). The simulated beam has the geometry and viscoelastic properties of that of test 4.

In figures 14, 15 we reported the results of a simulation carried out with 3 Tchebytchev polynomials for each of the curvature fields. The internal law (63) is applied over 2s with t i = 0s, t f = 0.25s, (T (2) , T (3) ) = (1, 0.25)s, (λ (2) , λ (3) ) = (1, 0.5)m, a 2 = a 3 = 0.2Nm. In figure 14, snapshots every 0.3s are drawn, while figures 15 display the two actuated curvature fields ξ a = (K 2 , K 3 ) T along the beam at the same instants. For assessing the consistency of the simulation, we computed the kinetic wrench (σ T tot , p T tot ) T = 1 0 Ad -T g (Mη)dX, with p tot and σ tot the total kinetic resultant and momentum of the beam in the inertial frame. As expected, owing to the absence of external forces, these 6 first integrals are nearly preserved along the simulation ( σ tot in kg.m 2 /s, and the norm of the center mass position in m, are both lower than 10 -3 ). Note that these conservation properties need to be fulfilled for future application to bio-inspired locomotion.

IX. CONCLUSION AND PERSPECTIVES

In this article we explored a new way of parameterizing an internally actuated Cosserat beam. As the PCS approach of [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF], the configuration of the beam is not parameterized by the field of its cross-sectional frame poses, but rather by its strain fields. However, in contrast to the PCS, the strains are not assumed to be piecewise constant but are rather reduced on a basis of strain modes. This reduction of the original fields resembles the modal reduction of small deformations theories as the floating frame approach (FFA) of multibody systems. However, it has two main advantages over the FFA. First, the strain modes do not require to fulfill any specific boundary conditions (a constraint which is known to be one of the major drawbacks of the FFA [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF]). Second, in contrast to the FFA, it is naturally suited to systems subject to deformations of high magnitudes. Based on this alternative parametrization, we propose a new geometrically exact formulation of the statics and dynamics of an internally actuated Cosserat beam which is expected to be helpful for continuum and soft robotics. Exploiting the basic fact that the strains are the continuous counterparts of the joint angles of a discrete rigid articulated robot, this formulation has been set in the standard Lagrangian matrix form of rigid robotics. In particular, it provides a minimal set of ordinary differential equations that can be used for control design purposes, or for simulations with usual explicit time-integration schemes. Due to the relative nature of the strain-based parametrization, all the geometric nonlinearities of the Cosserat theory, are shifted from the forces related to the material (e.g. the restoring forces), to those related to space (e.g. the inertial ones), and the strain-based static formulation is much more simple than those based on the usual absolute cross sectional poses. In particular the stiffness matrix is now linear and shape independent. Thus, in statics, the approach is not only a method for constructing a Lagrange model of these systems, but also a powerful solution for their fast simulation. On the other hand, regarding dynamics, the mass matrix becomes full, nonlinearly shape dependent and difficult to compute. This issue is here tackled with a new algorithm that uses the Newton-Euler inverse dynamic model (NE-IDM) of a fully actuated beam as a computational machine for reconstructing all the matrices of its Lagrangian model. Using this algorithm for calculating inertial and external forces, the approach can be applied to bioinspired, continuum and soft robotics as well as to nonlinear structural dynamics. In this latter regard, it can be considered as an alternative to the geometrically exact method of finite elements (GE-FEM). To assess the feasibility of the approach, we numerically applied it to several benches used in the past to validate the GE-FEM in statics, and to others in dynamics or in statics, but with a beam internally actuated by pulling cables routed in its body. These first results show that with a few strain modes, one can reproduce easily the results of the GE-FEM while opening promising perspectives for statics of soft robots. In spite of these first encouraging results, this work opens several perspectives of improvement. In statics, the convergence can be improved by partitioning the beam into several sections on which the Ritz reduction would be applied. In the dynamic case, though in principle, the algorithm can be easily extended to multi-beam systems, its computational burden would practically compromise such perspectives. Several improvements could then be imagined for accelerating its computational efficiency in future. First, linearizing the NE-IDM provides a tangent algorithm that could be used (through the computation of a Jacobian) to integrate the dynamics with faster unconditionally stable implicit time-integrators. Second, while we here only considered polynomial strain basis, one could adapt such a basis to the specific systems of continuum and soft robotics (especially to their actuation), in order to reduce further the dynamic model.

X. APPENDIX 1: LIE GROUPS NOTATIONS

In the paper a hat "∧" covering a vector Υ defines a matrix Υ = Υ ∧ whose definitions depends on the dimension of Υ. Namely, if Υ ∈ R 3 , Υ denotes the 3 × 3 skewsymmetric matrix such that ΥV = Υ × V for any V ∈ R 3 . If Υ = (W T , U T ) T ∈ R 6 , with W and U two any vectors of R 3 , then Υ is the 4 × 4 matrix defined by:

Υ = Ŵ U 0 0 . ( 66
)
Reciprocally the superscript "∨" is such that Υ∨ = Υ for any Υ in R 3 or R 6 . We also used some of the geometric notations of rigid body mechanics [START_REF] Murray | Robotic Manipulation[END_REF]. They essentially consist of the two operators ad and Ad, defined for any g = (R, r) ∈ SE(3) and Υ = (W T , U T ) T , with W, U ∈ R 3 , by the two 6 × 6 matrices:

Ad g = R 0 rR R , ad Υ = Ŵ 0 Û Ŵ . (67)
Since in the paper all the vectors and tensors are related to the mobile cross-sectional frames, ad is used to denote the contribution of the time (Υ = η) and space (Υ = ξ) derivatives of the cross sectional basis vectors, when deriving a vector with respect to the fixed frame. More occasionally, we used the finite operator Ad which allows changing a twist (velocity) vector from one frame to another, their relative pose being defined by a transformation g = (R, r) ∈ SE(3).

Note that 1 6×6 + ad is the infinitesimal counterpart of Ad, and that for any s-parameterized curve of frames defined by the field of poses s ∈ R + → g(s) ∈ SE(3), we have dAd g /ds = Ad g ad Υ with Υ = (g -1 (dg/s)) ∨ the twist (srate of variation) of the frame in s, related to itself. Duality between twists and wrenches provide other useful relations with transposal matrices.

XI. APPENDIX 2: DISCRETE IDM OF A FREE-FREE INTERNALLY ACTUATED COSSERAT BEAM

The forward dynamics [START_REF] Ekeberg | A combined neuronal and mechanical model of fish swimming[END_REF] are time-integrated in statespace form with ODE45, which calls the IDM input-output map of the form [START_REF] Tytell | Interactions between internal forces, body stiffness, and fluid environment in a neuromechanical model of lamprey swimming[END_REF]. This map, which is a numerical function, is detailed below. As introduced in section 7, it consists in space-integrating in cascade 2 sets of ODEs. The first is a forward ODE. It is fed by the kinematic inputs (g 0 = (Q 0 , r 0 ), η 0 , η0 , q, q, q) and calculates through a Xintegration from 0 to 1, the kinematic fields along the beam. The second is a backward ODE initialized by the kinematics in 1. It calculates through a X-integration from 1 to 0, the field of internal wrench Λ, and finally the root wrench F 0 , and the generalized internal forces Q a , which are the outputs of the algorithm.

Forward ODE:

Inputs: Q 0 , r 0 , η 0 , η0 , q, q, q. ξ a = ξ ao + Φq, ξ = Bξ a + Bξ c ,

(K T , Γ T ) T = ξ, ξ = BΦ q, ξ = BΦq,     Q r η η     (0) =     Q 0 r 0 η 0 η0     .
Forward integrate:

d dX     Q r η η     =     (2 Q ) -1 A(R(Q)K)Q R(Q)Γ -ad ξ η + ξ -ad ξ η -ad ξ η + ξ     , Outputs: Q 1 = Q(1)
, r 1 = r(1), η 1 = η(1), η1 = η(1).

Backward ODE:

Inputs: Q 1 , r 1 , η 1 , η1 , q, q, q. ξ = B(ξ ao + Φq) + Bξ c ,

(K T , Γ T ) T = ξ, ξ = BΦ q, ξ = BΦq,         Q r η η Λ Q a         (1) =         Q 1 r 1 η 1 η1 F 1 0        
, Backward integrate:

d dX         Q r η η Λ Q a         =         (2 Q ) -1 A(R(Q)K)Q R(Q)Γ -ad ξ η + ξ -ad ξ η -ad ξ η + ξ ad T ξ Λ + M η -ad T η Mη - F -Φ T B T Λ        
, Outputs: (F 0 , Q a ) = (-Λ(0), Q a (0)).

XII. ACKNOWLEDGMENT

This publication is partly supported by the Khalifa University of Science and Technology under Award No. FSU-2018-08 and Research Excellence (AARE) 2018-05.

Fig. 1 :

 1 Fig. 1: Beam subject to a concentrated tip torque: Equilibrium configurations for a few loading steps in 3D.

Fig. 2 :

 2 Fig. 2: Beam subject to a concentrated tip torque: Time evolution of K2 , K3 for several steps of loading.

Fig. 3 :

 3 Fig. 3: Beam subject to a concentrated tip torque: Beam's tip position vs tip torque magnitude.

Fig. 4 :

 4 Fig. 4: Beam subject to a follower tip force: Equilibrium configurations for a few loading steps.

Fig. 5 :

 5 Fig. 5: Beam subject to a follower tip force: Beam's tip position vs tip force magnitude.

Fig. 6 :

 6 Fig. 6: Flexion-torsion test: Evolution of the beam's configuration for several values of g 1 (t): (top-left) in 3D, (topright, bottom left and right) in planes x -y, x -z, y -z respectively.

Figures 7 (

 7 a,b,c,d) display several beam configurations obtained with T = 1, 2, 3, 4, 5N for each routine, except the parallel routine for which one obtains the same deformations as in figure 1 with C(t) = (R b /2)T (t). The simulated beam is that of test 1 but one hundred times softer (E = 10 6 Pa). All the deformations are planar except those of the spiral routine whose figure 7.c illustrates one perspective view. All simulations are carried out with 5 modes for each strain field, except in figure 7.d, i.e. for the parallel-truncated routine, for which we represented the results obtained with 5 and 10 modes. In fact except for the parallel-truncated routine, the Ritz discretisation has converged before 5 modes (the deformations do not change significantly by adding Y = R b /2, X ∈ [0, 1/2] 0 7.d truncated D Y = 0, X ∈]1/2, 1]

Fig. 7 :

 7 Fig. 7: Beam internally actuated with one cable. Configurations obtained by increasing the tension T for different routines: (a) Convergent. (b) Cross. (c) Spiral. (d) Paralleltruncated with 5 (dotted lines), and 10 modes (solid lines).

Fig. 8 :

 8 Fig. 8: Clamped beam released after bending with a tip torque C: Snapshots every 0.12s over the first second of simulation, when C = 0.8 (a), 2.3 (b) and 5.4Nm (c).

Fig. 9 :

 9 Fig. 9: Clamped beam bent with a tip torque C = 0.8,2.3,5.4Nm and released: Time evolution of its tip abscissa (a), and tip ordinate (b).

Fig. 10 :

 10 Fig. 10: Beam bent with C = 2.3Nm and released: Snapshots at t = 1s (a) and 1.5s (b) with n = 1, 2, 3, 4, 5 modes.

Fig. 11 :

 11 Fig. 11: Pendulum swinging in gravity: Snapshot every 0.15s.

Fig. 12 :

 12 Fig. 12: Pendulum swinging in gravity: Trajectories of the two ends in the y -z plane.

Fig. 13 :

 13 Fig. 13: Pendulum swinging in gravity: Time evolution of the tip position.

Fig. 14 :

 14 Fig. 14: Free-free internally actuated beam in vacuum: Snapshots every 0.3s in perspective view.

Fig. 15 :

 15 Fig. 15: Free-free internally actuated beam in vacuum: K 2 (left) and K 3 (right) every 0.3s.


Herve Parmentier

Dominique Andrieu

Sylvie Servain-Courant

Anne Riviere-Honegger

Gilles Armani

Integration of Blue Corridors in the urban requalification projects of riverside conurbations (Lyon, Nantes, France)