A total Lagrangian, objective and intrinsically locking-free Petrov-Galerkin SE(3) Cosserat rod finite element formulation

Based on more than three decades of rod finite element theory, this publication unifies all the successful contributions found in literature and eradicates the arising drawbacks like loss of objectivity, locking, path-dependence and redundant coordinates. Specifically, the idea of interpolating the nodal orientations using relative rotation vectors, proposed by Crisfield and Jeleni\'c in 1999, is extended to the interpolation of nodal Euclidean transformation matrices with the aid of relative twists; a strategy that arises from the SE(3)-structure of the Cosserat rod kinematics. Applying a Petrov-Galerkin projection method, we propose a novel rod finite element formulation where the virtual displacements and rotations as well as the translational and angular velocities are interpolated instead of using the consistent variations and time-derivatives of the introduced interpolation formula. Properties such as the intrinsic absence of locking, preservation of objectivity after discretization and parametrization in terms of a minimal number of nodal unknowns are demonstrated by conclusive numerical examples in both statics and dynamics.


Introduction
The theory of shear-deformable (spatial) rod formulations dates back to the pioneer works of Cosserat [1], Timoshenko [2], Reissner [3] and Simo [4]. Thus, depending on the chosen literature, shear-deformable rods are called (special) Cosserat rods [5], Simo-Reissner beams, spatial Timoshenko beams, geometrically exact beams [6], etc. Due to this ambiguity of naming conventions [7], it is best to describe the used rod finite element formulation by the underlying kinematics only [8,9]. Nonetheless, there is a wide agreement in English literature preferring the name special Cosserat rod. Thus, we adopt this notation but for simplicity suppress the prefix "special" in the subsequent treatment. In addition, there is a vast amount of rod finite element formulations found in literature. Introducing all of them concerning their historical development and their tiny differences or improvements is out of the scope of this publication and the interested reader is referred to the exhaustive literature survey of both shear-deformable and shear-rigid rods given by Meier et al. [10]. However, we want to introduce the most important developments found in literature that can be seen as individual milestones leading to the present formulation.
Cardona and Geradin [11] introduced the first total Lagrangian shear-deformable spatial rod finite element formulation where the nodal rotations are parametrized in terms of total rotation vectors. To prevent the well-known singularity of this parametrization, they introduced a strategy of also using the complement rotation vector. By application of a Bubnov-Galerkin projection, the virtual work contributions are discretized in space. Focusing on details of the rotation vector parametrization, the identical rod formulation is found in Ibrahimbegović et al. [12]. In 1999, Crisfield and Jelenić [13] made a groundbreaking discovery in the theory of shear-deformable rod finite element formulations. All previously existing rod finite element formulations violate the objectivity requirement in the discrete approximation. Moreover, they presented a solution to the described problem by interpolating the nodal orientations using relative rotation vectors. The follow-up publication [14] introduces an objective rod finite element formulation using incremental rotation vectors for the description of the nodal rotations and thus working in an updated Lagrangian setting. To simplify the linearization of the internal virtual work functional, the authors applied a Petrov-Galerkin [15] projection by introducing the virtual rotation expressed in the inertial frame. This results in a non-symmetric and configuration-dependent mass matrix. A remedy for this is to use instead the virtual rotation vector expressed in the cross-section-fixed basis together with the corresponding angular velocity (again expressed in the cross-section-fixed basis). As will be discussed in the course of the subsequent treatment, this combination yields a symmetric and possibly constant mass matrix, typically met in rigid body dynamics [16], however, at the cost of a non-symmetric stiffness matrix.
The previously mentioned rod finite element formulations rely on the uncoupled composition approximations of the centerline points as elements of R 3 and the cross-sectional orientations as elements of the special orthogonal group SO(3) denoted by R 3 × SO (3). In contrast, helicoidal approximations [17] and strain-based approaches [18,19] use a coupled interpolation of these two fields. Based on these developments, Sonneville et al. [20] extended the idea of interpolating the nodal orientations using relative rotation vectors [13] to the interpolation of nodal Euclidean transformation matrices using relative twists. That is, an interpolation strategy where positions and orientations are intrinsically coupled. Exclusively working on the Lie group SE (3), the authors [20] applied a Bubnov-Galerkin projection, i.e., the virtual displacements are given by the variation of the nodal values.
Inspired by this interpolation, the present paper introduces a novel rod finite element formulation using only the coupled interpolation strategy of SE (3) while retaining all other favorable properties of the classical formulations. In this way, the present rod finite element combines all the important properties of the formulations in the above literature and tries to eliminate all their drawbacks. In particular, the main contributions of this paper are the following: • We present a novel total Lagrangian (thus path-independent), intrinsically locking-free and objective rod finite element formulation parametrized by the nodal total rotation vectors and centerline points. Using the complement rotation vector in combination with a relative interpolation strategy circumvents possible occurring singularities for element deformations in which the relative rotation angle remains below π. Therefore, a minimal number of nodal unknowns is obtained in contrast to formulations relying on redundant coordinates [6,21,9].
• Application of a Petrov-Galerkin projection results in a very simple inertial virtual work functional that ultimately leads to a symmetric mass matrix that is constant for most applications. In contrast to Sonneville et al. [20], this approach neither requires the involved expressions of the SE(3)-tangent operator and its inverse nor their derivatives.
• Depending on the specific application, a system of ordinary differential equations or the nonlinear generalized force equilibrium is obtained. These can respectively be solved using a standard ODE solver (e.g. Runge-Kutta family, generalized-α) or root finding algorithms (e.g. Newton-Raphson, Riks) instead of requiring highly specialized Lie group versions of them [20]. Hence, the formulation can easily be integrated into existing flexible multibody dynamic software packages.
• Besides the presentation of the linearized internal forces, static and dynamic benchmark examples demonstrate the power of the new formulation. Using a two-node element -that can be integrated by a two-point Gauss-Legendre quadrature rule -second-order spatial convergence is achieved for both centerline and orientation.
The remainder of this paper is organized as follows. In Section 2, the Cosserat rod theory is briefly recapitulated. The section closes with the formulation of the internal, external and inertial virtual work functionals of the sheardeformable spatial rod. Section 3 introduces a novel Petrov-Galerkin rod finite element formulation based on a very efficient two-node SE(3)-interpolation strategy that preserves objectivity of the discretized strain measures. Introducing some bookkeeping enables the precise formulation of the discrete virtual work functionals that result in the tuples of internal, external and gyroscopic forces as well as the symmetric and possibly constant mass matrix. Further, the linearization of internal forces is given. To demonstrate the performance of the presented finite element formulation carefully selected benchmark examples are studied in Section 4, meaning that each of this minimal set of examples demonstrates at least one affirmed property of the proposed formulation. To close the paper, Section 5 presents concluding remarks.
While most parts of the paper can be read with a rudimentary knowledge of matrix Lie groups, the interested reader is referred to Appendix A for a deeper understanding of some manipulations. Therein, a concise overview of matrix Lie groups is given with all equations and references relevant to this work. In particular, the special orthogonal group SO(3) and the special Euclidean group SE(3) are examined as subgroups of GL (3) and GL(4), respectively. Similar introductions are found in literature [22,23,24,20,25]. For didactic reasons, the variation of the rod's strain measures is shown step by step in Appendix B. Since some of the derivatives of the used Lie group formulas cannot be found in literature, all required derivatives of the proposed SE(3)-interpolation strategy are presented in Appendix C. Finally, Appendix D discusses the discrete conservation properties of the proposed rod finite element formulation. That is, the conservation of total energy, linear and angular momentum. For this purpose, another Petrov-Galerkin formulation is presented, where the virtual nodal rotations and the nodal angular velocities are expressed in the inertial frame.

