The equations of motion for a rigid body using non-redundant unified local velocity coordinates

A novel formulation for the description of spatial rigid body motion using six non-redundant, homogeneous local velocity coordinates is presented. In contrast to common practice, the formulation proposed here does not distinguish between a translational and rotational motion in the sense that only translational velocity coordinates are used to describe the spatial motion of a rigid body. We obtain these new velocity coordinates by using the body-fixed translational velocity vectors of six properly selected points on the rigid body. These vectors are projected into six local directions and thus give six scalar velocities. Importantly, the equations of motion are derived without the aid of the rotation matrix or the angular velocity vector. The position coordinates and orientation of the body are obtained using the exponential map on the special Euclidean group SE(3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathit{SE}(3)$\end{document}. Furthermore, we introduce the appropriate inverse tangent operator on SE(3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathit{SE}(3)$\end{document} in order to be able to solve the incremental motion vector differential equation. In addition, we present a modified version of a recently introduced a fourth-order Runge–Kutta Lie-group time integration scheme such that it can be used directly in our formulation. To demonstrate the applicability of our approach, we simulate the unstable rotation of a rigid body.


Introduction
In this paper, we present a formulation for spatial rigid body motion based on a set of nonredundant, 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 SO (3) with the linear space R 3 either with the Cartesian product leading to  (3)) [3]. Note, both the translational velocity and the angular velocity are represented within the 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 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'tangentialdisplacement 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.
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.

Lie-group framework
Since the six unified velocity coordinates are defined in a body-fixed frame, standard integrators such as the generalized-α 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 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ā 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, 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 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 x P ∈ R 3 of any point P in the current configuration on the rigid body is expressed by the affine mapping where x ∈ 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 y P = Rȳ P represents the position of an arbitrary point P of the body with respect to the fixed reference frame I. Whereas the vectorȳ P ∈ R 3 represents the position of the point P in the body-fixed frame. The rotation matrix R maps the vectorȳ 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 The set of 4 × 4-matrices, which are formed by the semi-direct product SO(3) R 3 , where R ∈ SO(3) and x ∈ R 3 forms together with the matrix product (composition operation) a matrix Lie group, which is called the special Euclidean group SE(3) [30]. The velocity vectorẋ P ∈ R 3 of any point P on the rigid body can be expressed as linear relation betweenḢ and [ȳ T P 1] T in the form The time derivativeḢ readṡ denotes any element of se(3), the Lie algebra of SE (3), which is defined by with = R TṘ . The Lie algebra is isomorphic to a finite dimensional linear space R k with an invertible linear mapping where the tilde operator (•) depend on the specific Lie-group setting [3,7]. In the case of the Lie algebra se(3), the inverse mapping (9) maps a 4 × 4-matrix v L ∈ se(3) to a sixdimensional vector and vice versa. In the case of Lie algebra so(3), the inverse mapping (9) maps a 3 × 3-matrix ∈ so(3) to a three-dimensional vector and vice versa. In the vector representation v L of v L , the sub-vector represents the translational velocity in a body-fixed frame and the sub-vector 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 H andḢ [18,30]. This differential equation is not solved directly. Instead, the incremental motion vector differential equationq with q 0 being initial condition is solved with a time integration scheme for where While the subindex T in q T stands for translation, the subindex R in 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 exp SE(3) on SE(3) denotes the inverse tangent operator corresponding to the exponential map on SE(3) shown in Eq. (17). The exponential map and the inverse tangent operator on SE(3) can be written in a closed form where exp SO(3) is the exponential map and T SO(3) the tangent operator on the special orthogonal group SO(3) [3,6,33], which are given by using the Lie bracket , ; see e.g. [22]. The matrices q R and q T are skew symmetric matrices of the form (8), which are computed from the vectors q R and q T .

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.

Non-redundant unified local velocity coordinates
The unification of the local translational velocity coordinates (U 1 , U 2 , U 3 ) and the local angular velocity coordinates (Ω 1 , Ω 2 , Ω 3 ) to six local velocity coordinates (w 1 , w 2 , w 3 , w 4 , w 5 , 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 directionb i and the velocity vector v i (ȳ i ), see Fig. 3, which is related to the velocity vectorẋ P by the relation generates the ith unified local velocity coordinate w i ∈ R, The ith local velocity vector of the ith point P i is indicated by the vector v i ∈ R 3 . The local position of point P i on the rigid body is marked with vectorȳ i ∈ R 3 . The vectorb i ∈ R 3 represents the ith vector used for the projection. Therefore, the vector of unified local velocity coordinates reads The term'unified' emphasizes the fact that all velocity coordinates 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 e i , y i , v i and y i represent body-fixed (physical) vectors. Whereas the vectors e i ,b i , v i andȳ i are column vectors whose components are represented in a local frame. The set of projection points P i and direction vectorsb i to generate the six unified velocity coordinates w i is limited by the condition that there must be a bijective mapping between the velocity coordinates w i and the vector form of v 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 v 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 b i and the position vectorsȳ i read for the hexahedron as depicted in Fig. 4a, with a quantity α ∈ R to be chosen. Note that other directions ofb i would be possible in the hexahedron as well. In Sect. 4.3, it will be shown that choosing α = 1/ √ 2 is advantageous, which is employed throughout the following description.

Velocity field
In order to relate the velocity field v(ȳ p ) to the unified velocities w i , we use Eq.
with the interpolation matrix N(ȳ p ). In order to determine the interpolation matrix N(ȳ p ), the componentv 1 is investigated, exemplary. The velocity componentv 1 is composed of three parts, which result from different forms of motion of the rigid body. Namely a translational velocityv transl 1 , which occurs due to a pure translational motion along e 1 , a velocityv according to Fig. 5b. The negative slope inv Summing up,v 1 can be represented in vector notation as in which N 1×6 1 (ȳ 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 N(ȳ p ) follow from similar considerations and read as follows:

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 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 v(ȳ p ). In order to do so, the velocity field v in Eq. (32) is transformed into the inertial frame, v(ȳ p ) =v 1 e 1 +v 2 e 2 +v 3 e 3 = N 1 (ȳ p )w e 1 + N 2 (ȳ p )w e 2 + N 3 (ȳ p )w e 3 , where (e 1 , e 2 , 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: where (ė 1 ,ė 2 ,ė 3 ) represent the time derivatives of the three body-fixed orthogonal basis vectors. Formulated in matrix-vector notation, the acceleration vector (41) takes the forṁ The components of the acceleration vector represented in the local frame read as follows: where (e 1 , e 2 , e 3 ) represent the three body-fixed basis vectors, with their (time dependent) coordinates represented in the local frame and (ė 1 ,ė 2 ,ė 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 motioṅ e i vanishes. But in the case of a rotational motion it follows thatė i = 0. Sinceė i can be interpreted as difference of two velocities, Eq. (32) is used to obtain, e.g., The velocity term v(0) in Eq. (44) corresponds to a pure translational motion atȳ p = 0 3×1 . Equation (44) can be evaluated using Eq. (32). By inserting the Cartesian basis vectors into Eq. (44), it follows that the time derivatives of the basis vectors reaḋ To obtain a compact vector-matrix notation of Eq. (43), the time derivatives of the basis vectors are collected in the matrix 1Ė Inserting Eqs. (47)-(48) and Eq. (45) into Eq. (43), the acceleration vectorv(ȳ P ) of any point P follows asv where the matrix P(ȳ p , w) is given by and the time derivative of the unified coordinates readṡ

Equations of motion of a free rigid body
In this section, we determine the governing equations of motion in terms of w andẇ. 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ẇ and equate the derivatives with the corresponding generalized forcesQ 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, Since the vectorẇ does not appear in the term G 3 , the partial derivative of G 3 with respect toẇ vanishes. Therefore, we only make use of the first two terms G 1 and G 2 , since their partial derivative with respect toẇ is obviously nonzero. Therefore, G 1 and G 2 follow to A further evaluation of the integrals in G 1 and G 2 lead us to the definition of a mass matrix M and a velocity dependent matrix Since the interpolation matrix N(ȳ p ) depends only on the chosen reference pointȳ p , which is constant in the body-fixed frame, it follows that the mass matrix M is also constant with respect to the body-fixed frame. The partial derivatives of G 1 and G 2 follow to Adding up the partial derivatives of G 1 and G 2 , the equations of motion follow to Note that Eqs. (62) have been obtained without the use of rotation matrices or the angular velocity vector.

Generalized forces
Corresponding to the equations of motion (62), we define generalized forces Q ∈ 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, For the case of a force f ∈ R 3 acting at a point P at the rigid body, the vector of generalized forces Q thus reads

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.

Relation between w and v L
In this section, a relationship between the unified local velocity coordinates w and the vector form of v 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 v(ȳ P ) = U + ×ȳ P , which can be rewritten in terms of v L as follows: v(ȳ P ) = U + ×ȳ P = U + ȳ By inserting Eq. (66) into Eq. (24) it follows that the ith unified velocity coordinate w i is related to v L by From the selection of six projection points P i and direction vectorsb i it follows that the vector of the unified velocity coordinates w is related to v L via the relation According to Eq. (68) the matrix D ∈ R 6×6 maps v L to the vector of the unified velocity coordinates. Therefore, the inverse of D maps the vector of the unified velocity coordinates to v L according to

Condition for the choice ofb i andȳ i
In order to obtain a valid set of unified local velocity coordinates w i , the following conditions need to be fulfilled for the transformation quantitiesb i andȳ i : 1. there is a linear mapping between the vector of unified local velocity coordinates w and the vector v L according to Eq. (68), and 2. the matrix D of the linear mapping is a regular matrix.

Relation between w and v 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 w and v L reads for the hexahedral case where D turns out to be an orthogonal matrix, i.e. D Finally, v L is given in terms of w, reading (73)

Relation between w and v L for the tetrahedral case
In this subsection we present a possible selection of direction vectorsb i and position vectors y i , the resulting matrix D tet , as well as their inverse D −1 tet for the velocity projection along the edges of a tetrahedron. For simplicity we choose a tetrahedron as shown in Fig. 7. The vectors (e 1 , e 2 , e 3 ) represent the orthogonal basis as shown in Eq. (45). We select the coordinates for each node (i, j, k) as le 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 inb with α being a coefficient, which scales the direction vectorsb i . The matrix D tet and its inverse D −1 tet therefore take the form (80) In contrast to the projection along the side surfaces of a hexahedron, the matrix D tet is in general not orthogonal for the tetrahedral case.
Looking at the determinant of the matrix D 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 α and the length l, due to the fact that det(D tet ) = α 6 l 3 . (81)

Mass matrix
In order to obtain a more suitable form of the mass matrix M, the integral in Eq.
Using Eq. (83), the mass matrix Eq. (59) is therefore obtained from with M being defined by where the matrix J 1 is computed by the integral, the matrix J 2 by the integral, and the matrix J 3 by the integral From the definition of the skew symmetric matrixỹ p one concludes that if the origin of the body axes is attached to the center of mass of the body, then the matrix J 2 is the null matrix [32], which is employed throughout the further description. The mass matrix thus follows: where the diagonal elements of the 3 × 3-matrix J 1 contain the rigid body's mass m and the 3 × 3-matrix J 3 represents the body's inertia tensor represented in the body-fixed frame. If we use a diagonal inertia tensor the mass matrix M exhibits in the hexahedral case using Eqs. (26)-(31) and setting l 1 = l 2 = l 3 = 1 a simple block diagonal structure,

Quadratic velocity vector
The quadratic velocity vector 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 results in Since Eq. (94) is integrated over the volume and the matrixĖ(w) does not depend on the volume, we obtain By comparing the integrals in Eq. (96) with the integrals in Eq. (84) it follows that the matrix reads with the hereby defined matrix Ψ , As described in Sect. 4.5, J 2 = 0 3×3 , therefore the matrix is computed as Using the diagonal inertia tensor of Eq. (92), the vector w reads for the hexahedral case using Eqs. (26)-(31) and setting l 1 = l 2 = l 3 = 1 as follows:

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 D where τ = ȳ p f ∈ R 3 represents the torque acting on the rigid body due to the local force f ∈ R 3 acting at point P . The local force vector f is connected to a physical force f ∈ R 3 by Since the mass matrix M and the velocity dependent matrix 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 D, Multiplication from the left hand side with D T gives Sinceẇ = Dv L and w = D v L the equations of motion can be written as 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],

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 w, we need to find a proper representation of the Lie algebra v L and furthermore a proper representation of the inverse tangent operator T  with respect to 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.

Lie algebra corresponding to w
Since the Lie algebra v L in the kinematic differential Eq. (6) is a function of the translational velocity vector U and the angular velocity vector , we need to find a corresponding representation of the Lie algebra with respect to the velocity coordinates 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 w using Eq. (32) and Eq. (71), The vector therefore follows in the hexahedral case: After applying Eq. (9) to Eq. (110), it follows from comparison ofĖ in Eq. (49) and that Since the local angular velocity tensor can be expressed with matrixĖ, the Lie algebra v L (w) corresponding to H ∈ SE(3) is given by the matrix

Inverse tangent operator
The differential equation for the incremental motion vector q is computed as already shown in Eq. (13).
As shown in Sect. 4.3 and Sect. 4.4, the matrix 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 T −1

Integration scheme
We focus next on a numerical time integration scheme, which can be used to solve the differential equationsẇ with q 0 being initial condition. Once a solution of Eq. (116) is obtained for a specific time step (j → j + 1), the corresponding SE(3) group element update is obtained from the update formula with q j = q j +1 − 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 definė in which f (H, 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, w i , t j , H j , q j , the slope estimations are obtained by The four coefficients K 1 , . . . , 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 k 1 , . . . , k 4 contain the slopes within one time step in the solution process of Eq. (115). The quantity t represents the time step size in the integration process. The equations update the quantities w and H.

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 All quantities are expressed in standard units. An unstable rotation about the third axis is expected, since J 11 > J 33 > J 22 . Throughout the numerical example, the unified local velocity coordinates 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], where U 0 = 0 1×3 represents the velocity vector at the center of mass and 0 = [0.01 0 100] T represents the body's local angular velocity vector. Therefore, the initial values for the unified local velocity coordinates w i are set according to Eq. (70) The initial position of the center of mass is set to x 0 = 0 3×1 and the body's initial orientation is given by the rotation matrix R 0 = I 3×3 . Therefore, the initial conditions in terms of an element of SE(3) follow: Since no external forces or torques are applied, the equations of motion follow from Eq. (62), 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 t = 10 −5 , results obtained with that time step can be considered as a reference.
The trajectory ofp, 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 pointp is changing between [ 0 0 1] T and [ 0 0 − 1 ] T . The components of the body's angular velocity vector are shown in Fig. 10. In order to obtain the results shown  in Figs. 8-10, we used a time step size of 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 det(R) − 1. In order to compare the formulation of Terze et al. and the formulation presented here with respect to det(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 − converged 2 versus the integration step size t , in which the reference solution is obtained with a time step t = 10 −5 . Figure 13 shows the error of the determinant of the rotation matrix, det(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 × 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 SO(3) to obtain the rotation matrix, whereas in the method of Terze et al. the rotation matrix is computed from unit quaternions.

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 rota-tional 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 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.