1 Introduction

In this paper, we present a formulation for spatial rigid body motion based on a set of non-redundant, homogeneous local (translational) velocity coordinates. The initial motivation for the present work was the investigation of local velocity coordinates of a first order H(curl)-conforming finite element proposed by Nédélec [25], which turned out to be relevant for rigid bodies as well, because the latter finite element offers six non-redundant and—regarding their physical interpretation—equivalent (= unified) displacement coordinates as degrees of freedom. Since the spatial motion of a rigid body is known, if the velocity (three translational velocity coordinates) of one of its points as well as its angular velocity (three angular velocity coordinates) are known [14], the question arises if it is possible to describe a rigid body motion using six unified local velocity coordinates, which correspond to the displacement coordinates of this first order H(curl)-conforming triangular or tetrahedral finite element and whether and how the resulting equations can be solved. Furthermore, the description of the spatial motion of a body in space on velocity level based on purely translational velocity coordinates could also be relevant for other formulations such as the absolute nodal coordinate formulation (ANCF), because no rotation parameters are used in ANCF [15]. The description of the motion of a rigid body in space with pure translational velocity coordinates could also be advantageous in model reduction methods such as the so-called generalized component mode synthesis [28]. Here, however, we will only focus on the description of rigid body motion.

Modeling mechanical systems in such a way that they can be efficiently simulated plays an important role in many areas, e.g. in real-time simulation and control [19]. The number, complexity and computational efficiency of the equations describing the system strongly depend on the choice of coordinates used to model the mechanical system [19, 36]. A common approach to describing the kinematics of a multibody system is to divide the rigid body motion into a translational motion and a rotational motion [4, 14, 26, 32], as is the case with the Newton–Euler equations, for example.

Six coordinates are at least required to define the configuration of a rigid body in space [26]. Three Cartesian coordinates usually define the origin of the body’s reference frame. In engineering applications, e.g. rotor dynamics, the body’s orientation is often characterized by three angles of rotation like Euler angles or Bryant angles [3]. However, in the case of large rotations, the representation of the body’s orientation in terms of three angles of rotation suffers from singularities [32]. Alternative representations, which provide a singularity-free description of the body’s orientation either use for example the components of a rotation matrix or Euler parameters as degrees of freedom [26].

An alternative set of coordinates has been presented by Garcia de Jalon and co-workers in the early 1980s [13]. They created and further developed the natural coordinates to describe the 2-D and 3-D motion of multibody systems [12]. The natural coordinates approach uses Cartesian coordinates of three (or more) non-aligned points to define the position of a solid in space [11]. Natural coordinates require neither angles nor any form of angle parametrization to describe the orientation of a body in space, which leads to constant inertia matrices and a simple form of the equations of motion [11]. Natural coordinates can be seen as a set of dependent parameters used to describe rotations [11]. In contrast to the coordinates used in [26], the description of the motion of a multibody system with the natural coordinates yields a unified form of the descriptive equations of motion but at the expense of the use of non-redundant position coordinates.

Another approach to describing the configuration of a rigid body is to connect the special orthogonal group \(\mathit{SO}(3)\) with the linear space \(\mathbb{R}^{3}\) either with the Cartesian product leading to \(\mathit{SO}(3)\times \mathbb{R}^{3}\) or with the semi-direct product leading to \(\mathit{SO}(3)\ltimes \mathbb{R}^{3}\) (special Euclidean group \(\mathit{SE}(3)\)) [3]. Note, both the translational velocity and the angular velocity are represented within the \(\mathit{SE}(3)\) formulation in a body fixed coordinate system. As a result, the equations of motion are differential equations on a manifold with tangent spaces parameterized by the corresponding Lie algebra [3]. These equations can be solved using matrix Lie-group time integration schemes [34].

Nédélec introduced two new families of mixed finite elements in \(\mathbb{R}^{3}\), which are conforming on the Sobolev spaces H(div) and H(curl) [25]. Those element families are commonly used to approximate solutions of Maxwell’s equations [20]. Nédélec elements are also used in solid mechanics, e.g., Pechstein and Schöberl [27] introduced the so-called’tangential-displacement normal–normal-stress’ (TDNNS) formulation for elasticity. Especially the H(curl) conforming element of lowest order, which is also called the Nédélec element of the first-kind motivates this work, due to the fact that the degrees of freedom are the average value of the tangential component of the vector field on each edge. Elements of this type are generally referred to as edge elements because the degrees of freedom are associated with the edges of the mesh at the lowest order [20]. A graphical representation of the degrees of freedom of the first order H(curl) conforming element is shown in Fig. 1. The number of degrees of freedom of the Nédélec element of the first kind is exactly the same as the minimum number of velocity coordinates required to describe the motion of a rigid body at velocity level. The idea is therefore to express a local angular velocity by two local translational velocities, whereby the translational velocities are assigned to properly selected directions similar to the degrees of freedom of the Nédélec element of the first kind. That means, the motion of a rigid body is described by using only local translational velocity coordinates. Moreover, we want to point out that these translational velocity coordinates are fully equivalent (= unified) regarding their physical interpretation, which means that no partition into translational and rotational parts is introduced. Since these translational velocity coordinates are represented in a body-fixed frame, Lie-group time integrators based on the exponential map such as e.g. the ones presented by Crouch and Grossmann [10] or Munthe-Kaas [22, 23] could be applied to obtain information about the position and orientation of a rigid body.

Fig. 1
figure 1

Degrees of freedom of a Nédélec element of the first kind

The authors would like to emphasize that the idea to replace a rotation by two opposed displacements (or vice versa), an angular velocity by two translational velocities, or a torque by a force couple is not new at all. However, we are not aware of any formulation of the equations of motion for a rigid body that is purely based on non-redundant translational velocities. Furthermore, we show how the resulting equations of motion can be solved numerically and how to obtain position and orientation in a stable way.

The remaining paper is organized as follows. In Sect. 2 we introduce the Lie-group formalism for rigid bodies, which is only used in the following to integrate the local unified velocity coordinates. We have avoided relying on the general theory of Lie groups in Sect. 2. Of course, a reader who is already familiar with this general theory will find no new information in this section. However, for a detailed description of the multibody kinematics and dynamics with Lie groups we would like to refer the reader to literature [3, 8]. Section 3 is devoted to the fundamentals of the rigid body formulation with non-redundant unified local velocity coordinates as well as to the derivation of the equations of motion and their generalized forces. In Sect. 4, we show the relation of the proposed formulation to the Newton–Euler equations and present a further analysis of the proposed rigid body formulation. The proper time integration of the equations of motion and the non-redundant unified local velocity coordinates is shown in Sect. 5. Lastly, the proposed formulation is validated through a numerical example in Sect. 6.

The present paper is an extension of [17] regarding a significantly improved notation and additional results of the velocity projection along the edges of a tetrahedron.

2 Lie-group framework

Since the six unified velocity coordinates are defined in a body-fixed frame, standard integrators such as the generalized-\(\alpha \) method [9] are inappropriate for obtaining (drift-free) information about the current position and orientation of a rigid body via direct numerical integration. However, Lie-group time integrators based on the exponential map such as e.g. the ones presented by Crouch and Grossmann [10], Munthe-Kaas [22, 23], Brüls and Cardona [5], Terze et al. [35] or by Arnold et al. [3] can be used in an adapted form, which is why we introduce the Lie-group representation of rigid body kinematics on the special Euclidean group \(\mathit{SE}(3)\) in this section.

In the following, capital boldface letters represent matrices and lower case boldface letters represent vectors except the local translational velocity vector and the local angular velocity vector are written in capital boldface letters according to the notation used in [3]. The horizontal line above a vector \(\bar{ \mathbf{a}}\) means that the vector is a column vector whose components are represented in a body-fixed (local) frame. The components of a vector are written in non-bold letters. If a horizontal bar is drawn above the components of a vector, this means that these are the coordinates of the vector represented in the body-fixed frame. The same convention is used for matrices and their components.

To obtain a singularity-free representation of the orientation of the rigid body with respect to a reference frame, for example either a rotation matrix,

$$ \mathbf{R}\in \mathit{SO}(3):= \bigl\{ \mathbf{R}\in \mathbb{R}^{3 \times 3} \, | \, \mathbf{R}^{T}\mathbf{R}= \mathbf{I}_{3\times 3}, \ \textrm{det}(\mathbf{R})=+1 \bigr\} $$
(1)

or Euler parameters are used as degrees of freedom [14]. Rigid body transformations are represented by frame transformations, which belong to the special Euclidean group \(\mathit{SE}(3)\) [33]. The action of a rigid body transformation on a body-fixed frame describes how the body-fixed frame rotates and displaces with respect to a inertial frame [24]. Rigid body transformations are such that the position vector \(\mathbf{x}_{P}\in \mathbb{R}^{3}\) of any point \(P\) in the current configuration on the rigid body is expressed by the affine mapping

$$ \mathbf{x}_{P} = \mathbf{x}+ \mathbf{y}_{P} = \mathbf{x}+ \mathbf{R}\bar{ \mathbf{y}}_{P} $$
(2)

where \(\mathbf{x}\in \mathbb{R}^{3}\) represents the current position vector of a reference point, which was located at the point \(O\) in the reference configuration [33]. This is illustrated in Fig. 2. The vector \(\mathbf{y}_{P} = \mathbf{R}\bar{\mathbf{y}}_{P}\) represents the position of an arbitrary point \(P\) of the body with respect to the fixed reference frame I. Whereas the vector \(\bar{\mathbf{y}}_{P}\in \mathbb{R} ^{3}\) represents the position of the point \(P\) in the body-fixed frame. The rotation matrix \(\mathbf{R}\) maps the vector \(\bar{\mathbf{y}} _{P}\) from the body-fixed frame into the reference frame and defines the orientation of the body-fixed frame with respect to the reference frame I. A rigid body transformation (2) can conveniently be represented in matrix form using the homogeneous representation of vectors as

$$ \begin{bmatrix} \mathbf{x}_{P} \\ 1 \end{bmatrix} = \begin{bmatrix} \mathbf{R}& \mathbf{x} \\ \mathbf{0}_{1\times 3} & 1 \end{bmatrix} \begin{bmatrix} \bar{\mathbf{y}}_{P} \\ 1 \end{bmatrix} = \mathbf{H}( \mathbf{R}, \,\mathbf{x}) \begin{bmatrix} \bar{\mathbf{y}}_{P} \\ 1 \end{bmatrix} . $$
(3)

The set of 4 × 4-matrices, which are formed by the semi-direct product \(\mathit{SO}(3)\ltimes \mathbb{R}^{3}\),

$$ \mathbf{H}= \mathbf{H}( \mathbf{R}, \,\mathbf{x}) := \bigl\{ \mathbf{H}\in \mathbb{R}^{4\times 4} \,\, | \,\, \mathbf{H}\in \mathit{SO}(3)\ltimes \mathbb{R}^{3} \bigr\} , \qquad \mathbf{H}= \begin{bmatrix} \mathbf{R}& \mathbf{x} \\ \mathbf{0}_{1\times 3} & 1 \end{bmatrix} $$
(4)

where \(\mathbf{R}\in \mathit{SO}(3)\) and \(\mathbf{x}\in \mathbb{R} ^{3}\) forms together with the matrix product (composition operation) a matrix Lie group, which is called the special Euclidean group \(\mathit{SE}(3)\) [30]. The velocity vector \(\dot{\mathbf{x}}_{P}\in \mathbb{R}^{3}\) of any point \(P\) on the rigid body can be expressed as linear relation between \(\dot{\mathbf{H}}\) and \([\bar{\mathbf{y}}_{P}^{T} \ 1 ]^{T}\) in the form

$$ \begin{bmatrix} \dot{\mathbf{x}}_{P} \\ 0 \end{bmatrix} = \dot{ \mathbf{H}} \begin{bmatrix} \bar{\mathbf{y}}_{P} \\ 1 \end{bmatrix} . $$
(5)

The time derivative \(\dot{\mathbf{H}}\) reads

$$ \dot{\mathbf{H}}= \begin{bmatrix} \mathbf{R}\widetilde{\overline{\boldsymbol{\Omega}}}& \mathbf{R}\overline{ \mathbf{U}} \\ \mathbf{0}_{1\times 3} & 0 \end{bmatrix} = \begin{bmatrix} \mathbf{R}& \mathbf{x} \\ \mathbf{0}_{1\times 3} & 1 \end{bmatrix} \begin{bmatrix} \widetilde{\overline{\boldsymbol{\Omega}}}& \overline{\mathbf{U}} \\ \mathbf{0}_{1\times 3} & 0 \end{bmatrix} = \mathbf{H}\,\widetilde{\overline{\mathbf{v}}}_{\textrm{L}}. $$
(6)

The 4 × 4-matrix

$$ \widetilde{\overline{\mathbf{v}}}_{\textrm{L}}= \begin{bmatrix} \widetilde{\overline{\boldsymbol{\Omega}}}& \overline{\mathbf{U}} \\ \mathbf{0}_{1\times 3} & 0 \end{bmatrix} $$
(7)

denotes any element of \(\mathfrak{se}(3)\), the Lie algebra of \(\mathit{SE}(3)\), which is defined by

$$ \widetilde{\overline{\boldsymbol{\Omega}}}\in \mathfrak{so}(3) := \bigl\{ \widetilde{\overline{\boldsymbol{\Omega}}}\in \mathbb{R}^{3\times 3} \, \, | \,\, \widetilde{\overline{\boldsymbol{\Omega}}}+ \widetilde{\overline{ \boldsymbol{\Omega}}}^{T} = \mathbf{0}_{3\times 3} \bigr\} , \qquad \widetilde{\overline{\boldsymbol{\Omega}}}= \begin{bmatrix} 0 & -\overline{\varOmega }_{3} & \overline{\varOmega }_{2} \\ \overline{\varOmega }_{3} & 0 & -\overline{\varOmega }_{1} \\ -\overline{\varOmega }_{2} & \overline{\varOmega }_{1} & 0 \end{bmatrix} , $$
(8)