Centerline and cross-section orientation
We introduce the Euclidean 3-space E 3 as an abstract 3-dimensional real inner-product space [5]. A basis for E 3 is a linearly independent set of three vectors e x , e y , e z ∈ E 3 . The basis is said to be right-handed if e x · (e y × e z ) > 0 and orthonormal if their base vectors are mutually orthogonal and have unit length. In this paper, only right-handed orthonormal bases are considered. For a given basis Thus, we carefully distinguish R 3 from the Euclidean 3-space E 3 . The rotation between two bases K 0 and K 1 is captured by the proper orthogonal transformation matrix A K0K1 ∈ SO (3), which relates the coordinate representations K0 a and K1 a in accordance with K0 a = A K0K1 K1 a. Let ξ ∈ J = [0, 1] ⊂ R denote the centerline parameter and (η, ζ) ∈ A(ξ) ⊂ R 2 the cross-section parameters of the cross-section area at ξ. Considering the rod as a three-dimensional continuum, a point Q of the rod can be addressed by Herein, I r OP are the Cartesian components with respect to the inertial basis I of the time-dependent centerline curve r OP , where the subscripts O and P refer to the origin O and the centerline points, see Figure 1. To not overload the notation in (1), we abstained from explicitly denoting the dependence of the points P and Q on ξ and ξ, respectively. Moreover, K r P Q denote the cross-section coordinates, which are transformed to the I-basis using the transformation matrices A IK . Note that the base vectors of the cross-section-fixed K-bases are functions of the centerline parameter ξ and time t, i.e., e K i = e K i (ξ, t), i ∈ {x, y, z}.

Variations, velocities and curvature
Let(•) and (•) ,ξ respectively denote the derivative with respect to time t and centerline parameter ξ. The variation is denoted by δ(•). The virtual displacement I δr P and the centerline velocity I v P are given by the variation as well as the time derivative of the centerline The angular velocity of the K-basis relative to the inertial I-basis, in components w.r.t. the K-basis, is defined by where j SO(3) : R 3 → so(3) = {B ∈ R 3×3 |B T = −B} is the linear and bijective map such that ωr = j SO (3) (ω)r = ω × r for all ω, r ∈ R 3 , see (73) from Appendix A. Analogously, we define the scaled curvature as and the virtual rotation as

Reference arc length
For the reference centerline curve I r 0 OP , the length of the rod's tangent vector is J = I r 0 OP,ξ . Thus, for a given centerline parameter ξ, the reference arc length s is defined by Following Harsch and Eugster [26], the derivative with respect to the reference arc length s of a function f : J ×R → R 3 can be computed by

Objective strain measures
The objective strain measures of a Cosserat rod, see Antman [5,Section 8.2], are given by which can be gathered in the six-dimensional tuple ε = ( K γ, K κ IK ) ∈ R 6 . Therein, the dilatation and shear is captured by K γ, while K κ IK measures torsion and bending, see Antman[5, Section 8.6].

Internal virtual work
As in Equation (2.10) of Simo and Vu-Quoc [4], we introduce the diagonal elasticity matrices C γ = diag(k e , k s , k s ) and C κ = diag(k t , k by , k bz ) with constant coefficients. In the following, the simple quadratic strain energy density is used, where the superscript 0 refers to the evaluation in the rod's reference configuration.

External virtual work
Assume the line distributed external forces I b : J × R → R 3 and moments K c : J × R → R 3 to be given as densities with respect to the reference arc length. Moreover, for i ∈ {0, 1}, point forces I b i : R → R 3 and point moments K c i : R → R 3 can be applied to the rod's boundaries at ξ 0 = 0 and ξ 1 = 1. By that, the corresponding external virtual work functional is given by

Inertial virtual work
Let ρ 0 : J → R denote the rod's mass density per unit reference volume and dA the cross-section surface element. It is convenient to define the following abbreviations Further, using the mass differential dm = ρ 0 JdAdξ and expressing the variation and second time derivative of r OQ in terms of the rod's kinematics (1), the inertial virtual work functional of the Cosserat rod can be written as (cf. Equation (9.54) of Eugster and Harsch [8] for a coordinate-free version) where we have introduced the two quantities Note that the quantity S ρ0 vanishes if the centerline points I r OP (ξ, t) coincide with the cross-sections' center of mass. For this case M and g get independent of the rod's configuration.

Rod kinematics in homogenous coordinates
The frame I = {O, e I x , e I y , e I z } is the set collecting the origin O together with the inertial base vectors e I i , i ∈ {x, y, z}. A generic and possibly non-inertial frame K = {P, e K x , e K y , e K z } is given by a point P together with a right-handed orthonormal basis K. The point Q relative to the K-frame is addressed by the triple K r P Q ∈ R 3 , which are the components of the vector r P Q ∈ E 3 between P and Q with respect to the K-basis. Rigid body motions between two frames K 0 and K 1 are captured by the Euclidean transformation matrix which relates the triples K1 r P1Q ∈ R 3 and K0 r P0Q ∈ R 3 in accordance with In fact, the Euclidean transformation matrix H K0K1 is an element of the special Euclidean group SE(3) considered here as a Lie subgroup of the general linear group GL (4). A brief introduction to the Lie group setting and particularly an overview of the herein required mappings is given in Appendix A. Direct computation readily verifies that the inverse of H K0K1 is Consequently, using (18), the motion of the rod (1) can be written in homogenous coordinates as The group structure of SE(3) then allows to decompose the Euclidean transformation matrix H IK into

