Johann Herault
email: ajohann.herault@imt-atlantique.fr

A geometrically exact approach for floating slender bodies with finite deformations

Keywords: Floating body, equilibrium configuration, finite deformation, geometrically exact approach, Cosserat beam

Fast prediction of the equilibrium configuration of a static floating rigid body is a challenging problem, especially if the body geometry is complex. To address this problem, a new model is developed to seek quickly the equilibrium configurations (position and orientation) of a continuous slender body with finite deformations relatively to the water surface. The non-linear description of the body geometry relies on the geometrically exact approach, where the orientations of the sections are parametrized by a continuous field of rotation matrices. After introducing the advantages of this approach, a new model giving the exact net wrench of buoyancy forces applied to the so-called Cosserat beam is reported. Then, an optimization method adapted to the Lie group manifold is developed to compute the equilibrium configurations. The numerical results are then compared to experimental data to validate the accuracy of the model. Finally, the inverse problem consisting of finding the relevant body deformation for a desired equilibrium configuration is introduced. After demonstrating that the inverse problem is ill-posed, a method is reported to deform continuously the body with a constant stable head configuration.

Introduction

Many maritime applications require to maintain a body in a stable position on a free surface [START_REF] Sørensen | A survey of dynamic positioning control systems[END_REF]. The stability of floating rigid bodies has been extensively studied thanks to the well-known stability criteria for single hulls [START_REF] Biran | Ship hydrostatics and stability[END_REF]. However, controlling or predicting the static stability of a strongly deformed body remains an open problem. Surprisingly, some animals are capable to solve an analog problem in a very elegant way. For instance, semi-aquatic snakes, like the Cottonmouth or grass snakes, deform significantly their body to maintain a static stable posture on a free surface, while their head remains above the surface [START_REF] Jayne | Swimming in constricting (elaphe g. guttata) and nonconstricting (nerodia fasciata pictiventris) colubrid snakes[END_REF]. These snakes use their body deformation to balance gravity and buoyancy forces and achieve (quasi-)static positioning. Understanding and reproducing the performance of these animals could lead to new applications in maritime engineering. In the field of bioinspired marine robotics, stability induced by body deformation could improve thrust production and energy efficiency of snake-like robots [START_REF] Ijspeert | From swimming to walking with a salamander robot driven by a spinal cord model[END_REF]. Moreover, determining the stable position of deformed cables on sea surface, which are widely used for offshore engineering, is also of great interest for maritime 25 applications [START_REF] Davidson | Mathematical modelling of mooring systems for wave energy converters-a review[END_REF]. Moreover, poly-articulated floating systems, like floating articulated arms or energy harvesting system, could adapt their shape to optimize the efficiency of their task thanks to a controlled stability.

30

Hence, one needs to understand the effect of the geometric deformation of a slender body on its equilibrium configurations. The objective of the paper is thus to introduce an efficient tool computing the equilibrium configurations of a slender deformed 35 body floating on flat sea surface. More precisely, the aim is to find the orientation and the position of a mobile frame attached to the body for a given body shape so that the net vertical force and the net torques cancel out. An approach based on ex-40 periments or DNS would be very time-consuming to perform a parametric study of the influence of the deformations. Consequently, a simplified model is derived to perform rapidly the calculations. The simplification of the problem consists of describing the slender body as a 3D beam plunged into a 3D Cartesian space with a flat surface. To avoid time consuming simulations, equilibrium configurations are sought via an optimization method finding the roots of the net wrench instead of solving the dynamics. To take into account large deformations, the model follows the approach introduced by J.C. Simo: the geometrically exact approach [START_REF] Simo | A finite strain beam formulation. the threedimensional dynamic problem. part i[END_REF][START_REF] Simo | A three-dimensional finitestrain rod model. part ii: Computational aspects[END_REF][START_REF] Simo | On the dynamics in space of rods undergoing large motions-a geometrically exact approach[END_REF]. In this framework, the parametrization of each section orientation of the beam is performed thanks to rotation matrices taking values in the Lie group SO [START_REF] Jayne | Swimming in constricting (elaphe g. guttata) and nonconstricting (nerodia fasciata pictiventris) colubrid snakes[END_REF]. The great advantage of this model is to preserve the geometric properties of the rotations and to avoid afterwards complexities or singularities of others parametrizations, like the Euler's angles. However, the price to be paid for using rotation matrices instead of linearised models (infinitesimal rotation) or other representations, is to perform differential calculus on a manifold rather than in a vector space. For instance, the standard optimization methods used to compute equilibrium in a vector space have to be adapted to the geometry of the Lie group, as it will be shown in section 4. Consequently, the present paper introduces in a pedagogical way some of the specific features of the geometrically exact approach.

In the context of ocean engineering, a geometrically exact approach has been introduced in [START_REF] Burgess | Bending stiffness in a simulation of undersea cable deployment[END_REF] to model the dynamics of underwater flexible cables. The geometrically exact approach has also been successfully used to describe the swimming motion of slender bodies in sub-marine robotic [START_REF] Boyer | Multibody system dynamics for bio-inspired locomotion: from geometric structures to computational aspects[END_REF][START_REF] Boyer | Fast dynamics of an eel-like robot-comparisons with navier-1090 stokes simulations[END_REF][START_REF] Boyer | Macrocontinuous dynamics for hyperredundant robots: application to kinematic locomotion bioinspired by elongated body animals[END_REF], and in fluid mechanics applied to biology [START_REF] Candelier | Three-dimensional extension of lighthill's large-amplitude elongated-body theory of fish locomotion[END_REF]. Indeed, this geometrical method is particularly suitable for problem requiring to transport the kinematics and forces from one frame to another [START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF]. This approach is thus useful for animal locomotion since the forces acting on each segment are usually expressed in the mobile frames, while the motion is commonly described in the Galilean laboratory frame [START_REF] Candelier | Three-dimensional extension of lighthill's large-amplitude elongated-body theory of fish locomotion[END_REF]. Therefore, the geometrically exact approach is also well-adapted to the present problem. Unlike the previously mentioned works [START_REF] Boyer | Multibody system dynamics for bio-inspired locomotion: from geometric structures to computational aspects[END_REF][START_REF] Candelier | Three-dimensional extension of lighthill's large-amplitude elongated-body theory of fish locomotion[END_REF] considering infinite volume, one of the originalities of the present paper is to introduce buoyancy forces breaking the material symmetry of the isotropic sub-marine space. Thanks to the exact geometrical parametrization of each section relative to the water surface, it is possible to compute exactly the local buoyancy forces acting on a segment in the local 95 body section frame, and then to express the resulting force in the Galilean reference frame. Hence, the reported model combines a geometrically exact description of the body with an exact model of buoyancy forces.

100

This paper is organized as follows. First, the geometrical model of the slender body is presented in section 2. The computation of the buoyancy forces is then developed in section 3. The exact calculation of the stiffness matrix is also reported. To 105 find the equilibrium configuration, an optimization method adapted to the manifold of the body configuration is reported in section 4. To benchmark the code, the equilibrium positions of a twisted elliptical cylinder and a helical rod are studied in sec-110 tion 5 thanks to this new model. The results are then compared to experimental data. Finally, the inverse problem consisting of finding the relevant body deformation for a desired equilibrium configuration is presented. It is demonstrated that the 115 inverse problem is ill-posed since it exists an infinite set of body deformations corresponding to the same equilibrium configuration. To illustrate this feature, a method is reported to deform continuously a slender body while maintaining the same 120 head configuration.

