Emmanuel Cottanceau

Olivier Thomas
email: olivier.thomas@ensam.eu

Philippe Véron

Marc Alochet

Renaud Deligny

A finite element/quaternion/asymptotic numerical method for the 3D simulation of flexible cables

Keywords:

In this paper, a method for the quasi-static simulation of flexible cables assembly in the context of automotive industry is presented. The cables geometry and behavior encourage to employ a geometrically exact beam model. The 3D kinematics is then based on the position of the centerline and on the orientation of the cross-sections, which is here represented by rotational quaternions. Their algebraic nature leads to a polynomial form of equilibrium equations. The continuous equations obtained are then discretized by the finite element method and easily recast under quadratic form by introducing additional slave variables. The asymptotic numerical method, a powerful solver for systems of quadratic equations, is then employed for the continuation of the branches of solution. The originality of this paper stands in the combination of all these methods which leads to a fast and accurate tool for the assembly process of cables. This is proved by running several classical validation tests and an industry-like example.

Introduction

During the last decades, the room available in car vehicles (e.g. in engine compartment) has plummeted because of the rapid development of on-board electronics. As a result, a need for very accurate numerical tools for design has appeared in automotive industry. In the meantime, a fast computation is necessary so that design duration remains suitable for industry. In this context, flexible pieces represent an outstanding challenge since, unlike most of car pieces, they cannot be modeled as rigid body solids in CAD software. This paper focuses on a specific type of flexible piece, namely electrical cables. Cables have a complex structure. A wire is made up of copper filaments wrapped in an elastomer duct. These wires are most of the time gathered in bundles which are themselves surrounded by various protections such as tape, PVC tube … Moreover, the full cable is often constituted of several drifted cable pieces forming a system with a complex geometry.

Due to its slender shape and its flexibility, one can consider a simple cable as a beam undergoing large displacements and large rotations. There exist several theories accounting for nonlinear beams, see Refs. [START_REF] Shabana | Application of the absolute nodal coordinate formulation to large rotation and large deformation problems[END_REF][START_REF] Hegemier | A nonlinear dynamical theory for heterogeneous, anisotropic, elastic rods[END_REF][START_REF] Hodges | A mixed variational formulation based on exact intrinsic equations for dynamics of moving beams[END_REF]. The most widely used theory is the geometrically exact beam model whose founding principles were established by Reissner [START_REF] Reissner | On one-dimensional finite-strain beam theory: the plane problem[END_REF][START_REF] Reissner | On finite deformations of space-curved beams[END_REF] and further generalized by Simo [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF]. Various finite element formula-tions (FEM) have been presented to solve these equations numerically. The notable works of [START_REF] Jelenić | Geometrically exact 3d beam theory: implementation of a strain-invariant finite element for statics and dynamics[END_REF][START_REF] Cardona | A beam finite element non-linear theory with finite rotations[END_REF][START_REF] Ibrahimbegović | Computational aspects of vector-like parametrization of three-dimensional finite rotations[END_REF][START_REF] Simo | A three-dimensional finite-strain rod model. part ii: computational aspects[END_REF][START_REF] Ritto-Corrêa | On the differentiation of the Rodrigues formula and its significance for the vector-like parameterization of Reissner-Simo beam theory[END_REF] and more recently [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF][START_REF] Ghosh | Consistent quaternion interpolation for objective finite element approximation of geometrically exact beam[END_REF] can, among others, be listed.

At the heart of all these formulations, the rotation parameterization is of paramount importance in the numerical models. The 3D rotations modeling indeed is not an easy task especially when computational efficiency is sought. The rotational vector-like parameterization used by many authors features only 3 parameters (minimum set in 3D). However, as no parameterization of less than 4 parameters can be singularity-free [START_REF] Stuelpnagel | On the parametrization of the three-dimensional rotation group[END_REF], this choice poses several numerical limitations and lacks robustness. A powerful alternative consists in using quaternions, a set of 4 singularity-free parameters. Firstly used only for storage in the numerical models, Zupan et al. [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF] have recently shown their utility when used as primary variables. They also have developed a model without rotation matrices exploiting the high potential of quaternion algebra, and very efficient for numerical purposes [START_REF] Zupan | On a virtual work consistent three-dimensional Reissner-Simo beam formulation using the quaternion algebra[END_REF]. In addition, quaternions offer an original description of rotations since they substitute the usual trigonometric functions by algebraic variables and lead to polynomial equilibrium equations.

The finite element method, very adapted to the assembly of complex geometries such as electrical cables ones, is applied to the continuous equilibrium equations. The algebraic system obtained is then generally solved by a classical predictor-corrector method (PCM) [START_REF] Criesfield | Non Linear Finite Element Analysis of Solids and Structures[END_REF]. Even if the appearance of arc-length methods [START_REF] Riks | An incremental approach to the solution of snapping and buckling problems[END_REF] has considerably enhanced the robustness of this type of method, they often require to choose a step size which may be very tricky for the user. A sufficiently small step size allows to compute even highly nonlinear part of equilibrium branches but in return may impractically increase the computation time, while a larger step size may spoil the convergence. Taking advantage of the polynomial form of the system of equations obtained when using quaternion parameters an alternative consists in replacing the PCM by the asymptotic numerical method (ANM) firstly presented by Damil and Potier-Ferry [START_REF] Damil | A new method to compute perturbed bifurcations: application to the buckling of imperfect elastic structures[END_REF] and by Cochelin [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF]. The ANM is a very powerful solver for quadratic problems and it overcomes all the drawbacks of the PCM. This technique indeed is very robust, does not require any tuning parameters and is thus well suited for an industrial use. In addition, Cochelin and Medale [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] have equipped the method with a bifurcation detector and improved its efficiency in the vicinity of bifurcation points.

Combining quaternions with the ANM has already been set up on a rod model discretized with a finite difference scheme by Lazarus et al. [START_REF] Lazarus | Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method[END_REF], which have got very promising results. We propose here to set up the technique on the finite-element based geometrically exact beam model, in what constitutes the main originality of this paper. Validations and illustrations of the method on very intricate problems are provided and discussed. A critical evaluation and future researches are presented by way of conclusion.

Governing equations

In this section, the classical quasi-static formulation of the geometrically exact beam model, based on the rotation vector, is firstly presented. It enables to explain all the main ingredients of the model and to discuss their physical meaning. Secondly, the equations are modified by using quaternions instead of the rotation vector. This leads to the formulation which serves our numerical model, presented in part 3.

The geometrically exact beam model

In the geometrically exact beam model, a beam is described by defining a family of cross-sections whose centroids form a curve called the centerline of the beam. The kinematic variables essential to this description are introduced following the notations used in Ref. [START_REF] Géradin | Flexible Multibody Dynamics: a Finite Element Approach[END_REF]. Let us consider a beam of initial length L, with an arbitrary cross-section Ω, as depicted Fig. 1. Let (u 1 , u 2 , u 3 ) be a fixed Cartesian (global) frame. The current position of the centerline in this frame is described by the vector field x 0 (s) which is a function of the curvilinear abscissa s ∈ [0, L] along the beam axis. A material (local) frame (e 1 (s), e 2 (s), e 3 (s)) is introduced to define the cross-section at abscissa s. The vectors e 2 and e 3 span the cross-section while the vector e 1 remains normal to the crosssection for every deformed configurations. It is essential to note that e 1 is not necessarily tangent to the centerline of the beam such that shear deformation is taken into account (Timoshenko model). The reference configuration is chosen such that the beam is unstressed. However, in this state, the reference position of the centerline is an arbitrary curve and not necessarily a straight line, thus accounting for the initial curvature of the beam. Its initial position is then a function of s denoted X 0 (s), such that the displacement of any point of the centerline is x 0 (s) -X 0 (s). Similarly, the material frame in the reference configuration depends on s and is denoted (E 1 (s), E 2 (s), E 3 (s)). To end up with the definition of the kinematic variables, let us introduce the two rotation operators R 0 (s) and R(s) which depict the orientation of the material frames in the global frame in the initial and in the current configuration respectively, that being E i (s) = R 0 (s)u i and e i (s) = R(s)u i , i ∈ {1, 2, 3}. As R 0 (s) defines the initial curvature of the beam it is constant through the deformation. With these notations, the position vector in the spatial frame of any point M ′ of the undeformed beam located in the section at s writes

X(s, X 2 , X 3 ) = X 0 (s) + R 0 (s)Y(X 2 , X 3 ), (1)
where Y(X 2 , X 3 ) = [0 X 2 X 3 ] T is the material position of M ′ in the cross-section. Under the assumption that cross-sections remain plane and do not undergo any deformations along the transformation, Y(X 2 , X 3 ) is constant and M ′ after deformation becomes M ″ whose position is given by

x(s, X 2 , X 3 ) = x 0 (s) + R(s)Y(X 2 , X 3 ). ( 2
)
The current configuration of the beam is then completely characterized by the position of the centerline x 0 (s) and the orientation of the crosssections R(s). One recovers here that the kinematic variables depend solely on the curvilinear abscissa s as for any beam model. As it is clearly implied, the dependency in s is dropped in the following. Furthermore, the derivative with respect to s is from now on denoted with the prime symbol ′ . For any rotation operator R, it is shown that the matrix R T R ′ is skew-symmetric [START_REF] Géradin | Flexible Multibody Dynamics: a Finite Element Approach[END_REF]. Any 3 × 3 skew-symmetric matrix A is made of only 3 independent components which can be gathered in a vector a = vect(A) = [a 1 a 2 a 3 ] T . The skew-symmetric matrix associated to a is denoted ã and is equal to:

ã = ⎡ ⎢ ⎢ ⎢ ⎣ 0 -a 3 a 2 a 3 0 -a 1 -a 2 a 1 0 ⎤ ⎥ ⎥ ⎥ ⎦ . ( 3
)
With these notations and hypothesis, the material translational and rotational strain measures are equal to [START_REF] Crisfield | Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation[END_REF].

𝚪 = [ Γ 1 Γ 2 Γ 3 ] T = R T x ′ 0 -𝚪 0 (4a) K = [ 𝜅 1 𝜅 2 𝜅 3 ] T = vect(R T R ′ ) -K 0 , ( 4b
)
where the intrinsic translational strain 𝚪 0 = R T 0 X ′ 0 = u 1 and the intrinsic rotational strain K 0 = vect(R T 0 R ′ 0 ) solely depend on the initial configuration of the beam. Let us point out that in the reference configuration, the strain measures 𝚪 and K are both equal to 0: the assumption of an initially unstressed configuration is thus well recovered by these measures definition.

To interpret these measures, let us consider the straightforward case in which the beam centerline is a straight line in the reference configuration. The intrinsic strains then obey 𝚪 0 = u 1 and K 0 = 0. x ′ 0 (s) is, by definition, a vector tangent to the centerline of the beam at abscissa s in the current configuration. In expression (4a), the operator R T maps back this vector projection from the material frame to the fixed frame. The measure 𝚪 then consists in comparing this quantity to its reference value u 1 . A fully equivalent interpretation is made by recalling that the scalar product of 2 vectors may be geometrically interpreted as the projection of one vector in the direction of the other. It follows that