SE(3)-interpolation
In order to introduce the idea of the SE(3)-interpolation strategy, we discretize the rod by a single two-node finite element. In contrast to Sonneville et al. [20], we choose a parametrization of H IK0 and H IK1 in terms of time-dependent generalized coordinates q(t) = I r OP0 (t), ψ 0 (t), I r OP1 (t), ψ 1 (t) ∈ R 12 given by the nodal positions I r OP0 (t) ∈ R 3 , I r OP1 (t) ∈ R 3 and the nodal rotation vectors ψ 0 (t) ∈ R 3 , ψ 1 (t) ∈ R 3 . Using the exponential map of SO (3), defined in (77), the parametrization is With the SE(3)-logarithm map from (85), we can compute the relative twist θ K0K1 corresponding to the relative Euclidean transformation H K0K1 = H −1 IK0 H IK1 as Herein ψ 01 corresponds to the relative rotation vector parameterizing the transformation between the K 0 -and K 1basis according to A K0K1 = Exp SO(3) (ψ 01 ). Due to the singularity of Log SO (3) for ψ 01 = π, see Appendix A.2, this interpolation strategy is restricted to applications in which the relative rotation angle satisfies ω = ψ 01 < π. A discretization with a higher number of elements always cures this problem. By linearly scaling the relative twist θ K0K1 and using the SE(3)-exponential map (85), a relative Euclidean transformation from the K 0 -frame to the K-frame can be constructed as originally proposed in Equation (55) of Sonneville et al. [20]. In order to see the interpolation for the position and orientation, we explicitly compute the exponential map in (25) resulting in Consequently, the rod's orientation is discretized by which corresponds to Equation (4.7) of Crisfield and Jelenić [13]. Since A IK (0, q) = A IK0 and A IK (1, q) = A IK1 , the discretization (27) is a highly nonlinear interpolation of the nodal transformation matrices. The rod's centerline is discretized by which also interpolates the nodal positions I r OP (0, q) = I r OP0 and I r OP (1, q) = I r OP1 in a highly nonlinear manner.

Objectivity, discretized strain measures and absence of locking of SE(3)-interpolation
Using the just introduced Euclidean transformation matrices, we can recognize the strain measures of the Cosserat rod theory (8) in the following equation With the proposed interpolation (25), the strain measures in depend only on the relative Euclidean transformation H K0K and its reference arc length derivative. Objectivity is the requirement that a quantity remains unaltered under a change of observer, i.e., a change of the I-frame. Whether we use an I-frame or an I + -frame, which are related by the time-dependent Euclidean transformation H I + I , the discretized strain measures in H −1 I + K H I + K,s correspond to the strain measures obtained by (30). This proofs objectivity of the interpolation (25).
Using a little bit of Lie algebra introduced in the Appendix A, the discretized strain measures can be drastically simplified. First, with the aid of (30), (7) and (82), they can be extracted from (29), by Inserting the relative interpolation (24) together with the definitions (85) and suppressing for a while the subscripts indicating SE (3), the discretized strain measures are where we have used the left-trivialized differential (64) 1 . Due to the linearity of j, we have d/dξ j(ξθ K0K1 ) = j(θ K0K1 ).
Since j(ξθ K0K1 ) and j(θ K0K1 ) commute, according to the comment below (68), the discretized strain measures simplify further to This means that the discretized strains are constant. As discussed in Balobanov and Niiranen [27], shear and membrane locking can occur if Kirchhoff and inextensibility constraints follow in the limit case of a parameter tending to zero. This appears for instance for the strain energy density (12), if the stiffnesses are computed in the sense of Saint-Venant by using the material's Young's and shear moduli E and G, respectively, together with the cross-section geometry. For the particular choice of a quadratic cross-section (width w, area A = w 2 , second moment of area I = w 4 /12), the stiffnesses would be given as k e = EA, k s = GA, k by = k bz = EI and k t = 2GI. Dividing the strain energy density (12) by w 4 , in the limit of w → 0, the scaled axial k e /w 4 and shear stiffnesses k s /w 4 tend to infinity. The strain energy remains only bounded if the dilatation remains one and the shear deformations zero, i.e., if K γ = (1, 0, 0); inextensibility K γ − 1 = 0 readily follows. Formulations that are prone to locking cannot fulfill these constraints exactly over the entire element and introduce parasitic dilatation and shear strains, see Figure 2 (e) for the two-node R 3 × SO(3)-interpolation [13] in a configuration where the nodal coordinates are given by the pure bending situation of a quarter circle. In contrast, according to (33), the proposed SE(3)-interpolation can exactly represent constant strain measures for all the individual contributions, can thus satisfy K γ = (1, 0, 0), and is hence intrinsically relieved from both shear as well as membrane locking. Obviously, the SE(3)-interpolation can represent the quarter circle exactly, see Figure 2

Kinematics of element nodes
The concepts of the previous paragraphs can be synthesized in a rod finite element with a piecewise two-node SE(3)interpolation that results in constant strains within an element. Let the rod be discretized by n el two-node elements. The nodal Euclidean transformation matrices H IKi , i = 0, . . . , n el are parameterized in terms of the time-dependent nodal generalized coordinates q i (t) = I r OPi (t), ψ i (t) ∈ R 6 given by the nodal centerline points I r OPi (t) ∈ R 3 and the nodal total rotation vectors ψ i (t) ∈ R 3 . Using the exponential map of SO(3), see comment before (22), the nodal parametrization is The rod's parameter space J is divided into n el linearly spaced element intervals J e = [ξ e , ξ e+1 ) via J = n el −1 e=0 J e . This means a totality of n el + 1 nodes and n q = 6(n el + 1) unknowns. The nodal quantities can be assembled in the global tuple of generalized coordinates q(t) = q 0 (t), q 1 (t), . . . , q N −1 (t) ∈ R nq . Introducing an appropriate Boolean connectivity matrix C e ∈ R 12×nq , the element generalized coordinates q e (t) = I r OPe (t), ψ e (t), I r OPe+1 (t), ψ e+1 (t) ∈ R 12 can be extracted from the global generalized coordinates q via q e = C e q. By introducing an additional set of Boolean connectivity matrices C r,0 , C r,1 , C ψ,0 , C ψ,1 ∈ R 3×12 , the centerline points I r OPe , I r OPe+1 and total rotation vectors ψ e , ψ e+1 can be extracted from the element generalized coordinates q e via I r OPe = C r,0 q e , I r OPe+1 = C r,1 q e , ψ e = C ψ,0 q e and ψ e+1 = C ψ,1 q e . Note, the Boolean connectivity matrices are only used for the mathematical description of these extraction procedures. During a numerical implementation a different and much more efficient approach is used.