Geometry of the Cosserat's beam

In this section, the framework of the exact geometrically approach is introduced in a pedagogically meaningful way. The reader interested in this 125 approach will find in [START_REF] Simo | A finite strain beam formulation. the threedimensional dynamic problem. part i[END_REF][START_REF] Simo | A three-dimensional finitestrain rod model. part ii: Computational aspects[END_REF][START_REF] Simo | On the dynamics in space of rods undergoing large motions-a geometrically exact approach[END_REF][START_REF] Boyer | Multibody system dynamics for bio-inspired locomotion: from geometric structures to computational aspects[END_REF][START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF][START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional eel-like robot[END_REF][START_REF] Marsden | Introduction to mechanics and symmetry: a basic exposition of classical mechanical systems[END_REF] more details concerning the geometrically exact approach and the Lie group theory applied to mechanics. The continuous buoyant slender body is modelled by an inextensible shear-less beam of length L with el-130 liptical cross section that is plunged into a Cartesian space (Fig. 1). A Galilean reference frame F G = (0, e 1 , e 2 , e 3 ) is pinned to the water interface so that the surface corresponds to the altitude z = 0. The centerline of the beam is defined by the 135 set of points {O s }, which are located by the position vectors p(s) = (x(s), y(s), z(s)), and expressed in F G , as a function of the curvilinear variable s ∈ I = [0, L]. The so-called Cosserat beam is modelled by a continuous stack of elliptical rigid sections of 140 infinitesimal thickness ds. The position and the orientation of these sections are defined from the set of material frames F s = (O s , t 1 (s), t 2 (s), t 3 (s)) with origin O s . The orthonormal vectors set {t i } i=1,2,3 Thanks to the rotation matrices and the position vectors, each point in a section referenced by the coordinate vector q (b) (s) in F s can be expressed in F G with coordinates q (G) (s) given by

q (G) (s) = p(s) + R(s)q (b) (s) (1)
Consequently the fields p(s) and R(s) fully characterize the configuration of the beam relatively to the reference frame F G . Equation 1 can be written in a condensed form via the rigid transformation

g s (p, R) so that q(G) (s) = g s q(b) (2)
with qT = (q T , 1) the homogeneous coordinates [START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF], and the matrix g ∈ R 4×4 defined by

g s = R(s) p(s) 0 1 (3)
The rigid transformations g s take values in the Euclidean group SE(3) = R 3 × SO(3). Hence, the configuration of the body in the reference frame F G corresponds to a continuous curve of transformation g s with s ∈ I in SE [START_REF] Jayne | Swimming in constricting (elaphe g. guttata) and nonconstricting (nerodia fasciata pictiventris) colubrid snakes[END_REF]. The transformation g s has an inverse transformation g -1 s given by

g -1 s = R T -R T p 0 1 (4)
To describe the deformation of the beam irrespectively of its configuration in the surface frame, we introduce the rigid transformation 175 fields 0 g s 0 p(s), 0 R(s) that characterizes the body shape relatively to the mobile head frame F 0 (Fig. 1). The transformations g s is given by the composition of the transformation 0 g s with g 0 (Fig. 1) so that 180 g s = g 0 0 g s (5) For an inextensible rod, the geometry of the beam is fully characterized by the field of vectors set 0 t i (s) forming the material frames expressed in the head frame F 0 . The deformation of the vectors vature vector κ(s) = (κ 1 , κ 2 , κ 3 ). This strain field defines the infinitesimal rotations of the sections expressed in the local material frame

F s κ(s) = i κ i (s) 0 t i (s). ( 6
)
with the parameters κ 1 , κ 2 and κ 3 quantifying respectively the torsion and both bending deforma-

190
tions of the rod. The vectors fields 0 t i are thus given by the evolution equation

d ds 0 t i =κ × 0 t i (7)
for i ∈ [1, 2, 3] with the initial condition given by 0 t i (s = 0) = u i , with u i a unitary vector with (u i ) i = 1, and 0 otherwise. Using the rotation matrix to express the transformation of the vector 0 t i so that 0 t i = 0 R s u i , equation 7 can be written into a matrix form with 0 R s

d ds 0 R s = 0 R s κ(s) (8)
and the initial condition 0 R 0 = Id. The skewsymmetry matrix κ is the matrix form of the vector 200 product κ× defined by

κ(s) =   0 -κ 3 κ 2 κ 3 0 -κ 1 -κ 2 κ 1 0   (9)
The vector κ is thus the twist coordinate of the matrix κ taking values in the linear space so(3). The right-hand side term of equation 8 is the composition of first an infinitesimal rotation κ(s) in the material frame, with a finite rotation 0 R(s) corresponding to the mapping from the material frame to the head frame [START_REF] Simo | On the dynamics in space of rods undergoing large motions-a geometrically exact approach[END_REF]. The solutions of equation 8 can be written under the following form 0 R s = exp( Ω(s)) [START_REF] Boyer | Multibody system dynamics for bio-inspired locomotion: from geometric structures to computational aspects[END_REF] with Ω(s) ∈ so(3) a skew-symmetric matrix, and exp(•) the matrix exponential defined by

exp( Ω) = I d + ∞ n=1 1 n! Ω n (11)
Equation 10 refers to the exponential map in the 210 Lie groups theory [START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF][START_REF] Marsden | Introduction to mechanics and symmetry: a basic exposition of classical mechanical systems[END_REF], mapping element from so(3) to SO(3) (Fig. 2). The exponential map is a convenient tool to parametrize finite rotation as a function of skew-symmetry matrices. If the strain κ is constant, the matrix Ω is equal to Ω(s) = κs.

215

In the present paper, the exponential matrix will be used for the computation of the head orientation R 0 in section 4. However when the strain field κ(s) is not constant, the skew-symmetric matrix Ω can only be approached by a Magnus expansion 220 [START_REF] Blanes | The magnus ex-1110 pansion and some of its applications[END_REF]. Instead of an approximation of this expansion, a parametrization of the rotation based on the quaternions will be used to solve equation 8 (for details see Appendix A.1.2). The main advantage of using the quaternion instead of the rotation matrix 225 is that it is easier to preserve the unit norm of Q than the orthonormal property of R. Then, the position vector 0 p s for an inextensible rod is given by

0 p s = s 0 0 t 1 (s)ds (12)
Finally, the transformation 0 g s is reconstructed 230 via the fields 0 p(s) and 0 R s thanks to equation 3. In the model, the field of transformation 0 g s defined by the strain field κ(s) is fixed, and the transformation g 0 will be the unknown of the problem. To take into account the possible anisotropy of a body section, the shape of each section is considered elliptical with semi-axis aligned with the vectors (t 2 , t 3 )(s). For a given section, all the points belonging to a body section indexed by the position vectors q (b) with components (q 1 , q 2 , q 3 ) in F s , satisfy the following relations

q 1 = 0, q 2 a 2 + q 3 b 2 < 1, (13)
with a and b the semi-radii. The body density ρ b 235 is considered constant, and smaller or equal to the water density ρ w . In the rest of the paper, all the lengths will be rescaled by the length L. Hence, the curvilinear variable s varies from 0 to 1 with I = [0, 1]. The parameters a and b become the 240 semi-radii normalized by the length L.

The numerical integration are performed thanks to a Chebyshev collocation method implemented in a MATLAB code, using differentiation matrices [START_REF] Yagcia | A spectral-tchebychev technique for solving linear and nonlinear beam equations[END_REF]. This spectral method provides a good con-245 servation of the quaternion norm. Moreover, this method exhibits a fast convergence of the solutions [START_REF] Trefethen | Spectral methods in MATLAB[END_REF], the so-called "spectral accuracy" with an exponential decrease of the numerical error as a function of the number of grid points. This modal decompo-250 sition will be also useful to solve the inverse problem reported in section 6.

There are two worth noting remarks. First, the present geometrical model of the body is continuous, unlike a lumped mass model commonly used for mooring systems [START_REF] Davidson | Mathematical modelling of mooring systems for wave energy converters-a review[END_REF]. We will show in the next section that the model of buoyancy forces relies also on an exact calculation of the immersed surface and of the buoyancy center for each section, thanks to the geometrically exact approach. Secondly, the model relies on a rod parametrization satisfying the Kirchhoff constraint, i.e. an inextensible beam without shear between the sections. Therefore, the sections remain orthonormal to the centerline but might be non-parallel between each other, if the fields κ 2 (s) and κ 3 (s) produce a local bending of the centerline. The helix reported in section 5 displays this feature. The use of the Kirchhoff model is justified for slender bodies submitted to moderate load, otherwise, the Ressner's model has to be considered [START_REF] Audoly | Elasticity and geometry[END_REF] without major modification of the present model.

Calculation of the buoyant forces and stiffness matrix

Buoyant forces exerted on a Cosserat beam

The objective of this section is to perform an exact calculation of the density of force and torque acting on each infinitesimal elliptical section. Assuming a mechanical equilibrium between the rod and the water surface at rest, the acceleration of the fluid and body are neglected. The density of wrench w a (s) of the buoyant force exerted on a section and expressed in the frame F s is defined by

w a (s) = π a (s) γ a (s) (14)
with π a (s) and γ a (s) defining respectively the density of force and torque. The density of buoyancy force expressed in the local frame F s is

π a = a g ρ w S im R(s) -1 e 3 ( 15
)
with a g the gravity acceleration, ρ w S im the linear density of water mass for a local immersed surface S im (s), and R -1 e 3 the vertical vector e 3 expressed in F s , via the inverse rotation R -1 . To calculate the torque associated with the density of buoyancy force, one needs to define the position of the buoyancy center q B (s) where the force is applied to estimate the moment arm. Thus, the density of torque is given by

γ a = a g ρ w S im q B × R -1 e 3 (16)
To calculate the net wrench, all the wrenches have to be expressed in the same reference frame. For convenience, the Galilean reference frame F G is chosen. The condensed relation w

(G) a (s) = Ad T
g -1 w a (s) (Eq.A.9, see Appendix A.1.4) will be 300 used in section 4 to express all the local wrenches w a (s) in their section frame into the Galilean frame F G . The net wrench is thus given by

W (G) a W (G) a = s∈I w (G) a (s)ds (17)
Finally, the net wrench W N (the index (G) is omitted) combining the buoyancy force and the 305 gravity force reads

           a g 1 0 (ρ w S im -ρ b S b ) ds e 3 a g 1 0 (ρ w S im ) Rq B + (ρ w S im -ρ b S b )pds × e 3 (18)
Due to the symmetry, only three components of W N can be non-zero. Indeed, the forces on the body are invariant by parallel body translation relatively to the sea surface, i.e. in the (0, e 1 , e 2 ) plane, 310 and also by solid body rotation along the vertical axis e 3 . The unknowns of Eq. 18 are the surface S im and the weighted barycenter S im q B that are calculated in the next section.

Immersed surface and weighted barycenter 315

For the sake of simplicity, the calculation of S im and S im q B is performed in the section frame F s . The waterline is parametrized by the equation F q (b) = 0, with q (b) the position vector expressed in the frame F s of the considered section.

320

The values of S im and q B in a body section are thus defined by

S im = Ds dq 2 dq 3 (19)
and

S im q B = Ds qdq 2 dq 3 (20)
with D s the domain of the immersed surface (detailed in Appendix A.3). To calculate the functions 325 S im and q B , we will follow the following steps consisting of : 1) parametrizing the waterline in the section frame F s , 2) mapping the elements from the ellipse to a circle to compute the integral, 3) defining the immersion conditions, and 4) integrating the equations. Since the water surface is defined by z = 0 in the frame F G , the position q b of the water-line in F s satisfies the following equality (Eq. 1)