Γ 1 = e T 1 x ′ 0 -1, Γ 2 = e T 2
x ′ 0 and Γ 3 = e T 3 x ′ 0 are the comparison of the projection of x ′ 0 (s) on the material frame to e 1 in the 3 material directions. As the beam is parameterized by the curvilinear abscissa in the initial configuration, the norm of x ′ 0 (s) is equal to 1 if the length of the centerline remains unchanged, is greater than 1 in case of elongation of the centerline and smaller than 1 in case of shortening. With regard to all these observations, we eventually infer that Γ 1 is the axial strain of the neutral axis, while Γ 2 and Γ 3 represent the shear strains in directions e 2 and e 3 .

In the case of a pure rotational transformation R, R T R ′ represents the gradient of the transformation along the beam axis. As this quantity is a skew-symmetric matrix, it is equivalent to write

R T R ′ Y and vect(R T R ′ ) × Y = K × Y.
It follows directly that in (4b) 𝜅 1 depicts the torsional strain while 𝜅 2 and 𝜅 3 represent the curvature strains around e 2 and e 3 respectively.

The expression of the aforementioned rotation operator R is generally obtained from a reduced representation. The choice of this representation is actually a crucial issue since various criteria compete against each other to get the best set of parameters: numerical efficiency (number of parameters), mathematical form, existence of singularities, geometric interpretation. Let us recall that any 3D rotation of a rigid body solid is made up of 3 independent parameters. The 3-parameters rotational vector 𝝑 presented in this paragraph thus constitutes the minimal set of parameters. We choose it here for its easy geometrical interpretation since it is defined by

𝝑 = 𝜗n = [ 𝜗 1 𝜗 2 𝜗 3 ] T , (5)
for a rotation of angle 𝜗 around the unit axis n. The rotation operator then writes in terms of this parameterization as [START_REF] Géradin | Flexible Multibody Dynamics: a Finite Element Approach[END_REF] R(𝝑

) = I 3 + sin(𝜗) 𝜗 θ + 1 -cos(𝜗) 𝜗 2 θθ . (6)

Virtual work principle

Using the quantities introduced in the previous section and the Cartesian rotational vector 𝝑 defined in [START_REF] Reissner | On finite deformations of space-curved beams[END_REF] as rotational parameters, the virtual work principle for a 3D shear elastic beam is stated here. Following Reissner's work [START_REF] Reissner | On one-dimensional finite-strain beam theory: the plane problem[END_REF][START_REF] Reissner | On finite deformations of space-curved beams[END_REF] and the generalization of Simo [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF], it writes

∫ L 0 (N • 𝛿𝚪 + M • 𝛿K) ds = ∫ L 0 ( n e • 𝛿x 0 + m e • 𝛿𝝑 ) ds + [ N e • 𝛿x 0 ] L 0 + [ M e • 𝛿𝝑 ] L 0 , (7)
with the generalized stress resultant in force and moment N and M, the external applied force and moment per unit length n e and m e , the external end force and moment N e and M e and the virtual positional and rotational vectors 𝛿x 0 and 𝛿𝝑. The left-hand of the equation is the opposite of the virtual work of internal forces -𝛿 i and the righthand is the virtual work of external forces

𝛿 e . N = [N T 2 T 3 ] T and M = [M t M 2 M 3 ]
T are work conjugates to the virtual translational strain 𝛿𝚪 and the virtual rotational strain 𝛿K respectively. Therefore, N stands for the axial force resultant, T 2 and T 3 the transverse force resultants in respective directions e 2 and e 3 , M t the torsional moment and M 2 and M 3 the bending moments around axis e 2 and e 3 .

In the scope of this paper, only linear elastic materials are considered. In this particular case, the stress-strain relationship states proportionality between stress resultant vector 𝚺 = [M T N T ] T and strain resultant vector  = [𝚪 T K T ] T that be

𝚺 = C = diag(C N , C M ). (8)
In [START_REF] Cardona | A beam finite element non-linear theory with finite rotations[END_REF], C N = diag(EA, GA 2 , GA 3 ) and C M = diag(GJ, EI 2 , EI 3 ) are the two 3 × 3 stress-strain submatrices corresponding solely to the force components and the moment components respectively. In these expressions, E and G are respectively the Young's modulus and the shear modulus of the material; A is the cross-sectional area; A 2 and A 3 are the shear areas in directions e 2 and e 3 ; I 2 , I 3 are the second moments of area with respect to axes e 2 , e 3 ; J is the polar moment of area.

Let us notice that this model can take into account the case of the origin of the material frame not being on the centerline effortlessly. In this case, out-of-diagonal terms simply appear in the stress-strain matrix C.

Finally, the virtual strains in [START_REF] Jelenić | Geometrically exact 3d beam theory: implementation of a strain-invariant finite element for statics and dynamics[END_REF] are given in terms of the kinematic variables by the equations [START_REF] Géradin | Flexible Multibody Dynamics: a Finite Element Approach[END_REF]:

𝛿𝚪 = R T 𝛿x ′ 0 + RT x ′ 0 T𝛿𝝑 = R T 𝛿x ′ 0 + ( Γ + Γ0 )T𝛿𝝑, 𝛿K = (R T R ′ T + T ′ )𝛿𝝑 + T𝛿𝝑 ′ = (( K + K0 )T + T ′ )𝛿𝝑 + T𝛿𝝑 ′ . (9)
In [START_REF] Ibrahimbegović | Computational aspects of vector-like parametrization of three-dimensional finite rotations[END_REF], T is a tangent operator which can be expressed as a function of 𝝑 only under the form:

T = I 3 + cos(𝜗) -1 𝜗 2 θ + 1 𝜗 2 ( 1 - sin(𝜗) 𝜗 ) θθ . ( 10
)

Strong form of equilibrium equations

Introducing expressions (9) into the virtual work principle [START_REF] Jelenić | Geometrically exact 3d beam theory: implementation of a strain-invariant finite element for statics and dynamics[END_REF], integrating by parts with some rearranging and using the fundamental lemma of calculus of variations lead to the following equilibrium equations:

