Nonunit quaternion parametrization of a Petrov–Galerkin Cosserat rod finite element

The application of the Petrov–Galerkin projection method in Cosserat rod finite element formulations offers significant advantages in simplifying the expressions within the discrete virtual work functionals. Moreover, it enables a straight‐forward and systematic exchange of the ansatz functions, specifically for centerline positions and cross‐section orientations. In this concise communication, we present a total Lagrangian finite element formulation for Cosserat rods that attempts to minimize the number of required theoretical concepts. The chosen discretization preserves objectivity and allows for large displacements/rotations and for large strains. The orientation parametrization with nonunit quaternions results in a singularity‐free formulation.


Cosserat rod theory
Let ξ ∈ J = [0, 1] ⊂ R be the centerline parameter and let t denote time.The motion of a Cosserat rod is captured by a time-dependent centerline curve represented in an inertial I-basis I r OP = I r OP (ξ, t) ∈ R 3 augmented by the cross-section orientations A IK = A IK (ξ, t) ∈ SO(3) = {A ∈ R 3×3 |A T A = 1 3×3 ∧ det(A) = 1}.The subscripts O and P in the centerline curve refer to the origin and the centerline point, respectively.The cross-section orientation A IK can also be interpreted as a transformation matrix that relates the representation of a vector in the cross-section-fixed K-basis to its representation in the inertial I-basis.
The derivatives with respect to time t and centerline parameter ξ are denoted by (•) and (•) ,ξ , respectively.The variation of a function is indicated by δ(•).With this, we can introduce the centerline velocity I v P = ( I r OP ) • and the virtual displacement The angular velocity of the cross-section-fixed K-basis relative to the inertial I-basis, in components with respect to the K-basis, is defined by K ω IK := j −1 A T IK (A IK ) • , where j : R 3 → so(3) = {B ∈ R 3×3 |B T = −B} is the linear and bijective map such that ωr = j(ω)r = ω × r for all ω, r ∈ R 3 .Analogously, the virtual rotations and the scaled curvature are defined as K δφ IK := j −1 A T IK δ (A IK ) and K κIK := j −1 A T IK (A IK ) ,ξ , respectively.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 increment is ds = Jdξ.The derivative with respect to the reference arc length s of a function f = f (ξ, t) ∈ R 3 can then be defined as f ,s (ξ, t) := f ,ξ (ξ, t)/J(ξ).The objective strain measures of a Cosserat rod are the curvature K κ IK = K κIK /J, which measures torsion and bending, together with the measures for dilatation and shear strains contained in K γ = K γ/J determined by K γ := (A IK ) T I r OP,ξ .
The internal virtual work of a Cosserat rod is defined as where K n and K m denote the resultant contact forces and moments, respectively.For hyperelastic material models with a strain energy density with respect to the reference arc length W = W ( K γ, K κ IK ; ξ), they can be determined by the constitutive relations Assume the line distributed external forces I b = I b(ξ, t) ∈ R 3 and moments K c = K c(ξ, t) ∈ R 3 to be given as densities with respect to the reference arc length.Moreover, for i ∈ {0, 1}, point forces I b i = I b i (t) ∈ R 3 and point moments be applied to the rod's boundaries at ξ 0 = 0 and ξ 1 = 1.The corresponding external virtual work functional is defined as In case I r OP is the line of centroids, the inertial virtual work functional of the Cosserat rod can be written as where A ρ0 is the cross-section mass density and K I ρ0 the constant cross-section inertia tensor represented in the cross-sectionfixed K-basis.