Element-wise SE(3)-interpolation
For a given element interval J e = [ξ e , ξ e+1 ) the corresponding linear Lagrange basis functions N e 0 , N e 1 and their first derivatives N e 0,ξ , N e 1,ξ are defined by Further, we introduce the characteristic function χ J e : J → {0, 1}, being one for ξ ∈ J e = [ξ e , ξ e+1 ) and zero elsewhere. By that, we can extend the global two-node SE(3)-interpolation (25) to a piecewise two-node SE(3)-interpolation given by For symmetry reasons, cf. Crisfield and Jelenić [13], it might be advantageous to use a modified SE(3)-interpolation given by Since this strategy requires twice as many SE (3)-exponential evaluations and three times as many SE (3)-logarithm evaluations, the interpolation (36) is preferred due to superior efficiency and simplicity.
Using similar arguments that resulted in (33), the interpolation (36) leads to piecewise constant strains Using the same reasoning following (33) and since the piecewise two-node SE(3)-interpolation can exactly represent constant strains within each element, neither membrane nor shear locking will appear with this discretization. As we need no further numerical strategies to avoid locking as for instance re-interpolation of strain measures [28,29] or mixed formulations [30,31,32], we call the finite element formulation intrinsically locking-free.

Petrov-Galerkin projection
Following a Petrov-Galerkin projection [15], we introduce the nodal virtual displacements δs i (t) = I δr Pi (t), Ki δφ IKi (t) ∈ R 6 given by the nodal centerline variation I δr Pi (t) ∈ R 3 and the nodal virtual rotation Ki δφ IKi (t) ∈ R 3 . They are related to the variations of the nodal generalized coordinates by I δr Pi = δ( I r OPi ) and Ki δφ IKi = T SO(3) (ψ i )δψ i , see Equation (49) of Cardona and Geradin [11]. Led by the observation that M in (16) is symmetric, it is advisable to introduce the nodal minimal velocities u i (t) = I v Pi (t), Ki ω IKi (t) ∈ R 6 given by the nodal centerline velocity I v Pi (t) ∈ R 3 and the nodal angular velocity Ki ω IKi (t) ∈ R 3 . Similar to the virtual displacements, the nodal minimal velocities are related to the time derivative of the nodal generalized coordinates by the nodal kinematic differential equationq Consequently, during a subsequent Galerkin projection, the symmetry of M is preserved and results in a symmetric mass matrix, see (44). Since the inverse tangent map T −1 SO(3) used in (39) exhibits singularities for ψ = k2π with k = 0, 1, 2, . . . , we apply the following strategy to avoid them. For k = 0, we use the first order approximation For k > 0, the concept of complement rotation vectors [11,12] is applied. Due to the Petrov-Galerkin projection, it is sufficient to introduce a nodal update that is performed after each successful time step (in statics after each successful Newton increment). This update, which corresponds to a change of coordinates for the orientation parametrization, is given by , there is no difference whether the nodal transformation matrix A IK is described by the rotation vector ψ or by its complement ψ C = 1 − 2π/ ψ ψ. Thus, for reasonable time steps (Newton increments) a both minimal and singularity free parametrization of SO(3) is obtained without changing the rod's virtual work formulation since all relevant nodal quantities are interpolated relatively 1 .
In analogy to the generalized coordinates, the nodal virtual displacements and minimal velocities are assembled in the global tuple of virtual displacements δs(t) = δs 0 (t), δs 1 (t), . . . , δs N −1 (t) ∈ R nq and global tuple of minimal velocities u(t) = u 0 (t), u 1 (t), . . . , u N −1 (t) ∈ R nq . Using again the boolean connectivity matrix C e , we can extract the element virtual displacements δs e (t) = I δr Pe (t), Ke δφ IKe (t), I δr Pe+1 (t), Ke+1 δφ IKe+1 (t) ∈ R 12 and the element minimal velocities u e (t) = I v Pe (t), Ke ω IKe (t), I v Pe+1 (t), Ke+1 ω IKe+1 (t) ∈ R 12 from the global quantities in agreement with δs e = C e δs and u e = C e u. Further, the centerline variations I δr Pe , I δr Pe+1 and centerline velocities I v Pe , I v Pe+1 can be extracted from the tuple of virtual displacements δs e and minimal velocities u e using I δr Pe = C r,0 δs e , I δr Pe+1 = C r,1 δs e and I v Pe = C r,0 u e , I v Pe+1 = C r,1 u e , respectively. Identical extraction operations can be defined for the nodal virtual rotations and angular velocities via C ψ,0 , C ψ, 1 .
In the sense of a Petrov-Galerkin projection [15], a different interpolation for the virtual displacements is chosen. We chose a simple interpolation strategy based on the previously introduced linear Lagrange basis function. As outlined in the previous paragraph, in order to obtain a symmetric and possibly constant mass matrix, the angular velocities are interpolated similarly. Formalizing that, we introduce the interpolations Choices like K δr P , i.e., centerline variations in the cross-section K-basis, introduce cumbersome changes of bases when for instance the generalized contact force directions should be computed, cf. Bosten et al. [34]. This is quite inconvenient especially when dealing with complex flexible multibody system connected by perfect bilateral constraints as treated in Géradin and Cardona [35].

Discrete virtual work formulations
Now we have all ingredients for discretizing the previously introduced continuous virtual work formulations in space. Inserting (41) into the internal virtual work (11), their discretization is obtained where we have introduced the internal forces f int together with their element contributions f int e . The discrete ansatzfunction of the interpolation of (36) is used in the computation of the appearing contributions. For the sake of readability, above and subsequently, we partly suppress the function arguments, which should be clear from the context. Similarly, the discretization of the external virtual work (13) takes the form where we have introduced the external forces f ext and their element contributions f ext e . Finally, the discretization of the inertial virtual work contributions (15) is given by where where we have introduced the symmetric mass matrix M and the gyroscopic forces f gyr together with their element contributions M e and f gyr e . Appendix D discusses the discrete conservation properties of the proposed rod finite element formulation. That is, the conservation of total energy, linear and angular momentum. Moreover, another SE(3) finite element formulation is presented, for which the nodal virtual rotations and angular velocities expressed in their cross-section-fixed bases are replaced by the nodal virtual rotations I δφ IKi (t) ∈ R 3 and nodal angular velocities I ω IKi (t) ∈ R 3 expressed in the inertial I-basis.
The linearization of the internal forces f int with respect to the generalized coordinates q is given by which requires the derivative of the rotation matrix A IK and the derivatives of the internal forces and moments Both are depending on the derivatives of the strain measures K γ and K κ IK introduced in (38). These derivatives together with the derivative of the interpolation formula (36), i.e., the position vector I r OP and the transformation matrix A IK are given in Appendix C. As already noted, the quantity S ρ0 vanishes if the centerline points I r OP coincide with the cross-sections' center of mass. For this case the gyroscopic forces f gyr and the term Mu get independent of the generalized coordinates q which is why we do not consider their linearizations here. The linearization of the external forces is omitted as well since a possible q dependence is defined by the specific form of the external forces and moments only. Arising element integrals of the form J e f (ξ)dξ encountered in the discretized internal, external and gyroscopic forces, as well as in the mass matrix, are subsequently computed using a two-point Gauss-Legendre quadrature rule.