with \(\widetilde{\overline{\boldsymbol{\Omega}}}= \mathbf{R}^{T} \dot{\mathbf{R}}\). The Lie algebra is isomorphic to a finite dimensional linear space \(\mathbb{R}^{k}\) with an invertible linear mapping

$$ \widetilde{ ( \bullet )} : \mathbb{R}^{k} \rightarrow \mathfrak{g}, \quad \mathbf{q} \mapsto \widetilde{ \mathbf{q}}, $$
(9)

where the tilde operator \(\widetilde{ ( \bullet )}\) depend on the specific Lie-group setting [3, 7]. In the case of the Lie algebra \(\mathfrak{se}(3)\), the inverse mapping (9) maps a 4 × 4-matrix \(\widetilde{\overline{\mathbf{v}}}_{\textrm{L}}\in \mathfrak{se}(3)\) to a six-dimensional vector

$$ \overline{\mathbf{v}}_{\textrm{L}}= \begin{bmatrix} \,\overline{\mathbf{U}}^{T} & \overline{\boldsymbol{\Omega}}^{T}\,\end{bmatrix} ^{T} $$
(10)

and vice versa. In the case of Lie algebra \(\mathfrak{so}(3)\), the inverse mapping (9) maps a 3 × 3-matrix \(\widetilde{\overline{\boldsymbol{\Omega }}}\in \mathfrak{so}(3)\) to a three-dimensional vector

$$ \overline{\boldsymbol{\Omega}}= \begin{bmatrix} \overline{\varOmega }_{1} & \overline{\varOmega }_{2} & \overline{\varOmega } _{3} \end{bmatrix} ^{T} $$
(11)

and vice versa. In the vector representation \(\overline{\mathbf{v}} _{\textrm{L}}\) of \(\widetilde{\overline{\mathbf{v}}}_{\textrm{L}}\), the sub-vector

$$ \overline{\mathbf{U}}= [\,\overline{U}_{1} \quad \overline{U} _{2} \quad \overline{U}_{3} ]^{T} $$
(12)

represents the translational velocity in a body-fixed frame and the sub-vector \(\overline{\boldsymbol{\Omega }}\) represents the angular velocity vector in a body-fixed frame [3, 6]. Equation (6) is a first order differential equation on a matrix Lie group, which produces a relation between \(\mathbf{H}\) and \(\dot{\mathbf{H}}\) [18, 30]. This differential equation is not solved directly. Instead, the incremental motion vector differential equation

$$ \dot{\mathbf{q}}(t) = \mathbf{T}_{\mathit{SE}(3)}^{-1} \bigl(\mathbf{q}(t)\bigr) \overline{\mathbf{v}}_{\textrm{L}}(t), \quad \mathbf{q}_{0}= \mathbf{0} , $$
(13)

with \(\mathbf{q}_{0}\) being initial condition is solved with a time integration scheme for

$$ \mathbf{q}= \begin{bmatrix} q_{1} & q_{2} & q_{3} & q_{4} & q_{5} & q_{6} \end{bmatrix} ^{T} = \begin{bmatrix} \mathbf{q}_{T}^{T} & \mathbf{q}_{R}^{T} \end{bmatrix} ^{T} $$
(14)

where

$$ \mathbf{q}_{T}:= \begin{bmatrix} q_{1} & q_{2} & q_{3} \end{bmatrix} ^{T} \quad \textrm{and}\quad \mathbf{q}_{R}:= \begin{bmatrix} q_{4} & q_{5} & q_{6} \end{bmatrix} ^{T}. $$
(15)

While the subindex \(T\) in \(\mathbf{q}_{T}\) stands for translation, the subindex \(R\) in \(\mathbf{q}_{R}\) stands for rotation. The solution of the differential equation (6) can then be determined for each incremental time step using the matrix exponential map \(\textrm{exp}_{\mathit{SE}(3)}\) on \(\mathit{SE}(3)\)

$$\begin{aligned} \mathbf{H}(t+h) = \mathbf{H}(t) \textrm{exp}_{\mathit{SE}(3)}(\Delta \mathbf{q}), \qquad \mathbf{H}(0)=\mathbf{H} _{0} , \end{aligned}$$
(16)

with \(\Delta \mathbf{q}= \mathbf{q}_{t+h} - \mathbf{q}_{t}\) where \(h\) is a (usually small) step forward in time. In Eq. (13), \(\mathbf{T}_{ \mathit{SE}(3)}^{-1}\) denotes the inverse tangent operator corresponding to the exponential map on \(\mathit{SE}(3)\) shown in Eq. (17). The exponential map and the inverse tangent operator on \(\mathit{SE}(3)\) can be written in a closed form

$$\begin{aligned} \textrm{exp}_{\mathit{SE}(3)}(\mathbf{q}) &= \begin{bmatrix} \textrm{exp}_{\mathit{SO}(3)}(\mathbf{q}_{R}) & \mathbf{T}_{\mathit{SO}(3)}^{T}( \mathbf{q}_{R})\mathbf{q}_{T}\\ \mathbf{0}_{1\times 3} & 1 \end{bmatrix} , \end{aligned}$$
(17)
$$\begin{aligned} \mathbf{T}_{\mathit{SE}(3)}^{-1}(\mathbf{q}) &= \begin{bmatrix} \mathbf{T}_{\mathit{SO}(3)}^{-1}(\mathbf{q}_{R}) & \mathbf{T}_{q_{T}q_{R}-}( \mathbf{q}_{T},\, \mathbf{q}_{R}) \\ \mathbf{0}_{3\times 3} & \mathbf{T}_{\mathit{SO}(3)}^{-1}(\mathbf{q}_{R}) \end{bmatrix} , \end{aligned}$$
(18)

where \(\textrm{exp}_{\mathit{SO}(3)}\) is the exponential map and \(\mathbf{T}_{\mathit{SO}(3)}\) the tangent operator on the special orthogonal group \(\mathit{SO}(3)\) [3, 6, 33], which are given by

$$\begin{aligned} \textrm{exp}_{\mathit{SO}(3)}(\mathbf{q}_{R}) &= \mathbf{I}_{3\times 3} + \frac{ \sin (\Vert \mathbf{q}_{R}\Vert )}{\Vert \mathbf{q}_{R}\Vert } \widetilde{ \mathbf{q}}_{R}+ \frac{1-\cos (\Vert \mathbf{q}_{R}\Vert )}{ \Vert \mathbf{q}_{R}\Vert ^{2}}\widetilde{ \mathbf{q}}_{R}^{2}, \end{aligned}$$
(19)
$$\begin{aligned} \mathbf{T}_{\mathit{SO}(3)}(\mathbf{q}_{R}) &= \mathbf{I}_{3\times 3} + \frac{ \cos (\Vert \mathbf{q}_{R}\Vert )-1}{\Vert \mathbf{q}_{R}\Vert ^{2}} \widetilde{ \mathbf{q}}_{R}+ \frac{1-\frac{\sin (\Vert \mathbf{q}_{R} \Vert )}{\Vert \mathbf{q}_{R}\Vert }}{\Vert \mathbf{q}_{R}\Vert ^{2}} \widetilde{ \mathbf{q}}_{R}^{2}, \end{aligned}$$
(20)
$$\begin{aligned} \mathbf{T}_{\mathit{SO}(3)}^{-1}(\mathbf{q}_{R}) &= \mathbf{I}_{3\times 3} + \frac{1}{2}\widetilde{\mathbf{q}}_{R}+ \frac{1-\frac{\Vert \mathbf{q} _{R}\Vert }{2}\cot (\frac{\Vert \mathbf{q}_{R}\Vert }{2})}{\Vert \mathbf{q}_{R}\Vert ^{2}}\widetilde{\mathbf{q}}_{R}^{2}, \end{aligned}$$
(21)
$$\begin{aligned} \mathbf{T}_{q_{T}q_{R}-}(\mathbf{q}_{T},\, \mathbf{q}_{R}) &= \frac{1}{2}\widetilde{\mathbf{q}}_{T}+ \frac{\frac{1}{2}(2-\cot (\frac{ \Vert \mathbf{q}_{R}\Vert }{2}))}{\Vert \mathbf{q}_{R}\Vert ^{2}} \lceil \widetilde{\mathbf{q}}_{T},\,\widetilde{ \mathbf{q}}_{R}\rceil + \cdots \\ &\quad {}+\frac{\frac{1}{2} ( \frac{\Vert \mathbf{q}_{R}\Vert ^{2}}{1- \cos (\Vert \mathbf{q}_{R}\Vert )} + \cot (\frac{\Vert \mathbf{q}_{R} \Vert }{2}) )}{\Vert \mathbf{q}_{R}\Vert ^{4}} \bigl(\mathbf{q} _{R}^{T} \mathbf{q}_{T} \bigr)\widetilde{\mathbf{q}}_{R}^{2}, \end{aligned}$$
(22)

using the Lie bracket \(\lceil ,\rceil \); see e.g. [22]. The matrices \(\widetilde{\mathbf{q}}_{R}\) and \(\widetilde{\mathbf{q}}_{T}\) are skew symmetric matrices of the form (8), which are computed from the vectors \(\mathbf{q}_{R}\) and \(\mathbf{q}_{T}\).

Fig. 2
figure 2

Rigid body transformation

3 Spatial rigid body motion using non-redundant unified local velocity coordinates

In many approaches [4, 14, 32], the kinematics of the rigid body and its equations of motion are divided into a translational part and a rotational part of motion. This subdivision of motion at velocity level is common in multibody dynamics; see e.g. [3, 6, 33]. A unification of the velocity coordinates within the description of the rigid body motion could be advantageous, which is shown in the following.

3.1 Non-redundant unified local velocity coordinates

The unification of the local translational velocity coordinates \((\overline{U}_{1}, \overline{U}_{2}, \overline{U}_{3} )\) and the local angular velocity coordinates \((\overline{\varOmega }_{1}, \overline{\varOmega }_{2}, \overline{\varOmega }_{3})\) to six local velocity coordinates \((\overline{w}_{1}, \overline{w}_{2}, \overline{w}_{3}, \overline{w}_{4}, \overline{w}_{5}, \overline{w}_{6})\) is achieved, by projecting the local translational velocity coordinates of six selected points on the rigid body along six selected directions. The scalar product of the direction \(\bar{ \mathbf{b}}_{i}\) and the velocity vector \(\overline{\mathbf{v}}_{i}(\bar{ \mathbf{y}}_{i})\), see Fig. 3, which is related to the velocity vector \(\dot{\mathbf{x}}_{P}\) by the relation

$$ \overline{\mathbf{v}}_{P}(\bar{\mathbf{y}}_{P}) = \mathbf{R}^{T} \dot{\mathbf{x}}_{P}(\bar{ \mathbf{y}}_{P}), $$
(23)

generates the \(i\)th unified local velocity coordinate \(\overline{w} _{i}\in \mathbb{R,}\)

$$ \overline{w}_{i} := \bar{\mathbf{b}}_{i}^{T} \overline{\mathbf{v}}_{i}(\bar{ \mathbf{y}}_{i}). $$
(24)
Fig. 3
figure 3

Basic projection idea to obtain unified local velocity coordinates. \(S\) denotes the center of mass of the rigid body. The six points, used for the velocity projection are denoted with \(P_{1} \ldots P_{6}\)

The \(i\)th local velocity vector of the \(i\)th point \(P_{i}\) is indicated by the vector \(\overline{\mathbf{v}}_{i}\in \mathbb{R}^{3}\). The local position of point \(P_{i}\) on the rigid body is marked with vector \(\bar{\mathbf{y}}_{i}\in \mathbb{R}^{3}\). The vector \(\bar{\mathbf{b}}_{i}\in \mathbb{R}^{3}\) represents the \(i\)th vector used for the projection. Therefore, the vector of unified local velocity coordinates reads

$$ \mathbf{\overline{w}}:= \begin{bmatrix} \overline{w}_{1} & \overline{w}_{2} & \overline{w}_{3} & \overline{w} _{4} & \overline{w}_{5} & \overline{w}_{6} \end{bmatrix} ^{T}. $$
(25)

The term’unified’ emphasizes the fact that all velocity coordinates \(\overline{w}_{i}\) have the same kinematic meaning and that no distinction is made between a translational motion and a rotational motion. The general idea of the projection is depicted in Fig. 3. The vectors \(\mathbf{e}_{i}\), \(\mathbf{y}_{i}\), \(\mathbf{v}_{i}\) and \(\mathbf{y}_{i}\) represent body-fixed (physical) vectors. Whereas the vectors \(\overline{\mathbf{e}}_{i}\), \(\bar{\mathbf{b}}_{i}\), \(\overline{\mathbf{v}}_{i}\) and \(\bar{\mathbf{y}}_{i}\) are column vectors whose components are represented in a local frame.

The set of projection points \(P_{i}\) and direction vectors \(\bar{ \mathbf{b}}_{i}\) to generate the six unified velocity coordinates \(\overline{w}_{i}\) is limited by the condition that there must be a bijective mapping between the velocity coordinates \(\overline{w}_{i}\) and the vector form of \(\widetilde{\overline{\mathbf{v}}}_{\textrm{L}}\). In Sect. 4.2, the condition is shown which is required to achieve this bijective mapping. This restriction to the set of permitted projection points and direction vectors is necessary, as otherwise the Lie algebra \(\widetilde{\overline{\mathbf{v}}}_{\textrm{L}}\) cannot be clearly reconstructed from the unified local velocity coordinates. Two permitted projection examples are the projection towards the lateral surfaces of a hexahedron and towards the edges of a tetrahedron; see Fig. 4. In order to examine the description of the rigid body motion with unified local velocity coordinates more closely in the following subsection, we limit ourselves for the moment to the velocity projection along the side faces of a hexahedron with side lengths (\(2l_{1}\), \(2l _{2}\), \(2l_{3}\)) as depicted in Fig. 4a. The projection vectors \(\bar{\mathbf{b}}_{i}\) and the position vectors \(\bar{\mathbf{y}}_{i}\) read for the hexahedron as depicted in Fig. 4a,