ē3 • g s q(b) = 0 ( 21
)
where ē3 corresponds to the homogeneous coordinates of the vertical vector e 3 (Eq. 3). The border of the elliptical section being characterized by q 1 = 0, the implicit function F (q 2 , q 3 ) = 0 of the waterline is given by

F (q 2 , q 3 ) = R 32 q 2 + R 33 q 3 + p 3 ( 22
)
with R ij the component of the rotation matrix R, and p 3 the third component of p. This equation becomes ill-posed when R 32 = R 33 = 0. This problem appears if the vertical components of the vectors t 2 and t 3 vanish, when the section is parallel to the surface with t 1 ∝ e 3 . In this case, the concept of waterline is meaningless since since the section is either fully wet or fully dry. To regularise the problem, we fix S im (s) = S b if p 3 < 0, and 0 otherwise. Therefore, this case is easily solved and implemented in our model while all the other configurations are treated in the further calculations.

To simplify the integration, we introduce a bijective linear transformation (Fig. 3) that maps the elements from an ellipse with semi-radii (a, b) to a disk of radius r 0 = √ ab defined by

L : q 2 q 3 → q2 q3 = b a 0 0 a b q 2 q 3 ( 23
)
The waterline corresponds now to a chord of the disk (Fig. 3). Due to the symmetry of the problem, the new barycenter qB belongs to the line directed by the radial vector e r perpendicular to the waterline, and can be written

qB = r B (h)e r ( 24
)
with a radius r B > 0 calculated thanks to the level of immersion h. The variable h is defined from the position vector qh = he r of the point of intersection between the waterline and the radius orthogonal to To obtain the level of immersion h (see Appendix A.3), the vector qh is projected on e r . Finally, the level of immersion h is given by

370 h = p 3 R 2 32 a b + R 2 33 b a (26)
This height h is algebraical in order to estimate the transition from a dry body to a fully immersed body. Four cases can be identified. First, if h is larger than the radius r 0 = √ ab, then the section is dry. When h = r 0 , the section has a punctual 375 contact with the water surface. When |h| < r 0 , the section is partially immersed. Finally, if h < -r 0 , then the section is fully immersed. To summarize, the surface is immersed (third and fourth case), if

