Frames, rotations and coordinate systems
In this section, we introduce frames, which provide reference points and rotations attached to the ground or to bodies, thus providing a bases for coordinate systems.
Reference points and reference frames
In contrast to the reference coordinates (which may include reference position and rotations), the reference point of a body, marker or joint is available in any configuration (current, initial, etc.). The reference point and orientation (or reference frame) of a rigid body coincides with the position and orientation of the underlying node. The local position of a point of the body is expressed relative to the body’s reference frame, which may be different from the body’s center of mass. The reference frame is also used for FFRF objects. In most cases, reference points are denoted by \(\pRef\).
Coordinate systems of frames
Many considerations and computations in kinematics can be formulated coordinate-free. However, ultimately local positions are given by the user in different coordinates from results given in global (world) coordinates.
For this reason, we use left upper indices in symbols to indicate underlying frames of points or vectors, e.g., \(\LU{0}{{\mathbf{u}}}\) represents a displacement vector in the global (world) coordinate system \(0\), while \(\LU{m1}{{\mathbf{u}}}\) represents the displacement vector in marker 1 (\(m1\)) coordinates.
Typical coordinate systems (frames) are:
\(\LU{0}{{\mathbf{u}}}\) \(\ldots\) global or world coordinates
\(\LU{b}{{\mathbf{u}}}\) \(\ldots\) body-fixed, local coordinates
\(\LU{m0}{{\mathbf{u}}}\) \(\ldots\) local coordinates of (the body or node of) marker \(m0\)
\(\LU{m1}{{\mathbf{u}}}\) \(\ldots\) local coordinates of (the body or node of) marker \(m1\)
\(\LU{J0}{{\mathbf{u}}}\) \(\ldots\) local coordinates of joint \(0\), related to marker \(m0\)
\(\LU{J1}{{\mathbf{u}}}\) \(\ldots\) local coordinates of joint \(1\), related to marker \(m1\)
To transform the local coordinates \(\LU{m0}{{\mathbf{u}}}\) of marker 0 into global coordinates \(\LU{0}{{\mathbf{x}}}\), we use the relation
in which \(\LU{0,m0}{\Rot}\) is the transformation matrix of (the body or node of) the underlying marker 0.
In particular, note that any transformation, e.g., of the angular velocity in body (b) and world coordinates,
only transforms the vector from one frame to the other, but does not change the vector \(\tomega\) itself.
Homogeneous transformations
In order to considere rotations and translations in frames, homogeneous transformations are used. An arbitrary vector \(\LU{i}{{\mathbf{r}}}\) in coordinates of frame \(\mathcal{F}_i\) can be expressed relative to \(\mathcal{F}_j\) with the vector \(\LU{j}{{\mathbf{p}}_i}\), pointing from \(O_j\) to \(O_i\) and a rotation matrix according to Eq. (13):
Eq. (14) can be written as a linear mapping,
in which the matrix
is referred to as the homogeneous transformation matrix, and \(\LU{j}{\hat {\mathbf{r}}} = [\LU{j}{{\mathbf{r}}} \;\; 1]^T\) is the homogeneous vector corresponding to \(\LU{j}{{\mathbf{r}}}\).
Analogous to the rotation matrix, \(\LU{ji}{{\mathbf{T}}}\) has an inverse, which follows as
Sequential application of homogeneous transformations simplifies to
We could also actively move vectors, by assuming that the homogeneous transformation matrix \(\LU{i0,i1}{{\mathbf{T}}}\) transforms coordinates within the same frame \(\mathcal{F}_i\), between two (time) steps 0 and 1,
which advances the vector \(\LU{i0}{\hat {\mathbf{r}}}\) in time.
Note that an efficient implementation would only include \({\mathbf{R}}\) and \({\mathbf{p}}\), without necessarily computing \(4 \times 4\) matrices.
All homogeneous transformations form the Special Euclidean Group SE(3) – the motion group – for describing rigid body movements in 3D. Therefore, homogeneous transformations also meet the criteria for groups (closure, associativity, existence of neutral and inverse elements).
Elementary rotations about the \({\mathbf{e}}_z\)-axis and translations along the \({\mathbf{e}}_x\)-axis, for example, are given by:
Arbitrary composite transformations can be constructed from elementary transformations. In general, the following applies:
Pre-multiplication (Pre-Multiplication): Transformation in the global coordinate system
Post-multiplication (Post-Multiplication): Transformation in the rotated coordinate system
By exploiting the tools of so-called Lie groups, rigid body movements can be interpolated using matrix logarithm (LogSE3(...)
) and matrix exponential function (ExpSE3(...)
).
Rotations
In this section, we discuss in particular rotations, as already introduced for coordinate systems and homogeneous transformations. Rotations can be represented by transformation matrices, leading to a linear transformation
which transforms the vector \(\LU{1}{{\mathbf{r}}}\) accordingly. However, rotations are inherently nonlinear. First, rotation parametrizations – except for the components of the rotation matrix itself – are highly nonlinear functions of rotation parameters. Second, subsequent rotations couple in a multiplicative (i.e., nonlinear) way. Therefore rotations have to be treated differently from translations.
Derivation of the transformation matrix from coordinate transformations
The vector \({\mathbf{a}}\) can be represented in coordinates of frame \(\mathcal{F}_1\) with basis vectors \(({\mathbf{e}}_{x1},\,{\mathbf{e}}_{y1},\,{\mathbf{e}}_{z1})\) and of frame \(\mathcal{F}_2\) with basis vectors \(({\mathbf{e}}_{x2},\,{\mathbf{e}}_{y2},\,{\mathbf{e}}_{z2})\),
Multiplying Eq. (15) with the basis vectors \(({\mathbf{e}}_{x1},\,{\mathbf{e}}_{y1},\,{\mathbf{e}}_{z1})\) consecutively, we obtain three relations,
where we apply the orthonormality relationships \({\mathbf{e}}^T_{x1} {\mathbf{e}}_{x1} = 1\), \(\ldots\) as well as \({\mathbf{e}}^T_{x1} {\mathbf{e}}_{y1} = 0\), \({\mathbf{e}}^T_{x1} {\mathbf{e}}_{z1} = 0\), \(\ldots\) througout.
We can represent the result as linear transformation
or in short
The Transformationsmatrix \(\LU{12}{\Rot}\) transforms coordinates of vector \({\mathbf{a}}\) from \(\mathcal{F}_2\) into \(\mathcal{F}_1\). We observe that the columns of \(\LU{12}{\Rot}\) contain unit vectors \({\mathbf{e}}_{x2}\), \({\mathbf{e}}_{y2}\), and \({\mathbf{e}}_{z2}\) represented in frame \(\mathcal{F}_1\), and that the rows can be identified as the unit vectors \({\mathbf{e}}_{x1}^T\), \({\mathbf{e}}_{y1}^T\), and \({\mathbf{e}}_{z1}^T\) represented in frame \(\mathcal{F}_2\):
Note that:
Columns and rows of \(\LU{12}{\Rot}\) represent an orthonormal right-hand system.
Because the determinant \(\det(\LU{12}{\Rot})= + 1\), we call \(\LU{12}{\Rot}\) to be proper-orthogonal.
Rotations conserve length, \(|{\mathbf{A}} {\mathbf{x}}| = |{\mathbf{x}}|\).
For the same reason, rotations conserve angles.
During transformations, we observe that reading left upper indices from left to right, indices are changed only by transformation matrices (which are not transposed): \(\LU{2}{{\mathbf{a}}} = \LU{21}{\Rot} \LU{1}{{\mathbf{a}}}\)
We note the following rules:
Elementary rotations
Elementary rotations are rotations about a single axis, using one of the orthogonal basis vectors. Consider basis \(({\mathbf{e}}_{x1},\,{\mathbf{e}}_{y1},\,{\mathbf{e}}_{z1})\) and a rotation around axis \({\mathbf{e}}_{x1}\) with angle \(\varphi_1\) in positive rotation sense, obtaining a rotated frame \(({\mathbf{e}}_{x2},\,{\mathbf{e}}_{y2},\,{\mathbf{e}}_{z2})\), see Fig. 19. The rotation matrix for this case reads:
In this case, the coordinates of a vector \({\mathbf{r}}\) are transformed according to
In a second example, a rotation with angle \(\varphi_2\) around \({\mathbf{e}}_{y2}\) is performed to transform from basis \(({\mathbf{e}}_{x2},\,{\mathbf{e}}_{y2},\,{\mathbf{e}}_{z2})\) into \(({\mathbf{e}}_{x3},\,{\mathbf{e}}_{y3},\,{\mathbf{e}}_{z3})\), see Fig. 20. The rotation matrix for this case reads:
Finally, a rotation around \({\mathbf{e}}_{z3}\) with angle \(\varphi_3\) would give similarly,
Linearized rotations
In this case we can interpret the small rotations as a constant angular velocity over small time \(\Delta t\), \(\boldVarPhi = \Delta t \cdot \tomega\), which results in
Leading to the linearized rotation matrix
Alternatively, using linearized rotations \(\boldVarPhi = [\varphi_1,\;\varphi_2,\;\varphi_3]\tp\), with \(|\boldVarPhi| \ll 1\), we can approximate \(\sin \varphi_1\) as:
Similarly, we can approximate \(\cos \varphi_1 \approx 1\). This immediately gives a linearization for elementary rotations – where we observe that contributions for each axis can be added.
The axis-angle representation (see later) leads to the same result by linearization of the Rodrigues formula,
in which \(\varphi\) represents the infinitesimal angle and \({\mathbf{u}}\) is the rotation axis.
Note: whenever you would like to work with linearized rotations, or if you do not know the sequence of small rotations, still do not use the linearized formula as it gives immediately rotation matrices without strict orthogonality and determinant of 1. In order to avoid computational problems, use for example RotationVector2RotationMatrix(phi)
to compute a consistent rotation matrix based on the matrix exponential.
Active and passive rotations
It is important to distinguish two fundamentally different concepts of rotations. There is the active rotation, which can be represented by a coordinate-free relation,
in which a vector \({\mathbf{v}}\) is actively rotated into a new configuration \({\mathbf{v}}^\prime\) using a rotation tensor \(\Rot\). Eq. (16) can be represented in any coordinate system, but keeping it fixed.
Furthermore, we denote a passive rotation as a coordinate transformation, as used many times before,
in which the vector \({\mathbf{v}}\) is not changed at all, but only its coordinates are represented in two different coordinate systems.
While passive rotations (and homogeneous transformations) are used, e.g., to compute current positions of body-attached points through several relative coordinate systems, active rotations can represent the rotations from one time step to the next one.
Successive rotations
The successive application of rotation matrices is non-commutative. An exception is the planar case, i.e., all rotations occur around the same axis. Specifically, we see that for successive rotations, the following must be considered:
As an example, we consider in Figure Successive rotations are not commutative. the different order of rotations of a block.
In Fig. 21a, the block shown is first rotated 90° about the \(\mathbf{e}_2\) axis and then 90° about the \(\mathbf{e}_3\) axis. In Figure Successive rotations are not commutative.b, the same rotations are applied in reverse order, i.e., the block is first rotated 90° about the \(\mathbf{e}_3\) axis and then 90° about the \(\mathbf{e}_2\) axis. It is immediately apparent that the resulting orientation of the block is different in both cases.
Finally, it should be noted that there are two different types of successive rotations. In the first variant, also known as the single-frame method, the same reference frame is chosen for each rotation. This variant is best understood in the active rotation of vectors, where these rotations always take place, for example, in the global reference system.
In the second variant, also known as the multi-frame method, a new reference frame created by the preceding rotation is used for each rotation. This variant can be illustrated by applications in robotics (e.g., articulated arm robots), where (passive) rotations of the reference systems of the robot’s respective arms occur, with each rotation taking place in the reference system of the preceding arm.
Rotation tensor: axis-angle representation
In the following, we will elaborate the rotation tensor, inherently linked to the axis-angle representation. Note that the rotation axis times the angle is denoted as rotation vector throughout.
Considering Euler’s theorem on rotations, we assume a transformation
which purely follows from a rotation. This transformation thus can be represented by a rotation axis \({\mathbf{u}}\) and an angle \(\varphi(t)\),
both of them given in a coordinate-free manner. We may therefore conclude
Using Fig. 22, we may now consider relations of the two frames \(({\mathbf{e}}_{x0},\,{\mathbf{e}}_{y0},\,{\mathbf{e}}_{z0})\) and \(({\mathbf{e}}_{x1},\,{\mathbf{e}}_{y1},\,{\mathbf{e}}_{z1})\), solely defined by the angle-axis \(({\mathbf{u}}(t), \, \varphi(t))\) relation.
According to Fig. 23, we may split the vector \({\mathbf{r}}\), which has the length \(r = \vert {\mathbf{r}} \vert\), into
with the unit vectors given as
The coefficients can be calculated as follows,
Putting everything together gives
From the relation \({\mathbf{r}} = \Rot({\mathbf{u}} ,\varphi) {\mathbf{r}}_0\), we can reinterpret the rotation tensor \(\Rot\)
With the additional relations
we finally obtain
which is also known as the Rodrigues formula for rotations. While this formula is coordinate-free, it may be represented in any coordinate system, such as
Rotation vector
Using the rotation vector \({\mathbf{v}}_\mathrm{rot} = \varphi\cdot {\mathbf{u}}\) and \(\varphi = |{\mathbf{v}}_\mathrm{rot}|\), another representation of the rotation tensor follows as
Note, that we inherently assume that \(\varphi \in [0,\pi]\). In Exudyn, there are the following functions and items related to the rotation vector and the axis-angle representation:
NodeRigidBodyRotVecLG
: A 3D rigid body node based on rotation vector and Lie group methods; can be used for explicit integration methods, not leading to singularities when integrating with Lie-group methodsRotationVector2RotationMatrix
: computes the rotation matrix from a given rotation vector \({\mathbf{v}}_\mathrm{rot}\)RotationMatrix2RotationVector
: computes the rotation vector \({\mathbf{v}}_\mathrm{rot}\) from a given rotation matrix, based on a reconstruction of Euler parametersComputeRotationAxisFromRotationVector
: computes the rotation axis \({\mathbf{u}}\) from the rotation vector by using \(\varphi = |{\mathbf{v}}_\mathrm{rot}|\)
Note the following properties of the rotation tensor:
The axis-angle representations \(({\mathbf{u}}, \varphi)\) and \((-{\mathbf{u}}, -\varphi)\) represent the same tensor \(\Rot({\mathbf{u}},\varphi)=\Rot(-{\mathbf{u}},-\varphi)\)
The rotation tensor \(\Rot\) is orthogonal, i.e., \(\Rot\tp \, \Rot = {\mathbf{E}}\).
The inverse rotation tensor \(\Rot^{-1} = \Rot\tp\) is represented by the identical axis-angle tuples \(({\mathbf{u}}, -\varphi)\) and \((-{\mathbf{u}}, \varphi)\), leading to \({\mathbf{r}}_0=\Rot({\mathbf{u}}, -\varphi) {\mathbf{r}}\).
We furthermore not the identities: \(\Rot({\mathbf{u}}, -\varphi) \equiv \Rot(-{\mathbf{u}}, \varphi) \equiv \Rot\tp ({\mathbf{u}}, \varphi)\)
The rotation axis is invariant to the rotation: \(\Rot({\mathbf{u}}, \varphi) {\mathbf{u}} = {\mathbf{u}}\)
Euler angles: Tait-Bryan angles
In the following section, we discuss rotation matrices built by Tait-Bryan angles, which are consecutive \(xyz\)-rotations.
Note that there are many other representations of Euler angles, however, they are not implemented or used in Exudyn.
The output variable Rotation
gives Tait-Bryan angles, which is why it is important to define their properties.
The transformation given by Tait-Bryan angles \((\alpha,\beta,\gamma)\) follows from three successive rotations,
The Tait-Bryan rotation matrix is thus defined as
While not fully obvious from the latter equations, we note that there is a singularity at \(|\beta| = \frac{\pi}{2}\), which causes any computations to stall whenever getting close enough to this value.
It is also possible to reconstruct Tait-Bryan angles from a given rotation matrix. Assume, we have the rotation matrix \(\LU{03}{\Rot}=(A_{ij})\). There are two solutions for \(\beta\) which follow from
For beta fulfilling \(|\beta| \ne \frac{\pi}{2}\), we can uniquely compute \(\gamma\) and \(\alpha\),
Note that improvements are possible using the \(\mathrm{atan2}()\) function. We can finally resolve non-uniqueness by restricting \(\beta\) within the range \(-\frac{\pi}{2} < \beta < \frac{\pi}{2}\).
Tait-Bryan angles: angular velocity vector
Angular velocities are essential for the formulation of equations of motion of rigid bodies, for constraints and for evaluation. Therefore, basic relations are shown for Tait-Bryan angles, in particular the velocity transformation matrix \({\mathbf{G}}_{TB}\).
The angular velocity vector \(\omega\) can either be computed from the basic relation \(\dot \Rot = \tilde \omega \Rot\), which however involves many trigonometric terms, or from partial angular velocities, being part of the co-rotated rotation axes. With the latter approach, we see that
We can now evaluate this equation in frame \(\mathcal{F}_0\) with \(\tbeta=[\alpha\;\beta\;\gamma]\tp\), written as
Evaluating the consecutive rotations, we find the respective unit vectors as
which results in the relations for the (global) angular velocity vector
The matrix \(\LUR{0}{{\mathbf{G}}}{TB}\) is used to compute partial derivatives of the angular velocity vector w.r.t. rotations. However, it is also used to project torques and angular velocities into a space given by \({\mathbf{G}}_{TB}^\mathrm{T}\) (which is indeed not equivalent to \({\mathbf{G}}_{TB}^{-1}\), but gives a symmetric mass matrix – see the equations of motion for the rigid body).
The inverse relation between \(\dot \tbeta\) and \(\LUR{0}{\tomega}{30}\) gives
or in compact form
We again see that for the case \(\vert \beta \vert = \pi/2\) the matrix\(\LUR{0}{{\mathbf{G}}}{TB}\) becomes singular, and Eq. (18) cannot be resolved for \(\dot{\tbeta}\)!
Finally, we also provide the relations of \(\tomega_{30}\) in body-fixed (local) coordinates \(\mathcal{F}_3\), which are frequently used in the implementation,
which reads in compact form (\(\LUR{3}{{\mathbf{G}}}{TB} = \LUR{b}{{\mathbf{G}}}{TB}\)),
In Exudyn, there are the following functions and items related to Tait-Bryan angles:
NodeRigidBodyRxyz
: A 3D rigid body node based on Tait-Bryan anglesRotXYZ2RotationMatrix
: computes the rotation matrix from given Tait-Bryan anglesRotationMatrix2RotXYZ
: computes Tait-Bryan angles from a given rotation matrixRotXYZ2G
: computes the matrix \(\LU{0}{{\mathbf{G}}}\) relating time derivatives of Tait-Bryan angles to the (global) angular velocity vectorRotXYZ2Glocal
: computes the matrix \(\LU{b}{{\mathbf{G}}}\) relating time derivatives of Tait-Bryan angles to the (local, body-fixed) angular velocity vector
Euler parameters and unit quaternions
As one of the most important forms of rotational parameters for multibody systems, the 4 Euler parameters are introduced. They can be mathematically represented as unit quaternions and provide particularly simple expressions and equations of motion, albeit at the expense of an additional constraint.
The Euler parameters are defined as
Here, \(p_s=\cos \frac\varphi 2\) is the scalar part and \({\mathbf{p}} = {\mathbf{u}} \, \sin \frac\varphi 2\) is the vector part of the Euler parameters \(\underline{{\mathbf{p}}}\). The four Euler parameters \(\underline{{\mathbf{p}}}\) are subject to the normalization condition
First, the following trigonometric transformations are introduced,
When these are substituted into the rotation tensor using Eq. (19), it follows
Euler parameters can initially be expressed in the global reference frame \(\mathcal{F}_0\),
Thus, the coordinate representation of the rotation tensor is given by
Note that the axis of rotation is invariant to the rotation, \(\LU{0}{{\mathbf{u}}} = \LU{1}{{\mathbf{u}}}\), thus \(\LU{0}{\underline{{\mathbf{p}}}} = \LU{1}{\underline{{\mathbf{p}}}}\)Due to the ambiguity \(\Rot(\underline{{\mathbf{p}}}) = \Rot(-\underline{{\mathbf{p}}})\), the scalar part of the Euler parameters is usually normalized, i.e., \(p_s \ge 0\).
The careful reader may observe that matrix \(\Rot\) in Eq. (20) is not identical to the implementation in Exudyn. This is true, because there is always an alternative formulation for the diagonal terms by adding the normalization condition times a factor.
Quaternion operations
Calculations with Euler parameters are efficiently performed through the rules of quaternion algebra. Active perspective as vector rotation: The rotation of vector \({\mathbf{r}}_0\) into \({\mathbf{r}}(t)\),
is formulated using the quaternion \(\LU{0}{\underline{{\mathbf{p}}}}=\underline{{\mathbf{p}}}(\LU{0}{{\mathbf{u}}}, \varphi)\) and its conjugate \(\myoverline{\underline{{\mathbf{p}}}}\) with the help of the double quaternion product (operator \(\circ\)),
or in short form,
With the multiplication rule for quaternions, the scalar part yields \(0=0\) and the vector part the vector rotation (21). For multiple rotations, we have
Given two quaternions
we provide two important rules of quaternion algebra, which is the sum,
and the multiplication
Angular velocity vector
To derive the angular velocity from Euler parameters, there are both the possibilities of deriving the rotation matrix with respect to time and directly deriving the angular velocity from Euler parameters.
For an infinitesimal rotation \(\dd \psi\), the quaternion follows
from which we conclude that from \(\underline{{\mathbf{p}}}({\mathbf{e}}, \dd \psi) \, \circ \, \underline{{\mathbf{p}}}(t) = \underline{{\mathbf{p}}}(t) + \dd\underline{{\mathbf{p}}}\) and thus \(\dd\underline{{\mathbf{p}}} = \vp{0}{{\mathbf{e}} \frac{\dd\psi}{2}}\). Dividing by \(\dd t\) and with the angular velocity \(\tomega_{10} = {\mathbf{e}} \, \dot \psi\), it follows
or in short form
With the quaternion product in matrix notation, it follows
or with the matrix \(\underline{{\mathbf{P}}}\) it can be written as
This form can be represented in global coordinates due to \({\underline{{\mathbf{P}}}}^{-1} = {\underline{{\mathbf{P}}}}\tp\) as
Thus, the velocity transformation follows to
as a matrix with 3 rows and 4 columns. With \(\LUR{0}{{\mathbf{G}}}{EP}\), the angular velocity vector can be conveniently calculated,
In Exudyn, there are the following functions and items related to Euler parameters:
NodeRigidBodyEP
: A 3D rigid body node based on Euler parametersEulerParameters2RotationMatrix
: computes the rotation matrix from given Euler parametersRotationMatrix2EulerParameters
: computes Euler parameters from a given rotation matrixEulerParameters2G
: computes the matrix \(\LUR{0}{{\mathbf{G}}}{EP}\) relating time derivatives of Euler parameters to the (global) angular velocity vectorEulerParameters2Glocal
: computes the matrix \(\LUR{b}{{\mathbf{G}}}{EP}\) relating time derivatives of Euler parameters to the (local, body-fixed) angular velocity vector