b¯1=[α00]T,y¯1=[0l20]T,
(26)
b¯2=[α00]T,y¯2=[0l20]T,
(27)
b¯3=[0α0]T,y¯3=[00l3]T,
(28)
b¯4=[0α0]T,y¯4=[00l3]T,
(29)
b¯5=[00α]T,y¯5=[l100]T,
(30)
b¯6=[00α]T,y¯6=[l100]T,
(31)

with a quantity \(\alpha \in \mathbb{R}\) to be chosen. Note that other directions of \(\bar{\mathbf{b}}_{i}\) would be possible in the hexahedron as well. In Sect. 4.3, it will be shown that choosing \(\alpha = 1/\sqrt{2}\) is advantageous, which is employed throughout the following description.

Fig. 4
figure 4

Basic geometries for unified local velocity coordinates. \(S\) denotes the center of mass of the rigid body. The six points, used for the velocity projection are denoted with 1 … 6. The corresponding six unified velocities are denoted by \(\overline{w}_{1} \ldots \overline{w}_{6}\), which represent the local velocities projected into the indicated directions. Figure (a) shows the hexahedral case, Fig. (b) shows the tetrahedral case

3.1.1 Velocity field

In order to relate the velocity field \(\overline{\mathbf{v}}(\bar{ \mathbf{y}}_{p})\) to the unified velocities \(\overline{w}_{i}\), we use Eq. (24) with the quantities defined in Eqs. (26)–(31) and assume a linear interpolation of the velocity field in the rigid body. The local velocity vector \(\overline{\mathbf{v}}(\bar{ \mathbf{y}}_{p})= \begin{bmatrix} \bar{v}_{1} & \bar{v}_{2} & \bar{v}_{3} \end{bmatrix} ^{T}\) of each point \(P\) on the rigid body is thus written as

$$ \overline{\mathbf{v}}(\bar{\mathbf{y}}_{p}) = \overline{\mathbf{N}}(\bar{ \mathbf{y}}_{p})\mathbf{\overline{w}}= \begin{bmatrix} \overline{\mathbf{N}}_{1}^{1\times 6}(\bar{\mathbf{y}}_{p}) \\ \overline{\mathbf{N}}_{2}^{1\times 6}(\bar{\mathbf{y}}_{p}) \\ \overline{\mathbf{N}}_{3}^{1\times 6}(\bar{\mathbf{y}}_{p}) \end{bmatrix} \begin{bmatrix} \overline{w}_{1} & \overline{w}_{2} & \overline{w}_{3} & \overline{w} _{4} & \overline{w}_{5} & \overline{w}_{6} \end{bmatrix} ^{T} $$
(32)

with the interpolation matrix \(\overline{\mathbf{N}}(\bar{\mathbf{y}} _{p})\). In order to determine the interpolation matrix \(\overline{ \mathbf{N}}(\bar{\mathbf{y}}_{p})\), the component \(\bar{v}_{1}\) is investigated, exemplary.

The velocity component \(\bar{v}_{1}\) is composed of three parts, which result from different forms of motion of the rigid body. Namely a translational velocity \(\bar{v}_{1}^{\mathrm{transl}}\), which occurs due to a pure translational motion along \(\mathbf{e}_{1}\), a velocity \(\bar{v}_{1}^{\mathbf{e}_{2}}\), which results from a pure rotary motion around the axis \(\mathbf{e}_{2}\), as well as a velocity \(\bar{v}_{1} ^{\mathbf{e}_{3}}\), which results from a rotational motion around the axis given by \(\mathbf{e}_{3}\); see Fig. 5. In the case of a pure translational motion, the velocity of any reference point on the rigid body reads according to Fig. 5a

$$ \bar{v}_{1}^{\mathrm{transl}}= \frac{1}{\sqrt{2}} (\overline{w} _{1}-\overline{w}_{2} ). $$
(33)

The velocity component \(\bar{v}_{1}^{\mathbf{e}_{3}}\) follows from linear interpolation,

$$ \bar{v}_{1}^{\mathbf{e}_{3}}(\bar{y}_{2}) = - \frac{1}{\sqrt{2}} \biggl(\frac{\overline{w}_{1}}{l_{2}} + \frac{\overline{w}_{2}}{l_{2}} \biggr) \bar{y}_{2} $$
(34)

according to Fig. 5b. The negative slope in \(\bar{v}_{1}^{\mathbf{e}_{3}}\) follows from the positive sense of rotation as indicated in Fig. 5b. In the same sense the velocity component \(\bar{v}_{1}^{\mathbf{e}_{2}}\) follows from linear interpolation,

$$ \bar{v}_{1}^{\mathbf{e}_{2}}(\bar{y}_{3}) = \frac{1}{\sqrt{2}} \biggl(\frac{ \overline{w}_{5}}{l_{1}} + \frac{\overline{w}_{6}}{l_{1}} \biggr) \bar{y}_{3} $$
(35)

according to Fig. 6b. Summing up the terms \(\bar{v}_{1}^{\mathrm{transl}}\), \(\bar{v}_{1} ^{\mathbf{e}_{2}}\) and \(\bar{v}_{1}^{\mathbf{e}_{3}}\) results in the first line of Eq. (32)

$$ \bar{v}_{1} = \bar{v}_{1}^{\mathrm{transl}}+ \bar{v}_{1}^{\mathbf{e} _{2}}+ \bar{v}_{1}^{\mathbf{e}_{3}}= \frac{1}{\sqrt{2}} (\overline{w} _{1}-\overline{w}_{2} ) + \frac{1}{\sqrt{2}} \biggl(\frac{ \overline{w}_{5}}{l_{1}} + \frac{\overline{w}_{6}}{l_{1}} \biggr) \bar{y}_{3} - \frac{1}{\sqrt{2}} \biggl( \frac{\overline{w}_{1}}{l _{2}} + \frac{\overline{w}_{2}}{l_{2}} \biggr)\bar{y}_{2}. $$
(36)

Summing up, \(\bar{v}_{1}\) can be represented in vector notation as

$$\begin{aligned} \bar{v}_{1} &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1-\frac{\bar{y}_{2}}{l_{2}} & -1-\frac{\bar{y}_{2}}{l_{2}} & 0 & 0 & \frac{ \bar{y}_{3}}{l_{1}} & \frac{\bar{y}_{3}}{l_{1}} \end{bmatrix} \begin{bmatrix} \overline{w}_{1} & \overline{w}_{2} & \overline{w}_{3} & \overline{w} _{4} & \overline{w}_{5} & \overline{w}_{6} \end{bmatrix} ^{T} \end{aligned}$$
(37)
$$\begin{aligned} &= \overline{\mathbf{N}}_{1}^{1\times 6}(\bar{ \mathbf{y}}_{p}) \begin{bmatrix} \overline{w}_{1} & \overline{w}_{2} & \overline{w}_{3} & \overline{w} _{4} & \overline{w}_{5} & \overline{w}_{6} \end{bmatrix} ^{T} \end{aligned}$$
(38)

in which \(\overline{\mathbf{N}}_{1}^{1\times 6}(\bar{\mathbf{y}}_{p})\) is defined in Eq. (32). Of course, these relations would also follow from basic kinematic formulas using the angular velocity vector, which is not needed here. The remaining components of \(\overline{\mathbf{N}}(\bar{\mathbf{y}}_{p})\) follow from similar considerations and read as follows:

$$\begin{aligned} \overline{\mathbf{N}}(\bar{\mathbf{y}}_{p}) = \frac{1}{\sqrt{2}} \begin{bmatrix} 1-\frac{\bar{y}_{2}}{l_{2}} & -1-\frac{\bar{y}_{2}}{l_{2}} & 0 & 0 & \frac{ \bar{y}_{3}}{l_{1}} & \frac{\bar{y}_{3}}{l_{1}} \\ \frac{\bar{y}_{1}}{l _{2}} & \frac{\bar{y}_{1}}{l_{2}} & 1-\frac{\bar{y}_{3}}{l_{3}} & -1-\frac{ \bar{y}_{3}}{l_{3}} & 0 & 0 \\ 0 & 0 & \frac{\bar{y}_{2}}{l_{3}} & \frac{ \bar{y}_{2}}{l_{3}} & 1-\frac{\bar{y}_{1}}{l_{1}} & -1-\frac{\bar{y} _{1}}{l_{1}} \\ \end{bmatrix} . \end{aligned}$$
(39)
Fig. 5
figure 5

Representation of the motion of the cross-section spanned by the transformation points (1,2,5,6); see also Fig. 4a. Since only vectors in the global frame can be drawn in an image, in (a) the vectors \(\mathbf{b}_{1}\) and \(\mathbf{b}_{2}\) are represented in the global frame, which are related to the vectors \(\bar{\mathbf{b}}_{1}\) and \(\bar{\mathbf{b}}_{2}\) represented in the local frame by \(\mathbf{b} _{i}=\mathbf{R}\bar{\mathbf{b}}_{i}\)

Fig. 6
figure 6

Representation of the motion of the cross-section spanned by the transformation points \((3,4,5,6)\); see also Fig. 4a. Since only vectors in the global frame can be drawn in an image. In (a) the vectors \(\mathbf{b}_{5}\) and \(\mathbf{b}_{6}\) are represented in the global frame, which are related to the vectors \(\bar{\mathbf{b}}_{5}\) and \(\bar{\mathbf{b}}_{6}\) represented in the local frame by \(\mathbf{b} _{i}=\mathbf{R}\bar{\mathbf{b}}_{i}\)

3.2 Representation of accelerations

A unification of the kinematic description of the rigid body motion at velocity level is possible with the help of the six local velocity coordinates \(\overline{w}_{i}\) as shown in the previous section. The accelerations of the rigid body are determined based on the time derivative of the velocity field \(\overline{\mathbf{v}}(\bar{ \mathbf{y}}_{p})\). In order to do so, the velocity field \(\overline{ \mathbf{v}}\) in Eq. (32) is transformed into the inertial frame,

$$ \mathbf{v}(\bar{\mathbf{y}}_{p}) = \bar{v}_{1}\mathbf{e}_{1}+ \bar{v} _{2} \mathbf{e}_{2}+ \bar{v}_{3}\mathbf{e}_{3}= \bigl(\overline{ \mathbf{N}}_{1}(\bar{\mathbf{y}}_{p}) \mathbf{\overline{w}} \bigr) \mathbf{e}_{1}+ \bigl(\overline{ \mathbf{N}}_{2}(\bar{\mathbf{y}}_{p}) \mathbf{ \overline{w}} \bigr)\mathbf{e}_{2}+ \bigl(\overline{ \mathbf{N}}_{3}(\bar{\mathbf{y}}_{p})\mathbf{ \overline{w}} \bigr) \mathbf{e}_{3}, $$
(40)

where \((\mathbf{e}_{1}, \mathbf{e}_{2}, \mathbf{e}_{3})\) represent the three body-fixed orthogonal basis vectors, with their (time dependent) coordinates represented in the inertial frame. The acceleration vector of any point \(P\) in the current configuration can be determined from the time derivative of Eq. (40). Therefore, the acceleration vector follows:

$$\begin{aligned} \dot{\mathbf{v}}(\bar{\mathbf{y}}_{P})= &\,\bigl( \overline{\mathbf{N}} _{1}(\bar{\mathbf{y}}_{p})\dot{ \mathbf{\overline{w}}} \bigr) \mathbf{e}_{1}+ \bigl(\overline{ \mathbf{N}}_{2}(\bar{\mathbf{y}}_{p}) \dot{\mathbf{ \overline{w}}} \bigr)\mathbf{e}_{2}+ \bigl(\overline{ \mathbf{N}}_{3}(\bar{\mathbf{y}}_{p})\dot{\mathbf{ \overline{w}}} \bigr) \mathbf{e}_{3} \\ &{}+ \bigl(\overline{\mathbf{N}}_{1}(\bar{\mathbf{y}}_{p}) \mathbf{\overline{w}} \bigr)\dot{\mathbf{e}}_{1}+ \bigl(\overline{ \mathbf{N}}_{2}(\bar{\mathbf{y}}_{p})\mathbf{ \overline{w}} \bigr) \dot{\mathbf{e}}_{2}+ \bigl(\overline{ \mathbf{N}}_{3}(\bar{ \mathbf{y}}_{p})\mathbf{ \overline{w}} \bigr)\dot{\mathbf{e}}_{3}, \end{aligned}$$
(41)

where \((\dot{\mathbf{e}}_{1}, \,\, \dot{\mathbf{e}}_{2}, \,\, \dot{\mathbf{e}}_{3})\) represent the time derivatives of the three body-fixed orthogonal basis vectors. Formulated in matrix-vector notation, the acceleration vector (41) takes the form

