A FINITE ELEMENT/QUATERNION/ASYMPTOTIC NUMERICAL METHOD FOR THE 3D SIMULATION OF FLEXIBLE CABLES

To cite this version : Emmanuel COTTANCEAU, Olivier THOMAS, Philippe VERON, Marc ALOCHET, Renaud DELIGNY A finite element/quaternion/asymptotic numerical method for the 3D simulation of flexible cables In: International Congress of Theoretical and Applied Mechanics (Montréal; 24 ; 2016), Canada, 2016-08 Proceedings of the 24th International Congress of Theoretical and Applied Mechanics ICTAM 2016 2016


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. [1][2][3]. The most widely used theory is the geometrically exact beam model whose founding principles were established by Reissner [4,5] and further generalized by Simo [6]. Various finite element formula-tions (FEM) have been presented to solve these equations numerically. The notable works of [7][8][9][10][11] and more recently [12,13] 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 [14], 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. [12] 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 [15]. 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) [16]. Even if the appearance of arc-length methods [17] 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 [18] and by Cochelin [19]. 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 [20] 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. [21], 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. [22]. 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 where 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 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 [22]. 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: With these notations and hypothesis, the material translational and rotational strain measures are equal to [23].
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 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 for a rotation of angle around the unit axis n. The rotation operator then writes in terms of this parameterization as [22]

Virtual work principle
Using the quantities introduced in the previous section and the Cartesian rotational vector defined in (5) as rotational parameters, the virtual work principle for a 3D shear elastic beam is stated here. Following Reissner's work [4,5] and the generalization of Simo [6], it writes 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 In (8), 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 (7) are given in terms of the kinematic variables by the equations [22]: In (9), T is a tangent operator which can be expressed as a function of only under the form:

Strong form of equilibrium equations
Introducing expressions (9) into the virtual work principle (7), integrating by parts with some rearranging and using the fundamental lemma of calculus of variations lead to the following equilibrium equations: or the essential boundary conditions 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 (10) that the tangent operator T indeed becomes rank-deficient around 2 -angles and leads to a singularity. Procedures exist to circumvent this flaw [9], however, we find it more expedient to use an other parameterization not requiring any special procedure. As demonstrated in Ref. [14], 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 [24]. 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 [25] for more details. The notations of [12] on quaternion algebra are used here. A quaternionâ is defined bŷ where a 0 , a 1 , a 2 , a 3 ∈ ℝ and the imaginary numbers i, j and k are linked by the identities 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 , To distinguish quaternions from the three-dimensional vectors of ℝ 3 , the hat symbol is used: a ∈ ℝ 3 , a ∈ ℝ 4 . Let us introduce the handy notation a = 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 quaternionp 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 (14) and from the relations (15): 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 quaternionsp andq, it is defined aŝ 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 bŷ The rotation is then applied using the following formula: which is quadratic with respect toq 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 quaternionsq 1 andq 2 , the rotational quaternionq 2 ⚬q 1 , evidencing the multiplicative nature of quaternions, needs to be used in (20), that be: Let us notice for practical purposes that the norm of the rotational quaternion defined in equation (19) is 1 (the so-called unit quaternions) obeyinĝ Besides, deriving this relation with respect to s leads to the following useful relationŝ In order to manipulate equations in quaternion space, let us also introduce the following helpful relations, stemming from (17), (18) and the commutativity of inner product: • for arbitrary quaternionsp,q andr: for several values of n and Ne on the left; deformed shape during loading on the right.

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) with the intrinsic translational strain̂0 =û 1 and rotational strain̂0 = 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 (25) (see Appendix B or [12] for the full demonstration) leading tô with the virtual positional (pure) quaternionx 0 and the virtual rotational quaternionq. The resultant of the external moment is conjugated to the virtual rotational vector in the virtual work principle (7). 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 havê = 2q * ⚬q.
Substituting the quaternion expressions (26)- (28) to their vectorial counterpart in the virtual work principle (7) and using the relations (24) to obtain the resultant of the external moment, the new virtual work principle extended in quaternion space writes As noticed in paragraph 2.4, the four components of a rotational quaternionq are not independent but constrained by the condition of unitŷ q ·q − 1 = 0. This constraint is taken into account by the method of Lagrangian multipliers. Following the approach of Zupan et al. [15], the Lagrangian multiplier appended to this constraint and denoted is considered as a variable of the problem. The expression intro-  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: First, introducing the expressions of̂andK (27a) and (27b), then rearranging (29) 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): along with the natural boundary conditions or the essential boundary conditions whereq d is the prescribed rotation under quaternion form (19). It ensues from (31) that (s) = 0 is solution of the continuous problem while the remaining equations are strictly equivalent to the original form (11). 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 (33). The scalar component is indeed automatically known from the fulfillment of the unity constraint (31d). Another equivalent possibility is to replace (33b) byq =q d 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 Table 2 Free-end position of the bent cantilever beam for F = 600.

