#  A total Lagrangian, objective and intrinsically locking-free
Petrov–Galerkin ${S\\!E}(3)$ Cosserat rod finite element formulation

Jonas Harsch1, Simon Sailer1, Simon R. Eugster1
1Institute for Nonlinear Mechanics
University of Stuttgart
Pfaffenwaldring 9, 70569 Stuttgart, Germany

###### Abstract

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ć in 1999, is extended to the interpolation of nodal Euclidean
transformation matrices with the aid of relative twists; a strategy that
arises from the ${S\\!E}(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.

Keywords: rod finite elements, ${S\\!E}(3)$, objectivity, locking,
Petrov–Galerkin, large deformations

##  1 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
$\mathbb{R}^{3}$ and the cross-sectional orientations as elements of the
special orthogonal group $S\\!O(3)$ denoted by $\mathbb{R}^{3}\times
S\\!O(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 ${S\\!E}(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
${S\\!E}(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 $\pi$. 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 ${S\\!E}(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-$\alpha$) 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 recaptured. The section closes with the
formulation of the internal, external and inertial virtual work functionals of
the shear-deformable spatial rod. Section 3 introduces a novel Petrov–Galerkin
rod finite element formulation based on a very efficient two-node
${S\\!E}(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 $S\\!O(3)$ and the special
Euclidean group ${S\\!E}(3)$ are examined as subgroups of $G\\!L(3)$ and
$G\\!L(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 ${S\\!E}(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.

##  2 Cosserat rod theory

Figure 1: Kinematics of the centerline curve with their attached orthonormal
basis vectors.

###  2.1 Centerline and cross-section orientation

We introduce the Euclidean 3-space $\mathbb{E}^{3}$ as an abstract
3-dimensional real inner-product space[5]. A basis for $\mathbb{E}^{3}$ is a
linearly independent set of three vectors
$\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z}\in\mathbb{E}^{3}$. The basis is
said to be right-handed if
$\mathbf{e}_{x}\cdot(\mathbf{e}_{y}\times\mathbf{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
$K=\\{\mathbf{e}_{x}^{K},\ \mathbf{e}_{y}^{K},\ \mathbf{e}_{z}^{K}\\}$, the
respective components $a_{i}^{K}$, $i\in\\{x,y,z\\}$ of a vector
$\mathbf{a}=a_{x}^{K}\mathbf{e}_{x}^{K}+a_{y}^{K}\mathbf{e}_{y}^{K}+a_{z}^{K}\mathbf{e}_{z}^{K}\in\mathbb{E}^{3}$
can be collected in the triple ${}_{K}\mathbf{a}=(a_{x}^{K},\ a_{y}^{K},\
a_{z}^{K})\in\mathbb{R}^{3}$. Thus, we carefully distinguish $\mathbb{R}^{3}$
from the Euclidean 3-space $\mathbb{E}^{3}$. The rotation between two bases
$K_{0}$ and $K_{1}$ is captured by the proper orthogonal transformation matrix
$\mathbf{A}_{K_{0}K_{1}}\in S\\!O(3)$, which relates the coordinate
representations ${}_{K_{0}}\mathbf{a}$ and ${}_{K_{1}}\mathbf{a}$ in
accordance with
${}_{K_{0}}\mathbf{a}=\mathbf{A}_{K_{0}K_{1}}\,{}_{K_{1}}\mathbf{a}$. Let
$\xi\in\mathcal{J}=[0,1]\subset\mathbb{R}$ denote the centerline parameter and
$(\eta,\zeta)\in\mathcal{A}(\xi)\subset\mathbb{R}^{2}$ the cross-section
parameters of the cross-section area at $\xi$. Considering the rod as a three-
dimensional continuum, a point $Q$ of the rod can be addressed by

${}_{I}\mathbf{r}_{OQ}(\bm{\xi},t)={}_{I}\mathbf{r}_{OP}(\xi,t)+\mathbf{A}_{IK}(\xi,t)\,{}_{K}\mathbf{r}_{PQ}(\xi,\eta,\zeta)\,,\quad\bm{\xi}=(\xi,\
\eta,\
\zeta)\in\mathcal{B}=\mathcal{J}\times\mathcal{A}(\mathcal{J})\subset\mathbb{R}^{3}\,.$
(1)

Herein, ${}_{I}\mathbf{r}_{OP}$ are the Cartesian components with respect to
the inertial basis $I$ of the time-dependent centerline curve
$\mathbf{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 $\xi$ and $\bm{\xi}$, respectively. Moreover, ${}_{K}\mathbf{r}_{PQ}$
denote the cross-section coordinates, which are transformed to the $I$-basis
using the transformation matrices $\mathbf{A}_{IK}$. Note that the base
vectors of the cross-section-fixed $K$-bases are functions of the centerline
parameter $\xi$ and time $t$, i.e.,
$\mathbf{e}_{i}^{K}=\mathbf{e}_{i}^{K}(\xi,t),\,i\in\\{x,y,z\\}$.

###  2.2 Variations, velocities and curvature

Let $\dot{(\bullet)}$ and $(\bullet)_{,\xi}$ respectively denote the
derivative with respect to time $t$ and centerline parameter $\xi$. The
variation is denoted by $\delta(\bullet)$. The virtual displacement
${}_{I}\delta\mathbf{r}_{P}$ and the centerline velocity
${}_{I}\mathbf{v}_{P}$ are given by the variation as well as the time
derivative of the centerline

${}_{I}\delta\mathbf{r}_{P}=\delta\left({}_{I}\mathbf{r}_{OP}\right)\,,\quad{}_{I}\mathbf{v}_{P}={}_{I}\dot{\mathbf{r}}_{OP}\,.$
(2)

The angular velocity of the $K$-basis relative to the inertial $I$-basis, in
components w.r.t. the $K$-basis, is defined by

${}_{K}\bm{\omega}_{IK}\coloneqq
j^{-1}_{S\\!O(3)}\left({}_{K}\widetilde{\bm{\omega}}_{IK}\right)\,,\quad\text{with}\quad{}_{K}\widetilde{\bm{\omega}}_{IK}\coloneqq\mathbf{A}_{IK}^{\\!\mathrm{T}}\dot{\mathbf{A}}_{IK}\,,$
(3)

where
$j_{S\\!O(3)}\colon\mathbb{R}^{3}\to\mathfrak{so}(3)=\\{\mathbf{B}\in\mathbb{R}^{3\times
3}|\mathbf{B}^{\\!\mathrm{T}}=-\mathbf{B}\\}$ is the linear and bijective map
such that
$\widetilde{\bm{\omega}}\mathbf{r}=j_{S\\!O(3)}(\bm{\omega})\mathbf{r}=\bm{\omega}\times\mathbf{r}$
for all $\bm{\omega},\mathbf{r}\in\mathbb{R}^{3}$, see (72) from Appendix A.
Analogously, we define the scaled curvature as

${}_{K}\bar{\bm{\kappa}}_{IK}\coloneqq
j^{-1}_{S\\!O(3)}\left({}_{K}\widetilde{\bar{\bm{\kappa}}}_{IK}\right)\,,\quad\text{with}\quad{}_{K}\widetilde{\bar{\bm{\kappa}}}_{IK}\coloneqq\mathbf{A}_{IK}^{\\!\mathrm{T}}\mathbf{A}_{IK,\xi}$
(4)

and the virtual rotation as

${}_{K}\delta\bm{\phi}_{IK}\coloneqq
j^{-1}_{S\\!O(3)}\left({}_{K}\delta\widetilde{\bm{\phi}}_{IK}\right)\,,\quad\text{with}\quad{}_{K}\delta\widetilde{\bm{\phi}}_{IK}\coloneqq\mathbf{A}_{IK}^{\\!\mathrm{T}}\delta\mathbf{A}_{IK}\,.$
(5)

###  2.3 Reference arc length

For the reference centerline curve ${}_{I}\mathbf{r}_{OP}^{0}$, the length of
the rod’s tangent vector is $J=\|{}_{I}\mathbf{r}_{OP,\xi}^{0}\|$. Thus, for a
given centerline parameter $\xi$, the reference arc length $s$ is defined by

$s(\xi)\coloneqq\int_{0}^{\xi}J(\bar{\xi})\,\mathrm{d}\bar{\xi}\,.$ (6)

Following Harsch and Eugster[26], the derivative with respect to the reference
arc length $s$ of a function
$\mathbf{f}\colon\mathcal{J}\times\mathbb{R}\to\mathbb{R}^{3}$ can be computed
by

$\mathbf{f}_{,s}(\xi,t)\coloneqq\mathbf{f}_{,\xi}(\xi,t)\frac{1}{J(\xi)}\,.$
(7)

###  2.4 Objective strain measures

The objective strain measures of a Cosserat rod, see Antman[5, Section 8.2],
are given by

${}_{K}\bm{\kappa}_{IK}={}_{K}\bar{\bm{\kappa}}_{IK}/J\quad\text{and}\quad{}_{K}\bm{\gamma}={}_{K}\bar{\bm{\gamma}}/J\,,\quad\text{with}\quad{}_{K}\bar{\bm{\gamma}}\coloneqq\mathbf{A}_{IK}^{\\!\mathrm{T}}{}_{I}\mathbf{r}_{OP,\xi}\,,$
(8)

which can be gathered in the six-dimensional tuple
$\bm{\varepsilon}=({}_{K}\bm{\gamma},\
{}_{K}\bm{\kappa}_{IK})\in\mathbb{R}^{6}$. Therein, the dilatation and shear
is captured by ${}_{K}\bm{\gamma}$, while ${}_{K}\bm{\kappa}_{IK}$ measures
torsion and bending, see Antman[5, Section 8.6].

###  2.5 Internal virtual work

Without loss of generality, we restrict ourselves to hyperelastic material
models where the strain energy density with respect to the reference arc
length $W=W({}_{K}\bm{\gamma},{}_{K}\bm{\kappa}_{IK};\xi)$ depends on the
strain measures (8) and possibly explicitly on the centerline parameter $\xi$.
By that, the internal virtual work functional is defined as

$\displaystyle\delta W^{\mathrm{int}}\coloneqq$
$\displaystyle-\int_{\mathcal{J}}\delta
WJ\mathrm{d}\xi=-\int_{\mathcal{J}}\left\\{\delta({}_{K}\bar{\bm{\gamma}})^{\\!\mathrm{T}}{}_{K}\mathbf{n}+\delta({}_{K}\bar{\bm{\kappa}}_{IK})^{\\!\mathrm{T}}{}_{K}\mathbf{m}\right\\}\mathrm{d}\xi\,,$
(9)

where we have introduced the constitutive equations

${}_{K}\mathbf{n}\coloneqq\left(\frac{\partial
W}{\partial{}_{K}\bm{\gamma}}\right)^{\\!\mathrm{T}}\,,\quad{}_{K}\mathbf{m}\coloneqq\left(\frac{\partial
W}{\partial{}_{K}\bm{\kappa}_{IK}}\right)^{\\!\mathrm{T}}\,.$ (10)

Note that even in the inelastic case, where no strain energy density $W$ is
available, the internal virtual work (9) can be used, with internal forces and
moments ${}_{K}\mathbf{n}$ and ${}_{K}\mathbf{m}$ given by different
constitutive laws[8]. Using $\delta({}_{K}\bm{\gamma})$ and
$\delta({}_{K}\bm{\kappa}_{IK})$ derived in Appendix B, the internal virtual
work (9) takes the form

$\delta
W^{\mathrm{int}}=-\int_{\mathcal{J}}\Big{\\{}({}_{I}\delta\mathbf{r}_{P})_{,\xi}^{\\!\mathrm{T}}\,\mathbf{A}_{IK}\,{}_{K}\mathbf{n}+({}_{K}\delta\bm{\phi}_{IK})_{,\xi}^{\\!\mathrm{T}}\,{}_{K}\mathbf{m}-{}_{K}\delta\bm{\phi}_{IK}^{\\!\mathrm{T}}\left[{}_{K}\bar{\bm{\gamma}}\times{}_{K}\mathbf{n}+{}_{K}\bar{\bm{\kappa}}_{IK}\times{}_{K}\mathbf{m}\right]\Big{\\}}\mathrm{d}\xi\,.$
(11)

As in Equation (2.10) of Simo and Vu-Quoc[4], we introduce the diagonal
elasticity matrices
$\mathbf{C}_{\bm{\gamma}}=\operatorname{diag}(k_{\mathrm{e}},k_{\mathrm{s}},k_{\mathrm{s}})$
and
$\mathbf{C}_{\bm{\kappa}}=\operatorname{diag}(k_{\mathrm{t}},k_{\mathrm{b}_{y}},k_{\mathrm{b}_{z}})$
with constant coefficients. In the following, the simple quadratic strain
energy density

$W({}_{K}\bm{\gamma},{}_{K}\bm{\kappa}_{IK};\xi)=\frac{1}{2}\left({}_{K}\bm{\gamma}-{}_{K}\bm{\gamma}^{0}\right)^{\\!\mathrm{T}}\mathbf{C}_{\bm{\gamma}}\left({}_{K}\bm{\gamma}-{}_{K}\bm{\gamma}^{0}\right)+\frac{1}{2}\left({}_{K}\bm{\kappa}_{IK}-{}_{K}\bm{\kappa}_{IK}^{0}\right)^{\\!\mathrm{T}}\mathbf{C}_{\bm{\kappa}}\left({}_{K}\bm{\kappa}_{IK}-{}_{K}\bm{\kappa}_{IK}^{0}\right)$
(12)

is used, where the superscript $0$ refers to the evaluation in the rod’s
reference configuration.

###  2.6 External virtual work

Assume the line distributed external forces
${}_{I}\mathbf{b}\colon\mathcal{J}\times\mathbb{R}\to\mathbb{R}^{3}$ and
moments ${}_{K}\mathbf{c}\colon\mathcal{J}\times\mathbb{R}\to\mathbb{R}^{3}$
to be given as densities with respect to the reference arc length. Moreover,
for $i\in\\{0,1\\}$, point forces
${}_{I}\mathbf{b}_{i}\colon\mathbb{R}\to\mathbb{R}^{3}$ and point moments
${}_{K}\mathbf{c}_{i}\colon\mathbb{R}\to\mathbb{R}^{3}$ can be applied to the
rod’s boundaries at $\xi_{0}=0$ and $\xi_{1}=1$. By that, the corresponding
external virtual work functional is given by

$\delta
W^{\mathrm{ext}}=\int_{\mathcal{J}}\left\\{{}_{I}\delta\mathbf{r}_{P}^{\\!\mathrm{T}}{}_{I}\mathbf{b}+{}_{K}\delta\bm{\phi}_{IK}^{\\!\mathrm{T}}{}_{K}\mathbf{c}\right\\}J\mathrm{d}\xi+\sum_{i=0}^{1}\left[{}_{I}\delta\mathbf{r}_{P}^{\\!\mathrm{T}}{}_{I}\mathbf{b}_{i}+{}_{K}\delta\bm{\phi}_{IK}^{\\!\mathrm{T}}{}_{K}\mathbf{c}_{i}\right]_{\xi_{i}}\,.$
(13)

###  2.7 Inertial virtual work:

Let $\rho_{0}\colon\mathcal{J}\to\mathbb{R}$ denote the rod’s mass density per
unit reference volume and $\mathrm{d}A$ the cross-section surface element. It
is convenient to define the following abbreviations

$A_{\rho_{0}}(\xi)\coloneqq\int_{\mathcal{A}(\xi)}\rho_{0}\,\mathrm{d}A\,,\quad{}_{K}\mathbf{S}_{\rho_{0}}(\xi)\coloneqq\int_{\mathcal{A}(\xi)}\rho_{0}{}_{K}\widetilde{\mathbf{r}}_{PQ}\,\mathrm{d}A\,,\quad{}_{K}\mathbf{I}_{\rho_{0}}(\xi)\coloneqq\int_{\mathcal{A}(\xi)}\rho_{0}{}_{K}\widetilde{\mathbf{r}}_{PQ}\,{}_{K}\widetilde{\mathbf{r}}_{PQ}^{\,\mathrm{T}}\mathrm{d}A\,.$
(14)

Further, using the mass differential
$\mathrm{d}m=\rho_{0}J\mathrm{d}A\mathrm{d}\xi$ and expressing the variation
and second time derivative of $\mathbf{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)

$\displaystyle\delta W^{\mathrm{dyn}}$
$\displaystyle=-\int_{\mathcal{B}}{}_{I}\delta\mathbf{r}_{OQ}^{\\!\mathrm{T}}{}_{I}\ddot{\mathbf{r}}_{OQ}\,\mathrm{d}m=-\int_{\mathcal{J}}\begin{pmatrix}{}_{I}\delta\mathbf{r}_{P}\\\
{}_{K}\delta\bm{\phi}_{IK}\end{pmatrix}^{\\!\mathrm{T}}\left\\{\overline{\mathbf{M}}\begin{pmatrix}{}_{I}\mathbf{v}_{P}\\\
{}_{K}\bm{\omega}_{IK}\end{pmatrix}^{{\cdot}}+\mathbf{g}\right\\}J\mathrm{d}\xi\,,$
(15)

where we have introduced the two quantities

$\overline{\mathbf{M}}=\begin{pmatrix}A_{\rho_{0}}\mathbf{1}_{3\times
3}&\mathbf{A}_{IK}{}_{K}\mathbf{S}_{\rho_{0}}^{\\!\mathrm{T}}\\\
{}_{K}\mathbf{S}_{\rho_{0}}\mathbf{A}_{IK}^{\\!\mathrm{T}}&{}_{K}\mathbf{I}_{\rho_{0}}\end{pmatrix}\,,\quad\mathbf{g}=\begin{pmatrix}\mathbf{A}_{IK}{}_{K}\widetilde{\bm{\omega}}_{IK}{}_{K}\mathbf{S}_{\rho_{0}}^{\\!\mathrm{T}}{}_{K}\bm{\omega}_{IK}\\\
{}_{K}\widetilde{\bm{\omega}}_{IK}{}_{K}\mathbf{I}_{\rho_{0}}{}_{K}\bm{\omega}_{IK}\end{pmatrix}\,.$
(16)

Note that the quantity $\mathbf{S}_{\rho_{0}}$ vanishes if the centerline
points ${}_{I}\mathbf{r}_{OP}(\xi,t)$ coincide with the cross-sections’ center
of mass. For this case $\overline{\mathbf{M}}$ and $\mathbf{g}$ get
independent of the rod’s configuration.

##  3 Finite element formulation

###  3.1 Rod kinematics in homogenous coordinates

The frame
$\mathcal{I}=\\{O,\mathbf{e}_{x}^{I},\mathbf{e}_{y}^{I},\mathbf{e}_{z}^{I}\\}$
is the set collecting the origin $O$ together with the inertial base vectors
$\mathbf{e}_{i}^{I}$, $i\in\\{x,y,z\\}$. A generic and possibly non-inertial
frame
$\mathcal{K}=\\{P,\mathbf{e}_{x}^{K},\mathbf{e}_{y}^{K},\mathbf{e}_{z}^{K}\\}$
is given by a point $P$ together with a right-handed orthonormal basis $K$.
The point $Q$ relative to the $\mathcal{K}$-frame is addressed by the triple
${}_{K}\mathbf{r}_{PQ}\in\mathbb{R}^{3}$, which are the components of the
vector $\mathbf{r}_{PQ}\in\mathbb{E}^{3}$ between $P$ and $Q$ with respect to
the $K$-basis. Rigid body motions between two frames $\mathcal{K}_{0}$ and
$\mathcal{K}_{1}$ are captured by the Euclidean transformation matrix

$\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}_{1}}=\begin{pmatrix}\mathbf{A}_{K_{0}K_{1}}&{}_{K_{0}}\mathbf{r}_{P_{0}P_{1}}\\\
\mathbf{0}_{1\times 3}&1\end{pmatrix}\,,$ (17)

which relates the triples ${}_{K_{1}}\mathbf{r}_{P_{1}Q}\in\mathbb{R}^{3}$ and
${}_{K_{0}}\mathbf{r}_{P_{0}Q}\in\mathbb{R}^{3}$ in accordance with

$\begin{pmatrix}{}_{K_{0}}\mathbf{r}_{P_{0}Q}\\\
1\end{pmatrix}=\begin{pmatrix}{}_{K_{0}}\mathbf{r}_{P_{0}P_{1}}+\mathbf{A}_{K_{0}K_{1}}\,{}_{K_{1}}\mathbf{r}_{P_{1}Q}\\\
1\end{pmatrix}=\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}_{1}}\begin{pmatrix}{}_{K_{1}}\mathbf{r}_{P_{1}Q}\\\
1\end{pmatrix}\,.$ (18)

In fact, the Euclidean transformation matrix
$\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}_{1}}$ is an element of the special
Euclidean group ${S\\!E}(3)$ considered here as a Lie subgroup of the general
linear group $G\\!L(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
$\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}_{1}}$ is

$\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}_{1}}^{-1}=\begin{pmatrix}\mathbf{A}_{K_{0}K_{1}}^{\\!\mathrm{T}}&-\mathbf{A}_{K_{0}K_{1}}^{\\!\mathrm{T}}{}_{K_{0}}\mathbf{r}_{P_{0}P_{1}}\\\
\mathbf{0}_{1\times 3}&1\end{pmatrix}\,.$ (19)