$$\begin{aligned} \dot{\mathbf{v}}(\bar{\mathbf{y}}_{P})=&\, [ \mathbf{e}_{1} \quad \mathbf{e}_{2}\quad \mathbf{e}_{3} ] \begin{bmatrix} \overline{\mathbf{N}}_{1}^{T}(\bar{\mathbf{y}}_{p}) & \overline{ \mathbf{N}}_{2}^{T}(\bar{\mathbf{y}}_{p}) & \overline{\mathbf{N}}_{3} ^{T}(\bar{\mathbf{y}}_{p}) \end{bmatrix} ^{T}\dot{\mathbf{\overline{w}}} + \cdots \\ &{} + [\dot{\mathbf{e}}_{1}\quad \dot{ \mathbf{e}}_{2}\quad \dot{\mathbf{e}}_{3} ] \begin{bmatrix} \overline{\mathbf{N}}_{1}^{T}(\bar{\mathbf{y}}_{p}) & \overline{ \mathbf{N}}_{2}^{T}(\bar{\mathbf{y}}_{p}) & \overline{\mathbf{N}}_{3} ^{T}(\bar{\mathbf{y}}_{p}) \end{bmatrix} ^{T}\mathbf{\overline{w}}. \end{aligned}$$
(42)

The components of the acceleration vector represented in the local frame read as follows:

$$\begin{aligned} \overline{\dot{\mathbf{v}}}(\bar{\mathbf{y}}_{P})=&\, [ \overline{ \mathbf{e}}_{1}\quad \overline{\mathbf{e}}_{2} \quad \overline{ \mathbf{e}}_{3} ] \begin{bmatrix} \overline{\mathbf{N}}_{1}^{T}(\bar{\mathbf{y}}_{p}) & \overline{ \mathbf{N}}_{2}^{T}(\bar{\mathbf{y}}_{p}) & \overline{\mathbf{N}}_{3} ^{T}(\bar{\mathbf{y}}_{p}) \end{bmatrix} ^{T}\dot{\mathbf{\overline{w}}} + \cdots \\ &{} + [\overline{\dot{\mathbf{e}}}_{1}\quad \overline{ \dot{\mathbf{e}}}_{2}\quad \overline{\dot{\mathbf{e}}}_{3} ] \begin{bmatrix} \overline{\mathbf{N}}_{1}^{T}(\bar{\mathbf{y}}_{p}) & \overline{ \mathbf{N}}_{2}^{T}(\bar{\mathbf{y}}_{p}) & \overline{\mathbf{N}}_{3} ^{T}(\bar{\mathbf{y}}_{p}) \end{bmatrix} ^{T}\mathbf{ \overline{w}}, \end{aligned}$$
(43)

where \((\overline{\mathbf{e}}_{1}, \,\, \overline{\mathbf{e}}_{2}, \, \, \overline{\mathbf{e}}_{3})\) represent the three body-fixed basis vectors, with their (time dependent) coordinates represented in the local frame and \((\overline{\dot{\mathbf{e}}}_{1}, \,\, \overline{ \dot{\mathbf{e}}}_{2}, \,\, \overline{\dot{\mathbf{e}}}_{3})\) represent the time derivatives of the three body-fixed basis vectors, with their (time dependent) coordinates represented in the local frame. In the body-fixed frame it follows that in the case of a translational motion \(\overline{\dot{\mathbf{e}}}_{i}\) vanishes. But in the case of a rotational motion it follows that \(\overline{\dot{\mathbf{e}}}_{i} \neq \mathbf{0}\). Since \(\overline{\dot{\mathbf{e}}}_{i}\) can be interpreted as difference of two velocities, Eq. (32) is used to obtain, e.g.,

$$ \overline{\dot{\mathbf{e}}}_{i} = \overline{ \mathbf{v}}(\overline{ \mathbf{e}}_{i}) - \overline{\mathbf{v}}( \mathbf{0}). $$
(44)

The velocity term \(\overline{\mathbf{v}}(\mathbf{0})\) in Eq. (44) corresponds to a pure translational motion at \(\bar{\mathbf{y}}_{p}=\mathbf{0}_{3\times 1}\). Equation (44) can be evaluated using Eq. (32). By inserting the Cartesian basis vectors

$$ \overline{\mathbf{e}}_{1}= \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} ^{T}, \qquad \overline{ \mathbf{e}}_{2}= \begin{bmatrix} 0 & 1 & 0 \end{bmatrix} ^{T}, \qquad \overline{\mathbf{e}}_{3}= \begin{bmatrix} 0 & 0 & 1 \end{bmatrix} ^{T}, $$
(45)

into Eq. (44), it follows that the time derivatives of the basis vectors read

$$\begin{aligned} \overline{\dot{\mathbf{e}}}_{1} &= \begin{bmatrix} 0 & \frac{\overline{w}_{1} + \overline{w}_{2}}{\sqrt{2}l_{2}} & \frac{- \overline{w}_{5} - \overline{w}_{6}}{\sqrt{2}l_{1}} \end{bmatrix} ^{T}, \end{aligned}$$
(46)
$$\begin{aligned} \overline{\dot{\mathbf{e}}}_{2} &= \begin{bmatrix} \frac{-\overline{w}_{1} - \overline{w}_{2}}{\sqrt{2}l_{2}} & 0 & \frac{ \overline{w}_{3} + \overline{w}_{4}}{\sqrt{2}l_{3}} \end{bmatrix} ^{T}, \end{aligned}$$
(47)
$$\begin{aligned} \overline{\dot{\mathbf{e}}}_{3} &= \begin{bmatrix} \frac{\overline{w}_{5} + \overline{w}_{6}}{\sqrt{2}l_{1}} & \frac{- \overline{w}_{3} - \overline{w}_{4}}{\sqrt{2}l_{3}} & 0 \end{bmatrix} ^{T}. \end{aligned}$$
(48)

To obtain a compact vector-matrix notation of Eq. (43), the time derivatives of the basis vectors are collected in the matrixFootnote 1\(\overline{ \dot{\mathbf{E}}}\)

$$ \overline{\dot{\mathbf{E}}}:= \begin{bmatrix} \overline{\dot{\mathbf{e}}}_{1}& \overline{\dot{\mathbf{e}}}_{2}& \overline{ \dot{\mathbf{e}}}_{3} \end{bmatrix} = \renewcommand{\arraystretch}{1.2} \begin{bmatrix} 0 & \frac{-\overline{w}_{1} - \overline{w}_{2}}{\sqrt{2}l_{2}} & \frac{ \overline{w}_{5} + \overline{w}_{6}}{\sqrt{2}l_{1}} \\ \frac{\overline{w}_{1} + \overline{w}_{2}}{\sqrt{2}l_{2}} & 0 & \frac{- \overline{w}_{3} - \overline{w}_{4}}{\sqrt{2}l_{3}} \\ \frac{-\overline{w}_{5} - \overline{w}_{6}}{\sqrt{2}l_{1}} & \frac{ \overline{w}_{3} + \overline{w}_{4}}{\sqrt{2}l_{3}} & 0 \end{bmatrix} . $$
(49)

Inserting Eqs. (47)–(48) and Eq. (45) into Eq. (43), the acceleration vector \(\overline{\dot{\mathbf{v}}}(\bar{\mathbf{y}}_{P})\) of any point \(P\) follows as

$$ \overline{\dot{\mathbf{v}}}(\bar{\mathbf{y}}_{P})= \overline{ \mathbf{N}}(\bar{\mathbf{y}}_{p})\dot{\mathbf{ \overline{w}}} + \overline{ \mathbf{P}}(\bar{\mathbf{y}}_{p}, \mathbf{\overline{w}}) \mathbf{\overline{w}} $$
(50)

where the matrix \(\overline{\mathbf{P}}(\bar{\mathbf{y}}_{p}, \mathbf{\overline{w}})\) is given by

$$ \overline{\mathbf{P}}(\bar{\mathbf{y}}_{p}, \mathbf{\overline{w}})=\overline{ \dot{\mathbf{E}}}(\mathbf{\overline{w}}) \overline{\mathbf{N}}(\bar{ \mathbf{y}}_{p}) $$
(51)

and the time derivative of the unified coordinates reads

$$ \dot{\mathbf{\overline{w}}} := \begin{bmatrix} \dot{\overline{w}}_{1} & \dot{\overline{w}}_{2} & \dot{\overline{w}} _{3} & \dot{\overline{w}}_{4} & \dot{\overline{w}}_{5} & \dot{\overline{w}}_{6} \end{bmatrix} ^{T}. $$
(52)

3.3 Equations of motion of a free rigid body

In this section, we determine the governing equations of motion in terms of \(\mathbf{\overline{w}}\) and \(\dot{\mathbf{\overline{w}}}\). To obtain the governing equations of motion we use the Gibbs–Appel equations [1, 16, 31]. Since we are searching for equations of motion depicted in a body-fixed frame, we have to calculate the partial derivatives of the Gibbs function \(G\) with respect to the time derivative of the velocity coordinates \(\dot{\mathbf{\overline{w}}}\) and equate the derivatives with the corresponding generalized forces \(\bar{\mathbf{Q}}\)

$$ \frac{\partial G}{\partial \dot{\mathbf{\overline{w}}}} = \bar{ \mathbf{Q}} \quad \text{with} \ G = \frac{1}{2} \int _{V}\rho \overline{\dot{\mathbf{v}}}(\bar{ \mathbf{y}}_{P})^{T}\,\overline{\dot{\mathbf{v}}}(\bar{ \mathbf{y}} _{P})\,dV $$
(53)

which are defined in the body-fixed frame. Using Eq. (50), it turns out that the Gibbs function \(G\) can be written as a sum of three integrals,

$$\begin{aligned} G &= \frac{1}{2} \int _{V}\rho \overline{\dot{\mathbf{v}}}(\bar{ \mathbf{y}}_{P})^{T}\,\overline{\dot{\mathbf{v}}}(\bar{ \mathbf{y}} _{P})\,dV \end{aligned}$$
(54)
$$\begin{aligned} &=\underbrace{\frac{1}{2} \int _{V}\rho \dot{\mathbf{\overline{w}}} ^{T} \overline{\mathbf{N}}^{T}\overline{\mathbf{N}} \dot{\mathbf{ \overline{w}}}dV}_{G_{1}} + \underbrace{ \int _{V}\rho \dot{\mathbf{\overline{w}}}^{T} \overline{\mathbf{N}}^{T}\overline{ \mathbf{P}}\mathbf{ \overline{w}}dV}_{G_{2}} + \underbrace{\frac{1}{2} \int _{V}\rho \mathbf{\overline{w}}^{T} \overline{\mathbf{P}}^{T}\overline{ \mathbf{P}}\mathbf{ \overline{w}}dV}_{G_{3}} \end{aligned}$$
(55)
$$\begin{aligned} &= G_{1} + G_{2} + G_{3}. \end{aligned}$$
(56)

Since the vector \(\dot{\mathbf{\overline{w}}}\) does not appear in the term \(G_{3}\), the partial derivative of \(G_{3}\) with respect to \(\dot{\mathbf{\overline{w}}}\) vanishes. Therefore, we only make use of the first two terms \(G_{1}\) and \(G_{2}\), since their partial derivative with respect to \(\dot{\mathbf{\overline{w}}}\) is obviously nonzero. Therefore, \(G_{1}\) and \(G_{2}\) follow to

$$\begin{aligned} G_{1} &= \frac{1}{2}\dot{\mathbf{ \overline{w}}}^{T} \biggl( \int _{V} \rho \overline{\mathbf{N}}^{T}(\bar{ \mathbf{y}}_{p})\overline{ \mathbf{N}}(\bar{\mathbf{y}}_{p}) \,dV \biggr) \dot{\mathbf{\overline{w}}} = \frac{1}{2}\dot{\mathbf{ \overline{w}}} ^{T}\overline{\mathbf{M}}\dot{\mathbf{\overline{w}}}, \end{aligned}$$
(57)
$$\begin{aligned} G_{2} &= \dot{\mathbf{\overline{w}}}^{T} \biggl( \int _{V}\rho \overline{ \mathbf{N}}^{T}(\bar{ \mathbf{y}}_{p})\overline{\mathbf{P}}(\bar{ \mathbf{y}}_{p}, \mathbf{\overline{w}}) \,dV \biggr) \mathbf{\overline{w}}= \dot{\mathbf{ \overline{w}}}^{T}\overline{ \boldsymbol{\Gamma }}\mathbf{ \overline{w}}. \end{aligned}$$
(58)

A further evaluation of the integrals in \(G_{1}\) and \(G_{2}\) lead us to the definition of a mass matrix \(\overline{\mathbf{M}}\) and a velocity dependent matrix \(\overline{\boldsymbol{\Gamma }}\)

$$\begin{aligned} \overline{\mathbf{M}} &:= \int _{V}\rho \overline{\mathbf{N}}^{T}(\bar{ \mathbf{y}}_{p})\overline{\mathbf{N}}(\bar{\mathbf{y}}_{p}) \,dV, \end{aligned}$$
(59)
$$\begin{aligned} \overline{\boldsymbol{\Gamma }} &:= \int _{V}\rho \overline{\mathbf{N}} ^{T}(\bar{ \mathbf{y}}_{p})\overline{\mathbf{P}}(\bar{\mathbf{y}}_{p}, \mathbf{\overline{w}}) \,dV. \end{aligned}$$
(60)

Since the interpolation matrix \(\overline{\mathbf{N}}( \bar{\mathbf{y}}_{p})\) depends only on the chosen reference point \(\bar{\mathbf{y}}_{p}\), which is constant in the body-fixed frame, it follows that the mass matrix \(\overline{\mathbf{M}}\) is also constant with respect to the body-fixed frame. The partial derivatives of \(G_{1}\) and \(G_{2}\) follow to

$$\begin{aligned} \frac{\partial G_{1}}{\partial \dot{\mathbf{\overline{w}}}} = \overline{ \mathbf{M}}\dot{\mathbf{ \overline{w}}}\quad \text{and}\quad \frac{ \partial G_{2}}{\partial \dot{\mathbf{\overline{w}}}} &= \overline{ \boldsymbol{\Gamma }}\mathbf{\overline{w}}. \end{aligned}$$
(61)