{ (RN) ′ + n e = 0, (RM) ′ + x ′ 0 × (RN) + m e = 0, ∀s ∈ [0 L] ( 11
)
along with the natural boundary conditions { RN -N e = 0 at s = 0, L,

RM -M e = 0 at s = 0, L, (12)
or the essential boundary conditions

{ x 0 = x d at s = 0, L, 𝝑 = 𝝑 d at s = 0, L, (13)
where x d and 𝝑 d are respectively imposed positions and rotations.

Quaternion parameterization of rotations

As explained at the end of paragraph 2.1, the choice of rotation parameters is of paramount importance in numerical models. In a total Lagrangian framework (i.e. using total rotations), the rotational vector 𝝑 prevents from having angles greater than 2𝜋. It is clear from equation [START_REF] Simo | A three-dimensional finite-strain rod model. part ii: computational aspects[END_REF] that the tangent operator T indeed becomes rank-deficient around 2𝜋-angles and leads to a singularity. Procedures exist to circumvent this flaw [START_REF] Ibrahimbegović | Computational aspects of vector-like parametrization of three-dimensional finite rotations[END_REF], however, we find it more expedient to use an other parameterization not requiring any special procedure. As demonstrated in Ref. [START_REF] Stuelpnagel | On the parametrization of the three-dimensional rotation group[END_REF], the minimal set of parameters avoiding all singularities is made up of 4 parameters. With regard to this principle, the 4-parameters quaternions appear as an optimal choice for the representation of rotations. In addition, quaternions feature a special algebra totally equivalent to the algebra of rotation operators but which proves computationally more efficient [START_REF] Salamin | Application of Quaternions to Computation with Rotations[END_REF]. Finally, they offer a polynomial representation of rotations in comparison to classical trigonometric ones. We will show in the following part how this property is used as a leverage in our numerical solution.

The next developments of this paragraph aim at giving the fundamentals on quaternions necessary to the comprehension of the method. The interested reader may consult for instance [START_REF] Vince | Quaternions for Computer Graphics[END_REF] for more details. The notations of [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF] on quaternion algebra are used here. A quaternion â is defined by

â = a 0 + ia 1 + ja 2 + ka 3 , ( 14
)
where a 0 , a 1 , a 2 , a 3 ∈ ℝ and the imaginary numbers i, j and k are linked by the identities

i 2 = j 2 = k 2 = ijk = -1. (15)
1, i, j and k form a basis of the set of quaternions, denoted ℍ. As a result, a quaternion may also be seen as an element of the fourdimensional Euclidean linear space ℝ 4 , whose components are a 0 , a 1 , a 2 and a 3 in the Euclidean basis made up of the 4 vectors [1 0 0 0] T , [0 1 0 0] T , [0 0 1 0] T , [0 0 0 1] T . In this paper, no distinctions are made between the quaternions and their projection into a basis of the Euclidean space so that we can smoothly write a quaternion under the vector form â = [a 0 a 1 a 2 a 3 ] T . To distinguish quaternions from the three-dimensional vectors of ℝ 3 , the hat symbol is used:

a ∈ ℝ 3 , â ∈ ℝ 4 .
Let us introduce the handy notation

â = a 0 + a. ( 16
)
In this expression, a 0 and a = ia 1 + ja 2 + ka 3 are called respectively the scalar and the vector part of the quaternion. This notation helps drawing a parallel between quaternions and complex numbers, the scalar part amounting to the real part and the vector part amounting to the imaginary part. It is also helpful when one wants to express the extension of a vector of ℝ 3 in the quaternion space ℍ. Indeed, any basis (i 1 , i 2 , i 3 ) of the three-dimensional Euclidean space can be extended into a basis of the four-dimensional Euclidean space ( î 0 , î1 , î2 , î3 ), with î0 = 1 + 𝟎 and îk = 0 + i k , k ∈ {1, 2, 3}. As a result, any vector v of ℝ 3 has a quaternion counterpart in ℍ, v = 0 + v. This quaternion with a zero scalar part is called a pure quaternion and depicts the extension of a vector into a quaternion. Conversely, a restriction operation may be defined, which transforms a pure quaternion into a vector. This operation is denoted Vec, so that for any pure quaternion p of ℍ, Vec(p) = p ∈ ℝ 3 . In the following, no distinctions will be made between pure quaternions and other quaternions in the notation. For a better comprehension, the reader just has to keep in mind that any vector extension in the quaternion space is a pure quaternion. Using the analogy with the complex numbers, a conjugate quaternion is defined as â * = a 0 -a. The set of quaternions is equipped with the inner product â • b = a 0 b 0 + a • b and the quaternion norm ‖â‖ = √ â • â. Furthermore, the three following elementary operations are also defined: the addition â + b = (a 0 + b 0 ) + (a + b), the scalar multiplication 𝜆â = 𝜆a 0 + 𝜆a for 𝜆 ∈ ℝ and the quaternion multiplication, denoted ⚬, which follows directly from the definition of a quaternion [START_REF] Stuelpnagel | On the parametrization of the three-dimensional rotation group[END_REF] and from the relations [START_REF] Zupan | On a virtual work consistent three-dimensional Reissner-Simo beam formulation using the quaternion algebra[END_REF]:

â ⚬ b = a 0 b 0 -a • b + (b 0 a + a 0 b + a × b), ( 17
)
where • and × are the scalar product and the cross product of vectors. This operation is associative but not commutative, so left and right multiplication must be distinguished. Finally, it is convenient for the further presented formulation to extend the cross product of vectors of ℝ 3 to pure quaternions. For two pure quaternions p and q, it is defined as

p × q = 0 + p × q = p × q. ( 18
)
To perform the rotation of a vector x into a vector y, that be y = Rx, a special set of quaternions, called Euler parameters, is employed. For a rotation of angle 𝜗 around an axis oriented by the unit vector n, Euler parameters are defined by

q = cos ( 𝜗 2 ) + n sin ( 𝜗 2 ) . ( 19
)
The rotation is then applied using the following formula:

ŷ = q ⚬ x ⚬ q * , ( 20
)
which is quadratic with respect to q and is totally equivalent to using a rotation matrix (see Appendix F). Furthermore, if one wants to apply a combination of two successive rotations of respective quaternions q1 and q2 , the rotational quaternion q2 ⚬ q1 , evidencing the multiplicative nature of quaternions, needs to be used in [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF], that be:

ŷ = (q 2 ⚬ q1 ) ⚬ x ⚬ (q 2 ⚬ q1 ) * . ( 21
)
Let us notice for practical purposes that the norm of the rotational quaternion defined in equation ( 19) is 1 (the so-called unit quaternions) obeying

q ⚬ q * = q * ⚬ q = 1 = 1 + 0. (22)
Besides, deriving this relation with respect to s leads to the following useful relations

q′ ⚬ q * = -(q ⚬ q′ * ), q′ * ⚬ q = -(q * ⚬ q′ ). (23)
In order to manipulate equations in quaternion space, let us also introduce the following helpful relations, stemming from ( 17), [START_REF] Damil | A new method to compute perturbed bifurcations: application to the buckling of imperfect elastic structures[END_REF] and the commutativity of inner product:

• for arbitrary quaternions p, q and r:

p • ( q ⚬ r) = q • ( p ⚬ r * ) = ( p ⚬ r * ) • q, p • ( r ⚬ q) = q • ( r * ⚬ p) = ( r * ⚬ p) • q, ( 24
)
Table 1 Free-end position and rotation of the cantilever beam under the end-moment

M L = 3𝜋 EI 3 L
for several values of n and Ne on the left; deformed shape during loading on the right.

• for pure quaternions â and b:

â ⚬ b + b * ⚬ â = 2â × b. ( 25
)

Rewriting of governing equations

The weak form of equilibrium equations ( 7) along with ( 4) and ( 9) are now expressed in terms of the quaternion parameters. The rewriting operation consists in replacing vectors by their quaternion equivalent x → x and the rotation operation by its expression in the quaternion space y = Rx → ŷ = q ⚬ x ⚬ q * . Most of the relations introduced in the paragraph 2.4 serve this change of rotation parameters.

First, the strain measures (4) are rewritten by combining equations ( 22), ( 23) and ( 25) to get the pure quaternion strains (see Appendix A for the full demonstration)

Γ = q * ⚬ x′ 0 ⚬ q -Γ 0 (26a) K = 2q * ⚬ q′ -K0 , ( 26b
)
with the intrinsic translational strain Γ0 = û1 and rotational strain K0 = 2q * 0 ⚬ q′ 0 expressed in quaternion space. The virtual strain expressions (9) in the quaternion space are obtained by differentiating equations (26a) and (26b) and rearranging them with [START_REF] Vince | Quaternions for Computer Graphics[END_REF] (see Appendix B or [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF] for the full demonstration) leading to

𝛿 Γ = q * ⚬ 𝛿x ′ 0 ⚬ q + 2q * ⚬ ( x′ 0 × ( 𝛿q ⚬ q * ) ) ⚬ q, ( 27a
) 𝛿 K = 2q * ⚬ ( 𝛿q ⚬ q * ) ′ ⚬ q, ( 27b
)
with the virtual positional (pure) quaternion 𝛿x 0 and the virtual rotational quaternion 𝛿q. The resultant of the external moment is conjugated to the virtual rotational vector in the virtual work principle [START_REF] Jelenić | Geometrically exact 3d beam theory: implementation of a strain-invariant finite element for statics and dynamics[END_REF]. To express the resultant of the external moment in the quaternion space, a relationship between the virtual rotational vector and the virtual rotational quaternion must be established. One demonstrates (see Appendix D) that we have

𝛿 θ = 2q * ⚬ 𝛿q. ( 28
)
Substituting the quaternion expressions ( 26)-( 28) to their vectorial counterpart in the virtual work principle (7) and using the relations [START_REF] Salamin | Application of Quaternions to Computation with Rotations[END_REF] to obtain the resultant of the external moment, the new virtual work principle extended in quaternion space writes

∫ L 0 ( N • 𝛿 Γ + M • 𝛿 K) ds = ∫ L 0 ( ne • 𝛿x 0 + 2 ( q ⚬ me ) • 𝛿q ) ds + [ Ne • 𝛿x 0 ] L 0 + [ 2 ( q ⚬ Me ) • 𝛿q ] L 0 . ( 29
)
As noticed in paragraph 2.4, the four components of a rotational quaternion q are not independent but constrained by the condition of unity q • q -1 = 0. This constraint is taken into account by the method of Lagrangian multipliers. Following the approach of Zupan et al. [START_REF] Zupan | On a virtual work consistent three-dimensional Reissner-Simo beam formulation using the quaternion algebra[END_REF], the Lagrangian multiplier appended to this constraint and denoted 𝜇 is considered as a variable of the problem. The expression intro-Fig. 2. Influence of element order on the algorithm accuracy for the example of paragraph 5.1 (number of elements circled on second curve). duced in the new virtual work principle (29) hence is the term 𝜇(q • q -1) differentiated and integrated along the length of the beam, that be:

𝛿 𝜇 = 𝛿 ∫ L 0 𝜇(q • q -1) ds = ∫ L 0 ( (q • q -1)𝛿𝜇 + 2𝜇q • 𝛿q ) ds. ( 30
)
First, introducing the expressions of 𝛿 Γ and 𝛿 K (27a) and (27b), then rearranging [START_REF] Shabana | Dynamics of Multibody Systems[END_REF] with equations ( 24)-( 25), integrating by parts and applying the fundamental lemma of calculus of variations yield the strong form of equilibrium equations (see Appendix C for the full demonstration):

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Vec ( (q ⚬ N ⚬ q * ) ′ ) + n e = 𝟎, Vec ( (q ⚬ M ⚬ q * ) ′ + x′ 0 × (q ⚬ N ⚬ q * ) ) + me = 𝟎, 𝜇 = 0, q • q -1 = 0, ∀s ∈ [0, L] (31)
along with the natural boundary conditions

⎧ ⎪ ⎨ ⎪ ⎩ Vec ( q ⚬ N ⚬ q * ) -N e = 𝟎 at s = 0, L, Vec ( q ⚬ M ⚬ q * ) -Me = 𝟎 at s = 0, L, (32)
or the essential boundary conditions

{ x 0 = x d , at s = 0, L, Vec(q) = Vec(q d ) at s = 0, L, ( 33
)
where qd is the prescribed rotation under quaternion form [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF]. It ensues from (31) that 𝜇(s) = 0 is solution of the continuous problem while the remaining equations are strictly equivalent to the original form [START_REF] Ritto-Corrêa | On the differentiation of the Rodrigues formula and its significance for the vector-like parameterization of Reissner-Simo beam theory[END_REF]. Introducing the Lagrangian multiplier in the virtual work principle thus do not alter the equilibrium equations. Besides, it should be noticed that only the vector part of the quaternion is prescribed in the boundary conditions [START_REF] Zupan | On higher order integration of angular velocities using quaternions[END_REF]. The scalar component is indeed automatically known from the fulfillment of the unity constraint (31d). Another equivalent possibility is to replace (33b) by q = qd with |q d | = 1 and 𝜇 = 0 at s = 0, L. At this point, the continuous form of equilibrium equations of the problem is established. The equations are expressed in both weak form through equations ( 29)-( 30) and strong form through equations [START_REF] Romero | The interpolation of rotations and its application to finite element models of geometrically exact rods[END_REF] along with [START_REF] Bauchau | Interpolation of finite rotations in flexible multi-body dynamics simulations[END_REF], knowing relations (26a)-(27b). A vector equation is made up of 3 scalar equations. As a result, the set of governing equations in strong form ( 31) is composed of 8 (3 + 3 + 1 + 1) scalar equations. Besides, the unknowns of the problem are the 8 following fields of s: the 4 components of the rotational quaternion q = [q 0 (s) q 1 (s) q 2 (s) q 3 (s)] T depicting the rotation of the sections, the 3 components of the current position vector of the cross-sections centroids x 0 (s) = [x 1 (s) x 2 (s) x 3 (s)] T and the Lagrangian multiplier 𝜇(s). As a result, our problem is a well-posed 8 equations -8 unknowns problem and can henceforth be solved numerically.

For convenience, let us gather the 8 scalar unknown fields under a single vector field

z(s) = [ x 0 (s) T q(s) T 𝜇(s) ] T . ( 34
)
This allows to write the weak form of equilibrium equations ( 29)-( 30) under the compact form

∫ L 0 𝚺 • 𝛿 ds = ∫ L 0 (f e + f 𝜇 ) • 𝛿z ds + [ F e • 𝛿z ] L 0 , (35)
with

f e = [ n T e (2q ⚬ me ) T 0 ] T , F e = [ N T e (2q ⚬ Me ) T 0 ] T , f 𝜇 = [ 𝟎 (2𝜇q) T q • q -1 ] T . ( 36
)
In equation [START_REF] Schulz | Non-linear spatial Timoshenko beam element with curvature interpolation[END_REF], the strain and the stress resultant vector are still

equal to  = [ 𝚪 T K T ] T and 𝚺 = [ N T M T ] T but from this point for- ward are expressed in quaternion space, that be 𝚪 = Vec( Γ), K = Vec( K), N = Vec( N) = C N Vec( Γ) and M = Vec( M) = C M Vec( K)
, with Γ and K expressed as in (26a) and (26b). Equations (26a) and (26b) show clearly that Γ, K then N, M through (8) and as a consequence , 𝚺 are polynomial functions of z of degree 3. Similarly, equations (27a) and (27b) expose that 𝛿 is a polynomial function of z of degree 4. This in addition to the quadratic form of expressions [START_REF] Cochelin | Méthode Asymptotique Numérique[END_REF] evidences that the governing equation ( 35) is a polynomial equation of z of degree 7. Thus, using quaternion parameters instead of the rotational vector has allowed to transform the nonpolynomial system of 6 equations/unknowns (11) (it contains trigonometric functions due to the expression of the rotation and the tangent operators ( 6) and ( 10)) in a polynomial system of 8 equations -8 unknowns. How our method takes advantage of this specificity is shown in the following parts.

Finite element method

Discrete form of equations

The finite element method is now applied to the weak form of governing equations [START_REF] Schulz | Non-linear spatial Timoshenko beam element with curvature interpolation[END_REF] to obtain the discrete equilibrium equations. The beam is cut in N e elements, which for sake of simplicity are supposed of equal length L e = L/N e . Each element is composed of n + 2 equallyspaced nodes: n internal nodes and one node at each extremity of the element. In each element e, the unknown field z (e) value at node i ∈ {1, … , n + 2} is denoted z i .

The virtual work principle [START_REF] Schulz | Non-linear spatial Timoshenko beam element with curvature interpolation[END_REF] is discretized using a standard Galerkin approach in which the trial functions x 0 , q and 𝜇 and the test functions 𝛿x 0 , 𝛿q and 𝛿𝜇 are interpolated similarly. This choice is not obvious when dealing with rotational quantities and for this reason it is discussed briefly in section 3.4. The same interpolation order is used for all variables, since numerical studies have shown that a different order for each variable may lead to numerical divergence [START_REF] Bottasso | A non-linear beam space-time finite element formulation using quaternion algebra: interpolation of the Lagrange multipliers and the appearance of spurious modes[END_REF]. The interpolation then writes

z (e) (s) = n+2 ∑ i=1 N i (s)z i . ( 37
)
We use isoparametric elements so that the interpolation functions N i in equation ( 37) are the same as the shape functions and are chosen as standard Lagrange-type polynomials. Equation ( 37) may be written under the compact matrix form

z (e) (s) = P(s)u (e) , ( 38
)
where, if I 8 is the 8 × 8 identity matrix,

P(s) = [ N 1 (s)I 8 N 2 (s)I 8 … N n+2 (s)I 8 ] (39)
is an interpolation matrix of size 8 × 8(n + 2) and u (e) contains the unknowns nodal values of element e:

u (e) = [ (z 1 ) T (z 2 ) T … (z n+2 ) T ] T . ( 40
)
As the elementary fields all depend on s and as all the equations are written inside the element, except when needed for a better comprehension, the notations e and s are dropped in the following. As said earlier, the test functions are interpolated similarly to the trial functions so that we also have 𝛿z (e) (s) = P(s)𝛿u (e) . However, it is convenient for the writing to introduce a 2nd interpolation matrix Q of size 12 × 8(n + 2) for the test functions such that

[ (𝛿x ′ 0 ) T (𝛿q ′ ) T (𝛿q) T 𝛿𝜇 ] T = Q 𝛿u (e) , ( 41
)
with

Q = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ N ′ 1 I 3 𝟎 𝟎 … N ′ (n+2) I 3 𝟎 𝟎 𝟎 N ′ 1 I 4 𝟎 … 𝟎 N ′ (n+2) I 4 𝟎 𝟎 N 1 I 4 𝟎 … 𝟎 N (n+2) I 4 𝟎 𝟎 𝟎 N 1 … 𝟎 𝟎 N (n+2) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 42
)
We can henceforth write under discrete form the continuous equation [START_REF] Schulz | Non-linear spatial Timoshenko beam element with curvature interpolation[END_REF]. To get the discrete version of the left-hand side, we search for an expression of the virtual strains 𝛿 as a function of 𝛿u (e) . From the expression of the virtual strains (27a) and (27b) and from [START_REF] Bathe | Large displacement analysis of three-dimensional beam structures[END_REF], one shows (see Appendix E) that

[ N M ] • 𝛿 = Q T g(N, M) • 𝛿u (e) , ( 43
)
with the vector of 12 components

g(N, M) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Vec [ q * ⚬ N ⚬ q] 2q ⚬ M 2q ⚬ ( N × ( Γ + Γ0 ) + M × ( K + K0 ) ) + 2q ′ ⚬ M 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 44
)
Q T g(•,
•) may be interpreted as an elementary discrete gradient of the internal work, since it links 𝚺 • 𝛿 to the virtual nodal values.

The discrete version of the right-hand side of equation ( 35) in a given element stems directly from equations ( 38) and ( 41) and the virtual work principle may now be written in each element in a discrete manner. The elementary virtual work of internal forces, the elementary virtual constraint and the elementary virtual work of external forces respectively write

𝛿 (e) i = -∫ L e 0 𝚺 • 𝛿 ds = 𝛿u (e) • ( -∫ L e 0 Q T g(N, M) ds ) , 𝛿 (e) 𝜇 = -∫ L e 0 f 𝜇 • 𝛿z ds = 𝛿u (e) • ( -∫ L e 0 Q T c 𝜇 ds ) , 𝛿 (e) e = ∫ L e 0 f e • 𝛿z ds = 𝛿u (e) • ∫ L e 0 P T f e ds. ( 45
)
In (45), we have introduced the vector

c 𝜇 = [ 0 0 2𝜇q T q • q -1 ] T .
Besides, for the special cases of end elements, the loads on end-sections must be added to external work. They write respectively -F e (0) • 𝛿z(0) = 𝛿u (1) • (-P T (0)F e (0)) and F e (L) • 𝛿z(L) = 𝛿u (N e ) • (P T (L e )F e (L)) for the first and the last element. Summing these virtual quantities, the elementary virtual work principle takes the form

𝛿 (e) i + 𝛿 (e) e + 𝛿 (e) 𝜇 = 𝛿u (e) • (  (e) e - (e) i ) , ( 46
)
where  (e) i and  (e) e are respectively the elementary internal force vector (in which is included the unity constraint) and the elementary external force vector. They write

 (e) i = ∫ L e 0 Q T f (e) i ds,  (e) e = ∫ L e 0 P T f e ds. ( 47
)
As in (47), the vector  (e) e requires a special treatment at beam ends: the end forces -P T (0)F e (0) and P T (L e )F e (L) must be included in the first and the last element respectively. In this equation, the vector f (e) i is equal to

f (e) i = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Vec [ q ⚬ N ⚬ q * ] 2 q ⚬ M 2 q ⚬ ( N × ( Γ + Γ 0 ) + M × ( K + K0 ) ) + 2q ′ ⚬ M + 2𝜇 q q • q -1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 48
)
Few comments need to be made on this expression because it is at the heart of the computer code, since the internal force vector is directly linked to f (e) i . First, let us notice that it does not make use of the rotation matrix R which is replaced by computationally more efficient quaternion products (see for instance [START_REF] Salamin | Application of Quaternions to Computation with Rotations[END_REF]). This product is thus directly implemented in the code. Afterward, recalling that N = C N 𝚪, M = C M K and that the strains expressions are given by (26a) and (26b), one evidences that this expression is still polynomial of degree 7 but this time function of elementary fields. Finally, as quaternion expressions have 4 components and vectorial expression have 3, one observes that f (e) i is a vector of 12 (3 + 4 + 4 + 1) components. After being multiplied by Q T , it gives the vector  (e) i whose size is equal to the number of nodal values of the element, namely 8(n + 2).

In equations (47), the integrals over the element are evaluated numerically by a classical Gauss-Legendre quadrature of order N g . Integration is performed in the parent element characterized by the coordinate 𝜉 varying between -1 and 1. The Jacobian of the transformation from physical to parent space is the same for all elements and is equal to J = ds/d𝜉 = L e /2. One of the major problems of Timoshenko-like elements is the presence of shear locking [START_REF] Zienkiewicz | The Finite Element Method: Solid Mechanics[END_REF]. In order to avoid this numerical deficiency, reduced integration is used. It consists in choosing a quadrature order lower than the one required for exact numerical integration. In our case, for an element with n + 2 nodes, reduced integration is employed by using N g = n + 1 points in the quadrature formula [START_REF] Oak-Key | Reduced minimization theory in beam elements[END_REF].

At this point, all the necessary information has been given to compute elementary quantities. The last remaining step of the finite element method is to make an assembly operation on the elements to get the global discrete problem equations. For that purpose, let us define the global unknowns vector U containing all the degrees of freedom of the problem:

U = [ (x 1 0 ) T (q 1 ) T 𝜇 1 … (x N n 0 ) T (q N n ) T 𝜇 N n ] T , ( 49
)
where N n is the total number of nodes. Assembling expression [START_REF] Argyris | Nonlinear finite element analysis of elastic systems under nonconservative loading-natural formulation. part i. quasistatic problems[END_REF] for all elements and stating that the virtual work principle holds for any function 𝛿U, the global discrete problem simply writes

(U) =  i (U) - e (U) = 0, (50)
with the global internal force vector  i and the global external force vector  e . With eight degrees of freedom per node and N n = (n + 1)N e + 1 nodes, the finite element method has transformed the 8 continuous equilibrium equations [START_REF] Romero | The interpolation of rotations and its application to finite element models of geometrically exact rods[END_REF] in a set of N eq = 8(n + 1)N e + 8 equations forming the nonlinear global problem to solve. [START_REF] Zupan | Finite-element formulation of geometrically exact three-dimensional beam theories based on interpolation of strain measures[END_REF] n = 1 897.3 Cesarek et al. [START_REF] Češarek | Kinematically exact curved and twisted strain-based beam[END_REF] n = 1 906.57 899.69 Ibrahimbegovic et al. [START_REF] Ibrahimbegović | Computational aspects of vector-like parametrization of three-dimensional finite rotations[END_REF] n = 1 897.3 Simo and Vu-Quoc [START_REF] Simo | A three-dimensional finite-strain rod model. part ii: computational aspects[END_REF] n = 0 905.28

User control of a quasi-static problem

In the process of controlling equation (50) with a parameter 𝜆, we can make the distinction between 2 main types of control: load control (force/moment) and displacement/rotation control. In the former case, the control parameter is inserted into the known external force vector  e such that the problem (50) becomes

(U, 𝜆) =  i (U) - e (U, 𝜆) = 0. (51)
In the latter case, 𝜆 is included in the expression of the prescribed displacement/rotation degrees of freedom of vector U (a weight might be associated to 𝜆 when there are several prescribed displacements/rotations). As a degree of freedom is removed of the system for each imposed displacement/rotation, the corresponding equations are also removed of the system: the system would otherwise be overdetermined. These equations are not necessary to get the equilibrium configuration but are not hollow though since they give the reaction force at locked degrees of freedom. Mathematically, we can express control of the j-th degree of freedom with a prescribed displacement u d through

(U * , 𝜆) =  * i (U * ) - * e (U * ) = 0, U * = [ u 1 u 2 … u j-1 u d (𝜆) u j+1 … u N n ] T , ( 52
)
where  * i and  * e are the force vectors from which the j-th component has been removed.

Finally, let us comment the implementation of boundary conditions (clamped or pinned ends) in the presence of quaternions. Boundary conditions are applied classically in a Boolean manner: the equations corresponding to locked degrees of freedom are removed from the system. However, because the 4 components of a quaternion are constrained, the rotational degrees of freedom need to be dealt with care. In particular, it stems from the expression of essential boundary conditions [START_REF] Zupan | On higher order integration of angular velocities using quaternions[END_REF] that only the vector part of the quaternion needs to be prescribed: the fourth component is then determined by the unity constraint. It follows that to block all the rotations at one node, the three equations corresponding to q 1 , q 2 and q 3 at this node are to be removed. To partially block the rotations (two or less rotation directions), the equations corresponding to these quaternion components are classically removed. For instance, to block the rotations around e 2 and e 3 at one node (and supposing for sake of simplicity that the axes of the corresponding section are aligned with global axes), it suffices to remove the equations corresponding to q 2 and q 3 at this node.

Assembly of beams/rigid joints

With the purpose of simulating cables, it is necessary to be able to assemble several beams. Indeed, in automotive industry, cables often have a harness-type geometry. There exist several great references which explain how to model the different types of joints for flexible bodies, for instance [START_REF] Géradin | Flexible Multibody Dynamics: a Finite Element Approach[END_REF] or [START_REF] Shabana | Dynamics of Multibody Systems[END_REF]. However, it is less often expressed in terms of quaternions. We thus explain here how we implemented a rigid joint in our simulation as a basis for modeling other types of joints.

Let us assume that we want to bind the extremities of two beams that we will denote with roman numbers I and II. Each beam, discretized with the finite element method, has one node at its extremity. To create a rigid joint between the ends of the two beams, one must constrain the position of the 2 nodes to remain the same and the rotation between the cross-sections attached at each node to not vary along the deformation. If we denote respectively z I = [(x I 0 ) T (q I ) T 𝜇 I ] T and z II = [(x II 0 ) T (q II ) T 𝜇 II ] T the degrees of freedom of nodes I and II in the current state, the constraint writes for the position:

x I 0 -x II 0 = 0. ( 53
)
The rotation constraint is a bit more tricky to set up. Let us introduce the quaternions at each node in the reference configuration qI 0 and qII 0 . The variation of the rotation between the reference and the actual configuration for the nodes I and II may be represented by the two quaternions Δq I and Δq II which, following the rule of combination of rotations [START_REF] Lazarus | Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method[END_REF], then meet qI = Δq I ⚬ qI 0 , qII = Δq II ⚬ q II 0 .

(54)

As the variation of rotation should be the same for the bound nodes, that being Δq I = Δq II , and using relations (54) and ( 22) the constraint finally writes as a function of the variables of the problem

qI ⚬ q I * 0 -qII ⚬ q II * 0 = 0. ( 55
)
As recalled in the previous section, only three components of the quaternions need to be prescribed. Consequently, only three constraints have to be added to the equations and we choose here to use only the three components corresponding to the vector part of (55). The Lagrangian method is chosen to insert these constraints in the equations. Introducing the vector of constraints  jt for the joint which is then equal to

 jt = [ x I 0 -x II 0 Vec ( q I ⚬ q I * 0 -q II ⚬ qII * 0 ) ] , ( 56
)
and the associated Lagrangian multiplier 𝝁 jt , the constraint is included in the model by adding the quantity 𝝁 jt • 𝛿 jt +  jt • 𝛿𝝁 jt in the discrete form of the virtual work principle (29).

Comments on trial functions interpolation

As the virtual work principle (29) holds for any test functions 𝛿x 0 , 𝛿q and 𝛿𝜇, their interpolation can be carried out without precautions using standard schemes. On the contrary, more careful attention needs to be paid to the interpolation of trial functions [START_REF] Crisfield | Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation[END_REF]. The position x 0 and the Lagrangian multiplier 𝜇 both belong to a linear vector space, so the classical interpolation (37) is adequate. However, it is not the case for the rotations q which belong to a non-Euclidean space and thus are non-additive quantities.

For this reason, many articles focus on the subject [START_REF] Ibrahimbegovic | On the role of frame-invariance in structural mechanics models at finite rotations[END_REF][START_REF] Romero | The interpolation of rotations and its application to finite element models of geometrically exact rods[END_REF][START_REF] Bauchau | Interpolation of finite rotations in flexible multi-body dynamics simulations[END_REF][START_REF] Ghosh | Consistent quaternion interpolation for objective finite element approximation of geometrically exact beam[END_REF]]. Among the existing interpolation, the most promising is the Spherical Linear Interpolation (SLERP), widely used in computer graphics [START_REF] Romero | The interpolation of rotations and its application to finite element models of geometrically exact rods[END_REF]. However, this interpolation is thus far limited to linear elements (n = 0) and the attempt made by Ref. [START_REF] Ghosh | Consistent quaternion interpolation for objective finite element approximation of geometrically exact beam[END_REF] to extend it to superior orders has unfortunately been inconclusive. Recently, the work of [START_REF] Zupan | On higher order integration of angular velocities using quaternions[END_REF] on this matter looks interesting but it is not fully developed yet. Other nonadditive methods have been presented such as interpolating the incre- mental rotation or using the curvature as primary variable [START_REF] Zupan | Finite-element formulation of geometrically exact three-dimensional beam theories based on interpolation of strain measures[END_REF][START_REF] Schulz | Non-linear spatial Timoshenko beam element with curvature interpolation[END_REF] but are not adapted to the quadratic formalism of the ANM.