Equations of motion and static equilibrium
The principle of virtual work states that the totality of virtual work contributions has to vanish for arbitrary virtual displacements at all time instants t, see Chapter 8 of dell'Isola and Steigmann [36], i.e., Thus, the equations of motionu have to be fulfilled for each instant of time t. Further, the nodal minimal velocities u i are related to the time derivatives of the nodal generalized coordinatesq i via the nodal kinematic differential equation (39), which can be assembled to a global kinematic differential equation of the formq = B(q)u. Depending on the specific application, either the system of ordinary differential equationṡ q = B(q)u , or the nonlinear generalized force equilibrium is obtained. The system of ordinary differential equations can be solved using standard higher-order ODE solvers (e.g. family of explicit [37] and implicit [38,39] Runge-Kutta methods or structure-preserving algorithms [40,41]). This is a remarkable result since existing SE(3) rod formulations require highly customized Lie group solvers [42,20]. In order to apply well-established methods from structural dynamics like the Newmark-β [43] or the generalized-α method [44], a slightly modified update of the generalized coordinates has to be applied, see Equation (37a) and (37b) of Arnold et al. [25]. Alternatively, a generalized-α formulation for first-order differential equations [45] can be used without any modifications. The nonlinear generalized force equilibrium (51) is solved by any root finding algorithm (e.g. Newton-Raphson, Riks). Note that a system of linear equations with a non-symmetric matrix must be solved in each iteration.

Numerical experiments
This section demonstrates the power of the developed SE(3) finite element formulation by a variety of selected benchmark examples. The cantilever experiment slightly extends the investigations of Meier et al. [28] and gives numerical evidence for the absence of locking and further demonstrates a second-order spatial convergence of the proposed formulation. The helix experiment, which was used in Harsch et al. [9] to study locking for a director beam formulation, is used as a second example to show that the SE(3)-formulation is intrinsically locking-free. Objectivity after discretization is demonstrated in the example superimposed rotation of deformed cantilever. Large inhomogeneous deformations are studied in the example rod bent to a helical form [46]. Finally, dynamics is examined by the flexible heavy top example.
We used a Newton-Raphson method to solve all static experiments with an absolute tolerance atol in terms of the max error of the total generalized forces. Prescribed boundary conditions were incorporated into the principle of virtual work using perfect bilateral constraints [35]. The dynamic problem was solved using a standard fourth order Runge-Kutta method as well as the generalized-α method for first order differential equations [45]. For simplicity, boundary conditions were enforced on acceleration level only. Unless stated otherwise, arising element integrals were evaluated using a two-point Gauss-Legendre quadrature rule.

Cantilever experiment
We consider an initially straight cantilever rod of length L = 10 3 with a quadratic cross-section of width w subjected to a tip moment K c 1 = (0, 0, 0.5πk b /L) and an out-of-plane tip load I b 1 = A IK (0, 0, 0.5πk b /L 2 ). In order to show the absence of locking, different slenderness ratios ρ = L/w ∈ {10 1 , 10 2 , 10 3 , 10 4 } are considered, i.e., widths w ∈ {10 2 , 10 1 , 1, 10 −1 }. ranging from w = 10 −1 to w = 10 2 . Further, the elastic constants are given in terms of the Young's and shear moduli E = 1 and G = 0.5. That is, axial stiffness k e = EA, shear stiffness k s = GA, bending stiffnesses k b = k by = k bz = EI and torsional stiffness k t = 2GI, together with A = w 2 and I = w 4 /12. As there is no analytical solution for this load case, the numerical solutions obtained for a single SE(3)-element and for 64 SE(3)-elements are compared to a reference solution found by a finite element implementation similar to Eugster et al. [21] discretized by 256 linear elements. To reduce shear locking in the latter formulation, reduced integration was applied. The strain measures of the reference solution are plotted in Figure 4 (a) and (b). For the centerline, we define the error  where I r * OP denotes the centerline points of the reference solution and I r OP the same quantity of the rod in comparison. In order to investigate the error in orientations, we use the measure [24] Again, A * IK denotes the orthogonal transformation matrix of the reference solution and A IK the same quantity of the rod in comparison. Since the relative rotation vector ∆ψ should be tending to zero during spatial convergence, this is a well-suited error measure [47]. The final loads were applied for all slenderness ratios within 20 increments and absolute tolerances as documented in Table 1. It can be observed that the spatial convergence behavior of the rod is unaffected by the chosen slenderness ratio, although no reduced integration was applied, which numerically proofs the absence of locking. For the different slenderness ratios, the deformed configurations are visualized in Figure 3.
Due to the applied out-of-plane force, the resulting deformation is inhomogenous and cannot be described adequately by a small number of SE(3)-elements with constant strain measures, see Figure 4 (a) and (b). This makes the example suitable for the investigation of the spatial convergence of the proposed formulation. For that, we used a moderate slenderness ratio ρ = 10 3 , an absolute tolerance of atol = 10 −10 and 20 static increments. The numerical errors were computed with respect to the same reference solution as above. Figure 4 (c) visualizes the convergence behavior, which shows a second-order rate of spatial convergence in both centerline and rotation fields. This is in line with the observations made by Sonneville et al. [20] within a different experiment.

Helix experiment
In this example, an initially straight rod, is loaded such that it is deformed to a perfect helix [9] with n = 2 coils of radius R 0 = 10 and height h = 50, see Figure 5. Introducing the abbreviation c = h/(R 0 2πn) ≥ 1, the centerline of such a deformed rod is described by Depending on the slenderness ratio ρ, the rod has a circular cross-section of diameter d = L/ρ, radius r = d/2, area A = πr 2 and second moment of area I = π 4 r 4 . Using the Young's and shear moduli E = 1 and G = 0.5, respectively, the elastic constants are given by k e = EA, k s = GA, k b = k by = k bz = EI and k t = 2GI. As shown in Harsch et al. [9], the force boundary conditions for this specific example are found by a so called inverse procedure, see Section  slenderness tolerance force increments 10 1 10 −8 70 10 2 10 −9 100 10 3 10 −10 200 10 4 10 −11 500 5.2 of Ogden [48]. Evaluating the derivative of (54), the tangent vector of the curve is The rod should not be elongated during the deformation, i.e., J(ξ) = I r * OP,ξ (ξ) = L, from which the rod's total length L = √ 1 + c 2 R 0 2πn follows. Further, the Serret-Frenet frame of the deformed curve is given by Using (1 + c 2 ) − 1 2 = 2πnR 0 /L = α ,ξ R 0 /L together with the derivatives of (56), the rod's strain measures (8) compute as For the given strain energy density (12), evaluating the constitutive equations (10) results in the force boundary conditions that induce a pure torsion and bending deformation. Inserting (58) into the differential equation of the nonlinear Cosserat rod [8] results in the requirement k t = k b . This condition is satisfied for the given problem setup. For a numerical simulation, the rod has to be clamped at I r OP | ξ=0 = I r * OP | ξ=0 with an orientation given by A IK | ξ=0 = ( I e K x , I e K y , I e K z )| ξ=0 . Depending on the used slenderness ratio, a different number of force increments and error tolerances were chosen, see Table 2. As in the cantilever experiment, no locking is observed. In fact, the SE(3)-interpolation can exactly represent a helix, which is a configuration with constant strains (57). As mentioned in Section 3.2, a single SE(3)element is restricted to a relative rotation vector of absolute value bounded by π. Thus four elements can exactly represent the helix (54) with two coils. Nevertheless, during the Newton iterations locally the relative rotation vector might exceed an absolute value of π. Hence, we use five elements for the numerical investigation.