Petrov-Galerkin finite element formulation
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 .For a p-th order finite element, the closure of each of the intervals J e contains p + 1 evenly spaced points ξ e i ∈ cl(J e ) = [ξ e , ξ e+1 ] with i ∈ {0, . . ., p} such that ξ e 0 = ξ e < ξ e 1 < • • • < ξ e p = ξ e+1 .Note, for e ∈ {0, . . ., n el − 2}, the points ξ e p = ξ e+1 0 denote the same point ξ e+1 , which is the boundary point of the adjacent element intervals.It is convenient to use both indexations in the following.For a given element interval J e = [ξ e , ξ e+1 ), the p-th order Lagrange basis function and derivative of node i ∈ {0, . . ., p} are where ξ e i , ξ e j , and ξ e k are the points contained in the set {ξ e 0 = ξ e , ξ e 1 , . . ., ξ e p = ξ e+1 }.The centerline curve I r OP and the cross-section orientations A IK are approximated by interpolating nodal centerline points I r OP e i (t) ∈ R 3 and nodal transformation matrices A IK e i (t) ∈ SO(3).For each node i ∈ {0, . . ., p} within element e ∈ {0, . . ., n el − 1}, it will hold that I r OP e i (t) = I r OP (ξ e i , t) and A IK e i (t) = A IK (ξ e i , t).In contrast to [1,2], the nodal transformation matrices are parametrized by nodal non-unit quaternions P e i (t) = (p e 0,i (t), p e i (t)) ∈ R 4 with the scalar part p e 0,i (t) ∈ R and the vectorial part p e i (t) ∈ R 3 , see [7].Note that (5) is formulated in such a way to return orthogonal matrices also for non-unit quaternions.
Accordingly, the N = (pn el + 1) nodal generalized position coordinates q e i (t) = ( I r OP e i , P e i )(t) ∈ R 7 are given by the nodal centerline points I r OP e i and the nodal non-unit quaternions P e i resulting in n q = 7N positional degrees of freedom for the discretized rod.The nodal quantities can be assembled in the global tuple of generalized position coordinates q(t) = q 0 0 , . . ., q 0 p−1 , . . ., q e 0 , . . ., q e p−1 , . . ., For e ∈ {0, . . ., n el − 2}, the coordinates q e p = q e+1 0 refer to the same nodal coordinates.Introducing an appropriate Boolean connectivity matrix C e ∈ R 7(p+1)×nq , the element generalized position coordinates q e (t) = q e 0 , . . ., q e p (t) ∈ R 7(p+1) can be extracted from q via q e = C e q.Note that during a numerical implementation it is advisable to slice arrays instead of multiply them with Boolean matrices.
In the sense of [3,4], both the nodal centerline points and the cross-section orientations are interpolated by p-th order Lagrangian polynomials.Using the characteristic function χ J e : J → {0, 1}, which is one for ξ ∈ J e = [ξ e , ξ e+1 ) and zero elsewhere, together with the p-th order Lagrange basis functions (4), the ansatz functions for centerline and cross-section orientations are The discretized version of the curvature strain is computed as where the map Skw(M) = 1 2 (M − M T ) ∈ so(3) extracts the skew-symmetric part of the matrix M ∈ R 3×3 .Hence, the curvature can efficiently be computed using Copyright line will be provided by the publisher At the same N nodes as for the nodal generalized position coordinates, we introduce the nodal generalized virtual displacements δs e i (t) = ( I δr P e i , K e i δφ IK e i )(t) ∈ R 6 given by the nodal virtual centerline displacement I δr P e i (t) ∈ R 3 and the nodal virtual rotation K e i δφ IK e i (t) ∈ R 3 .In analogy to the generalized virtual displacements, we also introduce the nodal generalized velocities u e i (t) = ( I v P e i , K e i ω IK e i )(t) ∈ R 6 given by the nodal centerline velocity I v P e i (t) ∈ R 3 and the nodal angular velocity K e i ω IK e i (t) ∈ R 3 .Similar to the generalized position coordinates q, the nodal generalized virtual displacements and velocities are assembled in the global tuple of generalized virtual displacements δs(t) ∈ R nu and velocities u(t) ∈ R nu .In contrast to the nodal position coordinates, there are only six nodal generalized virtual displacements or velocity coordinates resulting in n u = 6N generalized virtual displacements or velocity degrees of freedom for the discretized rod.
Consequently, we require a new Boolean connectivity matrix C u,e ∈ R 6(p+1)×nu , which extracts the element generalized virtual displacements δs e (t) = (δs e 0 , . . ., δs e p )(t) ∈ R 6(p+1) and velocities u e (t) = (u e 0 , . . ., u e p )(t) ∈ R 6(p+1) from the global quantities via δs e = C u,e δs and u e = C u,e u.By further introducing the Boolean connectivity matrices C r,i ∈ R 3×6(p+1) , the nodal virtual centerline displacements I δr P e i and centerline velocities I v P e i can be extracted from the element generalized virtual displacements δs e and velocities u e via I δr P e i = C r,i δs e and I v P e i = C r,i u e , respectively.Identical extraction operations hold for the nodal virtual rotations K e i δφ IK e i = C φ,i δs e and angular velocities K e i ω IK e i = C φ,i u e , where C φ,i ∈ R 3×6(p+1) .The test functions are then given by interpolating the nodal generalized virtual displacements by p-th order Lagrangian basis functions (4) in agreement with χ J e (ξ) Note that the interpolation of the virtual rotations must be understood in the sense of a Petrov-Galerkin projection, where the virtual rotations are not obtained from a consistent variation of the ansatz functions (6).
To obtain a constant and symmetric mass matrix in the discretized formulation, see (13) below, the velocities are considered as independent fields and are interpolated with the same interpolation as the virtual displacements and rotations as The independent introduction of velocity fields ( 9) demands an additional relation defining the coupling between position coordinates q and velocity coordinates u.This coupling is given by the nodal kinematic differential equations i ω IK e i = F(q e i )u e i , where cf. [7].The nodal kinematic equations (10) can easily be assembled to a global kinematic differential equation of the form q = B(q)u.Note that the kinematic differential equation is linear in q too.This allows to write the relation also in the form q = D(u)q, see [7] for more details.
Inserting the test functions (8) together with the corresponding approximations for centerline, cross-section orientations (6) and strain measures into (1), the continuous internal virtual work is approximated by δW int (q; δs) = δs T f int (q), where the internal generalized forces are computed element-wise by Similarly, the external virtual work (2) is discretized by δW ext (t, q; δs) = δs T f ext (t, q) with Finally, inserting ( 8) and ( 9) into the inertial virtual work functional (3) yields the discrete counterpart δW dyn (u; δs) = −δs T M u + f gyr (u) , where we have introduced the symmetric and constant mass matrix and the gyroscopic forces Element integrals of the form J e f (ξ)dξ arising in the discretized external and gyroscopic forces, as well as in the mass matrix, are subsequently computed using a Gauss-Legendre quadrature rule with ceil[(p + 1) 2 /2] quadrature points.To alleviate locking, the internal generalized forces (11) are integrated by a reduced p-point quadrature rule.
Applying the principle of virtual work, which requires the total virtual work functional to vanish, we readily obtain the system dynamics in the form where the two lines correspond to the global kinematic differential equation and the equations of motion, respectively.Even though deviations from unit length of P e i do not affect the kinematic differential equation, to avoid numerical issues due to quaternion magnitudes near zero or floating point overflow, the nodal quaternions are normalized after each time-step, i.e., P e i = P e i / P i .For static problems, the n u = 6N nonlinear generalized force equilibrium equations must be augmented by the N constraint equations to ensure solvability.