p 3 < R 2 32 a 2 + R 2 33 b 2 . ( 27
)
Now the surface can be calculated thanks to the height h via a surface integration with a cartesian

S im = r 2 0 arccos h r 0 - h r 0 1 - h 2 r 2 0 . (28)
The variation of the immersed surface as a function of h/r 0 is reported in Fig. 4 (blue curve). The advantage of defining an algebraical height h is obvious since the surface increases smoothly as a function of h/r 0 ∈ [-1, 1]. In the same way, the distance r B introduced in Eq. 24 is given by

S im r B = 2r 3 0 3 1 - h r 0 2 3 2 . ( 29
)
The evolution of the rescaled distance r B /r 0 of the barycenter is reported in Fig. 4, varying from the center of the disk for h/r 0 < -1 to the radius r B = r 0 for a punctual contact at h = r 0 . Finally, the weighted position of the barycenter of the immersed portion of the ellipse given by q B = r B L -1 e r reads

S im q B = 2r 3 0 3 1 - h r 0 2 3 2 a b cos(θ) b a sin(θ) (30)
with θ given by Eq.25. Consequently, equations 29 and 30 giving the immersed surface and the weighted barycenter, can be used to compute the 395 net wrench exerted on the body via Eq. 18.

Stiffness matrix

In this section, the stiffness matrix quantifying the stability of a head configuration g 0 is performed in the framework of the geometrically exact ap-400 proach. The differential of W N written δW N is linearly related to the twist coordinate of a virtual head displacement δΨ (G) 0 ∈ R 3 in the spatial Galilean frame thanks to the stiffness matrix (D Ψ0 W N ) so that

405 δW N = D Ψ0 (G) W N δΨ 0 (G) (31)
with δΨ 0 (G) T = (δz h , δθ r , δθ p ) the three components corresponding respectively to a variation produced by heaving, rolling and pitching motions expressed in the Galilean frame. First, the variation δW N is expressed as a function of δg s with

410 δW N = s∈I (D gs w N ) δg s ds (32)
where (D gs w N ) is the jacobian matrix of the wrench density for a variation δg s , which is reported in Appendix A.4. Then, the variation δg s is expressed as a function of the virtual displacement δψ s ∈ R 6 in the body frame with

415 δg s = g s δ ψ s (33)
Finally, the adjoint operator Ad (g -1 s ) is used to transport the virtual displacement

δψ s = Ad (g -1 s ) δψ (G) 0 (34)
Combining Eqs. 31, 32, 33 and 34, the jacobian matrix reads

D Ψ0 (G) W N = L 0
(D gs w) g s T Ad (gs) P ds (35) with T a linear operator transforming the twist co-420 ordinates into the screw-symmetric tensor, and P an operator mapping δΨ

(G) 0 → δψ (G) 0 (R 3 → R 6 ), so that δψ (G) 0 = PδΨ (G)
0 . This jacobian matrix will be used in the Newton-Raphson method reported in the next section.

Optimization method

The purpose of this section is to derive an optimization algorithm seeking the equilibrium configuration for a given body shape defined by the field of configuration 0 g s . The choice of the mobile frame to perform the optimization method is of great importance, since the net wrench W N , and more precisely the net torque, depends on the position of the considered mobile frame where the resultant is estimated. If the head frame F 0 is useful to describe the geometry, it is less convenient for the optimization algorithm, and it shows poor convergence. The center of mass (COM) defined by equation A.6 will be preferred to the head frame. The configuration of the COM frame g com in F G is given by

g com = g 0 0 g com ( 36
)
with 0 g com a constant matrix defined in Eq. A.6. Thanks to Eq. A.7, the wrench W com N computed in the center of mass frame is written in a condensed form with the adjoint operator (cf. text after equation A.7 )

W com N = Ad T (g com ) -1 W N (37)
Hence, the algorithm will consist of finding the configuration g com canceling out W com N , but the final result will be expressed as a function of g 0 . To do so, the optimization method has to preserve the group property of the configuration g ∈ SE(3) (for the sake of clarity, we will omit the sub-index com). Indeed, the challenge consists of defining an incremental transformation g n → g n+1 at the n + 1 steps so that the wrench W N [g n+1 ] goes toward zero, while the transformation g n+1 remains in the group. Since the problem derived from a potential energy, the existence of equilibrium position is related to the extrema of the potential energy U defined by

U = 1 s=0 [ρ b S b z com (s) -ρ w S im (s)z B (s)] ds (38)
with z com the altitude of the center of mass of the section indexed by s, and z B (s) the altitude of the buoyancy center. A classical minimization method in R k will consist of finding the minima of U (x) as a function of an element x n ∈ R k via a linesearch method. This method is based on the update formula x n+1 = x n + t n ξ n , where ξ n ∈ R k is the search direction and t n ∈ R + the step size. In R k , the concept of moving in the direction of a vector is straightforward, since the tangent space is also R k . To keep the group structure, each new configura-470 tion g n+1 is deduced from g n thanks to the update formula with the transformation composition

g n+1 = g n ( n g n+1 ) (39)
where n g n+1 is the relative transformation mapping from the configuration g n to g n+1 (Fig. 5). Hence, at each incremental step from n to n + 1, the mate-

475
rial frame is transported and re-oriented thanks to a transformation n g n+1 ∈SE(3). The transformation n g n+1 is then parametrized by the step direction ξ n ∈ R 6 thanks to the exponential mapping given by Eq. 10

480 n g n+1 = exp(t n ξ n ) (40)
with ξ n the skew-symmetric matrix defined from the twist coordinate ξ n thanks to Eq. 9. Equation 40 defines a retraction, i.e. the projection of the step direction with a step size t n onto the manifold. The algorithm will be based either on a mod-485 ified gradient descent or a Newton-Raphson methods with an unitary step direction defined by

ξ n = 1 |DW n | DW n (41)
with W n the net wrench. For the gradient method, the matrix D is a diagonal matrix with positive elements, that are calculated to avoid unbalanced in 490 the weight of the components of the wrench. For the Newton-Raphson method, the matrix D is the inverse of the stiffness matrix D Ψ0 (G) W N -1 computed in section 3.3.

A satisfying condition for the algorithm will be to 495 impose a reduction (for a stable equilibrium) of the potential energy so that U n+1 < U n . However, this condition does not provide a sufficient rate of convergence to the solution. For a sufficient decrease, the Armijo condition states that the following inequality has to be satisfied

U n+1 < U n -c 1 t n ξ T n (W N ) n (42)
for c 1 ∈ [0, 1], (W N ) n and U n being the wrench and the energy estimated at the step n. This method provides a convergence to the solution (see Appendix A.5). The Armijo condition will determine the step size so that t n corresponds to the largest step size satisfying the inequality 42 . At a step n, we assume that U n and (g 0 ) n are known. We define cf < 1 a contraction factor reducing the step size at each iteration, and t the ini-510 tal time step. The algorithm determining t n combines equations 36, 37, 39, 40, and 42. It is reported in the algorithm 1. The computation of the wrench and the energy given by Eq. 18 and Eq. 38, are respectively referenced by the functions WrenchComputation and EnergyComputation.