Formulation
Order  For convenience, let us gather the 8 scalar unknown fields under a single vector field This allows to write the weak form of equilibrium equations (29)- (30) under the compact form with In equation (35) (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.

Discrete form of equations
The finite element method is now applied to the weak form of governing equations (35) 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 (35) 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 [26]. The interpolation then writes 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 where, if I 8 is the 8 × 8 identity matrix, is an interpolation matrix of size 8 × 8(n + 2) and u (e) contains the unknowns nodal values of element e: 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 We can henceforth write under discrete form the continuous equation (35). 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 (41), one with the vector of 12 components 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 In (45), we have introduced the vector c = 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 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 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 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 [24]). 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 [27]. 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 [28].
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: where N n is the total number of nodes. Assembling expression (44) for all elements and stating that the virtual work principle holds for any function U, the global discrete problem simply writes 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 (31) in a set of N eq = 8(n + 1)N e + 8 equations forming the nonlinear global problem to solve.

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 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 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 (33) 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 [22] or [29]. 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 T II ] T the degrees of freedom of nodes I and II in the current state, the constraint writes for the position: The rotation constraint is a bit more tricky to set up. Let us introduce the quaternions at each node in the reference configurationq I 0 and q II 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 (21), then meet 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 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 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 [23]. 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 rotationsq which belong to a non-Euclidean space and thus are non-additive quantities.
For this reason, many articles focus on the subject [30][31][32]13]. Among the existing interpolation, the most promising is the Spherical Linear Interpolation (SLERP), widely used in computer graphics [31]. However, this interpolation is thus far limited to linear elements (n = 0) and the attempt made by Ref. [13] to extend it to superior orders has unfortunately been inconclusive. Recently, the work of [33] 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 [34,35] 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 [31] 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 [19] and Damil and Potier-Ferry [18] 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 [19]. 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 (8) and the expression of f such that 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: 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 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 which allows to write the system (50) as the new quadratic problem  Q of N eq = 25(n + 1)N e + 8 equations: , (•) 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: 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 : 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): The tangent operator  t used in (66) is the same for all orders and is equal to 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: 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. [36] for the entire proof), the maximum step size is defined by A section of the branch of solution is thus defined as . 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 [19]. If ever necessary, for instance for non-polynomial nonlinearities, a Newton-Raphson type correction step can easily be coupled to the ANM [37].
• 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 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 (19) 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: 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 [38,21] to overcome this difficulty and to quadratically recast equation (73). Let us introduce the two new variables added to the global vector of unknowns W. Differentiating these expressions leads to which are quadratic equations. Let us now define the linear and bilinear vector-valued operators  h and  h such that (71) writes Karkar [38] noticed that if the vector W is written as the power series expansion (64), differentiating it reads 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 with the vector rot containing the 4 Lagrangian multipliers associated to the constraints.  [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.

Validation
The method presented above has been implemented in ManLab [39], 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 [40,10], is known analytically to give an arc of a circle. In particular, the analytical position and rotation of the right-end section 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 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.

Cantilever 45 • bend
This example first considered in Ref. [41] 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  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 [9,34,10,42] and to the reference solution of DaDeppo and Schmidt [43], 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 [34] and Ibrahimbegovic [9] though. This difference is mainly due to the different type of interpolation chosen (see paragraph 3.4).
Besides, as reported by Cesarek et al. [42], 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. [42] 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 [10] and Argyris and Symeonidis [44], 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. 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. [10] and the load of 18550 N obtained in Ref. [44]. For the following load case, we get a critical load of 35334 N against 35447 N in Ref. [10] and 36000 N in Ref. [44]. 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. [21] 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 [21]. 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. [21], 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.

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 tô Knowing the relation between the global and the local frame in quaternion spaceê i =q ⚬û i ⚬q * , deriving it and finally using the relation (25) we get Consequently, the curvature in global coordinates writeŝ As the rotational strain measure must equate 0 in the reference configuration, it is then obvious to define the strain measure aŝ witĥ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 (20) and thus replacing R T by it gives directlŷ

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 =q * ⚬x ′ 0 ⚬q +q * ⚬x ′ 0 ⚬q +q * ⚬x ′ 0 ⚬q. (B.1) On the other hand, using the relation (25) 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 havê = 2q * ⚬q ′ + 2q * ⚬q ⚬q ′ * ⚬q.

(B.4)
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)- (30). Firstly, let us notice that the differential form of equations (23) writeŝ q ⚬q * = −(q ⚬q * ), q * ⚬q = −(q * ⚬q).

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 [12] which means This expression is fully analogous to the definition of the global curvature. Expressed in terms of quaternions, it then writeŝ q * ⚬ê i ⚬q =̂×û i .

(E.2)
We thus recover the term factor ofq in equation (E.1). Moreover, stating thatq ⚬N ⚬q * andx ′ 0 are pure quaternions, we can write Finally, with the vector g(N, M) defined in (44), we can rewrite (E.1) aŝ and with equation (41) it results the final relation (44).

Appendix F. Expression of the rotation matrix with quaternions
Using the definition of the quaternion multiplication (17), one can express left and right quaternion multiplication under matrix form. They respectively writê a ⚬̂= L (â)̂= .

(F.2)
It is then trivial to express the rotation operation (20) under matrix form, that beinĝ y =q ⚬x ⚬q * = R (q * ) L (q)x = The matrix R in (F.3) is the classical rotation matrix expressed with quaternion parameters through expression The equivalence between quaternion operation and the classical rotation operator is then exhibited clearly under this form (see also Appendix B of [10] for an other proof).