Departing from these more elaborate interpolations, our use of a classical high order interpolation has been motivated by several arguments. The first reason is that this interpolation does not spoil the polynomial shape of the system (50) and allows to use the ANM. Then, the conclusions of the work of Romero [START_REF] Romero | The interpolation of rotations and its application to finite element models of geometrically exact rods[END_REF] show that no matter what interpolation is chosen, the inaccuracies are always cured by increasing the order of interpolation and the mesh refinement. Thus, only the computational cost is affected by a choice of interpolation less adapted to the rotation space. Finally, the classical additive interpolation benefits from an easy implementation while the interpolation devised on the above quoted article may reveal very tricky.

Asymptotic numerical method

To the best knowledge of the authors, up to now, solving the nonlinear systems (51) or (52) has always been carried out using classical predictor-corrector methods. In this part, the Asymptotic Numerical Method (ANM) introduced by Cochelin [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF] and Damil and Potier-Ferry [START_REF] Damil | A new method to compute perturbed bifurcations: application to the buckling of imperfect elastic structures[END_REF] is presented. It is a high order perturbation technique that allows to compute a large part of a nonlinear branch with only one matrix inversion. The method principle is, starting from a solution point, to seek the branch of solution as an asymptotic expansion of a path parameter. This leads to a semi-analytical solution, whose range of validity is automatically calculated. Applied step-by-step, it allows to describe the full branch of solution. Application of the method to the herein described mechanical problem is presented in the following.