Superimposed rotation of deformed cantilever
The objectivity of the SE(3)-formulation can be demonstrated numerically by using the cantilever rod from Section 4.1 with slenderness ratio ρ = 10 2 , discretized by a single element. Both tip force and moment were successively applied during 50 linearly spaced increments. During the subsequent 450 increments the total rotation vector of the clamped node was linearly increased up to a value of ψ 0 = (20π, 0, 0), i.e., the deformed rod performed 10 full rotations around the e I x -axis, see Figure 6 (a). Note, the very same problem can be solved using less than 100 static increments. However, to get a better resolution of the derived quantities, we used a totality of 500 static increments. It is well known that for non-objective finite element formulations, the total energy increases during a superimposed rigid body rotation [49]. This is not the case for the present formulation, see Figure 6 (b). After the final values of the tip force and moment are reached, the potential energy remains unaltered throughout the superimposed rotation of the clamped node. The tip displacement is shown in Figure 7 (a). Moreover, Figure 7 (b) visualizes that whenever the absolute value of the second nodal rotation vector would exceed π, its complement value is chosen. Thus, the trajectory of the individual components of the nodal rotation vector exhibit discontinuities.

Rod bent to a helical form
The next experiment investigates the capabilities of the presented formulation in describing large inhomogenous deformation. A well-suited example was introduced by Ibrahimbegovic [46] in which an initially straight cantilever rod of length L = 10, axial stiffness k e = 10 4 , shear stiffness k s = 10 4 , bending stiffness k b = k by = k bz = 10 2 and torsional stiffness k t = 10 2 is subjected to an incremental application of a tip moment K c 1 = A T IK (0, 0, 20πk b /L) and an out-of-plane tip load I b 1 = (0, 0, 50). By that, the rod is bend up to 10 full circles, while the e I z -component of the tip displacement "oscillates" with decreasing magnitude, see Figure 8. This observation is in perfect agreement with the results found in literature and is further confirmed by comparing the tip displacements of the e I z -component, shown in Figure 9 (a), with the results reported by Ibrahimbegovic [46]. Instead of the 100 or 200 linear finite elements used by Ibrahimbegovic [46] and Mäkinen [50], it was sufficient to discretize the rod with 30 elements of the present formulation.
Additionally, Figure 9 (b) and (c) visualize the strain measures of the final configuration. Not surprisingly, the dilatation γ 1 , both shear strains γ 2 , γ 3 and both curvatures κ 2 , κ 3 are step functions, since the used two-node element interpolation strategy (36) yields element-wise constant strain measures.

Flexible heavy top
Inspired by the investigation of Mäkinen [50], the dynamics of an elastic heavy top is studied here. On the one hand, this example demonstrates the capability of the presented formulation to be solved using standard ODE solvers. On the other hand, in the limit case of an infinite stiff rod, the rod shows the well-known behavior of a heavy top. The motion of the heavy top is described by Euler's equations, see Equation (1.83) and (3.35) of Magnus [51], whose solution is taken as reference solution in the subsequent investigation. For high stiffnesses of the rod, we thus expect the solution to be close to the one of a rigid body.
Let the top be given by a cylinder of radius r = 0.1 and length L = 0.5 with cross-section area A = πr 2 and second moment of area I = πr 4 /4. The cylinder is subjected to a constant distributed line force I b = ρ 0 A(0, 0, −9.81). The stiff rod should be made out of steel with a uniform density ρ 0 = 8000, Young's modulus E = 210 × 10 6 and shear modulus G = E/(2(1 + ν)), with a Poisson's ratio ν = 1/3. The stiff rod has consequently an axial stiffness k e = EA, shear stiffness k s = GA, bending stiffness k b = k by = k bz = EI and torsional stiffness k t = 2GI. We also considered a soft rod for which all stiffnesses were reduced by a factor 10 3 . The top was discretized using a single two-node SE(3)-element. The initial position is such that the top points from the origin in positive e I x -direction, i.e., q 0 = (0 3×1 , 0 3×1 , r 1 , 0 3×1 ) with r 1 = (L, 0, 0). Its initial velocities are chosen such that in the case of a rigid rod a perfect precession motion[51, Section 3.3.2 c)] is obtained, i.e., u 0 = (0 3×1 , Ω, v 1 , Ω) with the tip velocity v 1 = Ω × r 1 and the angular velocity Ω = (Ω, 0, Ω pr ), where Ω = 50π and Ω pr = gL/(r 2 Ω). Finally, the motion of the top is constrained such that the first node coincides with the origin for all times. This can either be guaranteed by removing the corresponding degrees of freedom from the set of unknowns, or by using the concept of perfect bilateral constraints, cf. Géradin [35].
Using two different numerical time integration schemes (a standard fourth order Runge-Kutta method with the absolute and relative tolerances atol = rtol = 1 × 10 −8 , see Hairer et al. [37] and the generalized-α method for first order differential equations proposed by Jansen et al. [45] with the spectral radius at infinity ρ ∞ = 0.9 and a step size h = 1 × 10 −5 ), the simulations were performed until a final time of t 1 = 2π/Ω pr was reached, i.e., the rigid top performed a full rotation. Since the numerical solution of both methods only differ in their numerical digits, we only show the results obtained by the Runge-Kutta method.
In Figure 10 (a), the spatial trajectories of the different tops' free ends are shown. When comparing the projections of the rod's free tip, it can be seen that the assertion is true that the solution of a stiff rod cannot be distinguished from the rigid body solution, while the soft rod's tip performs a fascinating oscillatory motion superimposed to the rigid body solution.