Consequently, using (18), the motion of the rod (1) can be written in
homogenous coordinates as

$\begin{pmatrix}{}_{I}\mathbf{r}_{OQ}\\\
1\end{pmatrix}=\mathbf{H}_{\mathcal{I}\mathcal{K}}\begin{pmatrix}{}_{K}\mathbf{r}_{PQ}\\\
1\end{pmatrix}\,,\quad\text{with}\quad\mathbf{H}_{\mathcal{I}\mathcal{K}}=\begin{pmatrix}\mathbf{A}_{IK}&{}_{I}\mathbf{r}_{OP}\\\
\mathbf{0}_{1\times 3}&1\end{pmatrix}\,.$ (20)

The group structure of ${S\\!E}(3)$ then allows to decompose the Euclidean
transformation matrix $\mathbf{H}_{\mathcal{I}\mathcal{K}}$ into

$\mathbf{H}_{\mathcal{I}\mathcal{K}}=\mathbf{H}_{\mathcal{I}\mathcal{K}_{0}}\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}\,.$
(21)

###  3.2 ${S\\!E}(3)$-interpolation

In order to introduce the idea of the ${S\\!E}(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
$\mathbf{H}_{\mathcal{I}\mathcal{K}_{0}}$ and
$\mathbf{H}_{\mathcal{I}\mathcal{K}_{1}}$ in terms of time-dependent
generalized coordinates
$\mathbf{q}(t)=\big{(}{}_{I}\mathbf{r}_{OP_{0}}(t),\bm{\psi}_{0}(t),{}_{I}\mathbf{r}_{OP_{1}}(t),\bm{\psi}_{1}(t)\big{)}\in\mathbb{R}^{12}$
given by the nodal positions ${}_{I}\mathbf{r}_{OP_{0}}(t)\in\mathbb{R}^{3}$,
${}_{I}\mathbf{r}_{OP_{1}}(t)\in\mathbb{R}^{3}$ and the nodal rotation vectors
$\bm{\psi}_{0}(t)\in\mathbb{R}^{3}$, $\bm{\psi}_{1}(t)\in\mathbb{R}^{3}$.
Using the exponential map of $S\\!O(3)$, defined in (76), the parametrization
is

$\mathbf{H}_{\mathcal{I}\mathcal{K}_{0}}(\mathbf{q})=\begin{pmatrix}\operatorname{Exp}_{S\\!O(3)}(\bm{\psi}_{0})&{}_{I}\mathbf{r}_{OP_{0}}\\\
\mathbf{0}_{3\times
1}&1\end{pmatrix}\quad\text{and}\quad\mathbf{H}_{\mathcal{I}\mathcal{K}_{1}}(\mathbf{q})=\begin{pmatrix}\operatorname{Exp}_{S\\!O(3)}(\bm{\psi}_{1})&{}_{I}\mathbf{r}_{OP_{1}}\\\
\mathbf{0}_{3\times 1}&1\end{pmatrix}\,.$ (22)

With the ${S\\!E}(3)$-logarithm map from (84), we can compute the relative
twist $\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}}$ corresponding to the
relative Euclidean transformation
$\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}_{1}}=\mathbf{H}_{\mathcal{I}\mathcal{K}_{0}}^{-1}\,\mathbf{H}_{\mathcal{I}\mathcal{K}_{1}}$
as

$\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}}=\operatorname{Log}_{{S\\!E}(3)}\big{(}\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}_{1}}\big{)}=\begin{pmatrix}\mathbf{T}_{S\\!O(3)}^{-\mathrm{T}}(\bm{\psi}_{01})\,{}_{K_{0}}\mathbf{r}_{P_{0}P_{1}}\\\
\bm{\psi}_{01}\end{pmatrix}\,.$ (23)

Herein $\bm{\psi}_{01}$ corresponds to the relative rotation vector
parameterizing the transformation between the $K_{0}$\- and $K_{1}$-basis
according to
$\mathbf{A}_{K_{0}K_{1}}=\operatorname{Exp}_{S\\!O(3)}(\bm{\psi}_{01})$. Due
to the singularity of $\operatorname{Log}_{S\\!O(3)}$ for
$\|\bm{\psi}_{01}\|=\pi$, see Appendix A.2, this interpolation strategy is
restricted to applications in which the relative rotation angle satisfies
$\omega=\|\bm{\psi}_{01}\|<\pi$. A discretization with a higher number of
elements always cures this problem. By linearly scaling the relative twist
$\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}}$ and using the
${S\\!E}(3)$-exponential map (84), a relative Euclidean transformation from
the $\mathcal{K}_{0}$-frame to the $\mathcal{K}$-frame can be constructed as

$\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}(\xi,\mathbf{q})=\operatorname{Exp}_{{S\\!E}(3)}\big{(}\xi\,\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}}(\mathbf{q})\big{)}\,.$
(24)

The Euclidean transformation $\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}$
satisfies
$\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}(0,\mathbf{q})=\mathbf{1}_{4\times 4}$
and
$\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}(1,\mathbf{q})=\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}_{1}}$.
To obtain the Euclidean transformation for a point within the finite element,
the expressions (21) and (24) motivate the following
${S\\!E}(3)$-interpolation

$\mathbf{H}_{\mathcal{I}\mathcal{K}}(\xi,\mathbf{q})\coloneqq\mathbf{H}_{\mathcal{I}\mathcal{K}_{0}}(\mathbf{q})\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}(\xi,\mathbf{q})=\mathbf{H}_{\mathcal{I}\mathcal{K}_{0}}(\mathbf{q})\operatorname{Exp}_{{S\\!E}(3)}\left(\xi\,\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}}(\mathbf{q})\right)\,,$
(25)

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

$\displaystyle\mathbf{H}_{\mathcal{I}\mathcal{K}}$
$\displaystyle=\begin{pmatrix}\mathbf{A}_{IK_{0}}&{}_{I}\mathbf{r}_{OP_{0}}\\\
\mathbf{0}_{1\times
3}&1\end{pmatrix}\begin{pmatrix}\operatorname{Exp}_{S\\!O(3)}(\xi\bm{\psi}_{01})&\quad\mathbf{T}_{S\\!O(3)}^{\\!\mathrm{T}}(\xi\bm{\psi}_{01})\,\xi\,\mathbf{T}_{S\\!O(3)}^{-\mathrm{T}}(\bm{\psi}_{01}){}_{K_{0}}\mathbf{r}_{P_{0}P_{1}}\\\
\mathbf{0}_{1\times 3}&1\end{pmatrix}$ (26)
$\displaystyle\begin{pmatrix}\mathbf{A}_{IK}&{}_{I}\mathbf{r}_{OP}\\\
\mathbf{0}_{1\times 3}&1\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\mathbf{A}_{IK_{0}}\operatorname{Exp}_{S\\!O(3)}(\xi\bm{\psi}_{01})&\quad{}_{I}\mathbf{r}_{OP_{0}}+\xi\,\mathbf{A}_{IK_{0}}\mathbf{T}_{S\\!O(3)}^{\\!\mathrm{T}}(\xi\bm{\psi}_{01})\mathbf{T}_{S\\!O(3)}^{-\mathrm{T}}(\bm{\psi}_{01}){}_{K_{0}}\mathbf{r}_{P_{0}P_{1}}\\\
\mathbf{0}_{1\times 3}&1\end{pmatrix}\,.$

Consequently, the rod’s orientation is discretized by

$\mathbf{A}_{IK}(\xi,\mathbf{q})=\mathbf{A}_{IK_{0}}(\mathbf{q})\operatorname{Exp}_{S\\!O(3)}\big{(}\xi\bm{\psi}_{01}(\mathbf{q})\big{)}\,,$
(27)

which corresponds to Equation (4.7) of Crisfield and Jelenić[13]. Since
$\mathbf{A}_{IK}(0,\mathbf{q})=\mathbf{A}_{IK_{0}}$ and
$\mathbf{A}_{IK}(1,\mathbf{q})=\mathbf{A}_{IK_{1}}$, the discretization (27)
is a highly nonlinear interpolation of the nodal transformation matrices. The
rod’s centerline is discretized by

${}_{I}\mathbf{r}_{OP}(\xi,\mathbf{q})={}_{I}\mathbf{r}_{OP_{0}}+\xi\,\mathbf{A}_{IK_{0}}(\mathbf{q})\mathbf{T}_{S\\!O(3)}^{\\!\mathrm{T}}\big{(}\xi\bm{\psi}_{01}(\mathbf{q})\big{)}\mathbf{T}_{S\\!O(3)}^{-\mathrm{T}}\big{(}\bm{\psi}_{01}(\mathbf{q})\big{)}{}_{K_{0}}\mathbf{r}_{P_{0}P_{1}}(\mathbf{q})\,,$
(28)

which also interpolates the nodal positions
${}_{I}\mathbf{r}_{OP}(0,\mathbf{q})={}_{I}\mathbf{r}_{OP_{0}}$ and
${}_{I}\mathbf{r}_{OP}(1,\mathbf{q})={}_{I}\mathbf{r}_{OP_{1}}$ in a highly
nonlinear manner.

###  3.3 Objectivity, discretized strain measures and absence of locking of
${S\\!E}(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

$\mathbf{H}_{\mathcal{I}\mathcal{K}}^{-1}\mathbf{H}_{\mathcal{I}\mathcal{K},s}=\begin{pmatrix}\mathbf{A}_{IK}^{\\!\mathrm{T}}&-\mathbf{A}_{IK}^{\\!\mathrm{T}}{}_{I}\mathbf{r}_{OP}\\\
\mathbf{0}_{1\times
3}&1\end{pmatrix}\begin{pmatrix}\mathbf{A}_{IK,s}&{}_{I}\mathbf{r}_{OP,s}\\\
\mathbf{0}_{1\times
3}&0\end{pmatrix}=\begin{pmatrix}\mathbf{A}_{IK}^{\\!\mathrm{T}}\mathbf{A}_{IK,s}&\mathbf{A}_{IK}^{\\!\mathrm{T}}{}_{I}\mathbf{r}_{OP,s}\\\
\mathbf{0}_{1\times
3}&0\end{pmatrix}=\begin{pmatrix}{}_{K}\widetilde{\bm{\kappa}}_{IK}&{}_{K}\bm{\gamma}\\\
\mathbf{0}_{1\times 3}&0\end{pmatrix}\,.$ (29)

With the proposed interpolation (25), the strain measures in

$\mathbf{H}_{\mathcal{I}\mathcal{K}}^{-1}\mathbf{H}_{\mathcal{I}\mathcal{K},s}=\big{(}\mathbf{H}_{\mathcal{I}\mathcal{K}_{0}}\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}\big{)}^{-1}\big{(}\mathbf{H}_{\mathcal{I}\mathcal{K}_{0}}\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}\big{)}_{,s}=\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}^{-1}\mathbf{H}_{\mathcal{K}_{0}\mathcal{K},s}$
(30)

depend only on the relative Euclidean transformation
$\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}$ 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 $\mathcal{I}$-frame. Whether
we use an $\mathcal{I}$-frame or an $\mathcal{I}^{+}$-frame, which are related
by the time-dependent Euclidean transformation
$\mathbf{H}_{\mathcal{I}^{+}\mathcal{I}}$, the discretized strain measures in
$\mathbf{H}_{\mathcal{I}^{+}\mathcal{K}}^{-1}\mathbf{H}_{\mathcal{I}^{+}\mathcal{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 (81), they can be extracted from (29), by

$\bm{\varepsilon}=\begin{pmatrix}{}_{K}\bm{\gamma}\\\
{}_{K}\bm{\kappa}_{IK}\end{pmatrix}=j_{{S\\!E}(3)}^{-1}\left(\mathbf{H}_{\mathcal{K}_{0}\mathcal{K}}^{-1}\mathbf{H}_{\mathcal{K}_{0}\mathcal{K},\xi}\right)\frac{1}{J}\,.$
(31)

Inserting the relative interpolation (24) together with the definitions (84)
and suppressing for a while the subscripts indicating ${S\\!E}(3)$, the
discretized strain measures are

$\displaystyle\bm{\varepsilon}(\xi,\cdot)$
$\displaystyle=j^{-1}\circ\left(\operatorname{Exp}(\xi\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}})^{-1}\frac{\mathrm{d}}{\mathrm{d}\xi}\operatorname{Exp}(\xi\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}})\right)\frac{1}{J}=j^{-1}\circ\left(\exp\big{(}j(\xi\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}})\big{)}^{-1}\frac{\mathrm{d}}{\mathrm{d}\xi}\exp\big{(}j(\xi\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}})\big{)}\right)\frac{1}{J}$
(32)
$\displaystyle=j^{-1}\circ\left(\operatorname{dexp}_{-j(\xi\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}})}\Big{(}\frac{\mathrm{d}}{\mathrm{d}\xi}j(\xi\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}})\Big{)}\right)\frac{1}{J}\,,$

where we have used the left-trivialized differential (63)1. Due to the
linearity of $j$, we have
$\mathrm{d}/\mathrm{d}\xi\big{(}j(\xi\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}})\big{)}=j(\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}})$.
Since $j(\xi\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}})$ and
$j(\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}})$ commute, according to the
comment below (67), the discretized strain measures simplify further to