Algorithm 1 Algorithm describing the update from the configuration (g 0 ) n to (g 0 ) n+1

(W N ) n ← WrenchComputation((g 0 ) n ) (W com N ) n ← Ad T (g com ) -1 n (W N ) n {Transport of the wrench to the COM} ξ n ← D (W com N ) n {Calculation of the step direc- tion} U n+1 ← U n {initiation of the energy} t n ← t {initiation of the step} while U n+1 > U n -c 1 t n ξ T n (W com N ) n do t n ← cf × t n {step size contraction} n (g com ) n+1 ← exp t n ξ n {Retraction} (g com ) n+1 ← (g com ) n n (g com ) n+1 {Update for- mula} (g 0 ) n+1 ← 0 g com -1 (g com ) n+1 {Head conFig. } U n+1 ← EnergyComputation((g 0 ) n+1 ) end while
After each loop, the configuration is updated with g n = g n+1 . To accelerate the code, the initial step size t at the iteration n is equal to the final step size t n-1 of the previous iteration for 520 n > 1. The algorithm 1 is then repeated into another loop with n iterations until the step size t n reaches a given threshold, usually chosen smaller than 10 -10 . Then, the algorithm is stopped because the energy minimum is almost reached. We 525 observe that the resulting order of magnitude of the wrench is around 10 -8 .

Note that the Newton-Raphson method could not converge if the initial state is not close enough to the solution. Therefore, the gradient method 530 could be preferred in some cases for its robustness.

Experimental benchmark

Body geometry

To illustrate and evaluate the performance of the code, we have chosen two kinds of body deforma-535 tion: (a) a helical rod with circular cross-sections, and (b) a piecewise-twisted elliptical cylinder.

The first body is a helix with a height h = 50mm, a pitch p = 35mm, a radius a = 15mm, and an angle α = 40 • parametrised by the free software 540 FreeCAD [START_REF] Freecad | [END_REF] (reported in Fig. 6). The volume estimated by the software is V = 2.48 × 10 -4 m 3 , and its length is L = 35.16 mm. The body has been 3D-printed with a mass of 64g so that the average body density is given by ρ s = 257.5 kg/m 3 .

545

After post-processing the mesh, the strain field κ is reported in Fig. 6.

The second body is made of 3D-printed elliptical sections that are stacked one next to the other. Each body section is elliptical with semi-550 radii a = 2cm, b = 3cm, and width e = 8mm. The angle of twist between each section is θ 1 = -π/10. The body reported in Fig. 7 is made of N seg = 17 sections that are stacked together by an axial rod. To model the piecewise rotation, we introduce the 555 following field of strain

κ 1 (s) = θ 1 Nseg-1 n=1 g (s -s n ) ( 43
)
where the sum is over the N seg -1 junctions, and κ 2 , κ 3 are set to zero. The function g is a Gaussian function centred in s n = n Nseg-1 of norm one ( gds = 1). The variance of the gaussian function 560 is small enough to model the localized twist so that we call this model the "quasi-discrete" model. The picture of the body and its numerical avatar are reported in Fig. 7.

Benchmark

Then, we compare the numerical results of the model with the experimental data for both geometrical configurations. We start by studying the stability of the helix with circular cross-sections. The picture of the floating body is reported in Fig. 8 so that its cross-section faces the lens. After a postprocessing, the contour of the upper part of the body is estimated, and reported in red in Fig. 8.

We now perform the calculation of the equilibrium configuration thanks to the geometry of the body introduced in section 5.1. It is worth noting that we use an average mass density ρ s estimated from the net mass and the net volume of the body, whereas the 3-D printing process could produce small inhomogeneous distribution of mass per cross section. The body is discretized on N = 100 points on the Chebyshev grid, which is sufficient to reproduce the strain field κ reported in Fig. 6. The efficiency of the Newton-Raphson (NR) and the gradient descent (GD) methods are compared.

The evolution of the norm of the vertical force and the torques are reported in Fig. 9 (red symbols for NR method and blue symbols for GD method). Both methods display a quasi exponential convergence to the equilibrium, but with different rates.

590

The GD method requires 50 iterations while the NR method needs only 30 iterations to reach a residue equal to 10 -8 . Blue symbols are the data from the gradient method while the red ones comes from the Newton-Raphson method.

To estimate the validity of the code, the body contour of the experimental data (red curve) is superimposed to the front picture of the numerical one (Fig. 10,top). Since the solution is invariant by translation parallel to the surface and rotation along the z-axis, the abscissae x and the angle of rotation θ z are computed to obtain the best match. 600 Hence, we observe a relative good qualitative agreement with only 2 free parameters. When the body is tilted with a slight inclination angle along the yaxis (θ y = -3 • ), the matching is even better (Fig. 10,bottom). Consequently, we conclude that the model gives a good estimation of the configuration with an incertitude of θ y = ±3 • , which could be easily explained by a small inhomogeneous mass distribution of the 3-D printing or a slight inclination of the camera.

For the twisted rod with elliptical segments, a quantitative study can be performed via the estimation of the inclination angle θ x along the x-axis of the first elliptical section. The results are reported in Fig. 12 with the experimental data (black circle) and the numerical estimations (blue square and dashed curve). The number of segments N seg has been varied from 12 to 17. The experimental angles are obtained by computing the inclination angle of the first segment (Fig. 12).

620

The numerical estimation of the equilibrium configuration is performed thanks to the quasi-discrete model given by Eq. 43 (blue squares in Fig. 11), and a second continuous model, which is parametrized by a constant torsion κ 1 = θ 1 (N seg -625 1) (dashed blue curve). The first model requires at least N = 500 points for a good convergence (to capture the sharp variation of twist angle between the sections), whereas only N = 100 points are necessary for the continuous model. We ob-630 serve that the quasi-discrete model is slower than the continuous model. For instance, the solution of the quasi-discrete model takes 40 ± 1.5 seconds, while the continuous model requires 10 ± 1s for a laptop with a CPU Intel-Core i7-6500U and 16Go continuous model can be reduced to N = 10 with a small change of accuracy (around 10 -4 ) in order to decrease the computation time to 1s.

The inclination angle θ x computed via the quasidiscrete model and the continuous one model are almost equal with a relative difference smaller than 10 -3 . Consequently, the beam can be modelled by a continuous field of twist. The evolution of the torques and the force as a function of the number of iterations is reported in Fig. 13. The convergence to the minimum is significantly slower than the one of the helix case, with some plateaus (see Fig. 13). We observe a good agreement between the experimental (black circles) and the numerical (blue) 650 angles θ x . The angle difference between the model and the experiments increases from θ x = 2 • for N seg = 12, to θ x = 8 • for N seg = 17. To evaluate a confidence interval, we have performed the calculation of θ x with a variation of plus or minus 655 5% on θ 1 , i.e. a variation of 1 • . The result is reported in grey in Fig. 11. The experimental data falls in the confidence interval. Thus, the angle difference between the experiment and the model can be explained by a possible slight variation of the 660 twist angle θ 1 between each segment.

Inverse problem and invariant manifold of

the head configuration

An ill-posed problem