Adding up the partial derivatives of \(G_{1}\) and \(G_{2}\), the equations of motion follow to

$$\begin{aligned} \overline{\mathbf{M}}\dot{\mathbf{\overline{w}}} + \overline{ \boldsymbol{\Gamma }}\mathbf{\overline{w}}= \overline{\mathbf{Q}}. \end{aligned}$$
(62)

Note that Eqs. (62) have been obtained without the use of rotation matrices or the angular velocity vector.

3.4 Generalized forces

Corresponding to the equations of motion (62), we define generalized forces \(\overline{\mathbf{Q}}\in \mathbb{R}^{6}\), which are obtained from the principle of virtual power [29]. Since the principle of virtual power is invariant under a coordinate transformation, it is written in the local (body-fixed) representation,

$$ \overline{\mathbf{Q}}^{T}\delta \mathbf{ \overline{w}}= \overline{ \mathbf{f}}^{T}\delta \overline{\mathbf{v}}( \bar{\mathbf{y}}_{p}) = \overline{ \mathbf{f}}^{T} \biggl(\frac{\partial \overline{\mathbf{v}}(\bar{ \mathbf{y}}_{p})}{\partial \mathbf{\overline{w}}}\delta \mathbf{\overline{w}} \biggr)= \overline{ \mathbf{f}}^{T} \biggl(\frac{ \partial }{\partial \mathbf{\overline{w}}} \bigl(\overline{ \mathbf{N}}(\bar{ \mathbf{y}}_{p})\mathbf{\overline{w}} \bigr)\delta \mathbf{\overline{w}} \biggr) = \overline{\mathbf{f}}^{T}\overline{ \mathbf{N}}(\bar{\mathbf{y}}_{p})\delta \mathbf{\overline{w}}. $$
(63)

For the case of a force \(\overline{\mathbf{f}}\in \mathbb{R}^{3}\) acting at a point \(P\) at the rigid body, the vector of generalized forces \(\overline{\mathbf{Q}}\) thus reads

$$ \overline{\mathbf{Q}} = \overline{\mathbf{N}}^{T}( \bar{\mathbf{y}} _{p})\,\overline{\mathbf{f}}. $$
(64)

4 Relation to rotation dependent rigid body motion

In this section we show the bijective mapping between the unified local velocity coordinates and the classical translational and angular velocity vectors. Furthermore, the relations between the equations of motion (62) and the Newton–Euler equations, as well as their respective mass matrices and inertia terms are presented.

4.1 Relation between \(\mathbf{\overline{w}}\) and \(\overline{\mathbf{v}}_{\textrm{L}}\)

In this section, a relationship between the unified local velocity coordinates \(\mathbf{\overline{w}}\) and the vector form of \(\widetilde{\overline{\mathbf{v}}}_{\textrm{L}}\) is derived.

It is well known that the local velocity vector of a point \(P\) on the rigid body can be computed using the basic kinematic formula

$$ \overline{\mathbf{v}}(\bar{\mathbf{y}}_{P}) = \overline{ \mathbf{U}}+ \overline{ \boldsymbol{\Omega }}\times \bar{\mathbf{y}}_{P} , $$
(65)

which can be rewritten in terms of \(\overline{\mathbf{v}}_{\textrm{L}}\) as follows:

$$ \overline{\mathbf{v}}(\bar{\mathbf{y}}_{P}) = \overline{\mathbf{U}}+ \overline{ \boldsymbol{\Omega }}\times \bar{ \mathbf{y}}_{P} = \overline{\mathbf{U}}+ \widetilde{\bar{ \mathbf{y}}}_{P}^{T}\overline{\boldsymbol{\Omega }}= \begin{bmatrix} \mathbf{I}_{3\times 3}& \widetilde{\bar{\mathbf{y}}}_{P}^{T} \end{bmatrix} \begin{bmatrix} \,\overline{\mathbf{U}} \\ \overline{\boldsymbol{\Omega }}\,\end{bmatrix} = \begin{bmatrix} \mathbf{I}_{3\times 3}& \widetilde{\bar{\mathbf{y}}}_{P}^{T} \end{bmatrix} \overline{ \mathbf{v}}_{\textrm{L}}. $$
(66)

By inserting Eq. (66) into Eq. (24) it follows that the \(i\)th unified velocity coordinate \(\overline{w}_{i}\) is related to \(\overline{\mathbf{v}} _{\textrm{L}}\) by

$$ \overline{w}_{i} = \bar{\mathbf{b}}_{i}^{T} \overline{\mathbf{v}}(\bar{ \mathbf{y}}_{i}) = \bar{ \mathbf{b}}_{i}^{T} \begin{bmatrix} \mathbf{I}_{3\times 3}& \widetilde{\bar{\mathbf{y}}}_{i}^{T} \end{bmatrix} \overline{\mathbf{v}}_{\textrm{L}}= \begin{bmatrix} \bar{\mathbf{b}}_{i}^{T} & (\tilde{\bar{\mathbf{y}}}_{i}\bar{ \mathbf{b}}_{i} )^{T} \end{bmatrix} \overline{\mathbf{v}}_{\textrm{L}}. $$
(67)

From the selection of six projection points \(P_{i}\) and direction vectors \(\bar{\mathbf{b}}_{i}\) it follows that the vector of the unified velocity coordinates \(\mathbf{\overline{w}}\) is related to \(\overline{ \mathbf{v}}_{\textrm{L}}\) via the relation

$$ \mathbf{\overline{w}}= \begin{bmatrix} \overline{w}_{1} \\ \overline{w}_{2} \\ \overline{w}_{3} \\ \overline{w}_{4} \\ \overline{w}_{5} \\ \overline{w}_{6} \end{bmatrix} = \begin{bmatrix} \bar{\mathbf{b}}_{1}^{T} & (\tilde{\bar{\mathbf{y}}}_{1}\bar{ \mathbf{b}}_{1} )^{T} \\ \bar{\mathbf{b}}_{2}^{T} & (\tilde{\bar{\mathbf{y}}}_{2}\bar{ \mathbf{b}}_{2} )^{T} \\ \bar{\mathbf{b}}_{3}^{T} & (\tilde{\bar{\mathbf{y}}}_{3}\bar{ \mathbf{b}}_{3} )^{T} \\ \bar{\mathbf{b}}_{4}^{T} & (\tilde{\bar{\mathbf{y}}}_{4}\bar{ \mathbf{b}}_{4} )^{T} \\ \bar{\mathbf{b}}_{5}^{T} & (\tilde{\bar{\mathbf{y}}}_{5}\bar{ \mathbf{b}}_{5} )^{T} \\ \bar{\mathbf{b}}_{6}^{T} & (\tilde{\bar{\mathbf{y}}}_{6}\bar{ \mathbf{b}}_{6} )^{T} \end{bmatrix} \overline{ \mathbf{v}}_{\textrm{L}}= \overline{\mathbf{D}}\,\overline{ \mathbf{v}}_{\textrm{L}}. $$
(68)

According to Eq. (68) the matrix \(\overline{ \mathbf{D}}\in \mathbb{R}^{6\times 6}\) maps \(\overline{\mathbf{v}} _{\textrm{L}}\) to the vector of the unified velocity coordinates. Therefore, the inverse of \(\overline{\mathbf{D}}\) maps the vector of the unified velocity coordinates to \(\overline{\mathbf{v}}_{\textrm{L}}\) according to

$$ \overline{\mathbf{v}}_{\textrm{L}}= \overline{ \mathbf{D}}^{-1}\, \mathbf{\overline{w}}. $$
(69)

4.2 Condition for the choice of \(\bar{\mathbf{b}}_{i}\) and \(\bar{\mathbf{y}}_{i}\)

In order to obtain a valid set of unified local velocity coordinates \(\overline{w}_{i}\), the following conditions need to be fulfilled for the transformation quantities \(\bar{\mathbf{b}}_{i}\) and \(\bar{ \mathbf{y}}_{i}\):

  1. 1.

    there is a linear mapping between the vector of unified local velocity coordinates \(\mathbf{\overline{w}}\) and the vector \(\overline{ \mathbf{v}}_{\textrm{L}}\) according to Eq. (68), and

  2. 2.

    the matrix \(\overline{\mathbf{D}}\) of the linear mapping is a regular matrix.

4.3 Relation between \(\mathbf{\overline{w}}\) and \(\overline{\mathbf{v}}_{\textrm{L}}\) for the hexahedral case

Using the direction and position vectors selected in Eqs. (26)–(31) and setting \(l_{1}=l_{2}= l_{3}=1\), which is employed throughout the following description, the relation between \(\mathbf{\overline{w}}\) and \(\overline{\mathbf{v}}_{\textrm{L}}\) reads for the hexahedral case

$$ \mathbf{\overline{w}}= \begin{bmatrix} \overline{w}_{1} \\ \overline{w}_{2} \\ \overline{w}_{3} \\ \overline{w}_{4} \\ \overline{w}_{5} \\ \overline{w}_{6} \end{bmatrix} = \begin{bmatrix} \bar{\mathbf{b}}_{1}^{T} & (\tilde{\bar{\mathbf{y}}}_{1}\bar{ \mathbf{b}}_{1} )^{T} \\ \bar{\mathbf{b}}_{2}^{T} & (\tilde{\bar{\mathbf{y}}}_{2}\bar{ \mathbf{b}}_{2} )^{T} \\ \bar{\mathbf{b}}_{3}^{T} & (\tilde{\bar{\mathbf{y}}}_{3}\bar{ \mathbf{b}}_{3} )^{T} \\ \bar{\mathbf{b}}_{4}^{T} & (\tilde{\bar{\mathbf{y}}}_{4}\bar{ \mathbf{b}}_{4} )^{T} \\ \bar{\mathbf{b}}_{5}^{T} & (\tilde{\bar{\mathbf{y}}}_{5}\bar{ \mathbf{b}}_{5} )^{T} \\ \bar{\mathbf{b}}_{6}^{T} & (\tilde{\bar{\mathbf{y}}}_{6}\bar{ \mathbf{b}}_{6} )^{T} \end{bmatrix} \overline{ \mathbf{v}}_{\textrm{L}}= \renewcommand{\arraystretch}{1.2} \begin{bmatrix} \frac{1}{\sqrt{2}} & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2}} \\ \frac{-1}{\sqrt{2}} & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2}} \\ 0 & \frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} & 0 & 0 \\ 0 & \frac{-1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} & 0 & 0 \\ 0 & 0 & \frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} & 0 \\ 0 & 0 & \frac{-1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} & 0 \\ \end{bmatrix} \overline{\mathbf{v}}_{\textrm{L}}= \overline{\mathbf{D}}\,\overline{ \mathbf{v}}_{\textrm{L}}, $$
(70)

where \(\overline{\mathbf{D}}\) turns out to be an orthogonal matrix, i.e. \(\overline{ \mathbf{D}}^{-1}=\overline{\mathbf{D}}^{T}\), which can be validated by direct evaluation of \(\overline{\mathbf{D}}^{T} \overline{\mathbf{D}}\). Thus, the conditions defined in Sect. 4.2 are obviously fulfilled. Finally, \(\overline{\mathbf{v}}_{\textrm{L}}\) is given in terms of \(\mathbf{\overline{w}}\), reading

$$ \overline{\mathbf{v}}_{\textrm{L}}:= \begin{bmatrix} \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \\ 0 & 0 & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2} } & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2} } & \frac{1}{\sqrt{2} } & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \overline{w}_{1} \\ \overline{w}_{2} \\ \overline{w}_{3} \\ \overline{w}_{4} \\ \overline{w}_{5} \\ \overline{w}_{6} \end{bmatrix} = \overline{\mathbf{D}}^{-1}\,\mathbf{\overline{w}}. $$
(71)

Note, the geometric interpretation of \(\overline{\mathbf{D}}\) is described in the components of \(\mathbf{\overline{w}}\); see Sect. 3.1.

Looking at the determinant of the matrix \(\overline{\mathbf{D}}\) for the general hexahedron case Eqs. (26)–(31), it follows that the conditions in section Sect. 3.1 are fulfilled for each non-trivial choice of the parameter \(\alpha \) and the lengths \((l_{1},\,\ l_{2},\,\,l_{3})\), due to the fact that

$$ \textrm{det}(\overline{\mathbf{D}}) = -8\, \alpha ^{6} l_{1} l_{2} l _{3}. $$
(72)

If the parameter \(\alpha =1 / \sqrt{2}\) is selected, the conditions in Sect. 3.1 are fulfilled for each non-trivial choice of lengths \((l_{1},\,\ l_{2},\,\,l_{3})\), due to the fact that

$$ \textrm{det}(\overline{\mathbf{D}}) = -l_{1} l_{2} l_{3}. $$
(73)

4.4 Relation between \(\mathbf{\overline{w}}\) and \(\overline{\mathbf{v}}_{\textrm{L}}\) for the tetrahedral case

In this subsection we present a possible selection of direction vectors \(\bar{\mathbf{b}}_{i}\) and position vectors \(\bar{\mathbf{y}}_{i}\), the resulting matrix \(\overline{\mathbf{D}}_{\mathit{tet}}\), as well as their inverse \(\overline{\mathbf{D}}_{\mathit{tet}}^{-1}\) for the velocity projection along the edges of a tetrahedron. For simplicity we choose a tetrahedron as shown in Fig. 7. The vectors \(( \mathbf{e}_{1},\,\, \mathbf{e}_{2},\,\, \mathbf{e}_{3})\) represent the orthogonal basis as shown in Eq. (45). We select the coordinates for each node \((i, \,\, j, \,\ k)\) as \(l\mathbf{e} _{i}\) with respect to the origin given by \(S\) with \(l\) being a coefficient, which scales the size of the tetrahedral. The position and direction vectors therefore result in