$\bm{\varepsilon}(\xi,\mathbf{q})=\left(j^{-1}\circ\operatorname{dexp}_{-j(\xi\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}}(\mathbf{q}))}\circ
j\right)\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}}(\mathbf{q})\frac{1}{J}=\bm{\theta}_{\mathcal{K}_{0}\mathcal{K}_{1}}(\mathbf{q})\frac{1}{J}\,.$
(33)

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_{\mathrm{e}}=EA$,
$k_{\mathrm{s}}=GA$, $k_{\mathrm{b}_{y}}=k_{\mathrm{b}_{z}}=EI$ and
$k_{\mathrm{t}}=2GI$. Dividing the strain energy density (12) by $w^{4}$, in
the limit of $w\to 0$, the scaled axial $k_{\mathrm{e}}/w^{4}$ and shear
stiffnesses $k_{\mathrm{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}\bm{\gamma}=(1,\ 0,\ 0)$; inextensibility
$\|{}_{K}\bm{\gamma}\|-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 $\mathbb{R}^{3}\times S\\!O(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
${S\\!E}(3)$-interpolation can exactly represent constant strain measures for
all the individual contributions, can thus satisfy ${}_{K}\bm{\gamma}=(1,\ 0,\
0)$, and is hence intrinsically relieved from both shear as well as membrane
locking. Obviously, the ${S\\!E}(3)$-interpolation can represent the quarter
circle exactly, see Figure 2 (a).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: Quarter circle in the
$\mathbf{e}_{x}^{I}$-$\mathbf{e}_{z}^{I}$-plane defined by $\bm{\psi}_{0}=(0,\
-\pi/2,\ 0)$, $\bm{\psi}_{1}=(0,\ 0,\ 0)$, ${}_{I}\mathbf{r}_{OP_{0}}=(0,\ 0,\
0)$ and ${}_{I}\mathbf{r}_{OP_{1}}=(1,\ 0,\ 1)$. Comparison between
${S\\!E}(3)$-interpolation and $\mathbb{R}^{3}\times S\\!O(3)$-interpolation.
(a) - (c) Visualization and strain measures of the two-node
${S\\!E}(3)$-element. (d) - (f) Visualization and strain measures of the two-
node $\mathbb{R}^{3}\times S\\!O(3)$-element.

###  3.4 Kinematics of element nodes

The concepts of the previous paragraphs can be synthesized in a rod finite
element with a piecewise two-node ${S\\!E}(3)$-interpolation that results in
constant strains within an element. Let the rod be discretized by
$n_{\mathrm{el}}$ two-node elements. The nodal Euclidean transformation
matrices $\mathbf{H}_{\mathcal{I}\mathcal{K}_{i}}$,
$i=0,\dots,n_{\mathrm{el}}$ are parameterized in terms of the time-dependent
nodal generalized coordinates
$\mathbf{q}_{i}(t)=\big{(}{}_{I}\mathbf{r}_{OP_{i}}(t),\bm{\psi}_{i}(t)\big{)}\in\mathbb{R}^{6}$
given by the nodal centerline points
${}_{I}\mathbf{r}_{OP_{i}}(t)\in\mathbb{R}^{3}$ and the nodal total rotation
vectors $\bm{\psi}_{i}(t)\in\mathbb{R}^{3}$. Using the exponential map of
$S\\!O(3)$, see comment before (22), the nodal parametrization is

$\mathbf{H}_{\mathcal{I}\mathcal{K}_{i}}(\mathbf{q}_{i})=\begin{pmatrix}\operatorname{Exp}_{S\\!O(3)}(\bm{\psi}_{i})&{}_{I}\mathbf{r}_{OP_{i}}\\\
\mathbf{0}_{3\times 1}&1\end{pmatrix}\,.$ (34)

The rod’s parameter space $\mathcal{J}$ is dived into $n_{\mathrm{el}}$
linearly spaced element intervals $\mathcal{J}^{e}=[\xi^{e},\xi^{e+1})$ via
$J=\bigcup_{e=0}^{n_{\mathrm{el}}-1}J^{e}$. This means a totality of
$n_{\mathrm{el}}+1$ nodes and $n_{\mathbf{q}}=6(n_{\mathrm{el}}+1)$ unknowns.
The nodal quantities can be assembled in the global tuple of generalized
coordinates
$\mathbf{q}(t)=\big{(}\mathbf{q}_{0}(t),\mathbf{q}_{1}(t),\dots,\mathbf{q}_{N-1}(t)\big{)}\in\mathbb{R}^{n_{\mathbf{q}}}$.
Introducing an appropriate Boolean connectivity matrix
$\mathbf{C}_{e}\in\mathbb{R}^{12\times n_{\mathbf{q}}}$, the element
generalized coordinates
$\mathbf{q}^{e}(t)=\big{(}{}_{I}\mathbf{r}_{OP_{e}}(t),\bm{\psi}_{e}(t),{}_{I}\mathbf{r}_{OP_{e+1}}(t),\bm{\psi}_{e+1}(t)\big{)}\in\mathbb{R}^{12}$
can be extracted from the global generalized coordinates $\mathbf{q}$ via
$\mathbf{q}^{e}=\mathbf{C}_{e}\mathbf{q}$. By introducing an additional set of
Boolean connectivity matrices
$\mathbf{C}_{\mathbf{r},0},\mathbf{C}_{\mathbf{r},1},\mathbf{C}_{\bm{\psi},0},\mathbf{C}_{\bm{\psi},1}\in\mathbb{R}^{3\times
12}$, the centerline points ${}_{I}\mathbf{r}_{OP_{e}}$,
${}_{I}\mathbf{r}_{OP_{e+1}}$ and total rotation vectors $\bm{\psi}_{IK_{e}}$,
$\bm{\psi}_{IK_{e+1}}$ can be extracted from the element generalized
coordinates $\mathbf{q}^{e}$ via
${}_{I}\mathbf{r}_{OP_{e}}=\mathbf{C}_{\mathbf{r},0}\mathbf{q}^{e}$,
${}_{I}\mathbf{r}_{OP_{e+1}}=\mathbf{C}_{\mathbf{r},1}\mathbf{q}^{e}$,
$\bm{\psi}_{IK_{e}}=\mathbf{C}_{\bm{\psi},0}\mathbf{q}^{e}$ and
$\bm{\psi}_{IK_{e+1}}=\mathbf{C}_{\bm{\psi},1}\mathbf{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.

###  3.5 Element-wise ${S\\!E}(3)$-interpolation

For a given element interval $\mathcal{J}^{e}=[\xi^{e},\xi^{e+1})$ the
corresponding linear Lagrange basis functions $N_{0}^{e}$, $N_{1}^{e}$ and
their first derivatives $N^{e}_{0,\xi}$, $N^{e}_{1,\xi}$ are defined by

$N^{e}_{0}(\xi)=\frac{\xi^{e+1}-\xi}{\xi^{e+1}-\xi^{e}}\,,\quad
N^{e}_{1}(\xi)=\frac{\xi-\xi^{e}}{\xi^{e+1}-\xi^{e}}\,,\quad
N^{e}_{0,\xi}(\xi)=\frac{-1}{\xi^{e+1}-\xi^{e}}\,,\quad
N^{e}_{1,\xi}(\xi)=\frac{1}{\xi^{e+1}-\xi^{e}}\,.$ (35)

Further, we introduce the characteristic function
$\chi_{\mathcal{J}^{e}}\colon\mathcal{J}\to\\{0,1\\}$, being one for
$\xi\in\mathcal{J}^{e}=[\xi^{e},\xi^{e+1})$ and zero elsewhere. By that, we
can extend the global two-node ${S\\!E}(3)$-interpolation (25) to a piecewise
two-node ${S\\!E}(3)$-interpolation given by

$\displaystyle\mathbf{H}_{\mathcal{I}\mathcal{K}}(\xi,\mathbf{q})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\chi_{\mathcal{J}^{e}}(\xi)\mathbf{H}_{\mathcal{I}\mathcal{K}_{e}}(\mathbf{q})\operatorname{Exp}_{{S\\!E}(3)}\left(N^{e}_{1}(\xi)\,\bm{\theta}_{\mathcal{K}_{e}\mathcal{K}_{e+1}}(\mathbf{q})\right)\,,$
(36) $\displaystyle\bm{\theta}_{\mathcal{K}_{e}\mathcal{K}_{e+1}}(\mathbf{q})$
$\displaystyle=\operatorname{Log}_{{S\\!E}(3)}\left[\mathbf{H}^{-1}_{\mathcal{I}\mathcal{K}_{e}}(\mathbf{q})\,\mathbf{H}_{\mathcal{I}\mathcal{K}_{e+1}}(\mathbf{q})\right]\,.$

For symmetry reasons, cf. Crisfield and Jelenić[13], it might be advantageous
to use a modified ${S\\!E}(3)$-interpolation given by

$\displaystyle\mathbf{H}_{\mathcal{I}\mathcal{K}}(\xi,\mathbf{q})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\chi_{\mathcal{J}^{e}}(\xi)\mathbf{H}_{\mathcal{I}\mathcal{R}_{e}}(\mathbf{q})\operatorname{Exp}_{{S\\!E}(3)}\left(N^{e}_{0}(\xi)\,\bm{\theta}_{\mathcal{R}_{e}\mathcal{K}_{e}}(\mathbf{q})+N^{e}_{1}(\xi)\,\bm{\theta}_{\mathcal{R}_{e}\mathcal{K}_{e+1}}(\mathbf{q})\right)\,,$
(37) $\displaystyle\mathbf{H}_{\mathcal{I}\mathcal{R}_{e}}(\mathbf{q})$
$\displaystyle=\mathbf{H}_{\mathcal{I}\mathcal{K}_{e}}(\mathbf{q})\operatorname{Exp}_{{S\\!E}(3)}\left(\tfrac{1}{2}\operatorname{Log}_{{S\\!E}(3)}\left[\mathbf{H}^{-1}_{\mathcal{I}\mathcal{K}_{e}}(\mathbf{q})\,\mathbf{H}_{\mathcal{I}\mathcal{K}_{e+1}}(\mathbf{q})\right]\right)\,,$
$\displaystyle\bm{\theta}_{\mathcal{R}_{e}\mathcal{K}_{e}}(\mathbf{q})$
$\displaystyle=\operatorname{Log}_{{S\\!E}(3)}\left[\mathbf{H}^{-1}_{\mathcal{I}\mathcal{R}_{e}}(\mathbf{q})\,\mathbf{H}_{\mathcal{I}\mathcal{K}_{e}}(\mathbf{q})\right]\,,\quad\bm{\theta}_{\mathcal{R}_{e}\mathcal{K}_{e+1}}(\mathbf{q})=\operatorname{Log}_{{S\\!E}(3)}\left[\mathbf{H}^{-1}_{\mathcal{I}\mathcal{R}_{e}}(\mathbf{q})\,\mathbf{H}_{\mathcal{I}\mathcal{K}_{e+1}}(\mathbf{q})\right]\,.$

Since this strategy requires twice as many ${S\\!E}(3)$-exponential
evaluations and three times as many ${S\\!E}(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

$\bm{\varepsilon}(\xi,\mathbf{q})=\sum_{e=0}^{n_{\mathrm{el}}-1}\chi_{\mathcal{J}^{e}}(\xi)\frac{\bm{\theta}_{\mathcal{K}_{e}\mathcal{K}_{e+1}}(\mathbf{q})}{\xi^{e+1}-\xi^{e}}\frac{1}{J}\,.$
(38)

Using the same reasoning following (33) and since the piecewise two-node
${S\\!E}(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.

###  3.6 Petrov–Galerkin projection

Following a Petrov–Galerkin projection[15], we introduce the nodal virtual
displacements
$\delta\mathbf{s}_{i}(t)=\big{(}{}_{I}\delta\mathbf{r}_{P_{i}}(t),{}_{K_{i}}\delta\bm{\phi}_{IK_{i}}(t)\big{)}\allowbreak\in\mathbb{R}^{6}$
given by the nodal centerline variation
${}_{I}\delta\mathbf{r}_{P_{i}}(t)\in\mathbb{R}^{3}$ and the nodal virtual
rotation ${}_{K_{i}}\delta\bm{\phi}_{IK_{i}}(t)\in\mathbb{R}^{3}$. They are
related to the variations of the nodal generalized coordinates by
${}_{I}\delta\mathbf{r}_{P_{i}}=\delta({}_{I}\mathbf{r}_{OP_{i}})$ and
${}_{K_{i}}\delta\bm{\phi}_{IK_{i}}=\mathbf{T}_{S\\!O(3)}(\bm{\psi}_{i})\delta\bm{\psi}_{i}$,
see Equation (49) of Cardona and Geradin [11]. Led by the observation that
$\overline{\mathbf{M}}$ in (16) is symmetric, it is advisable to introduce the
nodal minimal velocities
$\mathbf{u}_{i}(t)=\big{(}{}_{I}\mathbf{v}_{P_{i}}(t),{}_{K_{i}}\bm{\omega}_{IK_{i}}(t)\big{)}\in\mathbb{R}^{6}$
given by the nodal centerline velocity
${}_{I}\mathbf{v}_{P_{i}}(t)\in\mathbb{R}^{3}$ and the nodal angular velocity
${}_{K_{i}}\bm{\omega}_{IK_{i}}(t)\in\mathbb{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
equation

$\dot{\mathbf{q}}_{i}=\begin{pmatrix}{}_{I}\dot{\mathbf{r}}_{OP_{i}}\\\
\dot{\bm{\psi}}_{i}\end{pmatrix}=\begin{pmatrix}\mathbf{1}_{3\times
3}&\mathbf{0}_{3\times 3}\\\ \mathbf{0}_{3\times
3}&\mathbf{T}^{-1}_{S\\!O(3)}(\bm{\psi}_{i})\end{pmatrix}\begin{pmatrix}{}_{I}\mathbf{v}_{P_{i}}\\\
{}_{K_{i}}\bm{\omega}_{IK_{i}}\end{pmatrix}=\mathbf{B}_{i}(\mathbf{q}_{i})\mathbf{u}_{i}\,.$
(39)

Consequently, during a subsequent Galerkin projection, the symmetry of
$\overline{\mathbf{M}}$ is preserved and results in a symmetric mass matrix,
see (44). Since the inverse tangent map $\mathbf{T}^{-1}_{S\\!O(3)}$ used in
(39) exhibits singularities for $\|\bm{\psi}\|=k2\pi$ with $k=0,1,2,\dots$, we
apply the following strategy to avoid them. For $k=0$, we use the first order
approximation $\mathbf{T}^{-1}_{S\\!O(3)}(\bm{\omega})=\mathbf{1}_{3\times
3}+\frac{1}{2}\widetilde{\bm{\omega}}$. 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

$\bm{\psi}=\begin{cases}\bm{\psi}\,,&\|\bm{\psi}\|<\pi\,,\\\
\bm{\psi}^{C}=\big{(}1-2\pi/\|\bm{\psi}\|\big{)}\bm{\psi}\,,&\|\bm{\psi}\|\geq\pi\,,\end{cases}$
(40)

with
$\operatorname{Exp}_{S\\!O(3)}(\bm{\psi})=\operatorname{Exp}_{S\\!O(3)}(\bm{\psi}^{C})$,
i.e., there is no difference whether the nodal transformation matrix
$\mathbf{A}_{IK}$ is described by the rotation vector $\bm{\psi}$ or by its
complement $\bm{\psi}^{C}=\big{(}1-2\pi/\|\bm{\psi}\|\big{)}\bm{\psi}$. Thus,
for reasonable time steps (Newton increments) a both minimal and singularity
free parametrization of $S\\!O(3)$ is obtained without changing the rod’s
virtual work formulation since all relevant nodal quantities are interpolated
relatively111It should be mentioned that the proposed complement rotation
vector formalism is restricted to numerical single step methods, e.g.
Runge–Kutta family, generalized-$\alpha$ methods, etc., but cannot be used
straightforwardly for multi step algorithms like BDF. Alternatively, another
nodal rotation parametrization can be chosen, e.g., unit quaternions. Since in
numerics, we cannot deal with unit quaternions, a general quaternion is
constrained to be of unit length. While in statics these constraint equations
have to be taken into account, during a numerical time integration they can be
satisfied by applying an appropriate projection after each successful time
step. Likewise, a modified kinematic differential equation can be introduced,
see Section 6.9.6 of Egeland and Gravdahl[33]..

In analogy to the generalized coordinates, the nodal virtual displacements and
minimal velocities are assembled in the global tuple of virtual displacements
$\delta\mathbf{s}(t)=\big{(}\delta\mathbf{s}_{0}(t),\delta\mathbf{s}_{1}(t),\dots,\delta\mathbf{s}_{N-1}(t)\big{)}\in\mathbb{R}^{n_{\mathbf{q}}}$
and global tuple of minimal velocities
$\mathbf{u}(t)=\big{(}\mathbf{u}_{0}(t),\mathbf{u}_{1}(t),\dots,\mathbf{u}_{N-1}(t)\big{)}\in\mathbb{R}^{n_{\mathbf{q}}}$.
Using again the boolean connectivity matrix $\mathbf{C}_{e}$, we can extract
the element virtual displacements
$\delta\mathbf{s}^{e}(t)=\big{(}{}_{I}\delta\mathbf{r}_{P_{e}}(t),{}_{K_{e}}\\!\delta\bm{\phi}_{IK_{e}}(t),{}_{I}\delta\mathbf{r}_{P_{e+1}}(t),{}_{K_{e+1}}\\!\delta\bm{\phi}_{IK_{e+1}}(t)\big{)}\in\mathbb{R}^{12}$
and the element minimal velocities
$\mathbf{u}^{e}(t)=\big{(}{}_{I}\mathbf{v}_{P_{e}}(t),{}_{K_{e}}\\!\bm{\omega}_{IK_{e}}(t),{}_{I}\mathbf{v}_{P_{e+1}}(t),{}_{K_{e+1}}\\!\bm{\omega}_{IK_{e+1}}(t)\big{)}\in\mathbb{R}^{12}$
from the global quantities in agreement with
$\delta\mathbf{s}^{e}=\mathbf{C}_{e}\delta\mathbf{s}$ and
$\mathbf{u}^{e}=\mathbf{C}_{e}\mathbf{u}$. Further, the centerline variations
${}_{I}\delta\mathbf{r}_{P_{e}}$, ${}_{I}\delta\mathbf{r}_{P_{e+1}}$ and
centerline velocities ${}_{I}\mathbf{v}_{P_{e}}$, ${}_{I}\mathbf{v}_{P_{e+1}}$
can be extracted from the tuple of virtual displacements
$\delta\mathbf{s}^{e}$ and minimal velocities $\mathbf{u}^{e}$ using
${}_{I}\delta\mathbf{r}_{P_{e}}=\mathbf{C}_{\mathbf{r},0}\delta\mathbf{s}^{e}$,
${}_{I}\delta\mathbf{r}_{P_{e+1}}=\mathbf{C}_{\mathbf{r},1}\delta\mathbf{s}^{e}$
and ${}_{I}\mathbf{v}_{P_{e}}=\mathbf{C}_{\mathbf{r},0}\mathbf{u}^{e}$,
${}_{I}\mathbf{v}_{P_{e+1}}=\mathbf{C}_{\mathbf{r},1}\mathbf{u}^{e}$,
respectively. Identical extraction operations can be defined for the nodal
virtual rotations and angular velocities via
$\mathbf{C}_{\bm{\psi},0},\mathbf{C}_{\bm{\psi},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

$\displaystyle{}_{I}\delta\mathbf{r}_{P}(\xi,\delta\mathbf{s})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\chi_{\mathcal{J}^{e}}(\xi)\left(N^{e}_{0}(\xi){}_{I}\delta\mathbf{r}_{P_{e}}+N^{e}_{1}(\xi){}_{I}\delta\mathbf{r}_{P_{e+1}}\right)\,,$
(41) $\displaystyle{}_{I}\mathbf{v}_{P}(\xi,\mathbf{u})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\chi_{\mathcal{J}^{e}}(\xi)\left(N^{e}_{0}(\xi){}_{I}\mathbf{v}_{P_{e}}+N^{e}_{1}(\xi){}_{I}\mathbf{v}_{P_{e+1}}\right)\,,$
$\displaystyle{}_{K}\delta\bm{\phi}_{IK}(\xi,\delta\mathbf{s})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\chi_{\mathcal{J}^{e}}(\xi)\left(N^{e}_{0}(\xi)\,{}_{K_{e}}\\!\delta\bm{\phi}_{IK_{e}}+N^{e}_{1}(\xi)\,{}_{K_{e+1}}\\!\delta\bm{\phi}_{IK_{e+1}}\right)\,,$
$\displaystyle{}_{K}\bm{\omega}_{IK}(\xi,\mathbf{u})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\chi_{\mathcal{J}^{e}}(\xi)\left(N^{e}_{0}(\xi)\,{}_{K_{e}}\\!\bm{\omega}_{IK_{e}}+N^{e}_{1}(\xi)\,{}_{K_{e+1}}\\!\bm{\omega}_{IK_{e+1}}\right)\,.$

Choices like ${}_{K}\delta\mathbf{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].

###  3.7 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

$\displaystyle\delta W^{\mathrm{int}}(\delta\mathbf{s};\mathbf{q})$
$\displaystyle=\delta\mathbf{s}^{\\!\mathrm{T}}\mathbf{f}^{\mathrm{int}}(\mathbf{q})\,,\quad\mathbf{f}^{\mathrm{int}}(\mathbf{q})=\sum_{e=0}^{n_{\mathrm{el}}-1}\mathbf{C}_{e}^{\\!\mathrm{T}}\mathbf{f}^{\mathrm{int}}_{e}(\mathbf{C}_{e}\mathbf{q})\,,$
(42) $\displaystyle\mathbf{f}^{\mathrm{int}}_{e}(\mathbf{q}^{e})$
$\displaystyle=-\int_{\mathcal{J}^{e}}\sum_{i=0}^{1}\Big{\\{}N^{e}_{i,\xi}\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}\mathbf{A}_{IK}\,{}_{K}\mathbf{n}+N^{e}_{i,\xi}\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}{}_{K}\mathbf{m}-N^{e}_{i}\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}\left({}_{K}\bar{\bm{\gamma}}\times{}_{K}\mathbf{n}+{}_{K}\bar{\bm{\kappa}}_{IK}\times{}_{K}\mathbf{m}\right)\Big{\\}}\mathrm{d}\xi\,,$

where we have introduced the internal forces $\mathbf{f}^{\mathrm{int}}$
together with their element contributions $\mathbf{f}^{\mathrm{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

$\displaystyle\delta W^{\mathrm{ext}}(\delta\mathbf{s};\mathbf{q})$
$\displaystyle=\delta\mathbf{s}^{\\!\mathrm{T}}\mathbf{f}^{\mathrm{ext}}(\mathbf{q})\,,$
(43) $\displaystyle\mathbf{f}^{\mathrm{ext}}(\mathbf{q})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\mathbf{C}_{e}^{\\!\mathrm{T}}\mathbf{f}^{\mathrm{ext}}_{e}(\mathbf{C}_{e}\mathbf{q})+\mathbf{C}_{n_{\mathrm{el}}-1}^{\\!\mathrm{T}}\left[\mathbf{C}_{\mathbf{r},1}^{\\!\mathrm{T}}{}_{I}\mathbf{b}_{1}+\mathbf{C}_{\bm{\phi},1}^{\\!\mathrm{T}}{}_{K}\mathbf{c}_{1}\right]_{\xi=1}-\mathbf{C}_{0}^{\\!\mathrm{T}}\left[\mathbf{C}_{\mathbf{r},0}^{\\!\mathrm{T}}{}_{I}\mathbf{b}_{0}+\mathbf{C}_{\bm{\phi},0}^{\\!\mathrm{T}}{}_{K}\mathbf{c}_{0}\right]_{\xi=0}\,,$
$\displaystyle\mathbf{f}^{\mathrm{ext}}_{e}(\mathbf{q}^{e})$
$\displaystyle=\int_{\mathcal{J}^{e}}\sum_{i=0}^{1}\Big{\\{}N^{e}_{i}\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}{}_{I}\mathbf{b}+N^{e}_{i}\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}{}_{K}\mathbf{c}\Big{\\}}J\mathrm{d}\xi\,,$

where we have introduced the external forces $\mathbf{f}^{\mathrm{ext}}$ and
their element contributions $\mathbf{f}^{\mathrm{ext}}_{e}$. Finally, the
discretization of the inertial virtual work contributions (15) is given

$\displaystyle\delta W^{\mathrm{dyn}}(\delta\mathbf{s};\mathbf{q},\mathbf{u})$
$\displaystyle=-\delta\mathbf{s}^{\\!\mathrm{T}}\left\\{\mathbf{M}(\mathbf{q})\dot{\mathbf{u}}+\mathbf{f}^{\mathrm{gyr}}(\mathbf{q},\mathbf{u})\right\\}\,,$
(44) $\displaystyle\mathbf{M}(\mathbf{q})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\mathbf{C}_{e}^{\\!\mathrm{T}}\mathbf{M}_{e}(\mathbf{C}_{e}\mathbf{q})\mathbf{C}_{e}\,,\quad\mathbf{f}^{\mathrm{gyr}}(\mathbf{q},\mathbf{u})=\sum_{e=0}^{n_{\mathrm{el}}-1}\mathbf{C}_{e}^{\\!\mathrm{T}}\mathbf{f}^{\mathrm{gyr}}(\mathbf{C}_{e}\mathbf{q},\mathbf{C}_{e}\mathbf{u})\,,$
$\displaystyle\mathbf{M}_{e}(\mathbf{q}^{e})$
$\displaystyle=\begin{multlined}\int_{\mathcal{J}^{e}}\sum_{i=0}^{1}\sum_{k=0}^{1}N^{e}_{i}N^{e}_{k}\Big{\\{}\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}A_{\rho_{0}}\mathbf{1}_{3\times
3}\mathbf{C}_{\mathbf{r},k}+\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}\mathbf{A}_{IK}{}_{K}\mathbf{S}_{\rho_{0}}^{\\!\mathrm{T}}\mathbf{C}_{\bm{\psi},k}\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}{}_{K}\mathbf{S}_{\rho_{0}}\mathbf{A}_{IK}^{\\!\mathrm{T}}\mathbf{C}_{\mathbf{r},k}+\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}{}_{K}\mathbf{I}_{\rho_{0}}\mathbf{C}_{\bm{\psi},k}\Big{\\}}J\mathrm{d}\xi\,,\end{multlined}\int_{\mathcal{J}^{e}}\sum_{i=0}^{1}\sum_{k=0}^{1}N^{e}_{i}N^{e}_{k}\Big{\\{}\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}A_{\rho_{0}}\mathbf{1}_{3\times
3}\mathbf{C}_{\mathbf{r},k}+\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}\mathbf{A}_{IK}{}_{K}\mathbf{S}_{\rho_{0}}^{\\!\mathrm{T}}\mathbf{C}_{\bm{\psi},k}\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}{}_{K}\mathbf{S}_{\rho_{0}}\mathbf{A}_{IK}^{\\!\mathrm{T}}\mathbf{C}_{\mathbf{r},k}+\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}{}_{K}\mathbf{I}_{\rho_{0}}\mathbf{C}_{\bm{\psi},k}\Big{\\}}J\mathrm{d}\xi\,,$
$\displaystyle\mathbf{f}^{\mathrm{gyr}}_{e}(\mathbf{q}^{e},\mathbf{u}^{e})$
$\displaystyle=\int_{\mathcal{J}^{e}}\sum_{i=0}^{1}N^{e}_{i}\Big{\\{}\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}\mathbf{A}_{IK}{}_{K}\widetilde{\bm{\omega}}_{IK}{}_{K}\mathbf{S}_{\rho_{0}}^{\\!\mathrm{T}}{}_{K}\bm{\omega}_{IK}+\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}{}_{K}\widetilde{\bm{\omega}}_{IK}{}_{K}\mathbf{I}_{\rho_{0}}{}_{K}\bm{\omega}_{IK}\Big{\\}}J\mathrm{d}\xi\,,$