This section deals with the inverse problem con-665 sisting of seeking the body shape that provides a stable head configuration. This problem is motivated by applications like the static-positioning of deformed structures or the design of complex floaters to stabilize naval structures.

670

The direct problem treated in section 4 and 5 is well-posed since a countable number of solutions exists for a given strain field κ(s). The convergence to the solutions depends on the initial conditions. However, the inverse problem is of ill-posed kind, because the solutions are not unique. Indeed, the strain κ(s) takes value in a function space that is an infinite-dimensional vector space. In our case, this space is assumed to be a Hilbert space in order to apply the Chebyshev polynomial decomposition (cf.

The system is underdetermined for a truncation at order N p > 1 since there are 3N p degrees of freedom with only 3 conditions of stability for W N = 0. For an exact decomposition with N p → ∞, the dimension of the solution space is infinite, i.e. there is an infinite way to deform the body without modifying the head configuration. In the next section, we show how to compute the strain field that form a continuous manifold named M(g 0 ), the so-called invariant manifold defined by all all the strain κ(s) 695 corresponding to the same head configuration g 0 .

The invariant manifold of the head configuration

First, we introduce the relative variation of wrench δW N produced by an infinitesimal vector 700 δκ(s) defined by W N (g 0 , κ + δκ) = W N (g 0 , κ)+δW N (δκ)+o(|δκ(s)|) (44) with the following linear differential form

δW N (δκ) = L 0 (D κ w) δκds (45)
(D κ w) being the jacobian of the wrench density relative to the variation of twist δκ(s). The exact calculation of (D κ w) is reported in Appendix A.6.

It is worth noting that the calculation of (D κ W N ) requires to evaluate the variation of the local configuration δg s of a section as a function of δκ. In the geometrically exact framework, this transformation is given by a linear differential equation (see

710

A.29) that is relatively simple compared to other parametrizations (Euler's angles). This is one of the great advantage of the geometrically exact framework.

To pass from one equilibrium to another, the dif-715 ferential form has to cancel out with δW N = 0.

To produce no wrench variation, we seek the incremental twist δκ(s) that belongs to the null space of the linear operator defined by Eq. 45. Thanks to the Chebyshev polynomial decomposition, the field 720 δκ(•) can be parametrized by a truncated series

δκ(s) = Np-1 j=0 (δ κ) j T j (s) ( 46
)
This decomposition is then introduced in Eq. 45 to perform the integration. Thus, Eq. 45 can be rewritten as a function of the polynomial coefficients (δ κ) j ∈ R 3 stacked in a vector δ κ with

725 δW N = (D κW N ) δ κ ( 47
)
with the 3 × (3N T ) jacobian matrix (D κW N ). Thanks to a singular value decomposition, the basis vectors of the null space of (D κW N ) can be computed so that any strain can be written into the following form 730 δ κ = P δ(κ) + P ⊥ (δ κ) ⊥ (48)

with P a 3N T ×3(N T -1) projection matrix to the null space, P ⊥ a 3N T ×3 projection matrix with the three vectors producing a variation of the wrench locally. The sub-index refers to the local tangent space of M(g 0 ), and the sub-index ⊥ corresponds 735 to the vector space perpendicular to the tangent space.

The first right-hand side term in Eq. 48 produces a tangential variation to the manifold M(g 0 ). However, Eq. 44 is a linear approximation of the re-740 sulting wrench produced by a variation of the body shape, and the equilibrium is only satisfied at order |δκ|. Any strain δ κ produces a residual wrench (even small) that must be corrected to avoid a departure from the manifold M(g 0 ). Hence, we use 745 the component (δ κ) ⊥ to compensate the residual wrench. The second right-hand side term in Eq. 48 is thus used as a corrector to remain in the vicinity of the manifold.

Illustration of the invariant manifold 750

There are some practical reasons to use the invariant manifold of the head configuration for maritime applications, i.e. maintaining the configura- Initial (thick curves) and final (dashed line) strain fields (κ 1 , κ 2 , κ 3 )(s) (respectively the blue, black and red curves) resulting from the iterative body strain. Resulting variation of strain δκ between the initial and the final states is displayed in the insert.

tion of one extremity of the body while modifying its body shape. For instance, the local stability of a deformable vessel can be modified by a quasi-static modification of the body shape, either to improve the stability or to decrease it for better maneuvrability. To clean up macro plastic, a floating manipulator could grip an object while maintaining constant contact with it. The purpose of this section is thus to give an example of the process introduced in the previous section consisting of following the manifold M(g 0 ) by successive body deformations.

To illustrate this process, we consider the helical tangent space. The resulting vector δ κ is rescaled to one, and then multiplied by = 2 × 10 -3 to produce an infinitesimal variation. Finally, the resulting vector δ κ = P δ κ is projected from the modal space to the physical space, and then subtracted from the strain fields κ. The amplitude of is small enough to remain in the vicinity of M(g 0 ), but not too small for keeping a reasonable number of iterations. For each update of the body shape, the vector δ κ⊥ is then computed to produce a wrench compensating the one due to δ κ . During the process, the ratio |δκ ⊥ |/|δκ | is of or-der 10 -4 , and the body shape variation is mostly produced by the component δκ . This process has required 6000 iterations to reduce the norm of κ by 785 3. Every 300 iterations, the head configuration is recomputed thanks to the optimization method in order to check the equilibrium head configuration. The initial (thick curves) and final (dashed line) twist fields (κ 1 , κ 2 , κ 3 )(s) are reported in Fig. 14.

790

The configurations of the central-line in the Cartesian space are also reported in Fig. 15. The black curve corresponds to the initial state, and the progressive body shape evolution is illustrated by a variation from a blue curve (first iteration) to a 795 red curve (last iteration). The head is marked by the red circle. A rotation along the z axis and a shift to the origin (x, y) = (0, 0) are performed to align the bodies since the wrench is invariant under these transformations. At every step, the decrease 800 of the bending κ 3 and the torsion κ 1 is compensated by an increase of κ 2 . At the final stage, all three components are of same order of magnitude. We clearly observe on left Fig. 15 that the rods progressively uncoils. The oscillatory bending (see 805 Fig. 14 and15) produces ripples along the rods to maintain some sections fully immersed and other fully dry. The final state (red curve) of the central line is similar to a transversally curved bow with the tail and the center immersed.

810

For each iteration, the variation of the altitude, and the roll and pitch angles of the heads frame relative to the initial state are computed, and reported in Fig. 16. The relative deviation from the initial head configuration remains very small (around 815 10 -4 ), which confirms the efficiency of the correcting process.

Conclusion

This paper presents a continuous model of a buoyant slender body with finite deformations on 820 a flat free surface. The model relies on the geometrically exact approach, where the configuration of the slender body is parametrized by elements of the group SE(3). Elliptical sections have been considered, and an exact model of buoyancy wrench is re- good precision. Hence, this is the first model considering an exact description geometry of the body with an exact computation of the buoyant force act-835 ing on elliptical cross-section. Consequently, this new model is well suited for parametric investigations studying the effect of complex geometry on the stable configuration of continuous slender bodies floating on a flat surface.

Moreover, the parametrization of the body geometry thanks based on the infinitesimal rotation field κ allows for the consideration of elastic properties of the slender body, such as the bending or torsional We defined also the N × N forward transformation matrix Γ f = Γ -1 b . The matrix Γ f and Γ b are calculated explicitly thanks to the recurrence relation A.11 and the orthonormal properties of the polynomials [START_REF] Yagcia | A spectral-tchebychev technique for solving linear and nonlinear beam equations[END_REF][START_REF] Trefethen | Spectral methods in MATLAB[END_REF].

The differentiation of the vector y is obtained thanks to the differentiation matrix defined by D 1 with y = D 1 y. Once again, we use a recurrence relation to calculated each polynomial derivative as a function of the polynomials at lower order. These coefficient are stacked into a matrix D. Thus, any derivative at order n written D n is given by D n = Γ b D n Γ f . Reciprocally, the integral K(s) = s 0 κ(x)dx of the function y can be calculated thanks to matrix I containing the coefficients related to the integration of the polynomials [START_REF] Yagcia | A spectral-tchebychev technique for solving linear and nonlinear beam equations[END_REF]. Hence, integration is performed by a matrix product so that K = (Γ b IΓ f )κ. The resolution of the evolution equation of the quaternion Q(s) (Eq. A.3) is solved thanks to the recombination method. The quaternion are stacked in a vector Q T = (q 0 ) [1:N ] , . . . , (q 3 ) [1:N ] ∈ R 4N (with (q i ) j = q i (x j )) satisfying the linear system Appendix A. [START_REF] Jayne | Swimming in constricting (elaphe g. guttata) and nonconstricting (nerodia fasciata pictiventris) colubrid snakes[END_REF]

.3. Level of immersion h

The components of the vector qh satisfy the fol- Hence, the algebraical height h = -n•q h is finally variation of net buoyancy force density will be explicitly reported. Hence, the variation of net buoyancy force density as a function of (δp z , δR 32 , δR 33 ) reads The computation of the variation of net buoyancy force density requires more steps that can't be detailed here. For the sake of clarity, the explicit formula are not reported, and only the method is explained. First, the variation is applied to γ N so that δγ N = a g [ρ w Rδ (S im r B e r ) + (δR) (S im r B e r ) +ρ w δS im p + (ρ w S im -ρ b S b )δp] × e 3 (A.25) Then, we express the variation δ(S im r c ), δ(S im ), 1005 δe r as a function of δh and δθ, which are then expressed as function of δg s . The combination of the different terms leads to the jacobian (D δgs γ N ). It is worth noting that the variation δπ N depends only on (δp z , δR 32 , δR 33 ) defining the level of immersion 1010 while the variation δγ N depends on all the components of δg s since a variation is applied to R s and p s .

δπ N = -2a g ρ w r 0   1 - h r 0 2   ×              1 R 2
The term ad ξ is the adjoint action which may be identified with the derivation on the Lie-algebra.

The adjoint action ad A B can be expressed thanks to the screw-symmetric matrix with ad A B = [ Â, B] = Â B -B Â. To invert Eq. A.28, we substitute the spatial derivative by the linear differential operator D 1 (see Appendix Appendix A.2). Assuming that the linear operator is invertible, the variation of virtual body displacement δ( 0 ψ s ) reads

δ( 0 ψ s ) = [ad ξ + D 1 ] -1 δξ (A.29)
This equation is nothing but the resolution of the EDO A.28 with δ( 0 ψ s ) unknown. Finally, the variation of wrench as a function of a variation of body strain δξ reads (D κ w) δκ = (D gs w) g s T [ad ξ + D 1 ] -1 δξ (A.30) with T a linear operator transforming the twist coordinates into the screw-symmetric tensor.

Figure 1 :

 1 Figure 1: Parametrisation of the body considered as Cosserat beam

Figure 2 :

 2 Figure 2: Illustration of the geometry of the group SO(3). Each rotation R ∈ SO(3) can be parametrized by a skewsymmetric matrix Ω ∈so(3) thanks to the exponential map.

Figure 3 :

 3 Figure 3: Transformation from ellipse to circle

Figure 4 :

 4 Figure 4: Immersed surface (blue solid line) rescaled by the ellipse area πab as a function of the rescaled height h/ √ ab, and the rescaled distance r B /r 0 (dashed red curve) between the center of the ellipse and the barycenter of the immersed surface

Figure 5 :

 5 Figure 5: Schema illustrating the update in the incremental transformation to preserve the structure of the group SE(3).

3 Figure 6 :

 36 Figure 6:Illustration of the helix (BOT), and the parametrization of the infinitesimal rotation field with the torsion κ 1 (blue curve) and the bending κ 3 (red curve) (κ 2 = 0).

Figure 7 :Figure 8 :

 78 Figure 7: (top) Piecewise-twisted rod made of 17 3D-printed elliptical sections, and (bottom) its 3D reconstruction (we have only represented the lateral surface of the body) given by Eq. 43.

Figure 9 :

 9 Figure 9: Top: norm of the vertical force as a function of the iteration of the optimization process. Bottom: norm of both bending torque as a function of the iteration of the optimization process for the helical body. The circles (respectively squares) correspond to the torque τx (resp. τy).Blue symbols are the data from the gradient method while the red ones comes from the Newton-Raphson method.

Figure 10 :

 10 Figure 10: Superimposition in the x0z plane of the contour (red curve) of the body given by the experiment on the profile of the helix obtained by the model. The abscissae and the ordinate are dimensionless. Top: without correction. Bottom: with a slight tilt of θy = -3 • along the y-axis.

Figure 11 :Figure 12 :

 1112 Figure 11: Inclination angles θx in degrees along the x-axis by respect to the original vertical orientation for a rod made of N segments. The numerical data are reported in blue (square for the discrete model, dashed curve for the continuous model), and the experimental one in black (circle).

Figure 13 :

 13 Figure 13: Top: norm of the vertical force as a function of the iteration of the optimization process. Bottom: norm of both bending torque as a function of the iteration of the optimization process. The circles (respectively squares) correspond to the torque τx (resp. τy). Blue symbols are the data from the gradient method while the red ones comes from the Newton-Raphson method.

  Figure 14:Initial (thick curves) and final (dashed line) strain fields (κ 1 , κ 2 , κ 3 )(s) (respectively the blue, black and red curves) resulting from the iterative body strain. Resulting variation of strain δκ between the initial and the final states is displayed in the insert.

825ported.-Figure 15 :Figure 16 :

 1516 Figure 15: Iterative process uncoiling the helical body while maintaining the head configuration constant. The initial body shape corresponds to the black curve, and then the new configurations are coloured from blue (first iteration) to red (final iteration). Left: top view. Right: inclined view.

T 0

 0 stiffness. A step further would consist of study-845 ing the dynamics of this body based on the Euler-Poincaré equations, which will be investigated in a next paper.with θ i = iπ/(N -1). Hence, there is a linear relation between the amplitude of the mode κ k staked in the vector κ, and the value of the function κ i = κ(s i ) staked in the vector κ so that κ = Γ b κ with the N × N backward transformation matrix Γ b given by Γ (x N -1 ) . . . T N -1 (x N -1 )

  corresponds to the parametrization of the waterline, the second one to the radius. The solution of the system reads qh

of RAM. Note that the number of points of the

Acknowledgements

This work was supported by the Region Pays 850 de La Loire (NATRIX Project) and by the French National Research Agency ANR (grant no. ANR-2019-CE33-0004-01). The author is grateful to IMT Atlantique for hosting his research, and kindly acknowledges F. Boyer and V. Lebastard for their

set is parametrized thanks to the field of twist curappendix Appendix A.2). This method allows us to map elements from the vector space of the amplitudes of the polynomials (κ 1 , κ2 , κ3 ) [0:Np-1] to the functions of the Hilbert space. A bijective mapping is achieved when N p → ∞. body presented in section 5. The aim is to uncoil the helix by decreasing progressively the torsion κ 1 and the bending κ 3 . The method relies on the decomposition reported in Eq. 48. At each new step, the current strain fields κ are projected onto the 855 constructive remarks.

Appendix A. Appendices

Appendix A.1. Tools for the geometrically exact approach Appendix A.1.1. Lie group 860

The coordinate frame {t i } is right-handed and orthonormal, the determinant of R(s) is equal to 1. The inverse of an element R is given by its transpose R -1 = R T , and the neutral element is the 3×3 identity matrix I d . The Lie group SO(3) is a differ-865 entiable manifold of dimension 3 plunged into the vector space R 3×3 (Fig. 2). The non-linearity of the body geometry is related to the non-commutativity of the matrix multiplication, i.e. to the composition of rotations. This non-linearity implies that the manifold is globally curved (Fig. 2), even if the tangent space T R SO(3) at each element R has a vector space structure. Moreover, all the tangent spaces are isomorph to the tangent space at the neutral element called the Lie algebra so [START_REF] Jayne | Swimming in constricting (elaphe g. guttata) and nonconstricting (nerodia fasciata pictiventris) colubrid snakes[END_REF]. These properties will be important when the equilibrium configurations will be sought in section 4.

Appendix A.1.2. Unitary quaternion

Quaternions of norm one, also called unit quaternions, are in bijection with points of the real 3D sphere. Each unit quaternion Q are defined by the components Q = q 0 + q, with q = q 1 i + q 2 j + q 3 k and (q 0 , q 1 , q 2 , q 4 ) ∈ R 4 so that QQ * = 1. Each unit quaternion parametrizes a rotation R(s) of a vector u i by the following left and right multiplications

Using the properties of the quaternion, the rotation matrix can be then calculated by the following formula

Equation A.2 is analog to the Rodrigues' formula [START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF] used for computing the exponential map for constant skew-symmetric matrix. Hence, equation A.2 preserves the group structure of the rotation matrix like the exponential map. Combining equations 8 and A.1, the evolution equation for the quaternion Q(s) can be written into the following form

with the skew-symmetric matrix

Center of mass configuration

The center of mass 0 p com relative to the head frame F 0 is given by

with M the net mass of the body, L the length of the body, and S b = πab the surface of the section. The configuration of the center of mass reads

with I d the identity matrix. The wrench is transported to one frame to another via the transport of the density of buoyancy forces giving

the density of force and torque exerted on the section indexed by the variable s, 910 and expressed in the frame F G . These relation can be rewritten into a matrix form

with 0 3 a 3 × 3 null-matrix, and p the matrix form of the vector product p × (•) ( Eq. 9). The 6 × 6 matrix introduced in equation A.8 corresponds to 915 the transpose of the adjoint operator Ad g -1 . In the Lie group theory, the adjoint operator Ad g -1 transports twists from one coordinate frame to another [START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF], and by virtue of the conservation of the net work, the transposed operator Ad T g -1 allows for the 920 transport of wrenches. Therefore, the condensed form of Eq.A.8 is given by

. Chebyshev colocation method

In our code, each function is expressed using the Chebyshev series expansion at order N as

with κ k the amplitude of the Chebyshev polynomial T k (s) defined by recurrence relation

with T 0 = 1 and T 1 = x. The abscissae are discretized on a Chebyshev grid defined by s i = cos(θ i )

The boundary condition at x 0 can be written into a matrix form

with β a 4 × N matrix and (Q) 0 ∈ R 4 the values of the quaternion at x 0 . The matrix D being non invertible, the recombination method consists of writing the solution under the following form

with P the basis of the null space of the matrix β, z the component of the basis vectors of the null space, and R the pseudo-inverse of β. The unknowns are the components of z that satisfy the solvable linear system [START_REF] Marsden | Introduction to mechanics and symmetry: a basic exposition of classical mechanical systems[END_REF]) Finally, the left-hand side matrix is invert to ob-970 tain z that is used in Eq.A.15.