b¯1=[α00]T,y¯1=[l200]T,
(74)
b¯2=[0α0]T,y¯2=[0l20]T,
(75)
b¯3=[00α]T,y¯3=[00l2]T,
(76)
b¯4=[αα0]T,y¯4=[l2l20]T,
(77)
b¯5=[0αα]T,y¯5=[0l2l2]T,
(78)
b¯6=[α0α]T,y¯6=[l20l2]T,
(79)

with \(\alpha \) being a coefficient, which scales the direction vectors \(\bar{\mathbf{b}}_{i}\). The matrix \(\overline{\mathbf{D}}_{ \mathit{tet}}\) and its inverse \(\overline{\mathbf{D}}^{-1}_{ \mathit{tet}}\) therefore take the form

$$ \overline{\mathbf{D}}_{\mathit{tet}} = \renewcommand{\arraystretch}{1.2} \begin{bmatrix} \alpha & 0 & 0 & 0 & 0 & 0 \\ 0 & \alpha & 0 & 0 & 0 & 0 \\ 0 & 0 & \alpha & 0 & 0 & 0 \\ -\alpha & \alpha & 0 & 0 & 0 & l\alpha \\ 0 & -\alpha & \alpha & l\alpha & 0 & 0 \\ \alpha & 0 & -\alpha & 0 & l\alpha & 0 \\ \end{bmatrix} , \qquad \overline{ \mathbf{D}}_{\mathit{tet}}^{-1} = \renewcommand{\arraystretch}{1.2} \begin{bmatrix} \frac{1}{\alpha } & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{\alpha } & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{\alpha } & 0 & 0 & 0 \\ 0 & \frac{1}{l\alpha } & -\frac{1}{l\alpha } & 0 & \frac{1}{l\alpha } & 0 \\ -\frac{1}{l\alpha } & 0 & \frac{1}{l\alpha } & 0 & 0 & \frac{1}{l \alpha } \\ \frac{1}{l\alpha } & -\frac{1}{l\alpha } & 0 & \frac{1}{l\alpha } & 0 & 0 \\ \end{bmatrix} . $$
(80)

In contrast to the projection along the side surfaces of a hexahedron, the matrix \(\overline{\mathbf{D}}_{\mathit{tet}}\) is in general not orthogonal for the tetrahedral case.

Fig. 7
figure 7

Tetrahedral geometry for unified local velocity coordinates. \(S\) denotes the center of mass of the rigid body. The six points, used for the velocity projection are denoted with \(1 \ldots 6\). The corresponding six unified velocity coordinates are denoted by \(\overline{w}_{1} \ldots \overline{w}_{6}\)

Looking at the determinant of the matrix \(\overline{\mathbf{D}}_{ \mathit{tet}}\) for the tetrahedral case Eqs. (74)–(79), it follows that the conditions in section Sect. 3.1 are fulfilled for each non-trivial choice of the parameter \(\alpha \) and the length \(l\), due to the fact that

$$ \textrm{det}(\overline{\mathbf{D}}_{\mathit{tet}}) = \alpha ^{6} l^{3}. $$
(81)

4.5 Mass matrix

In order to obtain a more suitable form of the mass matrix \(\overline{ \mathbf{M}}\), the integral in Eq. (59) is investigated in this section.

By inserting equation (71) into equation (66), it follows that the velocity field \(\overline{\mathbf{v}}(\bar{\mathbf{y}}_{p})\) can be represented as

$$ \overline{\mathbf{v}}(\bar{\mathbf{y}}_{p}) = \begin{bmatrix} \mathbf{I}_{3\times 3}& \widetilde{\bar{\mathbf{y}}}_{p}^{T} \end{bmatrix} \overline{\mathbf{D}}^{-1} \mathbf{\overline{w}}. $$
(82)

By comparing Eq. (82) with Eq. (32), it follows that

$$ \overline{\mathbf{N}}(\bar{\mathbf{y}}_{p}) = \begin{bmatrix} \mathbf{I}_{3\times 3}& \widetilde{\bar{\mathbf{y}}}_{p}^{T} \end{bmatrix} \overline{\mathbf{D}}^{-1}. $$
(83)

Using Eq. (83), the mass matrix Eq. (59) is therefore obtained from

$$\begin{aligned} \overline{\mathbf{M}} &= \int _{V}\rho \,\overline{\mathbf{N}}^{T}(\bar{ \mathbf{y}}_{p})\overline{\mathbf{N}}(\bar{\mathbf{y}}_{p}) \,dV = \overline{ \mathbf{D}}^{-T}\left ( \int _{V}\rho \, \begin{bmatrix} \mathbf{I}_{3\times 3}\\ \widetilde{\bar{\mathbf{y}}}_{p} \end{bmatrix} \begin{bmatrix} \mathbf{I}_{3\times 3}& \widetilde{\bar{\mathbf{y}}}_{p}^{T} \end{bmatrix} \,dV\right )\overline{ \mathbf{D}}^{-1} \end{aligned}$$
(84)
$$\begin{aligned} &= \overline{\mathbf{D}}^{-T}\left ( \int _{V}\rho \, \begin{bmatrix} \mathbf{I}_{3\times 3}& \widetilde{\bar{\mathbf{y}}}_{p} \\ \widetilde{\bar{\mathbf{y}}}_{p} & & \widetilde{\bar{\mathbf{y}}}_{p} \widetilde{\bar{\mathbf{y}}}_{p}^{T} \end{bmatrix} \,dV\right )\overline{\mathbf{D}}^{-1} = \overline{ \mathbf{D}}^{-T} \begin{bmatrix} \overline{\mathbf{J}}_{1} & \overline{\mathbf{J}}_{2} \\ \overline{ \mathbf{J}}_{2} & \overline{\mathbf{J}}_{3} \end{bmatrix} \overline{\mathbf{D}}^{-1} \end{aligned}$$
(85)
$$\begin{aligned} &= \overline{\mathbf{D}}^{-T}\mathbf{M}\overline{ \mathbf{D}}^{-1} \end{aligned}$$
(86)

with \(\mathbf{M}\) being defined by

$$ \mathbf{M} := \begin{bmatrix} \,\overline{\mathbf{J}}_{1} & \overline{\mathbf{J}}_{2} \\ \overline{\mathbf{J}}_{2} & \overline{\mathbf{J}}_{3} \end{bmatrix} , $$
(87)

where the matrix \(\overline{\mathbf{J}}_{1}\) is computed by the integral,

$$ \overline{\mathbf{J}}_{1} = \int _{V}\rho \,\mathbf{I}_{3\times 3}\,dV = \begin{bmatrix} m & 0 & 0 \\ 0 & m & 0 \\ 0 & 0 & m \end{bmatrix} , $$
(88)

the matrix \(\overline{\mathbf{J}}_{2}\) by the integral,

$$\begin{aligned} \overline{\mathbf{J}}_{2}(\bar{ \mathbf{y}}_{p}) = \int _{V}\rho \, \tilde{\bar{\mathbf{y}}}_{p}^{T}\, dV, \end{aligned}$$
(89)

and the matrix \(\overline{\mathbf{J}}_{3}\) by the integral

$$\begin{aligned} \overline{\mathbf{J}}_{3}(\bar{ \mathbf{y}}_{p}) = \int _{V}\rho \, \tilde{\bar{\mathbf{y}}}_{p} \tilde{\bar{\mathbf{y}}}_{p}^{T}\, dV = \begin{bmatrix} \overline{J}_{11} & \overline{J}_{12} & \overline{J}_{13} \\ \overline{J}_{21} & \overline{J}_{22} & \overline{J}_{23} \\ \overline{J}_{31} & \overline{J}_{32} & \overline{J}_{33} \end{bmatrix} . \end{aligned}$$
(90)

From the definition of the skew symmetric matrix \(\tilde{\bar{\mathbf{y}}}_{p}\) one concludes that if the origin of the body axes is attached to the center of mass of the body, then the matrix \(\overline{\mathbf{J}}_{2}\) is the null matrix [32], which is employed throughout the further description. The mass matrix thus follows:

$$ \overline{\mathbf{M}}= \overline{\mathbf{D}}^{-T} \begin{bmatrix} \overline{\mathbf{J}}_{1} & \mathbf{0}_{3\times 3} \\ \mathbf{0}_{3\times 3}& \overline{\mathbf{J}}_{3} \end{bmatrix} \overline{\mathbf{D}}^{-1} = \overline{\mathbf{D}}^{-T}\mathbf{M}\overline{ \mathbf{D}}^{-1} $$
(91)

where the diagonal elements of the 3 × 3-matrix \(\overline{ \mathbf{J}}_{1}\) contain the rigid body’s mass \(m\) and the 3 × 3-matrix \(\overline{\mathbf{J}}_{3}\) represents the body’s inertia tensor represented in the body-fixed frame. If we use a diagonal inertia tensor

$$\begin{aligned} \overline{\mathbf{J}}_{3}(\bar{ \mathbf{y}}_{p}) = \begin{bmatrix} \overline{J}_{11} & 0 & 0 \\ 0 & \overline{J}_{22} & 0 \\ 0 & 0 & \overline{J}_{33} \end{bmatrix} , \end{aligned}$$
(92)

the mass matrix \(\overline{\mathbf{M}}\) exhibits in the hexahedral case using Eqs. (26)–(31) and setting \(l_{1}=l_{2}=l_{3}=1\) a simple block diagonal structure,

$$ \overline{\mathbf{M}} = \frac{1}{2} \begin{bmatrix} m+\overline{J}_{33} & \overline{J}_{33}-m & 0 & 0 & 0 & 0 \\ \overline{J}_{33}-m& m+\overline{J}_{33} & 0 & 0 & 0 & 0 \\ 0 & 0 & m+\overline{J}_{11} & \overline{J}_{11}-m & 0 & 0 \\ 0 & 0 & \overline{J}_{11}-m & m+\overline{J}_{11} & 0 & 0 \\ 0 & 0 & 0 & 0 & m+\overline{J}_{22} & \overline{J}_{22}-m \\ 0 & 0 & 0 & 0 & \overline{J}_{22}-m & m+\overline{J}_{22} \\ \end{bmatrix} . $$
(93)

4.6 Quadratic velocity vector

The quadratic velocity vector \(\overline{\boldsymbol{\Gamma }}\, \mathbf{\overline{w}}\) and the integral in Eq. (60) are related to classical inertia terms in the present section.

Inserting Eqs. (83) and (51), into Eq. (60), the matrix \(\overline{ \boldsymbol{\Gamma }}\) results in

$$\begin{aligned} \overline{\boldsymbol{\Gamma }} &= \int _{V}\rho \,\overline{\mathbf{N}} ^{T}(\bar{ \mathbf{y}}_{p})\overline{\dot{\mathbf{E}}}( \mathbf{\overline{w}}) \overline{\mathbf{N}}(\bar{\mathbf{y}}_{p})\, dV \end{aligned}$$
(94)
$$\begin{aligned} &= \overline{\mathbf{D}}^{-T}\left ( \int _{V}\rho \, \begin{bmatrix} \overline{\dot{\mathbf{E}}}(\mathbf{\overline{w}}) & \overline{ \dot{\mathbf{E}}}(\mathbf{\overline{w}})\widetilde{\bar{\mathbf{y}}} _{p} \\ \overline{\dot{\mathbf{E}}}(\mathbf{\overline{w}}) \widetilde{\bar{\mathbf{y}}}_{p} & \overline{\dot{\mathbf{E}}}( \mathbf{\overline{w}})\widetilde{\bar{\mathbf{y}}}_{p} \widetilde{\bar{\mathbf{y}}}_{p}^{T} \end{bmatrix}\, dV\right )\overline{\mathbf{D}}^{-1}. \end{aligned}$$
(95)

Since Eq. (94) is integrated over the volume and the matrix \(\overline{\dot{\mathbf{E}}}(\mathbf{\overline{w}})\) does not depend on the volume, we obtain

$$\begin{aligned} \overline{\boldsymbol{\Gamma }}= \overline{\mathbf{D}}^{-T}\, \overline{ \dot{\mathbf{E}}}(\mathbf{\overline{w}})\left ( \int _{V}\rho \, \begin{bmatrix} \mathbf{I}_{3\times 3}& \widetilde{\bar{\mathbf{y}}}_{p} \\ \widetilde{\bar{\mathbf{y}}}_{p} & \widetilde{\bar{\mathbf{y}}}_{p} \widetilde{\bar{\mathbf{y}}}_{p}^{T} \end{bmatrix}\, dV\right )\overline{\mathbf{D}}^{-1}. \end{aligned}$$
(96)

By comparing the integrals in Eq. (96) with the integrals in Eq. (84) it follows that the matrix \(\overline{\boldsymbol{\Gamma }}\) reads

$$\begin{aligned} \overline{\boldsymbol{\Gamma }}= \overline{\mathbf{D}}^{-T}\, \begin{bmatrix}\, \overline{\dot{\mathbf{E}}}(\mathbf{\overline{w}})\overline{ \mathbf{J}}_{1} & \overline{\dot{\mathbf{E}}}(\mathbf{\overline{w}})\overline{ \mathbf{J}}_{2} \\ \overline{\dot{\mathbf{E}}}(\mathbf{\overline{w}})\overline{ \mathbf{J}}_{2} & \overline{\dot{\mathbf{E}}}(\mathbf{\overline{w}})\overline{ \mathbf{J}}_{3} \end{bmatrix} \overline{\mathbf{D}}^{-1} = \overline{\mathbf{D}}^{-T}\mathbf{\varPsi } \mathbf{M}\overline{ \mathbf{D}}^{-1} \end{aligned}$$
(97)