Numerical experiments
In the following, the quadratic strain energy density is used.The superscript 0 refers to the evaluation in the rod's reference configuration.Moreover, K γ = diag(EA, GA, GA) and K κ = diag(G(I y + I z ), EI y , EI z ) denote the diagonal elasticity matrices with constant coefficients given by Saint-Venant's relations from linear elasticity.Therein, E and G, respectively denote the Young's and shear modulus.The crosssectional surface is denoted A and I y , I z are the respective second moments of area.
4.1 Helical spring Following [5], we investigate the elongation of an initially curved helical rod due to an applied external force at its tip, pointing in positive e I z -direction.The rod has a Young's modulus E = 10 11 N/m 2 and Poisson's ratio ν = 0.2, i.e., a shear modulus G = E/2(1 + ν).It has an undeformed shape of a perfect helix with n c = 10 coils, coil radius R = 10 mm, wire diameter d = 1 mm and unloaded pitch k = 5 mm, i.e., a total height of h = 50 mm.
In the simulation, the spring was discretized using 75 elements of the presented finite element formulation with p = 2. Reduced integration was performed with 2 quadrature points, while 5 points were used for all other integrals.The rod's curved initial configuration was obtained by solving the following minimization problem.Let ξ j = j m−1 ∈ [0, 1] for j ∈ {0, 1, . . ., m − 1} denote the m linearly spaced evaluation points of the reference helix curve Hence, the evaluation of the reference curve (19) at all ξ j 's leads to m target centerline points I r j = I r(ξ j ).Similarly, the corresponding cross-section orientations are given by evaluating the Serret-Frenet basis A IKj = ( I e for the individual ξ j 's.Following [1], the centerline positions and cross-section orientations can be assembled in the Euclidean transformations Copyright line will be provided by the publisher Using the SE(3)-logarithm map Log SE(3) introduced in [1], the optimal initial generalized position coordinates q 0 results from the nonlinear least squares problem in terms of the metric of relative twists.The minimization problem (21) can efficiently be solved using a Levenberg-Marquardt algorithm.The unity constraints of the nodal quaternions (17) can be incorporated into the optimization problem as equality constraints, albeit at the expense of employing a complex constrained nonlinear least squares solver.To simplify the process, we initially solved the unconstrained minimization problem and subsequently applied a projection step to normalize all nodal quaternions.Starting from q 0 , the maximal force of 100 N was applied within 500 linearly spaced force increments.During each iteration, the nonlinear equations ( 16) and (17) were solved up to an absolute error of 10 −8 .As can be seen in Fig. 1, the helical spring initially elongates proportional to the applied load.This is in line with classical helical spring theory [8], which assumes a linear force-displacement relation with linear equivalent stiffness Gd 4 /(64n c R 3 ) ≈ 65.1 N/m.When the elongation exceeds a certain value (approx.10 N), the linear theory does not longer agree with the numerically obtained nonlinear solution.This observation was also made by [5] and can be explained as follows.The helical spring unwinds gradually and approaches slowly a straight line with an altered linear stiffness EA.For comparison, we also solved the problem with the two-node SE(3)-interpolation strategy proposed in [1], using the same number of unknowns.As depicted in Fig. 1, the results are in line with the proposed quaternion formulation.

Wilberforce pendulum
More than 100 years ago, Lionel Robert Wilberforce did investigations On the Vibrations of a Loaded Spiral Spring [9].The experimental setup can be described as follows.While one end of a helical spring is clamped, at the other end a cylindrical bob is attached, see Fig. 2. When the cylinder in the gravitational field is displaced vertically, it starts oscillating up and down.Due to the coupling of bending and torsion of the deformed spring an additional torsional oscillation around the vertical axis of the cylinder is induced.When the cylinder's moment of inertia is properly adjusted, a beat phenomenon can be observed.In that case, the envelope of the vertical and torsional oscillations possess an almost perfect phase shift of π/2, i.e., the maximal amplitude of the vertical oscillations coincide with a zero torsional amplitude and vice versa.
To have a benchmark example that can be reproduced with reasonably computational effort, we introduce here a Wilberforce pendulum consisting of a spring with three coils modeled as a precurved rod.The rod has the properties of steel with mass density ρ 0 = 7850 kg/m 3 , shear modulus G = 81 • 10 9 N/m 2 and Poisson's ratio ν = 0.23, i.e., a Young's modulus E = 2G(1 + ν) = 199 • 10 9 N/m 2 .The undeformed shape is given by a perfect helix with n c = 3 coils, coil radius R = 16 mm, wire diameter d = 1 mm and an unloaded pitch of k = 1 mm.The bob is modeled as a cylindrical rigid body with radius r = 23 mm and height h = 36 mm also having the mass density of steel.
In the simulations, the rod was discretized using 18 elements of the presented Cosserat rod finite element with p = 2. Gravitational forces for the rod were neglected.Again, reduced integration was performed with 2 quadrature points, while for all other integrals 5 points were used.The bob was parameterized by the inertial position of the center of mass I r OS together with a non-unit quaternion P for the orientation.The bob was subjected to gravity with gravity constant g = 9.81 m/s 2 .For the governing equations describing such a parameterized rigid body under the influence of gravity, we refer to model 4 in [10].Cylinder and rod were rigidly connected by perfect bilateral constraints [11].Again, the optimal helical initial configuration q 0 was found by solving the minimization problem (21).The system was initialized at rest with initial velocity u 0 = 0.The resulting differential algebraic equations were solved using a first-order generalized-alpha method [12] for constrained mechanical systems of differential index 3, similar to the implementation found in [13].A constant step-size ∆t = 5 • 10 −3 s was chosen and the governing equations were solved up to a final time of t 1 = 8 s.Since the example includes high-frequency oscillations, we chose a spectral radius at infinity of ρ ∞ = 0.8.The internal Newton-Raphson method satisfied a tolerance of 10 −8 with respect to the maximum absolute error.In Fig. 2 the vertical position and the torsional angle of the rigid cylinder are plotted clearly showing the beat phenomenon of the Wilberforce pendulum.

Fig. 1 :
Fig. 1: Force displacement diagram and deformed configurations of the helical spring.

Fig. 2 :
Fig. 2: Vertical position z of the cylinder's center of mass and rotation angle α corresponding to the first Euler angle with sequence "zyx" over time.Snapshots of the Wilberforce pendulum are attached to the corresponding time instants.