where we have introduced the symmetric mass matrix $\mathbf{M}$ and the
gyroscopic forces $\mathbf{f}^{\mathrm{gyr}}$ together with their element
contributions $\mathbf{M}_{e}$ and $\mathbf{f}^{\mathrm{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 ${S\\!E}(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}\delta\bm{\phi}_{IK_{i}}(t)\in\mathbb{R}^{3}$
and nodal angular velocities ${}_{I}\bm{\omega}_{IK_{i}}(t)\in\mathbb{R}^{3}$
expressed in the inertial $I$-basis.

The linearization of the internal forces $\mathbf{f}^{\mathrm{int}}$ with
respect to the generalized coordinates $\mathbf{q}$ is given by

$\displaystyle\frac{\partial\mathbf{f}^{\mathrm{int}}}{\partial\mathbf{q}}(\mathbf{q})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\mathbf{C}_{e}^{\\!\mathrm{T}}\frac{\partial\mathbf{f}^{\mathrm{int}}}{\partial\mathbf{q}^{e}}(\mathbf{C}^{e}\mathbf{q})\,\mathbf{C}^{e}\,,$
(45)
$\displaystyle\frac{\partial\mathbf{f}^{\mathrm{int}}}{\partial\mathbf{q}^{e}}(\mathbf{q}^{e})$
$\displaystyle=-\int_{\mathcal{J}^{e}}\sum_{i=0}^{1}\begin{multlined}\bigg{\\{}N^{e}_{i,\xi}\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}\frac{\partial\left(\mathbf{A}_{IK}\,{}_{K}\mathbf{n}\right)}{\partial\mathbf{q}^{e}}+N^{e}_{i,\xi}\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}\frac{\partial{}_{K}\mathbf{m}}{\partial\mathbf{q}^{e}}\\\
-N^{e}_{i}\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}\left({}_{K}\widetilde{\bar{\bm{\gamma}}}\frac{\partial{}_{K}\mathbf{n}}{\partial\mathbf{q}^{e}}-{}_{K}\widetilde{\mathbf{n}}\frac{\partial{}_{K}\bar{\bm{\gamma}}}{\partial\mathbf{q}^{e}}+{}_{K}\widetilde{\bar{\bm{\kappa}}}_{IK}\frac{\partial{}_{K}\mathbf{m}}{\partial\mathbf{q}^{e}}-{}_{K}\widetilde{\mathbf{m}}\frac{\partial{}_{K}\bar{\bm{\kappa}}_{IK}}{\partial\mathbf{q}^{e}}\right)\bigg{\\}}\mathrm{d}\xi\,,\end{multlined}\bigg{\\{}N^{e}_{i,\xi}\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}\frac{\partial\left(\mathbf{A}_{IK}\,{}_{K}\mathbf{n}\right)}{\partial\mathbf{q}^{e}}+N^{e}_{i,\xi}\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}\frac{\partial{}_{K}\mathbf{m}}{\partial\mathbf{q}^{e}}\\\
-N^{e}_{i}\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}\left({}_{K}\widetilde{\bar{\bm{\gamma}}}\frac{\partial{}_{K}\mathbf{n}}{\partial\mathbf{q}^{e}}-{}_{K}\widetilde{\mathbf{n}}\frac{\partial{}_{K}\bar{\bm{\gamma}}}{\partial\mathbf{q}^{e}}+{}_{K}\widetilde{\bar{\bm{\kappa}}}_{IK}\frac{\partial{}_{K}\mathbf{m}}{\partial\mathbf{q}^{e}}-{}_{K}\widetilde{\mathbf{m}}\frac{\partial{}_{K}\bar{\bm{\kappa}}_{IK}}{\partial\mathbf{q}^{e}}\right)\bigg{\\}}\mathrm{d}\xi\,,$

which requires the derivative of the rotation matrix $\mathbf{A}_{IK}$ and the
derivatives of the internal forces and moments

$\displaystyle\frac{\partial{}_{K}\mathbf{n}}{\partial\mathbf{q}^{e}}=\frac{\partial{}_{K}\mathbf{n}}{\partial{}_{K}\bm{\gamma}}\frac{\partial{}_{K}\bm{\gamma}}{\partial\mathbf{q}^{e}}+\frac{\partial{}_{K}\mathbf{n}}{\partial{}_{K}\bm{\kappa}_{IK}}\frac{\partial{}_{K}\bm{\kappa}_{IK}}{\partial\mathbf{q}^{e}}\,,\quad\frac{\partial{}_{K}\mathbf{m}}{\partial\mathbf{q}^{e}}=\frac{\partial{}_{K}\mathbf{m}}{\partial{}_{K}\bm{\gamma}}\frac{\partial{}_{K}\bm{\gamma}}{\partial\mathbf{q}^{e}}+\frac{\partial{}_{K}\mathbf{m}}{\partial{}_{K}\bm{\kappa}_{IK}}\frac{\partial{}_{K}\bm{\kappa}_{IK}}{\partial\mathbf{q}^{e}}\,.$
(46)

Both are depending on the derivatives of the strain measures
${}_{K}\bm{\gamma}$ and ${}_{K}\bm{\kappa}_{IK}$ introduced in (38). These
derivatives together with the derivative of the interpolation formula (36),
i.e., the position vector ${}_{I}\mathbf{r}_{OP}$ and the transformation
matrix $\mathbf{A}_{IK}$ are given in Appendix C. As already noted, the
quantity $\mathbf{S}_{\rho_{0}}$ vanishes if the centerline points
${}_{I}\mathbf{r}_{OP}$ coincide with the cross-sections’ center of mass. For
this case the gyroscopic forces $\mathbf{f}^{\mathrm{gyr}}$ and the term
$\mathbf{M}\dot{\mathbf{u}}$ get independent of the generalized coordinates
$\mathbf{q}$ which is why we do not consider their linearizations here. The
linearization of the external forces is omitted as well since a possible
$\mathbf{q}$ dependence is defined by the specific form of the external forces
and moments only.

Arising element integrals of the form
$\int_{\mathcal{J}^{e}}f(\xi)\mathrm{d}\xi$ 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.

###  3.8 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.,

$\delta W^{\mathrm{tot}}=\delta W^{\mathrm{int}}+\delta
W^{\mathrm{ext}}+\delta
W^{\mathrm{dyn}}\stackrel{{\scriptstyle!}}{{=}}0\quad\forall\delta\mathbf{s}(t),\forall
t\,.$ (47)

Thus, the equations of motion

$\dot{\mathbf{u}}=\mathbf{M}^{-1}(\mathbf{q})\left(\mathbf{f}^{\mathrm{int}}(\mathbf{q})+\mathbf{f}^{\mathrm{ext}}(\mathbf{q})-\mathbf{f}^{\mathrm{gyr}}(\mathbf{q},\mathbf{u})\right)$
(48)

have to be fulfilled for each instant of time $t$. Further, the nodal minimal
velocities $\mathbf{u}_{i}$ are related to the time derivatives of the nodal
generalized coordinates $\dot{\mathbf{q}}_{i}$ via the nodal kinematic
differential equation (39), which can be assembled to a global kinematic
differential equation of the form
$\dot{\mathbf{q}}=\mathbf{B}(\mathbf{q})\mathbf{u}$. Depending on the specific
application, either the system of ordinary differential equations

$\displaystyle\dot{\mathbf{q}}$
$\displaystyle=\mathbf{B}(\mathbf{q})\mathbf{u}\,,$ (49)
$\displaystyle\dot{\mathbf{u}}$
$\displaystyle=\mathbf{M}(\mathbf{q})^{-1}\left(\mathbf{f}^{\mathrm{int}}(\mathbf{q})+\mathbf{f}^{\mathrm{ext}}(\mathbf{q})-\mathbf{f}^{\mathrm{gyr}}(\mathbf{q},\mathbf{u})\right)\,,$

or the nonlinear generalized force equilibrium

$\mathbf{f}^{\mathrm{int}}(\mathbf{q})+\mathbf{f}^{\mathrm{ext}}(\mathbf{q})=\mathbf{0}$
(50)

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 ${S\\!E}(3)$ rod formulations
require highly customized Lie group solvers[42, 20]. In order to apply well-
established methods from structural dynamics like the Newmark-$\beta$[43] or
the generalized-$\alpha$ 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-$\alpha$ formulation for
first-order differential equations[45] can be used without any modifications.
The nonlinear generalized force equilibrium (50) 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.

##  4 Numerical experiments