with the hereby defined matrix \(\mathbf{\varPsi }\),

$$ \mathbf{\varPsi } := \begin{bmatrix} \overline{\dot{\mathbf{E}}}(\mathbf{\overline{w}}) & \mathbf{0}_{3 \times 3} \\ \mathbf{0}_{3\times 3}& \overline{\dot{\mathbf{E}}}( \mathbf{\overline{w}}) \end{bmatrix} . $$
(98)

As described in Sect. 4.5, \(\overline{ \mathbf{J}}_{2}=\mathbf{0}_{3\times 3}\), therefore the matrix \(\overline{\boldsymbol{\Gamma }}\) is computed as

$$ \overline{\boldsymbol{\Gamma }}= \overline{\mathbf{D}}^{-T} \begin{bmatrix} \overline{\dot{\mathbf{E}}}\,\overline{\mathbf{J}}_{1} & \mathbf{0} _{3\times 3} \\ \mathbf{0}_{3\times 3}& \overline{\dot{\mathbf{E}}}\,\overline{ \mathbf{J}}_{3} \end{bmatrix} \overline{\mathbf{D}}^{-1} $$
(99)

Using the diagonal inertia tensor of Eq. (92), the vector \(\overline{\boldsymbol{\Gamma }}\,\mathbf{\overline{w}}\) reads for the hexahedral case using Eqs. (26)–(31) and setting \(l_{1}=l_{2}=l_{3}=1\) as follows:

$$ {\small \overline{\boldsymbol{\Gamma }}\mathbf{ \overline{w}}=\! \frac{1}{2 \sqrt{2}}\! \begin{bmatrix}\! m (\overline{w}_{1}+\overline{w}_{2} ) (\overline{w} _{4}-\overline{w}_{3} )\! + \! ( (-\overline{J}_{11}\!-\! \overline{J}_{22} ) (\overline{w}_{3}+\overline{w}_{4} )\!+\! m (\overline{w}_{5}\!-\!\overline{w}_{6} ) ) (\overline{w} _{5}+\overline{w}_{6} ) \\ m (\overline{w}_{1}+\overline{w}_{2} ) (\overline{w} _{3}-\overline{w}_{4} ) - ( (\overline{J}_{11}- \overline{J}_{22} ) (\overline{w}_{3}+\overline{w}_{4} )+m (\overline{w}_{5}-\overline{w}_{6} ) ) (\overline{w} _{5}+\overline{w}_{6} ) \\ m (\overline{w}_{1}^{2} - \overline{w}_{2}^{2} - (\overline{w} _{3}+\overline{w}_{4} ) (\overline{w}_{5}-\overline{w}_{6} ) ) - (\overline{J}_{22}-\overline{J}_{33} ) (\overline{w} _{1} +\overline{w}_{2} ) (\overline{w}_{5}+\overline{w} _{6} ) \\ -m (\overline{w}_{1}^{2} - \overline{w}_{2}^{2} - (\overline{w} _{3}+\overline{w}_{4} ) (\overline{w}_{5}-\overline{w}_{6} ) ) - (\overline{J}_{22}-\overline{J}_{33} ) (\overline{w} _{1} +\overline{w}_{2} ) (\overline{w}_{5}+\overline{w} _{6} ) \\ m (\overline{w}_{3}^{2}-\overline{w}_{4}^{2}- (\overline{w} _{1}-\overline{w}_{2} ) (\overline{w}_{5}+\overline{w} _{6} ) ) + (\overline{J}_{11}-\overline{J}_{33} ) (\overline{w}_{1}+\overline{w}_{2} ) (\overline{w} _{3}+\overline{w}_{4} ) \\ m (\overline{w}_{4}^{2}-\overline{w}_{3}^{2}+ (\overline{w} _{1}-\overline{w}_{2} ) (\overline{w}_{5}+\overline{w} _{6} ) ) + (\overline{J}_{11}-\overline{J}_{33} ) (\overline{w}_{1}+\overline{w}_{2} ) (\overline{w} _{3}+\overline{w}_{4} )\! \end{bmatrix}\! . } $$
(100)

4.7 Relation to the Newton–Euler equations of motion

In this section we present the relationship between the equations of motion (62) and the Newton–Euler equations of motion.

Inserting Eq. (83) in the equation for the generalized forces Eq. (64) leads to the representation of the generalized forces as a function of the matrix \(\overline{\mathbf{D}}\)

$$ \overline{\mathbf{Q}} = \overline{\mathbf{N}}^{T}(\bar{\mathbf{y}} _{p})\,\overline{\mathbf{f}} = \overline{\mathbf{D}}^{-T} \begin{bmatrix} \mathbf{I}_{3\times 3} \\ \widetilde{\bar{\mathbf{y}}}_{p} \end{bmatrix} \,\overline{\mathbf{f}} = \overline{\mathbf{D}}^{-T} \begin{bmatrix} \overline{\mathbf{f}} \\ \overline{\boldsymbol{\tau }} \end{bmatrix} , $$
(101)

where \(\overline{\boldsymbol{\tau }}=\widetilde{\bar{\mathbf{y}}}_{p}\overline{ \mathbf{f}}\in \mathbb{R}^{3}\) represents the torque acting on the rigid body due to the local force \(\overline{\mathbf{f}}\in \mathbb{R}^{3}\) acting at point \(P\). The local force vector \(\overline{\mathbf{f}}\) is connected to a physical force \(\mathbf{f}\in \mathbb{R}^{3}\) by

$$ \overline{\mathbf{f}} = \mathbf{R}^{T}\mathbf{f}. $$
(102)

Since the mass matrix \(\overline{\mathbf{M}}\) and the velocity dependent matrix \(\overline{\boldsymbol{\Gamma }}\) can be written as matrix product as presented in Eq. (91) and Eq. (99), the equations of motion (62) can be rewritten using the matrix \(\overline{\mathbf{D}}\),

$$ \overline{\mathbf{M}}\dot{\mathbf{\overline{w}}} + \overline{ \boldsymbol{\Gamma }}\mathbf{\overline{w}}= \overline{\mathbf{Q}} \quad \Longrightarrow \quad \overline{\mathbf{D}}^{-T}\mathbf{M}\overline{ \mathbf{D}}^{-1}\dot{\mathbf{\overline{w}}} + \overline{\mathbf{D}} ^{-T}\mathbf{\varPsi }\mathbf{M}\overline{\mathbf{D}}^{-1} \mathbf{\overline{w}}= \overline{\mathbf{D}}^{-T} \begin{bmatrix} \overline{\mathbf{f}} \\ \overline{\boldsymbol{\tau }}. \end{bmatrix} $$
(103)

Multiplication from the left hand side with \(\overline{\mathbf{D}} ^{T}\) gives

$$ \mathbf{M}\overline{\mathbf{D}}^{-1}\dot{\mathbf{ \overline{w}}} + \mathbf{\varPsi }\mathbf{M}\overline{\mathbf{D}}^{-1} \mathbf{\overline{w}}= \begin{bmatrix} \overline{\mathbf{f}} \\ \overline{\boldsymbol{\tau }} \end{bmatrix} . $$
(104)

Since \(\dot{\mathbf{\overline{w}}} =\overline{\mathbf{D}}\, \dot{\overline{\mathbf{v}}}_{L}\) and \(\mathbf{\overline{w}}=\overline{ \mathbf{D}}\,\overline{\mathbf{v}}_{\textrm{L}} \) the equations of motion can be written as

$$ \mathbf{M}\dot{\overline{\mathbf{v}}}_{\textrm{L}} + \mathbf{\varPsi } \mathbf{M}\overline{\mathbf{v}}_{\textrm{L}}= \begin{bmatrix} \overline{\mathbf{f}} \\ \overline{\boldsymbol{\tau }} \end{bmatrix} . $$
(105)

Splitting up the equations of motion (105) into a term for the translational motion and the rotational motion leads to the well known Newton–Euler equations see, e.g. the references [6, 33],

$$\begin{aligned} \overline{\mathbf{J}}_{1}\dot{\overline{ \mathbf{U}}} + \overline{ \mathbf{J}}_{1}\,\widetilde{\overline{ \boldsymbol{\Omega }}}\,\overline{ \mathbf{U}} &= \overline{\mathbf{f}} , \end{aligned}$$
(106)
$$\begin{aligned} \overline{\mathbf{J}}_{3}\dot{\overline{\boldsymbol{\Omega }}} + \widetilde{\overline{\boldsymbol{\Omega }}}\,\overline{\mathbf{J}}_{3} \,\overline{ \boldsymbol{\Omega }} &= \overline{\boldsymbol{\tau }}. \end{aligned}$$
(107)

5 Time integration of the unified local velocity coordinates

In order to obtain the current position and orientation of the rigid body from the unified velocity coordinates \(\mathbf{\overline{w}}\), we need to find a proper representation of the Lie algebra \(\widetilde{\overline{\mathbf{v}}}_{\textrm{L}}\) and furthermore a proper representation of the inverse tangent operator \(\mathbf{T}_{ \mathit{SE}(3)}^{-1}\) with respect to \(\overline{w}\). Hereafter, a proper time integration method can be applied to the incremental motion vector differential equation (13). The solution of the incremental motion is applied to Eq. (16) to obtain the current position and orientation of the rigid body.

5.1 Lie algebra corresponding to \(\mathbf{\overline{w}}\)

Since the Lie algebra \(\widetilde{\overline{\mathbf{v}}}_{\textrm{L}}\) in the kinematic differential Eq. (6) is a function of the translational velocity vector \(\overline{\mathbf{U}}\) and the angular velocity vector \(\overline{\boldsymbol{\Omega }}\), we need to find a corresponding representation of the Lie algebra with respect to the velocity coordinates \(\overline{w}_{i}\) in order to be able to apply the solution approach depicted in Eq. (16). This corresponding representation of the Lie algebra can easily be found, since the translational velocity of the reference point can be determined with respect to \(\mathbf{\overline{w}}\) using Eq. (32)

$$ \overline{\mathbf{U}}= \overline{\mathbf{v}}(\mathbf{0}) = \overline{ \mathbf{N}}(\mathbf{0})\mathbf{\overline{w}} $$
(108)

and Eq. (71),

$$ \overline{\mathbf{v}}_{\textrm{L}}= \renewcommand{\arraystretch}{1.2} \begin{bmatrix} \overline{\mathbf{U}} \\ \overline{\boldsymbol{\Omega }} \end{bmatrix} = \overline{\mathbf{D}}^{-1} \mathbf{\overline{w}}= \begin{bmatrix} \overline{\mathbf{D}}_{11} & \overline{\mathbf{D}}_{12} \\ \overline{\mathbf{D}}_{21} & \overline{\mathbf{D}}_{22} \end{bmatrix} \mathbf{ \overline{w}}\quad \Longrightarrow \quad \overline{ \boldsymbol{\Omega }}= \begin{bmatrix} \overline{\mathbf{D}}_{21} & \overline{\mathbf{D}}_{22} \end{bmatrix} \mathbf{\overline{w}}. $$
(109)

The vector \(\overline{\boldsymbol{\Omega }}\) therefore follows in the hexahedral case:

Ω=[Ω1Ω2Ω3]=[00121200000012121121120D21112112010100D2210100]w=[w3+w42w5+w62w1+w22].
(110)

After applying Eq. (9) to Eq. (110), it follows from comparison of \(\overline{\dot{\mathbf{E}}}\) in Eq. (49) and \(\widetilde{\overline{\boldsymbol{\Omega }}}\) that

$$ \widetilde{\overline{\boldsymbol{\Omega }}}= \overline{\dot{ \mathbf{E}}}( \mathbf{\overline{w}}). $$
(111)

Since the local angular velocity tensor \(\widetilde{\overline{\boldsymbol{\Omega }}}\) can be expressed with matrix \(\overline{\dot{\mathbf{E}}}\), the Lie algebra \(\widetilde{\overline{\mathbf{v}}}_{\textrm{L}}(\mathbf{\overline{w}})\) corresponding to \(\mathbf{H}\in \mathit{SE}(3)\) is given by the matrix

$$\begin{aligned} \widetilde{\overline{\mathbf{v}}}_{\textrm{L}}(\mathbf{ \overline{w}}) : &= \begin{bmatrix} \widetilde{\overline{\boldsymbol{\Omega }}}& \overline{\mathbf{U}}\\ \mathbf{0}_{1\times 3} & 0 \end{bmatrix} = \begin{bmatrix}\, \overline{\dot{\mathbf{E}}}(\mathbf{\overline{w}}) & \overline{ \mathbf{N}}(\mathbf{0}_{1\times 3})\mathbf{\overline{w}}\\ \mathbf{0} _{1\times 3} & 0 \end{bmatrix} \end{aligned}$$
(112)
$$\begin{aligned} &= \begin{bmatrix} 0 & \frac{-\overline{w}_{1} - \overline{w}_{2}}{\sqrt{2}} & \frac{ \overline{w}_{5} + \overline{w}_{6}}{\sqrt{2}} & \frac{\overline{w} _{1} - \overline{w}_{2}}{\sqrt{2}} \\ \frac{\overline{w}_{1} + \overline{w}_{2}}{\sqrt{2}} & 0 & \frac{-\overline{w}_{3} - \overline{w}_{4}}{\sqrt{2}} & \frac{\overline{w}_{3} - \overline{w} _{4}}{\sqrt{2}} \\ \frac{-\overline{w}_{5} - \overline{w}_{6}}{ \sqrt{2}} & \frac{\overline{w}_{3} + \overline{w}_{4}}{\sqrt{2}} & 0 & \frac{\overline{w}_{5} - \overline{w}_{6}}{\sqrt{2}} \\ 0 & 0 & 0 & 0 \\ \end{bmatrix} . \end{aligned}$$
(113)