Appendix A.3. Notes on the calculation of the immersed surface Appendix A.3.1. Integration domain

The domain D s is defined so that all the points of the ellipse plane formed by (0, t 2 , t 3 ) with q 1 = 0, satisfy both inequalities

The first inequality corresponds to the points in-975 side the ellipse, and the second one, to all the points below the water surface.

Appendix A.3.2. Linear mapping Since the jacobian determinant of transformation 23 is one, the calculation of the surface in the circle 980 domain reads

with the new domain of integration D defined by all the points (q 2 , q3 ) satisfying both inequalities

In the same way, the barycenter becomes

with qB the position of the barycenter of the im-985 mersed portion of the circle. Thanks to the inverse mapping L -1 , the weighted barycenter in the ellipse frame is given by

given by Eq.26.

Appendix A.4. Exact calculation of (D gs w)

The aim of this appendix is to calculate the relative variation of wrench density δw n as a function of the vector δg T s without approximation. Only the Appendix A.5. Convergence to the solution for the optimization method 1015 By construction, the quantity (t n ξ n ) T W n is homogeneous to a work, and it is always positive for t n > 0 ( Eq. 41), ensuring a convergence to the minimum. Consequently, this condition states that the energy reduction U n -U n+1 is smaller than a 1020 percentage c 1 of an extracted work during the iteration, with the linear estimation of the work given by (t n ξ n ) T W n . Thus, the step size t n corresponds here to a time scale, and the step direction ξ n is the unitary twist coordinate of a virtual body velocity 1025 of the mobile frame F 0 [START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF].

Appendix A.6. Jacobian matrix (D κ w N ) From Eq. 33 and 32, we already know how to calculate a variation of the wrench as a function of a virtual body shape displacement δ( 0 ψ s ). Conse-1030 quently, we have to find a linear relation between δ( 0 ψ s ) and δκ. First, the evolution equation of the matrix 8 has to be generalized to rigid body transformation by combining Eq. 8 and 12 [START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF] ( 0 g s ) = 0 g s ξ (A.26) with the body strain ξ T = (κ T , s T ) associated with 1035 the screw-symmetric matrix ξ defined by ξ = κ s 0 1,3 0 (A.27) with 0 1,3 a 1×3 null-vector, and s T = (1, 0, 0). The relation between δ( 0 ψ s ) and δξ reads [START_REF] Marsden | Introduction to mechanics and symmetry: a basic exposition of classical mechanical systems[END_REF] δξ = ad ξ δ( 0 ψ s ) + d ds δ( 0 ψ s ) (A.28)