Quadratic recast

As a first step, the nonlinear system of equilibrium equations ( 50) is recast into a quadratic form. The ANM framework allows a systematic solving of this particular form of problems thus providing a powerful solver method for a large class of physical problems [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF]. This recasting operation is set up by introducing additional variables, referred as the slave variables V, to the initial primary variables contained in the vector U. This operation is facilitated by the use of quaternions. As explained earlier, given the expressions of strains (26a) and (26b), the linear material law [START_REF] Cardona | A beam finite element non-linear theory with finite rotations[END_REF] and the expression of f (e) i (48), the system (50) is indeed a polynomial of degree 7. By introducing in each element e the

v (e) = [ v T 1 v T 2 v T 3 vT 4 vT 5 ] T , ( 57
) such that ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ v 1 = Vec( Γ) = Vec(q * ⚬ x ′ 0 ⚬ q -û1 ) = Vec(v 4 ⚬ q -û1 ), v 2 = Vec( K) = Vec(2q * ⚬ q′ -K0 ), v 3 = Vec( N × ( Γ + Γ 0 ) + M × ( K + K0 )), v4 = q * ⚬ x ′ 0 , v5 = q ⚬ N, ( 58
)
with N = ĈN 𝚪, M = ĈM K, and by introducing w (e) = [(u (e) ) T (v (e) ) T ] T , the vector f (e) i is now expressed under the quadratic form

(59)
where the symbols c, l(•) and q(•, •) highlight respectively the constant, linear and quadratic part of the vector. The equations defining the slave variables (58), themselves quadratic, are added to the global system under the form:

(60)

The 17 slave variables have only to be evaluated at Gauss points, which makes 17(n + 1) supplementary variables per element. They are all gathered in the vector V of size 17(n + 1)N e in sequence

V = [ ( 1 v (1) ) T ( 2 v (1) ) T … ( N g v (1) ) T … … ( 1 v (N e ) ) T … ( N g v (N e ) ) T ] T , ( 61
)
where the left exponent represents the Gauss point number and the right exponent the element number. The control parameter 𝜆 defined in paragraph 3.2 is included in the unknown vector so that it does not play a special role in the continuation process. This leads to a new vector of unknowns W defined by

W = [ U T V T 𝜆 ] T , ( 62
)
which allows to write the system (50) as the new quadratic problem  Q of N eq = 25(n + 1)N e + 8 equations:

 Q (W) =  + (W) + (W, W) = 0, (63)
, (•) and (•, •) being respectively constant, linear and bilinear vector valued operators.

The ANM continuation

Let us assume that we know a regular solution point W 0 . In the ANM process, the branch of solution passing through this point is sought in the form of a power series expansion of a path parameter a, truncated at order m:

W(a) = W 0 + aW 1 + a 2 W 2 + • • • + a m W m = W 0 + m ∑ p=1 a p W p . ( 64
)
The path parameter a is chosen as a pseudo arc-length measure: it is the projection of the vector of state variables increment W -W 0 on the normalized tangent vector W 1 :

a = ( W -W 0 ) T W 1 . ( 65
)
When the series ( 64) is inserted into the system (63) and the terms of each power of a are equated to zero, it yields to the m + 1 linear systems whose unknowns are the asymptotic expansion coefficients W p (p = 1, 2, … , m):

 + (W 0 ) + (W 0 , W 0 ) = 0 (order 0),  t W 1 = (W 1 ) + (W 0 , W 1 ) + (W 1 , W 0 ) = 0 (order 1),  t W p + p-1 ∑ q=1 (W q , W p-q ) = 0 (order 2 ≤ p ≤ m). ( 66
)
The tangent operator  t used in (66) is the same for all orders and is equal to

 t = 𝜕 Q 𝜕W (W 0 ) = (•) + 2(W 0 , •). (67)
The system at order 0 does not need to be solved since it stems directly from the definition of W 0 (solution of equation ( 50)). The m other systems which are to be solved are constituted of N eq -1 equations but have N eq unknowns, and are thus under-determined. The additional equations required to close the systems result from inserting (64) in the definition of the path parameter (65) which leads to:

W T 1 W 1 = 1 (order 1), ( 68
)
W T 1 W p = 0 (order 2 ≤ p ≤ m). ( 69
)
The m linear systems can then be solved successively since the p th system depends only on the solutions at previous orders.

Henceforth, an expression of the solution as a truncated series is known. However, this solution is generally admissible only in the vicinity of the starting point W 0 , limited by the finite radius of convergence a rc of the series. As the series ( 64) is truncated at order m, using a rc as the range of validity of the series may lead to erroneous solutions. Instead, a maximal utility range a max is defined such that a tolerance criterion 𝜖 (user-defined) is met: If  m+1 Q designates the coefficient of the power series of order m + 1, using the approximation ‖ Q (W(a))‖ ≈ a m+1 ‖ m+1 Q ‖ (see Ref. [START_REF] Cochelin | Méthode Asymptotique Numérique[END_REF] for the entire proof), the maximum step size is defined by

∀a ∈ [ 0, a max ] ‖ Q (W(a))‖ < 𝜖. ( 70
)
a max = ( 𝜖 ‖ m+1 Q ‖ ) 1 m+1 . ( 71
)
A section of the branch of solution is thus defined as

W(a), a ∈ [ 0 a max ]
. The next section is calculated by using W(a max ) as the new starting point. The full diagram can thereby be calculated iteratively.

Further comments

• The ANM presents the following advantages:

-The branch of solution is known analytically, section by section.

-The tangent matrix  t needs to be calculated only for order p = 1 and can then be reused for superior orders. It considerably outperforms the classical predictor-corrector methods which require an actualization of the Jacobian matrix at each iteration. -The calculation of an analytical expression of the tangent matrix  t is not required, since the quadratic form of the problem gives a direct access to its expression through the definition (67). -The range of validity is calculated automatically and without tuning parameter. Unlike predictor-corrector method, both parameters 𝜖 and m do not alter the convergence but solely the accuracy.