5.2 Inverse tangent operator

The differential equation for the incremental motion vector \(\mathbf{q}\) is computed as already shown in Eq. (13). The relationship between \(\overline{\mathbf{v}}_{\textrm{L}}\) and the vector of the unified coordinates \(\mathbf{\overline{w}}\) is established by the vectors \(\bar{\mathbf{b}}_{i}\) and the position vectors \(\bar{ \mathbf{y}}_{i}\) in the form of the inverse transformation matrix \(\overline{\mathbf{D}}^{-1}\) via the relationship \(\overline{ \mathbf{v}}_{\textrm{L}}= \overline{\mathbf{D}}^{-1} \mathbf{\overline{w}}\) as shown in Sect. 4.1. The inverse tangent operator \(\widehat{\mathbf{T}}_{\mathit{SE}(3)} ^{-1}\) matching the unified coordinates can thus be determined by matrix multiplication

$$ \widehat{\mathbf{T}}_{\mathit{SE}(3)}^{-1}( \mathbf{q}) = \mathbf{T} _{\mathit{SE}(3)}^{-1}(\mathbf{q})\,\overline{ \mathbf{D}}^{-1}. $$
(114)

As shown in Sect. 4.3 and Sect. 4.4, the matrix \(\overline{ \mathbf{D}}\) is regular in the cases we consider. Therefore, the tangent operator in Eq. (18) is decisive for the investigation of possible singularities in Eq. (114). As the problem of singularities in Eq. (18) also occurs in the Lie-group formulations reported in the literature, we use the method presented in [2] to evaluate \(\mathbf{T}_{ \mathit{SE}(3)}^{-1}\) robustly.

5.3 Integration scheme

We focus next on a numerical time integration scheme, which can be used to solve the differential equations

$$\begin{aligned} \dot{\mathbf{\overline{w}}} &= \overline{\mathbf{M}}^{-1} (\overline{ \mathbf{Q}} - \overline{\boldsymbol{\Gamma }}\mathbf{\overline{w}} ), \end{aligned}$$
(115)
$$\begin{aligned} \dot{\mathbf{q}} &= \widehat{\mathbf{T}}_{\mathit{SE}(3)}^{-1}( \mathbf{q}) \mathbf{\overline{w}}, \quad \mathbf{q}_{0}=\mathbf{0}, \end{aligned}$$
(116)

with \(\mathbf{q}_{0}\) being initial condition. Once a solution of Eq. (116) is obtained for a specific time step (\(j\)\(j+1\)), the corresponding \(\mathit{SE}(3)\) group element update is obtained from the update formula

$$ \mathbf{H}_{j+1} = \mathbf{H}_{j}\, \textrm{exp}_{\mathit{SE}(3)}( \Delta \mathbf{q}_{j}) $$
(117)

with \(\Delta \mathbf{q}_{j} = \mathbf{q}_{j+1} - \mathbf{q}_{j}\). For the purpose of numerical time integration, we adopt a recently presented fourth-order Runge–Kutta time integration scheme proposed by Terze et al. [35], such that it can be applied to our formulation. Following the notation used in [35], we define

$$ \dot{\mathbf{\overline{w}}} = f (\mathbf{H}, \, \mathbf{ \overline{w}}, \, t ) $$
(118)

in which \(f (\mathbf{H}, \,\mathbf{\overline{w}}, \, t )\) represents the function which provides the derivative of the unified velocity coordinates using the equations of motion (115). Starting with values at the previous step, \(\mathbf{\overline{w}}_{\textrm{i}}\), \(\Delta t_{j}\), \(\mathbf{H}_{j}\), \(\mathbf{q}_{j}\), the slope estimations are obtained by

K1=TˆSE(3)1(0)wj,k1=f(Hj,wj,tj),
(119)
K2=TˆSE(3)1(12K1)(wj+12k1),k2=f(Hj,wj+12k1,tj+Δt2),
(120)
K3=TˆSE(3)1(12K2)(wj+12k2),k3=f(Hj,wj+12k2,tj+Δt2),
(121)
K4=TˆSE(3)1(K3)(wj+k3),k4=f(Hj,wj+k3,tj+Δt).
(122)

The four coefficients \(\mathbf{K}_{1},\ldots ,\mathbf{K}_{4}\) represent the slopes within one time step (\(j\)\(j+1\)) in the solution process of the incremental motion vector differential equation (116) and the four coefficients \(\mathbf{k} _{1},\ldots ,\mathbf{k}_{4}\) contain the slopes within one time step in the solution process of Eq. (115). The quantity \(\Delta t\) represents the time step size in the integration process. The equations

$$\begin{aligned} \mathbf{H}_{j+1} &= \mathbf{H}_{j}\, \textrm{exp}_{\mathit{SE}(3)} (\Delta \mathbf{q}_{j} ), \end{aligned}$$
(123)
$$\begin{aligned} \Delta \mathbf{q}_{j} &= \frac{\Delta t}{6} ( \mathbf{K}_{1}+2 \mathbf{K}_{2}+2\mathbf{K}_{3}+ \mathbf{K}_{4} ), \end{aligned}$$
(124)
$$\begin{aligned} \mathbf{\overline{w}}_{j+1} &= \mathbf{\overline{w}}_{j} + \frac{ \Delta t}{6} (\mathbf{k}_{1} + 2\mathbf{k}_{2}+2 \mathbf{k}_{3}+ \mathbf{k}_{4} ), \end{aligned}$$
(125)

update the quantities \(\mathbf{\overline{w}}\) and \(\mathbf{H}\).

6 Numerical example

The rotation of a torque-free rigid body with an initial angular velocity vector oriented closely to its unstable axis is simulated according to [21, 35]. The inertia tensor of the rigid body with respect to the center of mass reads

$$ \overline{\mathbf{J}}_{3}=\textrm{diag}(5.2988, \, 1.1775, \, 4.3568). $$
(126)

All quantities are expressed in standard units. An unstable rotation about the third axis is expected, since \(\overline{J}_{11} > \overline{J}_{33} > \overline{J}_{22}\). Throughout the numerical example, the unified local velocity coordinates \(\overline{w}_{i}\) according to the hexahedral case with \(l_{1}=l_{2}=l_{3}=1\) are used. The initial local velocities are chosen according to [21, 35],

$$ {\overline{\mathbf{v}}_{\textrm{L}}}_{0} = \begin{bmatrix} \overline{\mathbf{U}}_{0}^{T} & \overline{\boldsymbol{\Omega }}_{0}^{T} \end{bmatrix} ^{T} = \begin{bmatrix} 0 & 0 & 0 & 0.01 & 0 & 100 \end{bmatrix} ^{T} $$
(127)

where \(\overline{\mathbf{U}}_{0}=\mathbf{0}_{1\times 3}\) represents the velocity vector at the center of mass and \(\overline{\boldsymbol{\Omega }}_{0}= [0.01\ 0\ 100] ^{T}\) represents the body’s local angular velocity vector.

Therefore, the initial values for the unified local velocity coordinates \(\overline{w}_{i}\) are set according to Eq. (70)

$$ \mathbf{\overline{w}}_{0} = \overline{\mathbf{D}}\,{\overline{ \mathbf{v}}_{\textrm{L}}}_{0} = \begin{bmatrix} \frac{100}{\sqrt{2}} & \frac{100}{\sqrt{2}} & \frac{0.01}{ \sqrt{2}} & \frac{0.01}{\sqrt{2}} & 0 & 0 \end{bmatrix} ^{T}. $$
(128)

The initial position of the center of mass is set to \(\mathbf{x}_{0} = \mathbf{0}_{3\times 1}\) and the body’s initial orientation is given by the rotation matrix \(\mathbf{R}_{0} = \mathbf{I}_{3\times 3}\). Therefore, the initial conditions in terms of an element of \(\mathit{SE}(3)\) follow:

$$ \mathbf{H}_{0}= \begin{bmatrix} \mathbf{R}_{0} & \mathbf{x}_{0} \\ \mathbf{0} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} . $$
(129)

Since no external forces or torques are applied, the equations of motion follow from Eq. (62),

$$ \overline{\mathbf{M}}\dot{\mathbf{\overline{w}}} + \overline{ \boldsymbol{\Gamma }}\mathbf{\overline{w}}= \mathbf{0}_{6\times 1} . $$
(130)

Despite the fact that no translational motion is studied in this example, we set the body’s mass \(m=1\) in Eq. (93) for the sake of completeness.

For a better understanding of the simulation results, the position of a point \(P\) given in the body-fixed frame by the vector \(\bar{ \mathbf{p}} = [\, 0 \,\,\, 0 \,\,\, 1\, ]^{T}\), is tracked during rotation. Numerical simulations are conducted in the time domain of one second. The reference solution, denoted as’quaternions’, is obtained following the approach of Terze et al. [35] based on Euler parameters and an appropriate fourth-order Runge–Kutta scheme which inherently respects the unit-length condition. Since no reference solution data from Terze et al.’s publication were available with which we can directly compare our results, the formulation of Terze et al. was implemented by ourselves. All formulations have been implemented in Matlab R2016a. Our unified local velocity coordinate formulation will be abbreviated by means of’ulvc (proposed)’ in the following figures. Since all methods converged to the same solution (up to machine precision) for a time step of \(\Delta t=10^{-5}\), results obtained with that time step can be considered as a reference.

The trajectory of \(\bar{\mathbf{p}}\), expressed in an inertial frame, is shown in the three-dimensional plot in Fig. 8, while its three components versus time are shown in Fig. 9. The position of the point \(\bar{\mathbf{p}}\) is changing between \([\, 0 \,\,\, 0 \, \,\, 1\, ]^{T}\) and \([\, 0 \,\,\, 0 \,\,\, -1\, ] ^{T}\). The components of the body’s angular velocity vector \(\overline{ \boldsymbol{\Omega }}\) are shown in Fig. 10. In order to obtain the results shown in Figs. 810, we used a time step size of \(\Delta t=10^{-4}\). As expected, the body’s rotation is unstable, the body flips over and continues its rotation as one can observe in Fig. 10. Figure 11 shows the deviation of the determinant of the rotation matrix from one, expressed by the formula \(\textrm{det}(\mathbf{R}) - 1\). In order to compare the formulation of Terze et al. and the formulation presented here with respect to \(\textrm{det}(\mathbf{R}) - 1\), we convert the unit quaternion computed with the formulation of Terze et al. into a rotation matrix, using Eq. (2.13) in [32] from which we calculate the determinant and use it for comparison. As depicted in Fig. 11, the proposed formulation shows a very good behavior in long term simulation as well. Finally, Fig. 12 illustrates the convergence order of the time integration scheme based on the norm of the error of the angular velocity vector \(\Vert \overline{ \boldsymbol{\Omega }}-\overline{\boldsymbol{\Omega }}_{\textrm{converged}} \Vert _{2}\) versus the integration step size \(\Delta t\), in which the reference solution is obtained with a time step \(\Delta t=10^{-5}\). Figure 13 shows the error of the determinant of the rotation matrix, \(\textrm{det}( \mathbf{R})-1\), versus the integration step size \(h\). As shown in Fig. 12 the unified local velocity coordinate formulation exhibits the same convergence and accuracy as the method proposed by Terze et al. [35], since both integration schemes are based on a fourth-order explicit Runge–Kutta method. As depicted in Fig. 13, the absolute error of the determinant of the rotation matrix is lower than \(2 \times 10^{-14}\) for both methods. The unified local velocity coordinate formulation shows a slightly better accuracy than the solution obtained with the method proposed by [35], whereas the method presented in [35] shows a slightly higher convergence rate. The better accuracy of the unified local velocity coordinate formulation compared with the one obtained with the method of Terze et al. may be explained by the fact that we use directly the exponential map on \(\mathit{SO}(3)\) to obtain the rotation matrix, whereas in the method of Terze et al. the rotation matrix is computed from unit quaternions.

Fig. 8
figure 8

Spatial trajectory of point \(\bar{\mathbf{p}}\)

Fig. 9
figure 9

Components of the point \(\bar{\mathbf{p}}\)

Fig. 10
figure 10

Components of the angular velocity vector in the body-fixed frame

Fig. 11
figure 11

Determinant of rotation matrix

Fig. 12
figure 12

Convergence order of the proposed time integration using the norm of the error of the angular velocity vector, shown on a double-logarithmic scale. As can be seen from the figure, both methods show fourth-order convergence

Fig. 13
figure 13

Absolute error of the determinant of the rotation matrix, \(\textrm{det}(\mathbf{R})-1\), versus time integration step size

7 Conclusions

Concluding, we found a formulation to describe spatial rigid body motion in terms of six non-redundant, homogeneous local velocity coordinates. This formulation is free of rotational parameters and allows the description of rigid body motion in space using purely translational velocity coordinates. The unified local velocity coordinates can be built using various projection geometries such as vectors within the faces of a hexahedron or vectors along the edges of a tetrahedron. Especially the velocity projection along the faces of a hexahedron proves to be a good choice, since in this case the transformation matrix may becomes orthogonal. We obtain equations of motion of similar structure as compared to the classic Newton–Euler equations, while not making use of the angular velocity vector or rotation matrices. The equations of motion given purely in terms of translational velocities can be integrated using Lie-group time integration schemes, since the Lie algebra can be represented w.r.t. the unified local velocity coordinates. Both the equations of motion and the incremental motion vector differential equation can be integrated directly using for example the Runge–Kutta formulas of the fourth order. The exponential map on \(\mathit{SE}(3)\) can be used to determine position and orientation of the rigid body in space for each time step. The presented formulation is numerically as accurate as existing formulations. However, the presented formulation promises to be useful in cases where all coordinates need to have an identical (unified) interpretation or when rotation matrix and angular velocity vector like to be avoided.