Conclusion
We have presented a rod finite element formulation based on a two-node SE(3)-interpolation, that is, the interpolation of relative Euclidean transformation matrices using relative twists. In contrast to the typical uncoupled interpolation of centerline points and relative rotations using relative rotation vectors, the proposed SE(3)-interpolation leads to element-wise constant strain measures. Thus, by construction, the resulting finite element formulation will show an absence of both membrane and shear locking and can be applied in scenarios where very high slenderness ratios are present. Objectivity of the discretized strain measures followed from the interpolation strategy with relative Euclidean transformation matrices.
With a Petrov-Galerkin projection method it is possible to discretize the arising virtual work functionals in space resulting in a set of ordinary differential equations (in dynamics) or a set of nonlinear equations (in statics) both of which can be solved using standard numerical schemes. Thus, the drawback of existing SE(3) rod formulations, which require highly specialized Lie group solvers is circumvented in an elegant way. Further, for typical applications, the arising mass matrix is symmetric and constant. Since the virtual displacements and rotations are interpolated instead of using the consistent variations of the proposed SE(3)-interpolation formula, a non-symmetric stiffness matrix is obtained, which can be disadvantageous in terms of performance and storage. The nodal kinematics of the proposed formulation, i.e., the nodal Euclidean transformation matrices, are parameterized using total centerline points and total rotation vectors together with the SO(3)-exponential map. Arising singularities are circumvented by introducing the concept of complement rotation vectors together with the proposed relative interpolation strategy. Thus, a parametrization with a minimal number of six nodal unknowns has been achieved. The paper is closed by demonstrating all the individual properties of the proposed formulation by numerical benchmark examples in statics and dynamics.

A.1 The general linear group GL(n)
In the course of this treatment, we restrict ourselves to matrix Lie groups over the real numbers, that is, Lie groups G that are closed subgroups of the general linear group GL(n) = {A ∈ R n×n | det(A) = 0}, which is the group of all n × n invertible matrices with real entries [52]. For A, B ∈ GL(n), the group operation is given by the usual matrix product (A, B) → AB. The inverse and transpose of A are respectively denoted by A −1 and A T . The identity is given by the n × n identity matrix 1 n×n . Analogously, we denote the n × n matrix containing only zeros by 0 n×n . Let gl(n) = T 1n×n GL(n) be the tangent space of GL(n) at the identity 1 n×n . In Section 3.3 of Baker [53], it is shown that gl(n) corresponds to the set of all real n × n matrices, which constitute an n 2 -dimensional vector space. The Lie algebra of GL(n) consists of the vector space gl(n) together with the bilinear map which is a Lie bracket as it is skew-symmetric and satisfies the Jacobi For a fixed argument X ∈ gl(n) the Lie bracket defines the linear operator As in Section 2.6 of Iserles et al. [54], for a smooth curve Z(t) ∈ gl(n), the derivative of the exponential map can be written as d dt exp (Z(t)) = dexp Z(t) Ż (t) exp (Z(t)) = exp (Z(t)) dexp −Z(t) Ż (t) , where, for X ∈ gl(n), the linear map dexp X is defined in terms of a power series as Rearranging the equations in (62) to it is reasonable that dexp −X(t) and dexp X(t) are called left-and right-trivialized differentials, respectively. According to Equation (2.46) of Iserles et al. [54], the inverse of dexp X can also be written by the power series where B i denotes the i-th Bernoulli number. Let the tuple (g, [·, ·]) be a subalgebra of (gl(n), [·, ·]), i.e., g is a subspace of gl(n) with dimension k ≤ n 2 and [X, Y] ∈ g for X, Y ∈ g. Then, we can always find a linear and bijective map Since j and dexp −j(x) are linear maps, it is convenient to construct the tangent operator T(x) ∈ R k×k , which is the linear map defined as where in the last equality we have introduced the mapping Often this mapping is applied in the form j x y = ad j(x) j(y) , for x, y ∈ R k . If j(x) and j(y) commute, then powers of the adjoint map ad i j(x) (j(y)) vanish for i > 0 resulting in T(x)y = j −1 • dexp −j(x) •j y = y. The inverse of the tangent operator is given by Since the matrix rank of x is at most k, the Cayley-Hamilton theorem guarantees the existence of functions a 0 ( x ), . . . , a k−1 ( x ) such that Hence, we can express the tangent operator (67) and its inverse (69) by a finite number of powers of x . Using similar arguments, also the exponential map (61) can be expressed by a finite number of matrix powers.

A.2 The special orthogonal group SO(3)
As a subgroup of GL (3), the special orthogonal group is the group of all proper orthogonal transformations. In Baker [53] where we have introduced the tilde symbol for compact notation. This map is related to the cross product by ωr = ω×r for ω, r ∈ R 3 . In order to minimize the number of transcendental function evaluations [23,20], we introduce the mappings with lim ω →0 α( ω) = lim ω →0 β( ω) = lim ω →0 γ( ω) = 1. Noting that ω 3 = − ω 2 ω and carefully separating even and odd parts of the power series of the exponential map (61), see Murray et al. [55] Section 2.2, the exponential map from so(3) to SO(3) has an analytical form known as Rodrigues' formula For the case Ω → 0, following the typical approach in computational methods, the first order approximation exp SO(3) (Ω) = 1 3×3 + Ω is applied. From explicitly computing (75) and denoting ω = ω , the diagonal terms reveal the identity cos ω = 1 2 (tr(A) − 1), which can be solved for ω = ω(A). The off diagonals directly lead to ω = Again, for ω → 0, we use the first order approximation Ω = (A − A T )/2. For notational convenience, we further introduce the capitalized exponential and logarithm maps which directly relate R 3 with SO(3) and vice versa. For a lighter notation, we shortly suppress the subscripts of j and • indicating SO (3). Let x, y, z ∈ R 3 , using the skew-symmetry and Jacobi's identity of the cross product, direct computation verifies j x y z = ad j(x) j(y) z = j(x)j(y) − j(y)j(x) z = x y − y x z = j( xy)z. Comparison of both sides of the equality, readily implies x SO(3) = j SO(3) (x) = x. Again using ω 3 = − ω 2 ω and carefully separating the terms in the power series (67) (see Iserles et al. [54] Equation (B.10) and Equation (18) of Park and Chung [23]) 2 , the SO(3)-tangent operator can be written in the form Similarly from (69) (see Iserles et al. [54] Equation (B.11), Park and Chung [23] Equation (19) and Equation (2.20) of Bullo [22]), the inverse SO(3)-tangent operator can efficiently be computed via For ω → 0, we use the first order approximations T SO (3) (ω) = 1 3×3 − 1 2 ω and T −1 SO(3) (ω) = 1 3×3 + 1 2 ω. Due to the skew-symmetry of ω, the tangent map and its inverse fulfill the additional properties T SO (3)

B Variation of strain measures
Using the skew-symmetry of K δ φ IK , the expression (5) 2 can be written as δA T IK = − K δ φ IK A T IK . Consequently, the variation of Kγ from (8) computes to In the last equality, we have recognized the virtual displacement (2) 1 since the variation and ξ-derivative commute for vector components with respect to a constant basis. The same property gives rise to the equality where I e K i = A IK K e K i , i ∈ {x, y, z}. Straightforward computation together with (4) and (5) yields as well as Inserting the latter two expressions in (87) and making use of the Jacobi identity and the skew-symmetry of the cross product, one readily sees that As this holds for all i ∈ {x, y, z}, the term in the bracket of (90) must vanish, resulting in the identity