From empirical results, setting 𝜖 = 1 × 10 -7 and 15 ≤ m ≤ 20 is found optimal to get the best trade-off between accuracy and speed and do not need any further tuning. -The ANM benefits from the same robustness as the arc-length method even for sharp turns, due to the use of a pseudo arc-length measure as path parameter. Moreover, for such smooth nonlinearities as in our equilibrium equations, a correction step is not necessary: in practice, it is rare to see a residue greater than 10𝜖 [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF]. If ever necessary, for instance for non-polynomial nonlinearities, a Newton-Raphson type correction step can easily be coupled to the ANM [START_REF] Karkar | A comparative study of the harmonic balance method and the orthogonal collocation method on stiff nonlinear systems[END_REF].

• The maximum step size defined in (71) in practice hardly differs from the radius of convergence a rc . • One of the main advantages of the ANM is its ability to follow a branch even when a bifurcation occurs. This remarkable property allows to describe a full bifurcation diagram. To switch from a branch to another, a special procedure is however required. One possible method is to perturb the initial problem with a constant vector of small norm cP such that the problem (63) becomes

 p (W) =  + (W) + (W, W) + cP = 0. ( 72
)
 p (W) is the new perturbed problem, P a normalized vector of constant random numbers and c is the intensity of the perturbation, assumed small.

The special case of imposed rotations by quaternions in the ANM

When using quaternions, imposing a rotation at one node proves to be a delicate task. Because of the redundancy of quaternions (q and -q represent the same rotation), applying a rotation in a classical way, as presented in paragraph 3.2, does not permit to impose rotations of angle greater than 𝜋. It is thus necessary to come back to the rotational quaternion definition [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF] to control the rotation. One way to perform a rotation of angle 𝜆 around a unit axis n at node i is under the form of a Lagrangian constraint:

 rot = qi - [ cos(𝜆∕2) sin(𝜆∕2) n ] = 0. (73)
However, an additional difficulty arises from using the ANM. This constraint is not in quadratic form and does not fall into the quadratic framework of the ANM. Fortunately, there exists a procedure [START_REF] Karkar | Méthodes numériques pour les systèmes dynamiques non linéaires: application aux instruments de musique auto-oscillants[END_REF][START_REF] Lazarus | Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method[END_REF] to overcome this difficulty and to quadratically recast equation (73). Let us introduce the two new variables

s = sin(𝜆∕2), c = cos(𝜆∕2), (74)
added to the global vector of unknowns W. Differentiating these expressions leads to

ds = c 2 d𝜆, dc = - s 2 d𝜆, (75)
which are quadratic equations. Let us now define the linear and bilinear vector-valued operators  h and  h such that (71) writes

 h (dW) =  h (W, dW), ( 76
) with  h (dW) = [ ds dc ] T and  h (W, dW) = [ c 2 d𝜆 -s 2 d𝜆
] T . Karkar [START_REF] Karkar | Méthodes numériques pour les systèmes dynamiques non linéaires: application aux instruments de musique auto-oscillants[END_REF] noticed that if the vector W is written as the power series expansion (64), differentiating it reads

= ( W 1 + 2aW 2 + • • • + ma m-1 W m ) da. ( 77
)
Inserting the power series expansions (64) and (77) in (76) and equating the terms of each order to 0 in a similar manner to paragraph 4.2, leads to the m linear equations

 h • W 1 =  h (W 1 ) - h (W 0 , W 1 ) = 0 (order 0),  h • W p+1 - p ∑ i=1 p + 1 -i p + 1  h (W i , W p+1-i ) = 0 (order 1 < p ≤ m -1), (78)
with the tangent operator

 h =  h (•) - h (W 0 , •)
. By adjoining the system of equations (78) to the main system of equations (66) and solving it, the values of s and c are known on the branch of solution. The prescribed rotations are then finally applied by adding the constraint (73) in the virtual work principle under its virtual form

𝝁 rot • 𝛿 rot +  rot • 𝛿𝝁 rot = 𝝁 rot • 𝛿q i + ( qi - [ c sn ]) • 𝛿𝝁 rot , (79)
with the vector 𝝁 rot containing the 4 Lagrangian multipliers associated to the constraints.

Validation

The method presented above has been implemented in ManLab [START_REF] Karkar | Manlab -an Interactive Path-following and Bifurcation Analysis Software[END_REF], an open-source Matlab package. In this section, we show the results of our code on several reference numerical experiments for which an analytical solution is known or which have been tested by other authors. These examples serve to demonstrate the validity and the accuracy of our model as well as its abilities. In these tests, the centerline is chosen such that it passes through the centroid of all the sections (no additional out-of-diagonal terms in equation ( 8)). For the comparison, both element order n and element number N e can be tuned.

In figures of deformed shapes, the color map represents the distribution of Von Mises stresses on the surface of the beam.

Cantilever beam under end-moment

We consider an initially straight beam, clamped at s = 0 and sub- jected to an end-moment M L = 3𝜋 EI 3 L around e 3 at s = L. This example, studied in numerous papers [START_REF] Geradin | Kinematics and dynamics of rigid and flexible mechanisms using finite elements and quaternion algebra[END_REF][START_REF] Simo | A three-dimensional finite-strain rod model. part ii: computational aspects[END_REF], is known analytically to give an arc of a circle. In particular, the analytical position and rotation of the right-end section 1. This example demonstrates a good agreement of our method with the analytical solution since it converges toward the exact solution when the number of elements is increased. Additionally, it shows that increasing element order allows to converge more rapidly. These results are confirmed by Fig. 2 which displays the absolute error on vertical position x 2 (L) with respect to both the number of elements and the computation time for n = 0 to n = 3. One also sees on this figure that, for our computer code and on this example, it is worthier to increase the element order than the number of elements to get a small error. This remark is not a general conclusion for our method and should be regarded with precautions since the computation time result is both code dependent and case dependent. However, it is a good indication on our code behavior.

x 1 (L) = EI 3 M L sin ( M L L EI 3 ) = 0, 𝜃 3 (L) = M L L EI 3 = 3𝜋, x 2 (L) = EI 3 M L ( 1 -cos ( M L L EI 3 )) = 2L 3𝜋 , (80)

Cantilever 45 • bend

This example first considered in Ref. [START_REF] Bathe | Large displacement analysis of three-dimensional beam structures[END_REF] has become a classical numerical experiment since it tests all together bending, shear, extensional and torsional deformations. An initially bent cantilever beam of radius R = 100 m and of angle 𝜙 = 𝜋/4 rad (thus of length L = 𝜙R), contained in the (u 1 , u 2 ) plane as pictured Fig. 3, is subjected to an end force F = 600 N normal to this plane. The cross-section is a unit square while the Young's and shear moduli are respectively E = 10 7 Pa and G = E/2. The mesh is constituted of 8 elements.

The experiment shows a good accordance with the numerical results of other authors, see Table 2. In addition, our algorithm proves to have a good convergence rate for increasing order since beyond n = 2 there is no more change on the 4th decimal.

Deep circular arch

In this well-known example, we consider an initially curved beam whose centerline forms an arc of a circle of angle 𝜙 = 215 ⚬ and of radius R = 100 m. The arch right-end is clamped while the left-end is hinged leading to a non-symmetrical structure. The remaining properties are EA = GA 2 = GA 3 = 10 6 N and GJ = EI 2 = EI 3 = 10 8 Nm 2 . A vertical load F 2 is applied at the middle of the beam.

The in-plane buckling and post-buckling of the beam are then studied. This experiment allows to test the ability of our code to detect accurately a buckling load. The results are compared to several results of the literature [START_REF] Ibrahimbegović | Computational aspects of vector-like parametrization of three-dimensional finite rotations[END_REF][START_REF] Zupan | Finite-element formulation of geometrically exact three-dimensional beam theories based on interpolation of strain measures[END_REF][START_REF] Simo | A three-dimensional finite-strain rod model. part ii: computational aspects[END_REF][START_REF] Češarek | Kinematically exact curved and twisted strain-based beam[END_REF] and to the reference solution of DaDeppo and Schmidt [START_REF] Dadeppo | Instability of clamped-hinged circular arches subjected to a point load[END_REF], F 1 2 = 897 N accurate up to 3 significant digits. The results shown on Table 3 totally concurs with other authors numerical findings. One notices that the numerical solution converges faster toward the analytical solution for Zupan and Saje [START_REF] Zupan | Finite-element formulation of geometrically exact three-dimensional beam theories based on interpolation of strain measures[END_REF] and Ibrahimbegovic [START_REF] Ibrahimbegović | Computational aspects of vector-like parametrization of three-dimensional finite rotations[END_REF] though. This difference is mainly due to the different type of interpolation chosen (see paragraph 3.4).

Besides, as reported by Cesarek et al. [START_REF] Češarek | Kinematically exact curved and twisted strain-based beam[END_REF], this example has been traditionally studied only as an in-plane problem. Nevertheless, when considered in 3D, this example is much more complex and the arch experiments out-of-plane buckling. Consequently, it appears as a perfect test for the branch-following feature of the ANM. As illustrated on Fig. 4, the ANM allows to describe the full diagram. To do so the perturbation coefficient c of equation ( 72) is set to 10 -6 . This perturbation is applied just before the bifurcation and then removed when the equilibrium point is "far enough" from the unstable branch. The corresponding deformed shapes are illustrated on Fig. 5.

On Fig. 4, the branch (1) is the in-plane branch of solution. Then, one can notice three "pairs" of pitchfork-type bifurcated branches. Each pair contains two symmetric branches, a u 3 positive and a u 3 negative, due to the symmetry of the problem. The branch (2) is very similar in shape to the one found in Ref. [START_REF] Češarek | Kinematically exact curved and twisted strain-based beam[END_REF] but the critical load differs since we find a value of F 3 2 = 78 N. In the absence of other comparisons, we will consider the exact value remains outstanding. The two other bifurcation points (branches (3) and ( 4)), to our best knowledge never mentioned in any other papers, are found at F 2 2 = 643 N before the in-plane turning point and at F 4 2 = 790 N right after this point.

In-plane hinged right-angle frame

In this example, for which a numerical and an analytical solution can be found in Simo and Vu-Quoc [START_REF] Simo | A three-dimensional finite-strain rod model. part ii: computational aspects[END_REF] and Argyris and Symeonidis [START_REF] Argyris | Nonlinear finite element analysis of elastic systems under nonconservative loading-natural formulation. part i. quasistatic problems[END_REF], the buckling and post-buckling behavior of a frame hinged at both ends is under study. The frame is made of two perpendicular legs linked by a rigid joint as presented in paragraph 3. and7). When running the example for a mesh of 10 quadratic elements per leg, a critical load F crit of 18552 N is obtained in the conservative case which matches with 3 significant digits the load of 18532 N obtained in Ref. [START_REF] Simo | A three-dimensional finite-strain rod model. part ii: computational aspects[END_REF] and the load of 18550 N obtained in Ref. [START_REF] Argyris | Nonlinear finite element analysis of elastic systems under nonconservative loading-natural formulation. part i. quasistatic problems[END_REF]. For the following load case, we get a critical load of 35334 N against 35447 N in Ref. [START_REF] Simo | A three-dimensional finite-strain rod model. part ii: computational aspects[END_REF] and 36000 N in Ref. [START_REF] Argyris | Nonlinear finite element analysis of elastic systems under nonconservative loading-natural formulation. part i. quasistatic problems[END_REF]. As our results concur with the results of other authors, it shows the validity of our method for a complex structure having joints. Once again, let us notice that the ANM follows the equilibrium effortlessly, even for the nonconservative load and despite the presence of a joint.