This section demonstrates the power of the developed ${S\\!E}(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
${S\\!E}(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-$\alpha$ 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.

###  4.1 Cantilever experiment

(a)

(b)

(c)

(d)

Figure 3: Visualization of deformed configuration of the _cantilever
experiment_ for (a) $\rho=${10}^{1}\text{\,}\frac{}{}$$, (b)
$\rho=${10}^{2}\text{\,}\frac{}{}$$, (c) $\rho=${10}^{3}\text{\,}\frac{}{}$$
and (d) $\rho=${10}^{4}\text{\,}\frac{}{}$$ for five elements of the proposed
${S\\!E}(3)$-rod element.

We consider an initially straight cantilever rod of length
$L=${10}^{3}\text{\,}\frac{}{}$$ with a quadratic cross-section of width $w$
subjected to a tip moment ${}_{K}\mathbf{c}_{1}=(0,\ 0,\ 0.5\pi
k_{\mathrm{b}}/L)$ and an out-of-plane tip load
${}_{I}\mathbf{b}_{1}=\mathbf{A}_{IK}\,(0,\ 0,\ 0.5\pi k_{\mathrm{b}}/L^{2})$.
In order to show the absence of locking, different slenderness ratios
$\rho=L/w\in\\{${10}^{1}\text{\,}\frac{}{}$,${10}^{2}\text{\,}\frac{}{}$,${10}^{3}\text{\,}\frac{}{}$,${10}^{4}\text{\,}\frac{}{}$\\}$
are considered, i.e., widths $w\in\\{${10}^{2}\text{\,}\frac{}{}$,\
${10}^{1}\text{\,}\frac{}{}$,\ $undef\text{\,}\frac{}{}$,\
${10}^{-1}\text{\,}\frac{}{}$\\}$. ranging from
$w=${10}^{-1}\text{\,}\frac{}{}$$ to $w=${10}^{2}\text{\,}\frac{}{}$$.
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_{\mathrm{e}}=EA$, shear
stiffness $k_{\mathrm{s}}=GA$, bending stiffnesses
$k_{\mathrm{b}}=k_{\mathrm{b}_{y}}=k_{\mathrm{b}_{z}}=EI$ and torsional
stiffness $k_{\mathrm{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 ${S\\!E}(3)$-element and for $64\text{\,}\frac{}{}$
${S\\!E}(3)$-elements are compared to a reference solution found by a finite
element implementation similar to Eugster et al.[21] discretized by
$256\text{\,}\frac{}{}$ 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

$e^{k}_{\mathbf{r}}=\frac{1}{k}\sqrt{\sum_{i=0}^{k-1}{}_{I}\Delta\mathbf{r}_{P}^{\\!\mathrm{T}}(\xi_{i})\,{}_{I}\Delta\mathbf{r}_{P}(\xi_{i})}\,,\quad{}_{I}\Delta\mathbf{r}_{P}(\xi_{i})={}_{I}\mathbf{r}_{OP}(\xi_{i})-{}_{I}\mathbf{r}_{OP}^{*}(\xi_{i})\,,\quad\xi_{i}=\frac{i}{k-1}\,,$
(51)

where ${}_{I}\mathbf{r}_{OP}^{*}$ denotes the centerline points of the
reference solution and ${}_{I}\mathbf{r}_{OP}$ the same quantity of the rod in
comparison. In order to investigate the error in orientations, we use the
measure[24]

$e^{k}_{\bm{\psi}}=\frac{1}{k}\sqrt{\sum_{i=0}^{k-1}\Delta\bm{\psi}^{\\!\mathrm{T}}(\xi_{i})\Delta\bm{\psi}(\xi_{i})}\,,\quad\Delta\bm{\psi}(\xi_{i})=\operatorname{Log}_{S\\!O(3)}\big{(}\mathbf{A}_{IK}^{\\!\mathrm{T}}(\xi_{i})\mathbf{A}_{IK}^{*}(\xi_{i})\big{)}\,,\quad\xi_{i}=\frac{i}{k-1}\,.$
(52)

Again, $\mathbf{A}_{IK}^{*}$ denotes the orthogonal transformation matrix of
the reference solution and $\mathbf{A}_{IK}$ the same quantity of the rod in
comparison. Since the relative rotation vector $\Delta\bm{\psi}$ 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\text{\,}\frac{}{}$ 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.

Table 1: Experimental parameters and errors of the _cantilever experiment_. slenderness | tolerance | $n_{\mathrm{el}}=1$ | $n_{\mathrm{el}}=64$
---|---|---|---
$e^{100}_{\mathbf{r}}$ | $e^{100}_{\bm{\psi}}$ | $e^{100}_{\mathbf{r}}$ | $e^{100}_{\bm{\psi}}$
${10}^{1}\text{\,}\frac{}{}$ | ${10}^{-8}\text{\,}\frac{}{}$ | $6.789\text{\,}\frac{}{}$ | $1.628\text{\,}\frac{}{}$ | $1.396\text{\times}{10}^{-03}\text{\,}\frac{}{}$ | $3.355\text{\times}{10}^{-04}\text{\,}\frac{}{}$
${10}^{2}\text{\,}\frac{}{}$ | ${10}^{-9}\text{\,}\frac{}{}$ | $6.792\text{\,}\frac{}{}$ | $1.629\text{\,}\frac{}{}$ | $1.397\text{\times}{10}^{-03}\text{\,}\frac{}{}$ | $3.357\text{\times}{10}^{-04}\text{\,}\frac{}{}$
${10}^{3}\text{\,}\frac{}{}$ | ${10}^{-10}\text{\,}\frac{}{}$ | $6.792\text{\,}\frac{}{}$ | $1.629\text{\,}\frac{}{}$ | $1.397\text{\times}{10}^{-03}\text{\,}\frac{}{}$ | $3.357\text{\times}{10}^{-04}\text{\,}\frac{}{}$
${10}^{4}\text{\,}\frac{}{}$ | ${10}^{-11}\text{\,}\frac{}{}$ | $6.792\text{\,}\frac{}{}$ | $1.629\text{\,}\frac{}{}$ | $1.397\text{\times}{10}^{-03}\text{\,}\frac{}{}$ | $3.356\text{\times}{10}^{-04}\text{\,}\frac{}{}$

Due to the applied out-of-plane force, the resulting deformation is
inhomogenous and cannot be described adequately by a small number of
${S\\!E}(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 $\rho=${10}^{3}\text{\,}\frac{}{}$$, an absolute tolerance
of $atol=${10}^{-10}\text{\,}\frac{}{}$$ and $20\text{\,}\frac{}{}$ 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.

(a)

(b)

(c)

Figure 4: (a), (b) Strain measures of the solution found by four
${S\\!E}(3)$-elements (solid) and those of the reference solution computed
with $n_{\mathrm{el}}=256$ elements (dashed). (c) Convergence behavior of the
${S\\!E}(3)$-formulation.

###  4.2 Helix experiment

(a)

(b)

(c)

(d)

Figure 5: Visualization of deformed configuration of the _helix experiment_
for (a) $\rho=${10}^{1}\text{\,}\frac{}{}$$, (b)
$\rho=${10}^{2}\text{\,}\frac{}{}$$, (c) $\rho=${10}^{3}\text{\,}\frac{}{}$$
and (d) $\rho=${10}^{4}\text{\,}\frac{}{}$$ for five elements of the proposed
${S\\!E}(3)$-rod element.

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\pi n)\geq 1$, the
centerline of such a deformed rod is described by

${}_{I}\mathbf{r}_{OP}^{*}(\xi)=R_{0}\begin{pmatrix}\sin\alpha(\xi)\\\
-\cos\alpha(\xi)\\\ c\alpha(\xi)\end{pmatrix}\,,\quad\alpha(\xi)=2\pi n\xi\,.$
(53)

Depending on the slenderness ratio $\rho$, the rod has a circular cross-
section of diameter $d=L/\rho$, radius $r=d/2$, area $A=\pi r^{2}$ and second
moment of area $I=\tfrac{\pi}{4}r^{4}$. Using the Young’s and shear moduli
$E=1$ and $G=0.5$, respectively, the elastic constants are given by
$k_{\mathrm{e}}=EA$, $k_{\mathrm{s}}=GA$,
$k_{\mathrm{b}}=k_{\mathrm{b}_{y}}=k_{\mathrm{b}_{z}}=EI$ and
$k_{\mathrm{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 5.2 of Ogden[48]. Evaluating the derivative of (53),
the tangent vector of the curve is

${}_{I}\mathbf{r}_{OP,\xi}^{*}(\xi)=R_{0}\alpha_{,\xi}\begin{pmatrix}\cos\alpha(\xi)\\\
\sin\alpha(\xi)\\\ c\end{pmatrix}\,,\quad\alpha_{,\xi}(\xi)=2\pi n\,.$ (54)

The rod should not be elongated during the deformation, i.e.,
$J(\xi)=\|{}_{I}\mathbf{r}_{OP,\xi}^{*}(\xi)\|=L$, from which the rod’s total
length $L=\sqrt{1+c^{2}}R_{0}2\pi n$ follows. Further, the Serret–Frenet frame
of the deformed curve is given by

$\displaystyle{}_{I}\mathbf{e}_{x}^{K}(\xi)=\frac{1}{\sqrt{1+c^{2}}}\begin{pmatrix}\cos\alpha(\xi)\\\
\sin\alpha(\xi)\\\
c\end{pmatrix}\,,\quad{}_{I}\mathbf{e}_{y}^{K}(\xi)=\begin{pmatrix}-\sin\alpha(\xi)\\\
\cos\alpha(\xi)\\\
0\end{pmatrix}\,,\quad{}_{I}\mathbf{e}_{z}^{K}(\xi)=\frac{1}{\sqrt{1+c^{2}}}\begin{pmatrix}-c\cos\alpha(\xi)\\\
-c\sin\alpha(\xi)\\\ 1\end{pmatrix}\,.$ (55)

Using $(1+c^{2})^{-\frac{1}{2}}=2\pi nR_{0}/L=\alpha_{,\xi}R_{0}/L$ together
with the derivatives of (55), the rod’s strain measures (8) compute as

${}_{K}\bm{\gamma}=(1,\ 0,\ 0)\ ,\quad{}_{K}\bm{\kappa}_{IK}=(c,\ 0,\
1)R_{0}\alpha^{2}_{,\xi}/L^{2}\,.$ (56)

For the given strain energy density (12), evaluating the constitutive
equations (10) results in the force boundary conditions

${}_{I}\mathbf{b}_{1}=(0,\ 0,\
0)\,,\quad{}_{K}\mathbf{c}_{1}=(k_{\mathrm{t}}c,\ 0,\
k_{\mathrm{b}})R_{0}\alpha^{2}_{,\xi}/L^{2}$ (57)

that induce a pure torsion and bending deformation. Inserting (57) into the
differential equation of the nonlinear Cosserat rod[8] results in the
requirement $k_{\mathrm{t}}=k_{\mathrm{b}}$. This condition is satisfied for
the given problem setup. For a numerical simulation, the rod has to be clamped
at ${}_{I}\mathbf{r}_{OP}|_{\xi=0}={}_{I}\mathbf{r}_{OP}^{*}|_{\xi=0}$ with an
orientation given by $\mathbf{A}_{IK}|_{\xi=0}=({}_{I}\mathbf{e}_{x}^{K},\
{}_{I}\mathbf{e}_{y}^{K},\ {}_{I}\mathbf{e}_{z}^{K})|_{\xi=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
${S\\!E}(3)$-interpolation can exactly represent a helix, which is a
configuration with constant strains (56). As mentioned in Section 3.2, a
single ${S\\!E}(3)$-element is restricted to a relative rotation vector of
absolute value bounded by $\pi$. Thus four elements can exactly represent the
helix (53) with two coils. Nevertheless, during the Newton iterations locally
the relative rotation vector might exceed an absolute value of $\pi$. Hence,
we use five elements for the numerical investigation.

Table 2: Experimental parameters used for the _helix example_. slenderness | tolerance | force increments
---|---|---
${10}^{1}\text{\,}\frac{}{}$ | ${10}^{-8}\text{\,}\frac{}{}$ | 70
${10}^{2}\text{\,}\frac{}{}$ | ${10}^{-9}\text{\,}\frac{}{}$ | 100
${10}^{3}\text{\,}\frac{}{}$ | ${10}^{-10}\text{\,}\frac{}{}$ | 200
${10}^{4}\text{\,}\frac{}{}$ | ${10}^{-11}\text{\,}\frac{}{}$ | 500

###  4.3 Superimposed rotation of deformed cantilever

The objectivity of the ${S\\!E}(3)$-formulation can be demonstrated
numerically by using the cantilever rod from Section 4.1 with slenderness
ratio $\rho=${10}^{2}\text{\,}\frac{}{}$$, 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
$\bm{\psi}_{0}=(20\pi,\ 0,\ 0)$, i.e., the deformed rod performed $10$ full
rotations around the $\mathbf{e}_{x}^{I}$-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 $\pi$, its complement value is chosen. Thus, the
trajectory of the individual components of the nodal rotation vector exhibit
discontinuities.

(a)

(b)

Figure 6: (a) Snapshots of deformed configurations of the rotated cantilever.
(b) Potential energy $U$ vs. static increments.

(a)

(b)

Figure 7: (a) Position vector of the second node vs. static increments. (b)
Absolute rotation vector of the second node vs. static increments.

###  4.4 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_{\mathrm{e}}=${10}^{4}\text{\,}\frac{}{}$$, shear stiffness
$k_{\mathrm{s}}=${10}^{4}\text{\,}\frac{}{}$$, bending stiffness
$k_{\mathrm{b}}=k_{\mathrm{b}_{y}}=k_{\mathrm{b}_{z}}=${10}^{2}\text{\,}\frac{}{}$$
and torsional stiffness $k_{\mathrm{t}}=${10}^{2}\text{\,}\frac{}{}$$ is
subjected to an incremental application of a tip moment
${}_{K}\mathbf{c}_{1}=\mathbf{A}_{IK}^{\\!\mathrm{T}}(0,\ 0,\ 20\pi
k_{\mathrm{b}}/L)$ and an out-of-plane tip load ${}_{I}\mathbf{b}_{1}=(0,\ 0,\
50)$. By that, the rod is bend up to 10 full circles, while the
$\mathbf{e}_{z}^{I}$-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 $\mathbf{e}_{z}^{I}$-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 $\gamma_{1}$, both shear
strains $\gamma_{2}$, $\gamma_{3}$ and both curvatures $\kappa_{2}$,
$\kappa_{3}$ are step functions, since the used two-node element interpolation
strategy (36) yields element-wise constant strain measures.

Figure 8: Linearly spaced snapshots of the helicoidal deformed rod centerline.

(a)

(b)

(c)

Figure 9: (a) Tip displacement vs. static increments. (b) and (c) Strain
measures of the final configuration.

###  4.5 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=\pi r^{2}$ and second moment of area $I=\pi r^{4}/4$.
The cylinder is subjected to a constant distributed line force
${}_{I}\mathbf{b}=\rho_{0}A(0,\ 0,\ -9.81)$. The stiff rod should be made out
of steel with a uniform density $\rho_{0}=8000$, Young’s modulus
$E=$210\text{\times}{10}^{6}\text{\,}\frac{}{}$$ and shear modulus
$G=E/(2(1+\nu))$, with a Poisson’s ratio $\nu=1/3$. The stiff rod has
consequently an axial stiffness $k_{\mathrm{e}}=EA$, shear stiffness
$k_{\mathrm{s}}=GA$, bending stiffness
$k_{\mathrm{b}}=k_{\mathrm{b}_{y}}=k_{\mathrm{b}_{x}}=EI$ and torsional
stiffness $k_{\mathrm{t}}=2GI$. We also considered a soft rod for which all
stiffnesses were reduced by a factor ${10}^{3}\text{\,}\frac{}{}$. The top was
discretized using a single two-node ${S\\!E}(3)$-element. The initial position
is such that the top points from the origin in positive
$\mathbf{e}_{x}^{I}$-direction, i.e., $\mathbf{q}^{0}=(\mathbf{0}_{3\times
1},\ \mathbf{0}_{3\times 1},\ \mathbf{r}_{1},\mathbf{0}_{3\times 1})$ with
$\mathbf{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., $\mathbf{u}^{0}=(\mathbf{0}_{3\times 1},\ \bm{\Omega},\
\mathbf{v}_{1},\bm{\Omega})$ with the tip velocity
$\mathbf{v}_{1}=\bm{\Omega}\times\mathbf{r}_{1}$ and the angular velocity
$\bm{\Omega}=(\Omega,\ 0,\ \Omega_{\mathrm{pr}})$, where $\Omega=50\pi$ and
$\Omega_{\mathrm{pr}}=gL/(r^{2}\Omega)$. 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\text{\times}{10}^{-8}\text{\,}\frac{}{}$$, see Hairer et al.[37]
and the generalized-$\alpha$ method for first order differential equations
proposed by Jansen et al.[45] with the spectral radius at infinity
$\rho_{\infty}=0.9$ and a step size
$h=$1\text{\times}{10}^{-5}\text{\,}\frac{}{}$$), the simulations were
performed until a final time of $t_{1}=2\pi/\Omega_{\mathrm{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.

(a)

(b)

Figure 10: Tip displacement of the rigid top vs. stiff and soft rod solutions.
The rigid top solution is drawn in black, the stiff rod in red and the soft
rod in blue. (a) Spatial tip trajectory. (b) Tip displacements vs. time.

##  5 Conclusion

We have presented a rod finite element formulation based on a two-node
${S\\!E}(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 ${S\\!E}(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 ${S\\!E}(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 ${S\\!E}(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 $S\\!O(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.

##  Appendix A Matrix Lie groups

###  A.1 The general linear group $G\\!L(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 $G\\!L(n)=\\{\mathbf{A}\in\mathbb{R}^{n\times
n}|\det(\mathbf{A})\neq 0\\}$, which is the group of all $n\times n$
invertible matrices with real entries[52]. For $\mathbf{A},\mathbf{B}\in
G\\!L(n)$, the group operation is given by the usual matrix product
$(\mathbf{A},\mathbf{B})\mapsto\mathbf{A}\mathbf{B}$. The inverse and
transpose of $\mathbf{A}$ are respectively denoted by $\mathbf{A}^{-1}$ and
$\mathbf{A}^{\\!\mathrm{T}}$. The identity is given by the $n\times n$
identity matrix $\mathbf{1}_{n\times n}$. Analogously, we denote the $n\times
n$ matrix containing only zeros by $\mathbf{0}_{n\times n}$.

Let $\mathfrak{gl}(n)=T_{\mathbf{1}_{n\times n}}G\\!L(n)$ be the tangent space
of $G\\!L(n)$ at the identity $\mathbf{1}_{n\times n}$. In Section 3.3 of
Baker[53], it is shown that $\mathfrak{gl}(n)$ corresponds to the set of all
real $n\times n$ matrices, which constitute an $n^{2}$-dimensional vector
space. The Lie algebra of $G\\!L(n)$ consists of the vector space
$\mathfrak{gl}(n)$ together with the bilinear map

$[\bullet,\bullet]\colon\mathfrak{gl}(n)\times\mathfrak{gl}(n)\to\mathfrak{gl}(n)\,,\quad(\mathbf{X},\mathbf{Y})\mapsto[\mathbf{X},\mathbf{Y}]=\mathbf{X}\mathbf{Y}-\mathbf{Y}\mathbf{X}\,,$
(58)

which is a Lie bracket as it is skew-symmetric and satisfies the Jacobi
identity
$[\mathbf{X},[\mathbf{Y},\mathbf{Z}]]+[\mathbf{Y},[\mathbf{Z},\mathbf{X}]]+[\mathbf{Z},[\mathbf{X},\mathbf{Y}]]=\mathbf{0}_{n\times
n}$ $\forall\mathbf{X},\mathbf{Y},\mathbf{Z}\in\mathfrak{gl}(n)$. For a fixed
argument $\mathbf{X}\in\mathfrak{gl}(n)$ the Lie bracket defines the linear
operator

$\operatorname{ad}_{\mathbf{X}}\colon\mathfrak{gl}(n)\to\mathfrak{gl}(n)\,,\quad\mathbf{Y}\mapsto\operatorname{ad}_{\mathbf{X}}(\mathbf{Y})=[\mathbf{X},\mathbf{Y}]\,,$
(59)

called the adjoint map. Thus, instead of
$[\mathbf{X},[\mathbf{X},[\mathbf{X},[\mathbf{X},\mathbf{Y}]]]]$, we can write
$\operatorname{ad}_{\mathbf{X}}^{4}(\mathbf{Y})$, see Hall[52] Section 3.1.
Note that the adjoint map to the power of zero is defined as
$\operatorname{ad}_{\mathbf{X}}^{0}(\mathbf{Y})\coloneqq\mathbf{Y}$. Further,
the identity
$\operatorname{ad}_{-\mathbf{X}}^{i}(\mathbf{Y})=(-1)^{i}\operatorname{ad}_{\mathbf{X}}^{i}(\mathbf{Y})$
readily follows from iteratively applying the bilinearity of the Lie bracket.
According to Baker[53] Section 2.1, the matrix exponential is defined by the
power series

$\exp\colon\mathfrak{gl}(n)\to
G\\!L(n)\,,\quad\mathbf{X}\mapsto\mathbf{A}=\exp(\mathbf{X})=\sum_{i=0}^{\infty}\frac{\mathbf{X}^{i}}{i!}\,.$
(60)

As in Section 2.6 of Iserles et al.[54], for a smooth curve
$\mathbf{Z}(t)\in\mathfrak{gl}(n)$, the derivative of the exponential map can
be written as

$\frac{\mathrm{d}}{\mathrm{d}t}\exp\left(\mathbf{Z}(t)\right)=\operatorname{dexp}_{\mathbf{Z}(t)}\left(\dot{\mathbf{Z}}(t)\right)\exp\left(\mathbf{Z}(t)\right)=\exp\left(\mathbf{Z}(t)\right)\operatorname{dexp}_{-\mathbf{Z}(t)}\left(\dot{\mathbf{Z}}(t)\right)\,,$
(61)

where, for $\mathbf{X}\in\mathfrak{gl}(n)$, the linear map
$\operatorname{dexp}_{\mathbf{X}}$ is defined in terms of a power series as

$\operatorname{dexp}_{\mathbf{X}}\colon\mathfrak{gl}(n)\to\mathfrak{gl}(n)\,,\quad\mathbf{Y}\mapsto\operatorname{dexp}_{\mathbf{X}}(\mathbf{Y})\coloneqq\sum_{i=0}^{\infty}\frac{1}{(i+1)!}\operatorname{ad}_{\mathbf{X}}^{i}(\mathbf{Y})\,.$
(62)

Rearranging the equations in (61) to

$\operatorname{dexp}_{-\mathbf{Z}(t)}\left(\dot{\mathbf{Z}}(t)\right)=\exp\left(\mathbf{Z}(t)\right)^{-1}\frac{\mathrm{d}}{\mathrm{d}t}\exp\left(\mathbf{Z}(t)\right)\quad\text{and}\quad\operatorname{dexp}_{\mathbf{Z}(t)}\left(\dot{\mathbf{Z}}(t)\right)=\frac{\mathrm{d}}{\mathrm{d}t}\exp\left(\mathbf{Z}(t)\right)\exp\left(\mathbf{Z}(t)\right)^{-1}\,,$
(63)

it is reasonable that $\operatorname{dexp}_{-\mathbf{X}(t)}$ and
$\operatorname{dexp}_{\mathbf{X}(t)}$ are called left- and right-trivialized
differentials, respectively. According to Equation (2.46) of Iserles et
al.[54], the inverse of $\operatorname{dexp}_{\mathbf{X}}$ can also be written
by the power series

$\displaystyle\operatorname{dexp}^{-1}_{\mathbf{X}}\colon\mathfrak{gl}(n)\to\mathfrak{gl}(n)\,,\quad\mathbf{Y}\mapsto\operatorname{dexp}^{-1}_{\mathbf{X}}(\mathbf{Y})\coloneqq\sum_{i=0}^{\infty}\frac{B_{i}}{i!}\operatorname{ad}_{\mathbf{X}}^{i}(\mathbf{Y})\,,$
(64)

where $B_{i}$ denotes the $i$-th Bernoulli number.

Let the tuple $(\mathfrak{g},[\cdot,\cdot])$ be a subalgebra of
$(\mathfrak{gl}(n),[\cdot,\cdot])$, i.e., $\mathfrak{g}$ is a subspace of
$\mathfrak{gl}(n)$ with dimension $k\leq n^{2}$ and
$[\mathbf{X},\mathbf{Y}]\in\mathfrak{g}$ for
$\mathbf{X},\mathbf{Y}\in\mathfrak{g}$. Then, we can always find a linear and
bijective map

$j\colon\mathbb{R}^{k}\to\mathfrak{g}\,,\quad\mathbf{x}\mapsto\mathbf{X}=j(\mathbf{x})\,.$
(65)

Since $j$ and $\operatorname{dexp}_{-j(\mathbf{x})}$ are linear maps, it is
convenient to construct the tangent operator
$\mathbf{T}(\mathbf{x})\in\mathbb{R}^{k\times k}$, which is the linear map
defined as

$\mathbf{T}(\mathbf{x})\coloneqq
j^{-1}\circ\operatorname{dexp}_{-j(\mathbf{x})}\circ
j=j^{-1}\circ\sum_{i=0}^{\infty}\frac{(-1)^{i}}{(i+1)!}\operatorname{ad}_{j(\mathbf{x})}^{i}\circ
j=\sum_{i=0}^{\infty}\frac{(-1)^{i}}{(i+1)!}\lceil\mathbf{x}\rfloor^{i}\,,$
(66)

where in the last equality we have introduced the mapping

$\lceil\bullet\rfloor\colon\mathbb{R}^{k}\to\mathbb{R}^{k\times
k}\,,\quad\mathbf{x}\mapsto\lceil\mathbf{x}\rfloor=j^{-1}\circ\operatorname{ad}_{j(\mathbf{x})}\circ
j\,.$ (67)

Often this mapping is applied in the form
$j\big{(}\lceil\mathbf{x}\rfloor\mathbf{y}\big{)}=\operatorname{ad}_{j(\mathbf{x})}\big{(}j(\mathbf{y})\big{)}$,
for $\mathbf{x},\mathbf{y}\in\mathbb{R}^{k}$. If $j(\mathbf{x})$ and
$j(\mathbf{y})$ commute, then powers of the adjoint map
$\operatorname{ad}_{j(\mathbf{x})}^{i}(j(\mathbf{y}))$ vanish for $i>0$
resulting in
$\mathbf{T}(\mathbf{x})\mathbf{y}=\left(j^{-1}\circ\operatorname{dexp}_{-j(\mathbf{x})}\circ
j\right)\mathbf{y}=\mathbf{y}$. The inverse of the tangent operator is given
by

$\mathbf{T}(\mathbf{x})^{-1}=\left(j^{-1}\circ\operatorname{dexp}_{-j(\mathbf{x})}\circ
j\right)^{-1}=j^{-1}\circ\operatorname{dexp}^{-1}_{-j(\mathbf{x})}\circ
j=j^{-1}\circ\sum_{i=0}^{\infty}\frac{(-1)^{i}B_{i}}{i!}\operatorname{ad}_{j(\mathbf{x})}^{i}\circ
j=\sum_{i=0}^{\infty}\frac{(-1)^{i}B_{i}}{i!}\lceil\mathbf{x}\rfloor^{i}\,.$
(68)

Since the matrix rank of $\lceil\mathbf{x}\rfloor$ is at most $k$, the
Cayley–Hamilton theorem guarantees the existence of functions
$a_{0}(\lceil\mathbf{x}\rfloor)$, $\dots$, $a_{k-1}(\lceil\mathbf{x}\rfloor)$
such that

$\lceil\mathbf{x}\rfloor^{k}=\sum_{i=0}^{k-1}a_{i}(\lceil\mathbf{x}\rfloor)\lceil\mathbf{x}\rfloor^{i}\,.$
(69)

Hence, we can express the tangent operator (66) and its inverse (68) by a
finite number of powers of $\lceil\mathbf{x}\rfloor$. Using similar arguments,
also the exponential map (60) can be expressed by a finite number of matrix
powers.

###  A.2 The special orthogonal group $S\\!O(3)$

As a subgroup of $G\\!L(3)$, the special orthogonal group

$S\\!O(3)=\\{\mathbf{A}\in
G\\!L(3)|\mathbf{A}^{\\!\mathrm{T}}\mathbf{A}=\mathbf{1}_{3\times
3},\,\operatorname{det}(\mathbf{A})=+1\\}$ (70)

is the group of all proper orthogonal transformations. In Baker[53] Section
3.3, it is shown that the tangent space at the identity

$\mathfrak{so}(3)=T_{\mathbf{1}_{3\times
3}}S\\!O(3)=\\{\bm{\Omega}\in\mathbb{R}^{3\times
3}|\bm{\Omega}+\bm{\Omega}^{\\!\mathrm{T}}=\mathbf{0}_{3\times 3}\\}$ (71)

corresponds to the three-dimensional space of skew-symmetric $3\times 3$
matrices, which in dynamics often corresponds to the space of angular
velocities. Since $\mathfrak{so}(3)$ is a subspace of $\mathfrak{gl}(3)$ and
the Lie bracket $[\bm{\Omega}_{1},\bm{\Omega}_{2}]\in\mathfrak{so}(3)$ for
$\bm{\Omega}_{1},\bm{\Omega}_{2}\in\mathfrak{so}(3)$, the tuple
$(\mathfrak{so}(3),[\cdot,\cdot])$ is a Lie subalgebra of
$(\mathfrak{gl}(3),[\cdot,\cdot])$.

For $\mathfrak{so}(3)$, the map (65) can be identified explicitly as

$j_{S\\!O(3)}(\bullet)=\widetilde{(\bullet)}\colon\mathbb{R}^{3}\to\mathfrak{so}(3)\,,\quad\bm{\omega}=\begin{pmatrix}\omega_{1}\\\
\omega_{2}\\\
\omega_{3}\end{pmatrix}\mapsto\bm{\Omega}=j_{S\\!O(3)}(\bm{\omega})=\widetilde{\bm{\omega}}=\begin{pmatrix}0&-\omega_{3}&\omega_{2}\\\
\omega_{3}&0&-\omega_{1}\\\ -\omega_{2}&\omega_{1}&0\end{pmatrix}\,,$ (72)

where we have introduced the tilde symbol for compact notation. This map is
related to the cross product by
$\widetilde{\bm{\omega}}\mathbf{r}=\bm{\omega}\times\mathbf{r}$ for
$\bm{\omega},\mathbf{r}\in\mathbb{R}^{3}$.

In order to minimize the number of transcendental function evaluations[23,
20], we introduce the mappings

$\alpha(\widetilde{\bm{\omega}})=\frac{\sin(\|\bm{\omega}\|)}{\|\bm{\omega}\|}\,,\quad\beta(\widetilde{\bm{\omega}})=2\frac{1-\cos(\|\bm{\omega}\|)}{\|\bm{\omega}\|^{2}}\,,\quad\gamma(\widetilde{\bm{\omega}})=\frac{\alpha(\widetilde{\bm{\omega}})}{\beta(\widetilde{\bm{\omega}})}=\frac{\|\bm{\omega}\|}{2}\cot(\|\bm{\omega}\|/2)\,,$
(73)

with $\lim_{\|\bm{\omega}\|\to
0}\alpha(\widetilde{\bm{\omega}})=\lim_{\|\bm{\omega}\|\to
0}\beta(\widetilde{\bm{\omega}})=\lim_{\|\bm{\omega}\|\to
0}\gamma(\widetilde{\bm{\omega}})=1$. Noting that
$\widetilde{\bm{\omega}}^{3}=-\|\bm{\omega}\|^{2}\widetilde{\bm{\omega}}$ and
carefully separating even and odd parts of the power series of the exponential
map (60), see Murray et al.[55] Section 2.2, the exponential map from
$\mathfrak{so}(3)$ to $S\\!O(3)$ has an analytical form known as Rodrigues’
formula

$\exp_{S\\!O(3)}\colon\mathfrak{so}(3)\to
S\\!O(3)\,,\quad\bm{\Omega}\mapsto\mathbf{A}=\exp_{S\\!O(3)}(\bm{\Omega})=\mathbf{1}_{3\times
3}+\alpha(\bm{\Omega})\bm{\Omega}+\frac{\beta(\bm{\Omega})}{2}\bm{\Omega}^{2}\,.$
(74)

For the case $\|\bm{\Omega}\|\to 0$, following the typical approach in
computational methods, the first order approximation
$\exp_{S\\!O(3)}(\bm{\Omega})=\mathbf{1}_{3\times 3}+\bm{\Omega}$ is applied.

From explicitly computing (74) and denoting $\omega=\|\bm{\omega}\|$, the
diagonal terms reveal the identity
$\cos\omega=\tfrac{1}{2}(\operatorname{tr}(\mathbf{A})-1)$, which can be
solved for $\omega=\omega(\mathbf{A})$. The off diagonals directly lead to
$\widetilde{\bm{\omega}}=\tfrac{\omega}{2\sin(\omega)}(\mathbf{A}-\mathbf{A}^{\\!\mathrm{T}})$.
Using
$\sin^{2}(\omega(\mathbf{A}))=1-(\operatorname{tr}(\mathbf{A})-1)^{2}/4$, an
analytic form of the $S\\!O(3)$-logarithm map is given by

$\log_{S\\!O(3)}\colon
S\\!O(3)\to\mathfrak{so}(3)\,,\quad\mathbf{A}\mapsto\bm{\Omega}=\log_{S\\!O(3)}(\mathbf{A})=\frac{\omega(\mathbf{A})}{2\sin\big{(}\omega(\mathbf{A})\big{)}}(\mathbf{A}-\mathbf{A}^{\\!\mathrm{T}})\,.$
(75)

Again, for $\omega\to 0$, we use the first order approximation
$\bm{\Omega}=(\mathbf{A}-\mathbf{A}^{\\!\mathrm{T}})/2$. For notational
convenience, we further introduce the capitalized exponential and logarithm
maps

$\displaystyle\operatorname{Exp}_{S\\!O(3)}\colon\mathbb{R}^{3}\to
S\\!O(3)\,,\quad$ $\displaystyle\bm{\omega}\mapsto\mathbf{A}$
$\displaystyle=\operatorname{Exp}_{S\\!O(3)}(\bm{\omega})=\exp_{S\\!O(3)}\circ
j_{S\\!O(3)}(\bm{\omega})\,,$ (76)
$\displaystyle\operatorname{Log}_{S\\!O(3)}\colon
S\\!O(3)\to\mathbb{R}^{3}\,,\quad$ $\displaystyle\mathbf{A}\mapsto\bm{\omega}$
$\displaystyle=\operatorname{Log}_{S\\!O(3)}(\mathbf{A})=j_{S\\!O(3)}^{-1}\circ\log_{S\\!O(3)}(\mathbf{A})$
$\displaystyle=\frac{\omega(\mathbf{A})}{2\sin\big{(}\omega(\mathbf{A})\big{)}}\big{(}A_{32}-A_{23},\,A_{13}-A_{31},\,A_{21}-A_{12}\big{)}\,,$

which directly relate $\mathbb{R}^{3}$ with $S\\!O(3)$ and vice versa.

For a lighter notation, we shortly suppress the subscripts of $j$ and
$\lceil\bullet\rfloor$ indicating $S\\!O(3)$. Let
$\mathbf{x},\mathbf{y},\mathbf{z}\in\mathbb{R}^{3}$, using the skew-symmetry
and Jacobi’s identity of the cross product, direct computation verifies
$j\big{(}\lceil\mathbf{x}\rfloor\mathbf{y}\big{)}\mathbf{z}=\operatorname{ad}_{j(\mathbf{x})}\big{(}j(\mathbf{y})\big{)}\mathbf{z}=\big{(}j(\mathbf{x})j(\mathbf{y})-j(\mathbf{y})j(\mathbf{x})\big{)}\mathbf{z}=\big{(}\widetilde{\mathbf{x}}\,\widetilde{\mathbf{y}}-\widetilde{\mathbf{y}}\,\widetilde{\mathbf{x}}\big{)}\mathbf{z}=j(\widetilde{\mathbf{x}}\mathbf{y})\mathbf{z}$.
Comparison of both sides of the equality, readily implies
$\lceil\mathbf{x}\rfloor_{S\\!O(3)}=j_{S\\!O(3)}(\mathbf{x})=\widetilde{\mathbf{x}}$.
Again using
$\widetilde{\bm{\omega}}^{3}=-\|\bm{\omega}\|^{2}\widetilde{\bm{\omega}}$ and
carefully separating the terms in the power series (66) (see Iserles et
al.[54] Equation (B.10) and Equation (18) of Park and Chung[23])222Note, both
references define the tangent operator and its inverse in terms of the right-
trivialized differential. Thus, their formulas differ in the sign of the skew-
symmetric part., the $S\\!O(3)$-tangent operator can be written in the form

$\mathbf{T}_{S\\!O(3)}(\bm{\omega})=\sum_{i=0}^{\infty}\frac{(-1)^{i}}{(i+1)!}\lceil\bm{\omega}\rfloor^{i}_{S\\!O(3)}=\mathbf{1}_{3\times
3}-\frac{\beta(\widetilde{\bm{\omega}})}{2}\widetilde{\bm{\omega}}+\frac{1-\alpha(\widetilde{\bm{\omega}})}{\|\bm{\omega}\|^{2}}\widetilde{\bm{\omega}}^{2}\,.$
(77)

Similarly from (68) (see Iserles et al.[54] Equation (B.11), Park and
Chung[23] Equation (19) and Equation (2.20) of Bullo[22]), the inverse
$S\\!O(3)$-tangent operator can efficiently be computed via

$\mathbf{T}^{-1}_{S\\!O(3)}(\bm{\omega})=\sum_{i=0}^{\infty}\frac{B_{i}}{i!}\lceil\bm{\omega}\rfloor^{i}_{S\\!O(3)}=\mathbf{1}_{3\times
3}+\frac{1}{2}\widetilde{\bm{\omega}}+\frac{1-\gamma(\widetilde{\bm{\omega}})}{\|\bm{\omega}\|^{2}}\widetilde{\bm{\omega}}^{2}\,.$
(78)

For $\|\bm{\omega}\|\to 0$, we use the first order approximations
$\mathbf{T}_{S\\!O(3)}(\bm{\omega})=\mathbf{1}_{3\times
3}-\tfrac{1}{2}\widetilde{\bm{\omega}}$ and
$\mathbf{T}^{-1}_{S\\!O(3)}(\bm{\omega})=\mathbf{1}_{3\times
3}+\frac{1}{2}\widetilde{\bm{\omega}}$. Due to the skew-symmetry of
$\widetilde{\bm{\omega}}$, the tangent map and its inverse fulfill the
additional properties
$\mathbf{T}_{S\\!O(3)}(-\bm{\omega})=\mathbf{T}_{S\\!O(3)}^{\\!\mathrm{T}}(\bm{\omega})$
and
$\mathbf{T}^{-1}_{S\\!O(3)}(-\bm{\omega})=\mathbf{T}_{S\\!O(3)}^{-\mathrm{T}}(\bm{\omega})$.

###  A.3 The special Euclidean group ${S\\!E}(3)$

As a subgroup of $G\\!L(4)$, the special Euclidean group

${S\\!E}(3)=S\\!O(3)\ltimes\mathbb{R}^{3}=\left\\{\mathbf{H}=\begin{pmatrix}\mathbf{A}&\mathbf{r}\\\
\mathbf{0}_{1\times 3}&1\end{pmatrix}\in G\\!L(4)\,\bigg{|}\,\mathbf{A}\in
S\\!O(3),\,\mathbf{r}\in\mathbb{R}^{3}\right\\}$ (79)

is the the group of all Euclidean transformations $\mathbf{H}$, which can be
expressed in terms of a translation $\mathbf{r}\in\mathbb{R}^{3}$ and a proper
orthogonal transformation $\mathbf{A}\in S\\!O(3)$. Its linearization at the
identity

$\mathfrak{se}(3)=T_{\mathbf{1}_{4\times
4}}{S\\!E}(3)=\mathfrak{so}(3)\ltimes\mathbb{R}^{3}=\left\\{\bm{\Theta}=\begin{pmatrix}\bm{\Omega}&\mathbf{v}\\\
\mathbf{0}_{1\times 1}&0\end{pmatrix}\in\mathbb{R}^{4\times
4}\bigg{|}\,\bm{\Omega}\in\mathfrak{so}(3),\,\mathbf{v}\in\mathbb{R}^{3}\right\\}$
(80)

is the space of all twists $\bm{\Theta}$, which is a $4\times 4$ matrix
composed of a translational and angular velocity $\mathbf{v}\in\mathbb{R}^{3}$
and $\bm{\Omega}\in\mathfrak{so}(3)$, respectively. Since $\mathfrak{se}(3)$
is a subgroup of $\mathfrak{gl}(4)$ and the Lie bracket
$[\bm{\Theta}_{1},\bm{\Theta}_{2}]\in\mathfrak{se}(3)$ for
$\bm{\Theta}_{1},\bm{\Theta}_{2}\in\mathfrak{se}(3)$, the tuple
$(\mathfrak{se}(3),[\cdot,\cdot])$ is a Lie subalgebra of
$(\mathfrak{gl}(4),[\cdot,\cdot])$.

For $\mathfrak{se}(3)$, the map (65) can be identified explicitly as

$j_{{S\\!E}(3)}\colon\mathbb{R}^{6}\to\mathfrak{se}(3)\,,\quad\bm{\theta}=\begin{pmatrix}\mathbf{v}\\\
\bm{\omega}\end{pmatrix}\mapsto\bm{\Theta}=j_{{S\\!E}(3)}(\bm{\theta})=\begin{pmatrix}j_{S\\!O(3)}(\bm{\omega})&\mathbf{v}\\\
\mathbf{0}_{1\times
3}&0\end{pmatrix}=\begin{pmatrix}\widetilde{\bm{\omega}}&\mathbf{v}\\\
\mathbf{0}_{1\times 3}&0\end{pmatrix}\,.$ (81)

The exponential map from $\mathfrak{se}(3)$ to ${S\\!E}(3)$ is given by (see
Example A.12 of Murray et al.[55])

$\exp_{{S\\!E}(3)}\colon\mathfrak{se}(3)\to{S\\!E}(3)\,,\quad\bm{\Theta}\mapsto\exp_{{S\\!E}(3)}(\bm{\Theta})=\begin{pmatrix}\exp_{S\\!O(3)}(\bm{\Omega})&\quad\mathbf{T}_{S\\!O(3)}^{\\!\mathrm{T}}(\bm{\omega})\mathbf{v}\\\
\mathbf{0}_{1\times 3}&1\end{pmatrix}\,,$ (82)

where $\bm{\omega}=j^{-1}_{S\\!O(3)}(\bm{\Omega})$. By computing
$\bm{\Theta}=\log_{{S\\!E}(3)}\circ\exp_{{S\\!E}(3)}(\bm{\Theta})$ it can be
verified that the ${S\\!E}(3)$-logarithm map is given in terms of $S\\!O(3)$
formulas only and reads as

$\log_{{S\\!E}(3)}\colon{S\\!E}(3)\to\mathfrak{se}(3)\,,\quad\mathbf{H}=\begin{pmatrix}\mathbf{A}&\mathbf{r}\\\
\mathbf{0}_{1\times
3}&1\end{pmatrix}\mapsto\bm{\Theta}=\begin{pmatrix}\bm{\Omega}&\quad\mathbf{T}_{S\\!O(3)}^{-\mathrm{T}}(\bm{\omega})\mathbf{r}\\\
\mathbf{0}_{1\times 3}&0\end{pmatrix}\,,$ (83)

where $\bm{\Omega}=\log_{S\\!O(3)}(\mathbf{A})$ and
$\bm{\omega}=j_{S\\!O(3)}^{-1}(\bm{\Omega})$. Again, the capitalized
exponential and logarithm maps relate $\mathbb{R}^{6}$ with ${S\\!E}(3)$ and
vice versa via

$\displaystyle\operatorname{Exp}_{{S\\!E}(3)}\colon\mathbb{R}^{6}\to{S\\!E}(3)\,,\quad$
$\displaystyle\bm{\theta}\mapsto\mathbf{H}$
$\displaystyle=\operatorname{Exp}_{{S\\!E}(3)}(\bm{\theta})=\exp_{{S\\!E}(3)}\circ
j_{{S\\!E}(3)}(\bm{\theta})\,,$ (84)
$\displaystyle\operatorname{Log}_{{S\\!E}(3)}\colon{S\\!E}(3)\to\mathbb{R}^{6}\,,\quad$
$\displaystyle\mathbf{H}\mapsto\bm{\theta}$
$\displaystyle=\operatorname{Log}_{{S\\!E}(3)}(\mathbf{H})=j_{{S\\!E}(3)}^{-1}\circ\log_{{S\\!E}(3)}(\mathbf{H})$
$\displaystyle=\begin{pmatrix}\mathbf{T}_{S\\!O(3)}^{-\mathrm{T}}\big{(}\operatorname{Log}_{S\\!O(3)}(\mathbf{A})\big{)}\mathbf{r}\\\
\operatorname{Log}_{S\\!O(3)}(\mathbf{A})\end{pmatrix}\,.$

##  Appendix B Variation of strain measures

Using the skew-symmetry of ${}_{K}\delta\widetilde{\bm{\phi}}_{IK}$, the
expression (5)2 can be written as
$\delta\mathbf{A}_{IK}^{\\!\mathrm{T}}=-{}_{K}\delta\widetilde{\bm{\phi}}_{IK}\mathbf{A}_{IK}^{\\!\mathrm{T}}$.
Consequently, the variation of ${}_{K}\bar{\bm{\gamma}}$ from (8) computes to

$\delta({}_{K}\bar{\bm{\gamma}})=\delta(\mathbf{A}_{IK}^{\\!\mathrm{T}}{}_{I}\mathbf{r}_{OP,\xi})=\mathbf{A}_{IK}^{\\!\mathrm{T}}\,\delta\left({}_{I}\mathbf{r}_{OP,\xi}\right)+\delta\mathbf{A}_{IK}^{\\!\mathrm{T}}{}_{I}\mathbf{r}_{OP,\xi}=\mathbf{A}_{IK}^{\\!\mathrm{T}}({}_{I}\delta\mathbf{r}_{P})_{,\xi}-{}_{K}\delta\bm{\phi}_{IK}\times{}_{K}\bar{\bm{\gamma}}\,.$
(85)

In the last equality, we have recognized the virtual displacement (2)1 since
the variation and $\xi$-derivative commute for vector components with respect
to a constant basis. The same property gives rise to the equality

$\mathbf{0}=\mathbf{A}_{IK}^{\\!\mathrm{T}}\left[\left(\delta\left({}_{I}\mathbf{e}_{i}^{K}\right)\right)_{,\xi}-\delta\left(\left({}_{I}\mathbf{e}_{i}^{K}\right)_{,\xi}\right)\right]\,,$
(86)

where ${}_{I}\mathbf{e}_{i}^{K}=\mathbf{A}_{IK}{}_{K}\mathbf{e}_{i}^{K}$,
$i\in\\{x,y,z\\}$. Straightforward computation together with (4) and (5)
yields

$\displaystyle\mathbf{A}_{IK}^{\\!\mathrm{T}}\left(\delta\left({}_{I}\mathbf{e}_{i}^{K}\right)\right)_{,\xi}$
$\displaystyle=\mathbf{A}_{IK}^{\\!\mathrm{T}}\left(\delta\mathbf{A}_{IK}{}_{K}\mathbf{e}_{i}^{K}\right)_{,\xi}$
(87)
$\displaystyle=\mathbf{A}_{IK}^{\\!\mathrm{T}}\left[\mathbf{A}_{IK}\left({}_{K}\delta\bm{\phi}_{IK}\times{}_{K}\mathbf{e}_{i}^{K}\right)\right]_{,\xi}$
$\displaystyle={}_{K}\bar{\bm{\kappa}}_{IK}\times\left({}_{K}\delta\bm{\phi}_{IK}\times{}_{K}\mathbf{e}_{i}^{K}\right)+\left({}_{K}\delta\bm{\phi}_{IK}\right)_{,\xi}\times{}_{K}\mathbf{e}_{i}^{K}$

as well as

$\displaystyle-\mathbf{A}_{IK}^{\\!\mathrm{T}}\delta\left(\left({}_{I}\mathbf{e}_{i}^{K}\right)_{,\xi}\right)$
$\displaystyle=-\mathbf{A}_{IK}^{\\!\mathrm{T}}\delta\left(\mathbf{A}_{IK,\xi}{}_{K}\mathbf{e}_{i}^{K}\right)$
(88)
$\displaystyle=-\mathbf{A}_{IK}^{\\!\mathrm{T}}\delta\left[\mathbf{A}_{IK}\left({}_{K}\bar{\bm{\kappa}}_{IK}\times{}_{K}\mathbf{e}_{i}^{K}\right)\right]$
$\displaystyle=-{}_{K}\delta\bm{\phi}_{IK}\times\left({}_{K}\bar{\bm{\kappa}}_{IK}\times{}_{K}\mathbf{e}_{i}^{K}\right)-\delta\left({}_{K}\bar{\bm{\kappa}}_{IK}\right)\times{}_{K}\mathbf{e}_{i}^{K}$
$\displaystyle={}_{K}\delta\bm{\phi}_{IK}\times\left({}_{K}\mathbf{e}_{i}^{K}\times{}_{K}\bar{\bm{\kappa}}_{IK}\right)-\delta\left({}_{K}\bar{\bm{\kappa}}_{IK}\right)\times{}_{K}\mathbf{e}_{i}^{K}\,.$

Inserting the latter two expressions in (86) and making use of the Jacobi
identity and the skew-symmetry of the cross product, one readily sees that

$\left(\left({}_{K}\delta\bm{\phi}_{IK}\right)_{,\xi}-{}_{K}\delta\bm{\phi}_{IK}\times{}_{K}\bar{\bm{\kappa}}_{IK}-\delta\left({}_{K}\bar{\bm{\kappa}}_{IK}\right)\right)\times{}_{K}\mathbf{e}_{i}^{K}=\mathbf{0}\,.$
(89)

As this holds for all $i\in\\{x,y,z\\}$, the term in the bracket of (89) must
vanish, resulting in the identity

$\delta\left({}_{K}\bar{\bm{\kappa}}_{IK}\right)=\left({}_{K}\delta\bm{\phi}_{IK}\right)_{,\xi}-{}_{K}\delta\bm{\phi}_{IK}\times{}_{K}\bar{\bm{\kappa}}_{IK}\,.$
(90)

##  Appendix C Linearization of internal forces

In order to evaluate the linearization of the internal force vector
$\mathbf{f}^{\mathrm{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 $\omega=\|\bm{\omega}\|\to
0$. Thus, during a numerical implementation a critical angle
$\omega_{\mathrm{crit}}=$1\text{\times}{10}^{-6}\text{\,}\frac{}{}$$ is
defined. Whenever the value of $\omega$ 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 $\mathbf{a}$ by
$a_{i}$, the $i,j$-th component of the matrix $\mathbf{A}$ by $A_{ij}$ and the
components of third order objects are indicated by three subscripts. Using the
identities $\widetilde{\omega}_{ij}=\omega_{k}\varepsilon_{kji}$,
$\partial\widetilde{\omega}_{ij}/\partial\omega_{k}=\varepsilon_{kji}$ and
$\partial(\widetilde{\omega}_{il}\widetilde{\omega}_{lj})/\partial\omega_{k}=\varepsilon_{kli}\widetilde{\omega}_{lj}+\widetilde{\omega}_{il}\varepsilon_{kjl}$,
where $\varepsilon_{ijk}$ denotes the Levi–Civita permutation symbol, we can
compute the derivatives of the auxiliary functions introduced in (73), with
respect to $\omega_{k}$ by

$\frac{\partial\alpha(\widetilde{\bm{\omega}})}{\partial\bm{\omega}}=\frac{\bm{\omega}}{\|\bm{\omega}\|^{2}}\left(\cos(\omega)-\alpha(\widetilde{\bm{\omega}})\right)\,,\quad\frac{\partial\beta(\widetilde{\bm{\omega}})}{\partial\bm{\omega}}=2\frac{\bm{\omega}}{\|\bm{\omega}\|^{2}}\left(\alpha(\widetilde{\bm{\omega}})-\beta(\widetilde{\bm{\omega}})\right)\,,\quad\frac{\partial\gamma(\widetilde{\bm{\omega}})}{\partial\bm{\omega}}=\bm{\omega}\left(\frac{\gamma(\widetilde{\bm{\omega}})}{\|\bm{\omega}\|^{2}}\big{(}1-\gamma(\widetilde{\bm{\omega}})\big{)}-1\right)\,.$
(91)

Thereby, the derivative of the rotation matrix components $A_{ij}$, obtained
from the capitalized $S\\!O(3)$ exponential map (76), with respect to
$\omega_{k}$ is

$\left[\frac{\partial\operatorname{Exp}_{S\\!O(3)}}{\partial\bm{\omega}}(\bm{\omega})\right]_{ijk}=\begin{cases}-\alpha\varepsilon_{ijk}+\frac{\cos\omega-\alpha}{\omega^{2}}\widetilde{\omega}_{ij}\omega_{k}+\frac{\alpha-\beta}{\omega^{2}}\widetilde{\omega}_{il}\widetilde{\omega}_{lj}\omega_{k}+\frac{\beta}{2}\left(\varepsilon_{kli}\widetilde{\omega}_{lj}+\widetilde{\omega}_{il}\varepsilon_{kjl}\right)\,,&\omega>\omega_{\mathrm{crit}}\,,\\\
-\varepsilon_{ijk}\,,&\omega\leq\omega_{\mathrm{crit}}\,.\\\ \end{cases}$ (92)

Similar relations are found in Gallego and Yezzi[56].

Multiplying $\widetilde{\omega}_{nm}=\varepsilon_{kmn}\omega_{k}$ with
$\varepsilon_{imn}$ and using
$\varepsilon_{kmn}\varepsilon_{imn}=2\delta_{ki}$ gives
$\omega_{i}=\frac{1}{2}\varepsilon_{imn}\widetilde{\omega}_{nm}$. Thus, we
have to compute

$\left[\frac{\partial\operatorname{Log}_{S\\!O(3)}}{\partial\mathbf{A}}(\mathbf{A})\right]_{ijk}=\frac{1}{2}\varepsilon_{imn}\left[\frac{\partial\log_{S\\!O(3)}}{\partial\mathbf{A}}(\mathbf{A})\right]_{nmjk}\,.$
(93)

Further, from $\cos\omega=\tfrac{1}{2}(\operatorname{tr}(\mathbf{A})-1)$ with
$\partial\arccos(x)/\partial x=-1/\sqrt{1-x^{2}}$ we find

$\left[\frac{\partial\omega(\mathbf{A})}{\partial\mathbf{A}}\right]_{jk}=-\frac{\delta_{jk}}{2\sin\omega(\mathbf{A})}\quad\text{and}\quad\left[\frac{\partial}{\partial\mathbf{A}}\left(\frac{\omega(\mathbf{A})}{2\sin\omega(\mathbf{A})}\right)\right]_{jk}=\frac{\omega(\mathbf{A})\cos\omega(\mathbf{A})-\sin\omega(\mathbf{A})}{2\sin^{3}\omega(\mathbf{A})}\delta_{jk}\,.$
(94)

Finally, using the identities $\partial A_{ab}/\partial
A_{cd}=\delta_{ac}\delta_{bd}$ and
$\varepsilon_{imn}(\delta_{mj}\delta_{nk}-\delta_{nj}\delta_{mk})=2\varepsilon_{ijk}$,
the derivative of $\omega_{i}$, obtained from the capitalized $S\\!O(3)$
logarithm map (76), with respect to the rotation matrix components $A_{jk}$ is
given by

$\left[\frac{\partial\operatorname{Log}_{S\\!O(3)}}{\partial\mathbf{A}}(\mathbf{A})\right]_{ijk}=\begin{cases}\frac{\omega\cos\omega-\sin\omega}{4\sin^{3}\omega}\varepsilon_{imn}(A_{nm}-A_{mn})\delta_{jk}+\frac{\omega}{2\sin\omega}\varepsilon_{ijk}\,,&\omega>\omega_{\mathrm{crit}}\,,\\\
\frac{1}{2}\varepsilon_{ijk}\,,&\omega\leq\omega_{\mathrm{crit}}\,.\end{cases}$
(95)

Further, the derivatives of the $S\\!O(3)$ tangent operators are required.
Using the identities introduced above they can be computed as

$\left[\frac{\partial\mathbf{T}_{S\\!O(3)}}{\partial\bm{\omega}}(\bm{\omega})\right]_{ijk}=\begin{cases}\frac{\beta}{2}\varepsilon_{ijk}+\widetilde{\omega}_{ij}\omega_{k}\frac{\beta-\alpha}{\omega^{2}}+\frac{1-\alpha}{\omega^{2}}\left(\varepsilon_{kli}\widetilde{\omega}_{lj}+\widetilde{\omega}_{il}\varepsilon_{kjl}\right)+\frac{3\alpha-2-\cos\omega}{\omega^{4}}\widetilde{\omega}_{il}\widetilde{\omega}_{lj}\omega_{k}\,,&\omega>\omega_{\mathrm{crit}}\,,\\\
\frac{1}{2}\varepsilon_{ijk}\,,&\omega\leq\omega_{\mathrm{crit}}\end{cases}$
(96)

and

$\left[\frac{\partial\mathbf{T}^{-1}_{S\\!O(3)}}{\partial\bm{\omega}}(\bm{\omega})\right]_{ijk}=\begin{cases}-\frac{1}{2}\varepsilon_{ijk}+\frac{1-\gamma}{\omega^{2}}\left(\varepsilon_{kli}\widetilde{\omega}_{lj}+\widetilde{\omega}_{il}\varepsilon_{kjl}\right)+\frac{1}{\omega^{2}}\left(2\frac{1-\gamma}{\omega^{2}}+\frac{\gamma}{\omega^{2}}(1-\gamma)-\frac{1}{4}\right)\widetilde{\omega}_{il}\widetilde{\omega}_{lj}\omega_{k}\,,&\omega>\omega_{\mathrm{crit}}\,,\\\
-\frac{1}{2}\varepsilon_{ijk}\,,&\omega\leq\omega_{\mathrm{crit}}\,.\end{cases}$
(97)

The derivatives of the capitalized ${S\\!E}(3)$ exponential and logarithm
introduced in (84) boils down to correctly assemble the just computed
derivatives. Recalling that $\bm{\theta}=(\mathbf{v},\ \bm{\omega})$, we get
the derivative of the exponential map

$\left[\frac{\partial\operatorname{Exp}_{{S\\!E}(3)}}{\partial\bm{\theta}}(\bm{\theta})\right]_{ijk}=\begin{cases}\left[\frac{\partial\operatorname{Exp}_{S\\!O(3)}}{\partial\bm{\omega}}(\bm{\omega})\right]_{ij(k-3)}\,,&\mathrm{for}\quad
1\leq i,j\leq 3\,,\quad 4\leq k\leq 6\,,\\\
\left[\mathbf{T}^{\\!\mathrm{T}}_{S\\!O(3)}(\bm{\omega})\right]_{ik}\,,&\mathrm{for}\quad
1\leq i\leq 3\,,\quad j=4\,,\quad 1\leq k\leq 3\,,\\\
\left[\mathbf{v}\right]_{l}\left[\frac{\partial\mathbf{T}_{S\\!O(3)}}{\partial\bm{\omega}}(\bm{\omega})\right]_{li(k-3)}\,,&\mathrm{for}\quad
1\leq i\leq 3\,,\quad j=4\,,\quad 4\leq k\leq 6\,,\\\
0\,,&\mathrm{else}\,.\end{cases}$ (98)

With

$\mathbf{H}=\begin{pmatrix}\mathbf{A}&\mathbf{r}\\\ \mathbf{0}_{1\times
3}&1\end{pmatrix}$ (99)

the derivative of the ${S\\!E}(3)$ logarithm map with respect to its argument
is given by

$\left[\frac{\partial\operatorname{Log}_{{S\\!E}(3)}}{\partial\mathbf{H}}(\mathbf{H})\right]_{ijk}=\begin{cases}\left[\mathbf{r}\right]_{l}\left[\frac{\partial\mathbf{T}^{-1}_{S\\!O(3)}}{\partial\bm{\omega}}(\bm{\omega})\right]_{lim}\left[\frac{\partial\operatorname{Log}_{S\\!O(3)}}{\partial\mathbf{A}}(\mathbf{A})\right]_{mjk}\,,&\mathrm{for}\quad
1\leq i,j,k\leq 3\,,\\\
\left[\mathbf{T}_{S\\!O(3)}^{-\mathrm{T}}(\bm{\omega})\right]_{ij}\,,&\mathrm{for}\quad\,1\leq
i,j\leq 3\,,\quad k=4,\\\
\left[\frac{\partial\operatorname{Log}_{S\\!O(3)}}{\partial\mathbf{A}}(\mathbf{A})\right]_{(i-3)jk}\,,&\mathrm{for}\quad\,4\leq
i\leq 6\,,\quad 1\leq j,k\leq 3,\\\ \end{cases}$ (100)

where internally $\bm{\omega}=\operatorname{Log}_{S\\!O(3)}(\mathbf{A})$ is
used. Note that again the first order approximation of the required $S\\!O(3)$
quantities have to be used whenever $\omega\leq\omega_{\mathrm{crit}}$.

Finally, the derivative of the element-wise ${S\\!E}(3)$-interpolation (36)
with respect to the generalized coordinates $\mathbf{q}$ can be computed as

$\left[\frac{\partial\mathbf{H}_{\mathcal{I}\mathcal{K}}}{\partial\mathbf{q}}(\xi,\mathbf{q})\right]_{ijk}=\begin{multlined}\sum_{e=0}^{n_{\mathrm{el}}-1}\chi_{\mathcal{J}^{e}}(\xi)\Bigg{\\{}\Bigg{[}\frac{\partial\mathbf{H}_{\mathcal{I}\mathcal{K}_{e}}}{\partial\mathbf{q}}\Bigg{]}_{ilk}\left[\operatorname{Exp}_{{S\\!E}(3)}\left(N^{e}_{1}(\xi)\,\bm{\theta}_{\mathcal{K}_{e}\mathcal{K}_{e+1}}\right)\right]_{lj}\\\
+\left[\mathbf{H}_{\mathcal{I}\mathcal{K}_{e}}\right]_{il}\left[\frac{\partial\operatorname{Exp}_{{S\\!E}(3)}}{\partial\bm{\theta}}\left(N^{e}_{1}(\xi)\,\bm{\theta}_{\mathcal{K}_{e}\mathcal{K}_{e+1}}\right)\right]_{ljm}N^{e}_{1}(\xi)\left[\frac{\partial\operatorname{Log}_{{S\\!E}(3)}}{\partial\mathbf{H}}(\mathbf{H}_{\mathcal{K}_{e}\mathcal{K}_{e+1}})\right]_{mno}\\\
\Bigg{(}\Bigg{[}\frac{\partial\mathbf{H}^{-1}_{\mathcal{I}\mathcal{K}_{e}}}{\partial\mathbf{q}}\Bigg{]}_{npk}\left[\mathbf{H}_{\mathcal{I}\mathcal{K}_{e+1}}\right]_{po}+\left[\mathbf{H}^{-1}_{\mathcal{I}\mathcal{K}_{e}}\right]_{np}\left[\frac{\partial\mathbf{H}_{\mathcal{I}\mathcal{K}_{e+1}}}{\partial\mathbf{q}}\right]_{pok}\Bigg{)}\Bigg{\\}}\,.\end{multlined}\sum_{e=0}^{n_{\mathrm{el}}-1}\chi_{\mathcal{J}^{e}}(\xi)\Bigg{\\{}\Bigg{[}\frac{\partial\mathbf{H}_{\mathcal{I}\mathcal{K}_{e}}}{\partial\mathbf{q}}\Bigg{]}_{ilk}\left[\operatorname{Exp}_{{S\\!E}(3)}\left(N^{e}_{1}(\xi)\,\bm{\theta}_{\mathcal{K}_{e}\mathcal{K}_{e+1}}\right)\right]_{lj}\\\
+\left[\mathbf{H}_{\mathcal{I}\mathcal{K}_{e}}\right]_{il}\left[\frac{\partial\operatorname{Exp}_{{S\\!E}(3)}}{\partial\bm{\theta}}\left(N^{e}_{1}(\xi)\,\bm{\theta}_{\mathcal{K}_{e}\mathcal{K}_{e+1}}\right)\right]_{ljm}N^{e}_{1}(\xi)\left[\frac{\partial\operatorname{Log}_{{S\\!E}(3)}}{\partial\mathbf{H}}(\mathbf{H}_{\mathcal{K}_{e}\mathcal{K}_{e+1}})\right]_{mno}\\\
\Bigg{(}\Bigg{[}\frac{\partial\mathbf{H}^{-1}_{\mathcal{I}\mathcal{K}_{e}}}{\partial\mathbf{q}}\Bigg{]}_{npk}\left[\mathbf{H}_{\mathcal{I}\mathcal{K}_{e+1}}\right]_{po}+\left[\mathbf{H}^{-1}_{\mathcal{I}\mathcal{K}_{e}}\right]_{np}\left[\frac{\partial\mathbf{H}_{\mathcal{I}\mathcal{K}_{e+1}}}{\partial\mathbf{q}}\right]_{pok}\Bigg{)}\Bigg{\\}}\,.$
(101)

##  Appendix 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}\mathbf{r}_{OP}(\xi,t)$ coincide with the
cross-sections’ center of mass and for which $\mathbf{S}_{\rho_{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

$T(\mathbf{q},\mathbf{u})=\frac{1}{2}\int_{\mathcal{B}}{}_{I}\dot{\mathbf{r}}_{OQ}^{\\!\mathrm{T}}{}_{I}\dot{\mathbf{r}}_{OQ}\,\mathrm{d}m=\frac{1}{2}\sum_{e=0}^{n_{\mathrm{el}}-1}\int_{\mathcal{J}^{e}}\left\\{{}_{I}\mathbf{v}_{P}^{\\!\mathrm{T}}A_{\rho_{0}}{}_{I}\mathbf{v}_{P}+{}_{K}\bm{\omega}_{IK}^{\\!\mathrm{T}}{}_{K}\mathbf{I}_{\rho_{0}}{}_{K}\bm{\omega}_{IK}\right\\}J\,\mathrm{d}\xi$
(102)

and the potential energy

$U(\mathbf{q})=\sum_{e=0}^{n_{\mathrm{el}}-1}\int_{\mathcal{J}^{e}}W({}_{K}\bm{\gamma},{}_{K}\bm{\kappa}_{IK};\xi)J\,\mathrm{d}\xi\,.$
(103)

Their respective time derivatives are

$\dot{T}(\mathbf{q},\mathbf{u})=\sum_{e=0}^{n_{\mathrm{el}}-1}\int_{\mathcal{J}^{e}}\Big{\\{}{}_{I}\mathbf{v}_{P}^{\\!\mathrm{T}}A_{\rho_{0}}{}_{I}\dot{\mathbf{v}}_{P}+{}_{K}\bm{\omega}_{IK}^{\\!\mathrm{T}}{}_{K}\mathbf{I}_{\rho_{0}}({}_{K}\bm{\omega}_{IK})^{{\cdot}}\Big{\\}}J\,\mathrm{d}\xi$
(104)

and

$\displaystyle\dot{U}(\mathbf{q})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\int_{\mathcal{J}^{e}}\left\\{{}_{K}\mathbf{n}^{\\!\mathrm{T}}({}_{K}\bar{\bm{\gamma}})^{{\cdot}}+{}_{K}\mathbf{m}^{\\!\mathrm{T}}({}_{K}\bar{\bm{\kappa}}_{IK})^{{\cdot}}\right\\}\mathrm{d}\xi$
(105)
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\int_{\mathcal{J}^{e}}\left\\{({}_{I}\mathbf{v}_{P})_{,\xi}^{\\!\mathrm{T}}\,\mathbf{A}_{IK}{}_{K}\mathbf{n}+({}_{K}\bm{\omega}_{IK})_{,\xi}^{\\!\mathrm{T}}\,{}_{K}\mathbf{m}-{}_{K}\bm{\omega}_{IK}^{\\!\mathrm{T}}({}_{K}\bar{\bm{\gamma}}\times{}_{K}\mathbf{n}+{}_{K}\bar{\bm{\kappa}}_{IK}\times{}_{K}\mathbf{m})\right\\}\mathrm{d}\xi\,,$

where we have used the identities
$({}_{K}\bar{\bm{\gamma}})^{{\cdot}}=\mathbf{A}_{IK}^{\\!\mathrm{T}}({}_{I}\mathbf{v}_{P})_{,\xi}-{}_{K}\bm{\omega}_{IK}\times{}_{K}\bar{\bm{\gamma}}$
and
$({}_{K}\bar{\bm{\kappa}}_{IK})^{{\cdot}}=({}_{K}\bm{\omega}_{IK})_{,\xi}-{}_{K}\bm{\omega}_{IK}\times{}_{K}\bar{\bm{\kappa}}_{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
$\delta\mathbf{s}^{i}=({}_{I}\mathbf{v}_{P_{i}},\
{}_{K_{i}}\bm{\omega}_{IK_{i}})$ 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=\delta W^{\mathrm{dyn}}(\delta\mathbf{s};\mathbf{q},\mathbf{u})+\delta
W^{\mathrm{int}}(\delta\mathbf{s};\mathbf{q})=-\big{(}T(\mathbf{q},\mathbf{u})+U(\mathbf{q})\big{)}^{{\cdot}}\,.$
(106)

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

${}_{I}\mathbf{L}(\mathbf{u})=\int_{\mathcal{B}}{}_{I}\dot{\mathbf{r}}_{OQ}\mathrm{d}m=\sum_{e=0}^{n_{\mathrm{el}}-1}\int_{\mathcal{J}^{e}}A_{\rho_{0}}{}_{I}\mathbf{v}_{P}J\,\mathrm{d}\xi\,.$
(107)

For an arbitrary ${}_{I}\mathbf{c}\in\mathbb{R}^{3}$, let the nodal virtual
displacements be given by

$\delta\mathbf{s}^{i}=\begin{pmatrix}{}_{I}\mathbf{c}\\\ \mathbf{0}_{3\times
1}\end{pmatrix}\,.$ (108)

This choice represents a virtual translation of the rod. Inserting (108) into
the discrete internal virtual work (42) and recognizing that
$N_{0,\xi}^{e}(\xi)+N_{1,\xi}^{e}(\xi)=0$, one readily sees that $\delta
W^{\mathrm{int}}(\delta\mathbf{s};\mathbf{q})=0$. Using (108) in the discrete
inertial virtual work (44) together with $N_{0}^{e}(\xi)+N_{1}^{e}(\xi)=1$
leads to

$\delta
W^{\mathrm{dyn}}(\delta\mathbf{s};\mathbf{q},\mathbf{u})=-{}_{I}\mathbf{c}^{\\!\mathrm{T}}\left(\sum_{e=0}^{n_{\mathrm{el}}-1}\int_{\mathcal{J}^{e}}A_{\rho_{0}}{}_{I}\mathbf{v}_{P}J\,\mathrm{d}\xi\right)^{{\cdot}}=-{}_{I}\mathbf{c}^{\\!\mathrm{T}}{}_{I}\dot{\mathbf{L}}(\mathbf{u})\stackrel{{\scriptstyle!}}{{=}}0\,,$
(109)

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 ${}_{I}\dot{\mathbf{L}}=\mathbf{0}_{3\times 1}$ vanishes
and the linear momentum ${}_{I}{\mathbf{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
${}_{K_{i}}\delta\bm{\phi}_{IK_{i}}(t)\in\mathbb{R}^{3}$. For an arbitrary
${}_{I}\bm{\eta}\in\mathbb{R}^{3}$, the nodal virtual displacements
$\delta\mathbf{s}^{i}=({}_{I}\bm{\eta}\times{}_{I}\mathbf{r}_{OP_{i}},\
\mathbf{A}_{IK_{i}}^{\\!\mathrm{T}}{}_{I}\bm{\eta})$ 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}\bm{\eta}$. 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 ${}_{K_{i}}\delta\bm{\phi}_{IK_{i}}(t)\in\mathbb{R}^{3}$
expressed in the cross-section-fixed basis $K_{i}$, the nodal virtual
rotations ${}_{I}\delta\bm{\phi}_{IK_{i}}(t)\in\mathbb{R}^{3}$ expressed in
the inertial basis $I$ are used. Analogously, the nodal angular velocities
${}_{K_{i}}\bm{\omega}_{IK_{i}}(t)\in\mathbb{R}^{3}$ are replaced by
${}_{I}\bm{\omega}_{IK_{i}}(t)\in\mathbb{R}^{3}$. Proceeding similar to the
original formulation, new discretized internal and inertial virtual work
contributions are obtained. They are

$\displaystyle\delta W^{\mathrm{int}}(\delta\mathbf{s};\mathbf{q})$
$\displaystyle=\delta\mathbf{s}^{\\!\mathrm{T}}\mathbf{f}^{\mathrm{int}}(\mathbf{q})\,,\quad\mathbf{f}^{\mathrm{int}}(\mathbf{q})=\sum_{e=0}^{n_{\mathrm{el}}-1}\mathbf{C}_{e}^{\\!\mathrm{T}}\mathbf{f}^{\mathrm{int}}_{e}(\mathbf{C}_{e}\mathbf{q})\,,$
(110) $\displaystyle\mathbf{f}^{\mathrm{int}}_{e}(\mathbf{q}^{e})$
$\displaystyle=-\int_{\mathcal{J}^{e}}\sum_{i=0}^{1}\Big{\\{}N^{e}_{i,\xi}\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}\mathbf{A}_{IK}{}_{K}\mathbf{n}+N^{e}_{i,\xi}\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}\mathbf{A}_{IK}{}_{K}\mathbf{m}-N^{e}_{i}\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}\left({}_{I}\mathbf{r}_{OP,\xi}\times\mathbf{A}_{IK}{}_{K}\mathbf{n}\right)\Big{\\}}\mathrm{d}\xi\,,$

and

$\displaystyle\delta W^{\mathrm{dyn}}(\delta\mathbf{s};\mathbf{q},\mathbf{u})$
$\displaystyle=-\delta\mathbf{s}^{\\!\mathrm{T}}\left\\{\mathbf{M}(\mathbf{q})\dot{\mathbf{u}}+\mathbf{f}^{\mathrm{gyr}}(\mathbf{q},\mathbf{u})\right\\}\,,$
(111) $\displaystyle\mathbf{M}(\mathbf{q})$
$\displaystyle=\sum_{e=0}^{n_{\mathrm{el}}-1}\mathbf{C}_{e}^{\\!\mathrm{T}}\mathbf{M}_{e}(\mathbf{C}_{e}\mathbf{q})\mathbf{C}_{e}\,,\quad\mathbf{f}^{\mathrm{gyr}}(\mathbf{q},\mathbf{u})=\sum_{e=0}^{n_{\mathrm{el}}-1}\mathbf{C}_{e}^{\\!\mathrm{T}}\mathbf{f}^{\mathrm{gyr}}(\mathbf{C}_{e}\mathbf{q},\mathbf{C}_{e}\mathbf{u})\,,$
$\displaystyle\mathbf{M}_{e}(\mathbf{q}^{e})$
$\displaystyle=\begin{multlined}\int_{\mathcal{J}^{e}}\sum_{i=0}^{1}\sum_{k=0}^{1}N^{e}_{i}N^{e}_{k}\Big{\\{}\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}A_{\rho_{0}}\mathbf{1}_{3\times
3}\mathbf{C}_{\mathbf{r},k}+\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}{}_{I}\mathbf{I}_{\rho_{0}}\mathbf{C}_{\bm{\psi},k}\Big{\\}}J\mathrm{d}\xi\,,\end{multlined}\int_{\mathcal{J}^{e}}\sum_{i=0}^{1}\sum_{k=0}^{1}N^{e}_{i}N^{e}_{k}\Big{\\{}\mathbf{C}_{\mathbf{r},i}^{\\!\mathrm{T}}A_{\rho_{0}}\mathbf{1}_{3\times
3}\mathbf{C}_{\mathbf{r},k}+\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}{}_{I}\mathbf{I}_{\rho_{0}}\mathbf{C}_{\bm{\psi},k}\Big{\\}}J\mathrm{d}\xi\,,$
$\displaystyle\mathbf{f}^{\mathrm{gyr}}_{e}(\mathbf{q}^{e},\mathbf{u}^{e})$
$\displaystyle=\int_{\mathcal{J}^{e}}\sum_{i=0}^{1}N^{e}_{i}\mathbf{C}_{\bm{\psi},i}^{\\!\mathrm{T}}{}_{I}\widetilde{\bm{\omega}}_{IK}{}_{I}\mathbf{I}_{\rho_{0}}{}_{I}\bm{\omega}_{IK}J\mathrm{d}\xi\,,$

with

${}_{I}\mathbf{I}_{\rho_{0}}(\mathbf{q})=\mathbf{A}_{IK}(\mathbf{q}){}_{K}\mathbf{I}_{\rho_{0}}\mathbf{A}_{IK}(\mathbf{q})^{\\!\mathrm{T}}\,.$
(112)

Note, with this reformulation the mass matrix $\mathbf{M}(\mathbf{q})$ and the
gyroscopic forces $\mathbf{f}^{\mathrm{gyr}}(\mathbf{q},\mathbf{u})$ depend on
the generalized coordinates $\mathbf{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

$\delta\mathbf{s}^{i}=\begin{pmatrix}{}_{I}\bm{\eta}\times{}_{I}\mathbf{r}_{OP_{i}}\\\
{}_{I}\bm{\eta}\end{pmatrix}\,,$ (113)