C Linearization of internal forces
In order to evaluate the linearization of the internal force vector f int , the derivatives of the exponential and logarithm maps of the underlying parametrization are required. Since most of these derivatives are not well established in literature and in order to be self-consistent, we briefly introduce them in the subsequent treatment. Most of the presented formulas have a removable singularity for ω = ω → 0. Thus, during a numerical implementation a critical angle ω crit = 1 × 10 −6 is defined. Whenever the value of ω falls below this critical angle, the first order approximations introduced in the first section are used. Thus, we have to introduce both, the derivative of the required formulas and the corresponding first order approximation. In the following, we identify the i-th component of a tuple a by a i , the i, j-th component of the matrix A by A ij and the components of third order objects are indicated by three subscripts. Using the identities ω ij = ω k ε kji , ∂ ω ij /∂ω k = ε kji and ∂( ω il ω lj )/∂ω k = ε kli ω lj + ω il ε kjl , where ε ijk denotes the Levi-Civita permutation symbol, we can compute the derivatives of the auxiliary functions introduced in (74), with respect to ω k by Thereby, the derivative of the rotation matrix components A ij , obtained from the capitalized SO(3) exponential map (77), with respect to ω k is Similar relations are found in Gallego and Yezzi [56].
Multiplying ω nm = ε kmn ω k with ε imn and using ε kmn ε imn = 2δ ki gives ω i = 1 2 ε imn ω nm . Thus, we have to compute Further, from cos ω = 1 2 (tr(A) − 1) with ∂ arccos(x)/∂x = −1/ Finally, using the identities ∂A ab /∂A cd = δ ac δ bd and ε imn (δ mj δ nk − δ nj δ mk ) = 2ε ijk , the derivative of ω i , obtained from the capitalized SO(3) logarithm map (77), with respect to the rotation matrix components A jk is given by Further, the derivatives of the SO(3) tangent operators are required. Using the identities introduced above they can be computed as and The derivatives of the capitalized SE(3) exponential and logarithm introduced in (85) boils down to correctly assemble the just computed derivatives. Recalling that θ = (v, ω), we get the derivative of the exponential map , With the derivative of the SE(3) logarithm map with respect to its argument is given by

D Discrete conservation properties
Following Romero and Armero [57], this section discusses the discrete conservation properties of the proposed rod finite element formulation. That is, the conservation of total energy, linear and angular momentum. For the sake of brevity, but without loss of generality, we consider the case for which the centerline points I r OP (ξ, t) coincide with the cross-sections' center of mass and for which S ρ0 of (14) vanishes.

D.1 Discrete conservation of total energy
With the discrete kinematics from (36) and (41), the total energy E = T + U of the discretized rod is given in terms of the kinetic energy and the potential energy Their respective time derivatives arė where we have used the identities ( Kγ ) · = A T IK ( I v P ) ,ξ − K ω IK × Kγ and ( KκIK ) · = ( K ω IK ) ,ξ − K ω IK × KκIK , which follow from the derivation given in Appendix B by replacing the variations with the time derivatives.
Since the principle of virtual work holds for arbitrary virtual displacements, we can choose the nodal virtual displacements δs i = ( I v Pi , Ki ω IKi ) to be composed of the nodal minimal velocities. Inserting these virtual displacements into (44) and (42), in the absence of any external force contributions, the virtual work principle leads to 0 = δW dyn (δs; q, u) + δW int (δs; q) = − T (q, u) + U (q) · .
Consequently, conservation of the total energy is guaranteed also for the discrete equations.

D.2 Discrete conservation of linear momentum
The discrete linear momentum of the Cosserat rod is For an arbitrary I c ∈ R 3 , let the nodal virtual displacements be given by This choice represents a virtual translation of the rod. Inserting (109) into the discrete internal virtual work (42) and recognizing that N e 0,ξ (ξ) + N e 1,ξ (ξ) = 0, one readily sees that δW int (δs; q) = 0. Using (109) in the discrete inertial virtual work (44) together with N e 0 (ξ) + N e 1 (ξ) = 1 leads to δW dyn (δs; q, u) = − I c T which, in agreement with the principle of virtual work, must vanish if no external forces are applied to the rod. Consequently, the rate of change of the linear momentum IL = 0 3×1 vanishes and the linear momentum I L is preserved.

D.3 Discrete conservation of angular momentum
The present rod finite element formulation uses nodal virtual rotations expressed in the individual cross-section-fixed bases Ki δφ IKi (t) ∈ R 3 . For an arbitrary I η ∈ R 3 , the nodal virtual displacements δs i = ( I η × I r OPi , A T IKi I η) lead to a virtual rotation around the origin of all nodal points. With respect to the inertial basis, by construction, the nodal cross-sections are all virtually rotated the same. However, due to the chosen interpolation (41) 3 , the cross-sections within the elements are not all rotated with the same constant virtual rotation I η. Hence, for this formulation a virtual rotation of the entire rod is not contained in the set of admissible virtual displacements. Consequently, there is no algorithmic access to the discrete angular momentum and it is not possible to construct a numerical time integration scheme that preserves the discrete total angular momentum.
If such a conservation property is crucial for a specific application, the present rod formulation can easily be modified. Instead of using the nodal virtual rotations Ki δφ IKi (t) ∈ R 3 expressed in the cross-section-fixed basis K i , the nodal virtual rotations I δφ IKi (t) ∈ R 3 expressed in the inertial basis I are used. Analogously, the nodal angular velocities Ki ω IKi (t) ∈ R 3 are replaced by I ω IKi (t) ∈ R 3 . Proceeding similar to the original formulation, new discretized internal and inertial virtual work contributions are obtained. They are δW int (δs; q) = δs T f int (q) , f int (q) = with I I ρ0 (q) = A IK (q) K I ρ0 A IK (q) T .
Note, with this reformulation the mass matrix M(q) and the gyroscopic forces f gyr (q, u) depend on the generalized coordinates q, which is clearly a drawback with respect to the original formulation. Conservation of energy and linear momentum follows in the same way as just shown and is not repeated here. As the virtual rotations are now formulated with respect to the inertial basis, the nodal virtual displacements induce a virtual rotation for the entire rod and conservation of total angular momentum can be shown. Indeed, for this discretization, the discrete angular momentum takes the form Inserting (114) into the discrete internal virtual work (111) and recognizing the identities already used for the linear momentum leads to a vanishing internal virtual work functional δW int (δs; q) = 0. Using (114) in the discrete inertial virtual work (112) results in which, according to the principle of virtual work, must be zero for vanishing external virtual work functionals. It readily follows that IJ = 0 3×1 and that the total angular momentum I J is conserved in these cases.