Besides, although it is not displayed in this paper for brevity, this example exhibits out-of-plane equilibrium branches which could be followed easily with the ANM.

Out-of-plane buckling of a clamped-clamped beam

The classical example illustrated on Fig. 9, is a buckling problem for which the beam is clamped at both ends and a force F applied in the axis of the beam on the right end. The beam measures L = 1 m, has a circular cross-section of radius R = 0.02 m and material properties E = 105 GPa and 𝜈 = 0.3. As the out-of-plane buckling is under study, a very small displacement 𝜖 = 10 -6 in u 3 direction is applied before loading. In addition, for practical purposes the buckling is favored in u 2 direction by applying a vertical force F 2 = 10 -6 F at the center of the beam during loading.

The critical force for this example is known analytically to be

F crit = 4EI𝜋 2
L 2 = 5.2 10 5 N (Euler load). The numerical experiment matches well with this expected value: on the left equilibrium curve of Fig. 8, it corresponds to the first change of slope. The post-buckling is usually studied in-plane and the blue equilibrium curve obtained. However, there exists an out-of-plane solution more stable than the in-plane one. The ability of the ANM to switch from one branch of solution to an other allows to get this solution effortlessly. The bifurcation occurs at F = 7.3 10 5 N which is visible on both equilibrium curves. The corresponding deformed shape right after bifurcation is exposed on Fig. 9. The three following deformed shapes show the evolution of this solution when pursuing the loading. An interesting phenomenon appears: a loop is formed and then rotates on itself, which physically translates as a coupling between torsion and bending. This illustrates the ability of the beam model to take into account complex phenomena.

Extremely twisted elastic rod

This example, more unusual, stems from an article of Lazarus et al. [START_REF] Lazarus | Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method[END_REF] for which an experimental procedure has been devised. For this experiment, two types of elastomeric rods have been fabricated under thoroughly controlled conditions. Thanks to this procedure, all the geometrical and material properties of the rods are accurately known. In particular, the two rods of same length 300 mm have been manufactured such that one is naturally straight and the other has a constant initial (natural) curvature 𝜅 3 = 44.84 m -1 , forming the rolls at the top of Fig. 10. The other parameters of the experiment are the following: a circular cross-section of radius R = 1.55 mm, a Young's modulus E = 1300 kPa and a shear modulus G = 433 kPa. The setup, illustrated on Fig. 10, consists in clamping both ends of the rod in two horizontally aligned chucks separated by a distance of 220 mm. Then, the rotation angle of the right-end is quasi-statically increased and the deformed shape observed by a camera.

Simulation-wise, the preparation of the setup requires preliminary loading steps (see Fig. 10). The curved rod is in a first step unrolled by applying a pure moment at one extremity while the other is clamped: as explained in paragraph 5.1, the analytical solution is known for this case and we can thus apply the predetermined moment such that the rod is straight at the end of the step. As a second step, a displacement in the axis of the rod u 1 = -80 mm is applied to the right-end of both rods. During this step, the gravity load is taken into account such that the beam buckles in the direction -u 2 and takes the same V-shape as in the real experiment. At this point, both rods are in the same configuration and solely the prestress differs between them because of the different initial curvature. The main step then consists in applying a rotation angle 𝜃 1 using the special procedure described in paragraph 4.4.

The deformed shapes obtained for the initially straight rod (Fig. 11) and for the initially curved rod (Fig. 12) match perfectly the experimental observations of [START_REF] Lazarus | Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method[END_REF]. One observes that both deformed shapes are very different. In the first case, a plectoneme (loop) forms at the center of the rod after 𝜃 = 11.3𝜋. In the second case, the rod takes a wavy shape with one twist per wave from 2.2𝜋 to 16𝜋. At the latter angle, a plectoneme also forms but this time on the right end. These qualitative observations illustrate the paramount importance of the initial curvature on the behavior of the beam.

This example also validates our method for rotation control and shows its ability to catch a torsion-bending coupling with ease. Besides, both equilibrium curves exhibited Fig. 13 are very intricate with a lot of turning points and unstable branches. As experimentally exposed by Lazarus et al. [START_REF] Lazarus | Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method[END_REF], the branch indeed becomes unstable after the turning point at 𝜃 = 11.3𝜋 in the first case and is unstable from 0 to 2.2𝜋 in the second case. In spite of that, thanks to the ANM, the continuation of the branch is realized automatically and without efforts for the user.

Finally, let us highlight the fact that this experiment as a whole confronts our numerical solution to a real test case (with nevertheless thoroughly controlled conditions) and shows a very good agreement between the simulation and the reality.

Application on a practical example

The examples of the previous section give a good overview of the robustness and accuracy of the herein presented method. To show the potential of the method in an industrial environment, we henceforth present an application to a real example. In this example pictured on Fig. 14, a cable harness made of a main fixed cable which first divides itself in two parts themselves fixed and then in 3 movable pieces is studied. The mounting operation consists in plugging these 3 cables to 3 connectors of known positions and orientations. The needs for the designers are to know accurately the final shape of the cable, the trajectory during the assembly, the stress inside the cable (to avoid damaging) and the feasibility of the operation.

The first step for the simulation is to build the geometry of the initial configuration. The 3 fixed parts of the cable harness are included in the geometry here only for illustration purposes but could be removed without affecting the results. Thanks to the choice of the FEM (see paragraph 3.3), the building of the complex geometry of the cable (6 pieces) is very easy and does not need any further efforts. For this simulation, the 3 movable cables are supposed to be initially straight. Then, for the end of the 3 cables, the final positions and orientations are set from CAD data. The loading thus consists in displacing and rotating linearly these ends from their initial (𝜆 = 0) to their final configuration (𝜆 = 1).

The shape of the cable during loading as well as the stress inside are displayed on Fig. 15. Of course, these data are known during the whole loading process and the trajectory of the cable can be retrieved from the code. Here, an arbitrary material is chosen but when the right material is defined, the simulation gives access to the right reaction force at the ends of the 3 cables. Finally, let us notice that the initial shape guessed in the CAD environment is not far from the "physical" shape obtained with our code. However, on the left cable of Fig. 14, the abrupt change of curvature at the joint seems unrealistic on the contrary of our results, and could lead to unpredicted problems. In addition, this type of mockup of course does not give access to the physical data previously mentioned, i.e. inside stress, reaction force …

Conclusion

In this paper, a finite element formulation of the geometrically exact beam model using rotational quaternions is presented. Their singularity-free property, their numerical efficiency and their algebraic nature are the main arguments for this choice. In particular, the latter trait leads to a polynomial formulation of the discrete equilibrium equations that we exploit through the use of the asymptotic numerical method. This method is a powerful solver for systems of quadratic nonlinear equations which, among a lot of benefits, is fully automatic and is very robust since it gives a semi-analytical representation of the equilibrium branches. It replaces the classical predictor-corrector methods and the combination of the finite element method, the quaternions and the asymptotic numerical method then forms a remarkable and innovative tool for the simulation of nonlinear beams. Clues for using this method as an industrial tool for the simulation of cables in automotive industry are given, namely the assembly of several pieces of cables and the displacement/rotation control of the algorithm. Application of the method on several numerical and experimental examples exhibits great robustness and accuracy, which makes it very adapted for industry.

Further developments will consist in testing the method on a great variety of cables. In particular, the determination of material parameters which are inputs of the model appears to be a great challenge, especially for complex bundle-type structures. Refined constitutive laws may need to be devised. Finally, adding a stability studying tool to the numerical code will be a necessary step, in order to better comprehend the real phenomena.

N • 𝛿 Γ = (q ⚬ N ⚬ q * ) • 𝛿x ′ 0 + (( 2(q ⚬ N ⚬ q * ) × x ′ 0 ) ⚬ q) • 𝛿q. (C.4)
Besides, still using relations [START_REF] Salamin | Application of Quaternions to Computation with Rotations[END_REF] several times and the expression [START_REF] Crisfield | Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation[END_REF] in the direct differentiation of (27b) 𝛿 K = 2𝛿q * ⚬ q′ + 2q * ⚬ 𝛿q ′ , the right part of the internal work develops as

M • 𝛿 K = 2(q ⚬ M) • 𝛿q ′ -2 M • (q * ⚬ 𝛿q ⚬ q * ⚬ q′ ), = 2(q ⚬ M) • 𝛿q ′ -2( M ⚬ q′ * ⚬ q) • (q * ⚬ 𝛿q), = 2(q ⚬ M) • 𝛿q ′ -2(q ⚬ M ⚬ q′ * ⚬ q) • 𝛿q, = 2(q ⚬ M) • 𝛿q ′ + 2(q ⚬ M ⚬ q * ⚬ q′ ) • 𝛿q. (C.5)
Then, integrating by parts the virtual work of internal forces in which expressions (C.4) and (C.5) have been inserted we get

∫ L 0 ( N • 𝛿 Γ + M • 𝛿 K) ds = [ (q ⚬ N ⚬ q * ) • 𝛿x 0 + 2(q ⚬ M) • 𝛿q ] L 0 -∫ L 0 (( (q ⚬ N ⚬ q * ) ′ • 𝛿x 0 + 2 ( (q ⚬ M ⚬ q * ) ′ + x ′ 0 × (q ⚬ N ⚬ q * ) ) ⚬ q)
• 𝛿q ) ds.

(C.6)

Finally, inserting this expression in the weak form ( 29)-( 30) and using the fundamental lemma of calculus of variations leads to

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Vec [ (q ⚬ N ⚬ q * ) ′ ]
+ n e = 0, ( (q ⚬ M ⚬ q * ) ′ + x′ 0 × (q ⚬ N ⚬ q * ) + me -𝜇 1) ⚬ q = 0, ∀s ∈ [0, L] q • q -1 = 0, (C.7)

along with

⎧ ⎪ ⎨ ⎪ ⎩ Vec [ (q ⚬ N ⚬ q * ) ′ ]
-N e = 0 at s = 0, L, ( (q ⚬ M ⚬ q * ) -Me ) ⚬ q = 0 at s = 0, L.

(C.8)

As |q| ≠ 0, multiplying the second equation of (C.7) on the right by q * leads to the fully equivalent equation:

(q ⚬ M ⚬ q * ) ′ + x′ 0 × (q ⚬ N ⚬ q * ) + me -𝜇 1 = 0. (C.9)

It is noteworthy that the term 𝜇 1 in (C.9) is a scalar quaternion while the other terms of this equation are pure quaternions. It results that this equation is fully equivalent to 𝜇(s) = 0 and Vec [ (q ⚬ M ⚬ q * ) ′ + x′ 0 × (q ⚬ N ⚬ q * ) + me ] = 0, leading to the strong form of equilibrium equations [START_REF] Romero | The interpolation of rotations and its application to finite element models of geometrically exact rods[END_REF]. Similarly, the second equation of (C.8) may be multiplied on the right by q * , leading to the final form of natural boundary conditions [START_REF] Bauchau | Interpolation of finite rotations in flexible multi-body dynamics simulations[END_REF].

Appendix D. Expression of the virtual rotational vector in terms of quaternions

The virtual rotational vector 𝛿𝝑 is defined as the rate of change of the material basis expressed in the global frame [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF] which means

R T 𝛿e i = 𝛿𝝑 × ûi . (D.1)
This expression is fully analogous to the definition of the global curvature. Expressed in terms of quaternions, it then writes q * ⚬ 𝛿ê i ⚬ q = 𝛿 θ × ûi . (D.2)

Using the relation between the material and the global frame as well as the equation ( 25), we demonstrate q * ⚬ 𝛿ê i ⚬ q = q * ⚬ 𝛿q ⚬ ûi + ûi ⚬ 𝛿q * ⚬ q = ûi × (2𝛿q * ⚬ q) = (2q * ⚬ 𝛿q) × ûi .

(D.3)
We thereby deduce the relation [START_REF] Oak-Key | Reduced minimization theory in beam elements[END_REF]: 𝛿 θ = 2q * ⚬ 𝛿q.

Appendix E. Discrete expression of the virtual strains

Using the results obtained in Appendix C, we have

N • 𝛿 Γ + M • 𝛿 K = (q ⚬ N ⚬ q * ) • 𝛿x ′ 0 + 2(q ⚬ M) • 𝛿q ′ + ( x ′ * 0 ⚬ q ⚬ N + q ⚬ N ⚬ q * ⚬ x ′ 0 ⚬ q + 2(q ⚬ M ⚬ q * ⚬ q′ ) ) • 𝛿q. (E.1)
Now developing the third component of [START_REF] Argyris | Nonlinear finite element analysis of elastic systems under nonconservative loading-natural formulation. part i. quasistatic problems[END_REF] with equations (26a), (26b) and ( 25), we show that 2q ⚬

( N × ( Γ + Γ 0 ) + M × ( K + K0 ) ) + 2q ′ ⚬ M = 2q ⚬ ( N × (q * ⚬ x ′ 0 ⚬ q) ) + 2q ⚬ ( M × (2q * ⚬ q′ ) ) + 2q ′ ⚬ M = q ⚬ N ⚬ q * ⚬ x ′ 0 ⚬ q + q ⚬ q * ⚬ x ′ * 0 ⚬ q ⚬ N + 2q ⚬ M ⚬ q * ⚬ q′ + 2q ⚬ q′ * ⚬ q ⚬ M + 2q ′ ⚬ M = q ⚬ N ⚬ q * ⚬ x ′ 0 ⚬ q + x ′ * 0 ⚬ q ⚬ N + 2q ⚬ M ⚬ q * ⚬ q′ -2q ′ ⚬ q * ⚬ q ⚬ M + 2q ′ ⚬ M = q ⚬ N ⚬ q * ⚬ x ′ 0 ⚬ q + x ′ * 0 ⚬ q ⚬ N + 2q ⚬ M ⚬ q * ⚬ q′ . (E.2)
We thus recover the term factor of 𝛿q in equation (E.1). Moreover, stating that q ⚬ N ⚬ q * and 𝛿x ′ 0 are pure quaternions, we can write

(q ⚬ N ⚬ q * ) • 𝛿x ′ 0 = Vec(q ⚬ N ⚬ q * ) • 𝛿x ′ 0 . (E.3)
Finally, with the vector g(N, M) defined in [START_REF] Argyris | Nonlinear finite element analysis of elastic systems under nonconservative loading-natural formulation. part i. quasistatic problems[END_REF], we can rewrite (E.1) as

N • 𝛿 Γ + M • 𝛿 K = g(N, M) • ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝛿x ′ 0 𝛿q ′ 𝛿q 𝛿𝜇 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (E.4)
and with equation ( 41) it results the final relation [START_REF] Argyris | Nonlinear finite element analysis of elastic systems under nonconservative loading-natural formulation. part i. quasistatic problems[END_REF].

Appendix F. Expression of the rotation matrix with quaternions

Using the definition of the quaternion multiplication ( 17 It is then trivial to express the rotation operation (20) under matrix form, that being ŷ = q ⚬ x ⚬ q * = 𝚽 R (q * )𝚽 L (q)x =

[ 1 0 T 0 R ] x. (F.3)
The matrix R in (F.3) is the classical rotation matrix expressed with quaternion parameters through expression R = q 2 0 I 3 + 2q 0 q + qq + qq T . (F.4)

The equivalence between quaternion operation and the classical rotation operator is then exhibited clearly under this form (see also Appendix B of [START_REF] Simo | A three-dimensional finite-strain rod model. part ii: computational aspects[END_REF] for an other proof).

Fig. 1 .

 1 Fig. 1. Kinematics of the geometrically exact beam model.

Fig. 3 .

 3 Fig. 3. Cantilever 45-degree bend subject to a terminal out-of-plane force.

Fig. 4 .

 4 Fig. 4. Bifurcation diagrams for the deep circular arch. Displacement of the force application point in u 2 direction (left) and u 3 direction (right) vs. the force applied. The 4 branches are numbered and their corresponding deformed shapes are found on Fig. 5.

Fig. 5 .

 5 Fig. 5. Deep circular arch example. Deformed shape at several solution points for each branch. Corresponding branch number in Fig. 4 in parenthesis.

Fig. 6 .

 6 Fig. 6. Hinged right-angle frame under conservative vertical load: bifurcation diagram and deformed shapes for several intensities of the force. On the bifurcation diagram, the displacement of the load application point is displayed.

Fig. 7 .

 7 Fig. 7. Hinged right-angle frame under following vertical load: bifurcation diagram and deformed shapes for several intensities of the force. On the bifurcation diagram, the displacement of the load application point is displayed.

Fig. 8 .

 8 Fig. 8. Equilibrium curve of the clamped-clamped beam for in-plane and out-of-plane (3D) buckling. Displacement of the center of the beam (s = L/2) in u 2 direction on the left and maximum displacement in u 3 direction on the right.

Fig. 9 .

 9 Fig. 9. Deformed shape of the clamped-clamped beam for an out-of-plane buckling.

  dW

Fig. 10 .

 10 Fig. 10. Initial configuration of the experiment of Lazarus et al. [21]. At the top, pictures of the bench from Ref. [21]. Below, results of our simulation for the 2 initialization steps of the curved rod: (1) unrolling (2) positioning of the right end.

Fig. 11 .

 11 Fig. 11. Deformed shape of the initially straight rod for the Lazarus et al. experiment. Corresponding point on the equilibrium curve are found on Fig. 13.

  are compared to our numerical solution. The geometrical and material properties chosen for the simulation are: length L = 1 m, circular crosssection of radius R = 4 × 10 -2 m, Young's modulus E = 105 × 10 9 Pa and shear modulus G = 40.3 × 10 9 Pa. The numerical results for various element orders n and number of elements N e are gathered in Table

Fig. 12 .

 12 Fig. 12. Deformed shape of the initially curved rod for the Lazarus et al. [21] experiment. Corresponding point on the equilibrium curve are found on Fig. 13.

Fig. 13 .

 13 Fig. 13. Equilibrium curve of straight (left) and curved (right) Lazarus et al. [21] rods.

3 .

 3 Both legs are identical of length L = 120 m. Their cross-section area is A = 2 m 2 , their quadratic moment I = 6 m 4 , their Young's modulus E = 7.2•10 6 Pa and their Poisson's ratio 𝜈 = 0.3. The loading consists in a vertical force applied at one fifth of the length of the horizontal leg. The example is looked at for both a conservative and a following load (see Figs.6

Fig. 14 .

 14 Fig. 14. Flexible cable designed shape in a CAD environment.

Fig. 15 .

 15 Fig. 15. Numerical simulation results on an industrial example. Evolution of the deformed shape during loading.

  ), one can express left and right quaternion multiplication under matrix form. They respectively write â ⚬ b = 𝚽 L (â) b =

Table 2

 2 Free-end position of the bent cantilever beam for F = 600.

	Formulation	Order	x 1 (L)	x 2 (L)	x 3 (L)
	Present	n = 0	47.3061	15.7009	53.3364
		n = 1	47.1507	15.6845	53.4751
		n = 2	47.1501	15.6847	53.4755
		n = 3	47.1501	15.6847	53.4755
	Bathe and Blorchi [41]	n = 0	47.2	15.9	53.4
	Simo and Vu-Quoc [10]	n = 0	47.23	15.79	53.37
	Zupan et al. [15]	n = 15	47.4159	15.2861	53.4725

Table 3

 3 Critical forces obtained numerically with our code and the code of several authors for in-plane buckling of the deep circular arch.

		N e	10	12	20	40
	Present	n = 1	905.3	901.3	897.8	897.3
	Zupan and Saje

Appendix A. Strain measures in quaternion algebra

Let us prove here the expressions of the measures of strain in quaternion space. For the rotational strain measure, we need to come back to its definition and transpose it to quaternion space. Let us first define the local curvature 𝜿 L which stems from Frenet-Serret formulas through

The expression of this curvature in the global frame, denoted 𝜿 G , is obtained by mapping back this expression from the material to the global frame and then expressing it with the global coordinates which writes

Extending this definition in the quaternion space leads to

Knowing the relation between the global and the local frame in quaternion space êi = q ⚬ ûi ⚬ q * , deriving it and finally using the relation [START_REF] Vince | Quaternions for Computer Graphics[END_REF] we get

Consequently, the curvature in global coordinates writes

As the rotational strain measure must equate 0 in the reference configuration, it is then obvious to define the strain measure as

with κ0 G = 2q * 0 ⚬ q′ 0 the intrinsic curvature of the beam. The expression of the translational strain measure is much more easier to obtain and derives simply from the classic expression (4). It suffices to recall that the rotation in terms of quaternions is defined by [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] and thus replacing R T by it gives directly

Appendix B. Virtual strains in quaternion space

Let us prove here the expressions in quaternion space of the virtual strains (27a) and (27b). For the translational strain, let us on one hand differentiate simply the expression of the strain (26a) from which we get

On the other hand, using the relation [START_REF] Vince | Quaternions for Computer Graphics[END_REF] allows to develop the final expression (27a) in

Finally using the relations ( 22) in (B.2), we recover the expression (B.1) and in that way prove the final expression (27a).

We proceed similarly to get the virtual rotational strain: we first differentiate the expression of the strain (26b) and get

Then, developing the final expression (27b), we have

Finally, using relations ( 23) twice for the second term, we get

and thus recover expression (B.3).

Appendix C. Strong form of equilibrium equations in quaternion space

Let us demonstrate here how to deduce the strong form of equilibrium equations (31) from the weak form ( 29)- [START_REF] Ibrahimbegovic | On the role of frame-invariance in structural mechanics models at finite rotations[END_REF]. Firstly, let us notice that the differential form of equations [START_REF] Crisfield | Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation[END_REF] writes 𝛿q ⚬ q * = -(q ⚬ 𝛿q * ), 𝛿q * ⚬ q = -(q * ⚬ 𝛿q).

Then, from the direct differentiation of the virtual translational strain 𝛿q * ⚬ x ′ 0 ⚬ q + q * ⚬ 𝛿x ′ 0 ⚬ q + q * ⚬ x ′ 0 ⚬ 𝛿q, and by using several times the relations ( 24) and (C.1) we show that the left part of the internal work is equal to

Then factorizing on the right by q and using that, as a pure quaternion, x ′ 0 meets x ′ * 0 = -x ′ 0 , we have