Coordinate representations for rigid parts in multibody dynamics

Abstract

The paper is concerned with coordinate representations for rigid parts in multibody dynamics. The discussion is based on the general theory of the dynamics of multibody systems under constraints. General configuration coordinates are introduced and the requirement of their regularity is discussed. The use of Euler angles, quaternions and linear coordinates, where quaternions and linear coordinates require constraint conditions, is analysed in detail. These coordinate systems are all shown to be regular. Equations of motion are formulated, using Lagrange’s as well as Euler’s equations, and they are supplemented by the appropriate constraint conditions in the cases of quaternions and linear coordinates. Mass matrices are derived, and in terms of the Euler angles the mass matrix components are products of trigonometric functions whereas in terms of quaternions the matrix components are quadratic polynomials. Using linear coordinates gives rise to a constant mass matrix. Thus, there is a decreasing degree of complexity, regarding mass matrix components, when going from Euler angles to linear coordinates. This is obtained at the expense of an increasing gross number of degrees of freedom and the necessary introduction of constraint conditions. The different equations of motion obtained are compared with respect to their structural complexity. In all representations the components of the angular velocity are explicitly calculated. This is not always the case in previous investigations of this subject. The paper also gives a new proof of the well-known relation between angular velocity and unit quaternions and their time derivative.

Keywords

Multibody dynamics rigid body rotation tensor configuration coordinates Euler angles quaternions linear coordinates equations of motion

1. Introduction

The rigid body is an important model concept for parts in multibody systems. The rotational motion of a rigid body may be represented by a trajectory in the special orthogonal group (the group of rotations). The differential equations governing this motion are traditionally obtained by using Lagrange’s or Euler’s equations of motion. These equations may be expressed in terms of coordinates representing the elements in the rotation group, a group which is a three-dimensional manifold. It is then sufficient to use three real-valued parameters for the coordinate representation. Examples of such coordinate patches are, for instance, Euler or Bryant angles and variants thereof. However, it is not possible to cover all of the rotation group by using only one coordinate patch. The rotation group is not in a one-to-one smooth correspondence with an open set in three-dimensional Euclidean space. Singularities will appear and there may be the need for more than one three-dimensional coordinate patch in order to cover the entire motion. One way to avoid these singularities is to use quaternions as coordinates for the rotation group. This is done at the expense of having to deal with a constraint condition guaranteeing that the quaternions used are so-called unit quaternions. Besides the quaternions there are other possibilities as well.

The use of Euler angles in rigid body dynamics originates with Euler [1]. A classical example, where Euler angles are traditionally employed, is the heavy symmetrical (Lagrange’s) top; see Goldstein et. al. [2]. Bryant angles are often used for representing vehicle and spacecraft attitudes in cases where only a part of the rotation group is embarked upon by the motion. An example of the use of Bryant angles, in connection with Gough–Stewart platforms, is given in Afzali-Far et. al. [3]. The innovation and use of quaternions in mechanics originates with Hamilton [4]. A predecessor of quaternions appeared in Euler [5] and Rodrigues [6] in the form of Euler–Rodrigues parameters. Quaternions are, with some preference, used when the entire rotation group is entered by the motion. By using quaternions one avoids the singularities which afflict the Euler angles. Many applications, using quaternions, are found in robotics (see Selig [7]) as well as in attitude control of space-crafts. There is abundant literature discussing the rotation group and its relation to other abstract groups. A full review of this will not be undertaken here and the reader is referred to the general survey of attitude representations given in Shuster [8]. For the quaternion representation one may consult Altmann [9]. Both these publications have extensive lists of references.

The motion of a rigid body $B$ takes place in the three-dimensional Euclidean point space $E$ with corresponding translation vector space $V$ . Let $SO (V)$ denote the rotation group on $V$ . Furthermore, let $ν \in N$ and let $q \in Ω \subseteq ℝ^{ν \times 1}$ , where $Ω$ is an open set, be coordinates representing an element $R \in SO (V)$ through a smooth mapping $Ω ∋ q \to R = R (q)$ . The rotational motion of a rigid body may then be represented by a mapping $[0, \infty [\mapsto R = R (q (t))$ where the function $q : [0, \infty [\to ℝ^{ν \times 1}$ is assumed to be differentiable. Whether or not Euler’s or Lagrange’s methods are used, the equations representing the rotational motion of the rigid body may, in general, be put on the form

A \ddot{q} = F (t, q, \overset{\cdot}{q})

(1)

where $A : Ω \to ℝ^{ν \times ν}$ and $F : [0, \infty [\times Ω \times ℝ^{ν \times 1} \to ℝ^{ν \times 1}$ . Here $A = A (q) \in ℝ^{ν \times ν}$ is a non-zero matrix which may, in some cases, be identified with the so-called mass matrix and $F (t, q, \overset{\cdot}{q}) \in ℝ^{ν \times 1}$ is a column vector representing the force sum of inertial and external forces and moments acting on the body. The character of $A$ and $F$ will depend on the coordinate system used as well as on whether the Euler or Lagrange method is used to formulate the equations. An obvious strategy, in order to obtain simple equations, would be to choose configuration coordinates so that these matrices become as simple functions of $(t, q, \overset{\cdot}{q})$ as possible. A major simplification is obtained if one may, for instance, choose coordinates which make the mass matrix constant, i.e. independent of $q$ . This will, however, usually lead to a more complex expression for the force sum in terms of $(t, q, \overset{\cdot}{q})$ . In general, associated with a choice of configuration coordinates, there will be a payoff between $A$ and $F$ concerning the complexity. It is, in any case, desirable to obtain a non-singular matrix $A$ . If one chooses coordinates $q \in Ω \subseteq ℝ^{ν \times 1}$ , where $ν > 3$ , then constraint conditions will have to be imposed on these coordinates since

R (q) R {(q)}^{T} = 1_{V}, q \in Ω

(2)

where $1_{V}$ is the unit tensor on $V$ .

In the case of Euler angles, $ν = 3$ and the condition expressed in equation (2) is identically satisfied. In the case of quaternions, $ν = 4$ and there will be one constraint condition which is equivalent to the requirement that $q$ is a unit quaternion. A third case is obtained by choosing as coordinates the nine components of the rotation tensor, in some orthonormal basis. These coordinates are, in the present paper, called linear coordinates and they will have to fulfil six constraint conditions imposed by equation (2). Euler angles and linear coordinates are opposite extremes in terms of the number of coordinates used while quaternions represent a middle course. Using Euler angles will result in a mass matrix with components which are non-linear functions of the coordinates, typically products of sines and cosines. Using quaternions, the mass matrix elements will be quadratic polynomials in the coordinates and using linear coordinates will result in a constant mass matrix. So, from this perspective, linear coordinates or quaternions seem to be a better choice than Euler angles. In addition to simpler mass matrices, the avoidance of singularities is obtained by making these choices. The drawback is, of course, that one has to handle a larger number of configuration coordinates and the even more serious problem that one has to consider the appropriate constraint conditions when formulating the equations of motion and in connection with their solution.

In this paper, we derive the equations of motion for a rigid body in terms of Euler angles, quaternions and linear coordinates, respectively. The derivations are obtained by using Lagrange’s equations as well as Euler’s equations. In the case of Lagrange the constraint equations will give rise to Lagrangian multipliers in the equations of motion; see Lidström [10]. Equations based on Euler angles are found in almost every book on advanced mechanics such as Goldstein et. al. [2] and it will here be used primarily as a reference. The use of quaternions is more rarely discussed in mechanics textbooks. Quaternions in combination with Euler’s equations are presented in the papers by Rapaport [11], Chou [12] and Koshlyakov [13]. Examples of more recent papers using Lagrange’s equations are, for instance, Betsch and Siebert [14], Udwadia and Schutte [15], Zoric et. al. [16] and Choi and Yang [17]. Linear coordinates are dealt with in Betch and Siebert [14] and Lidström [18] as well as in Krenk and Nielsen [19]. In the last paper mentioned, equations of motion are formulated using Hamilton’s equations. Some of these papers deal with numerical procedures illustrated by solving the equations of motion for specific problems, such as the Lagrange’s top.

The present paper is not concerned with numerical procedures. Instead we focus on the basic requirements on coordinate representations such as their regularity and the specific structure that the choice of coordinates gives to the equations of motion and constraint equations. It is demonstrated that Euler angles, quaternions and linear coordinates all constitute regular coordinate systems. Regularity is equivalent to positive definiteness of the mass matrix, a property which is of theoretical as well as numerical importance. The equations of motion are formulated using Lagrange’s method as well as Euler’s and it turns out that in the case of Euler angles the two methods give equations of similar complexity, however with a small advantage for Lagrange’s method since it delivers $A$ as the mass matrix which is symmetric. Using quaternions does favour Euler’s equations through a simpler and more direct incorporation of the constraint condition. Furthermore the matrix $A$ is, in this case, orthonormal. Linear coordinates seem, however, to be more convenient when using Lagrange’s equations. This is partly due to the fact that the mass matrix is constant in this case.

Using Lagrange’s method requires a computation of the generalized external force. It is shown that this force may be calculated by invoking the astatic load tensor

A = \int_{\partial B_{t}} (x - x_{c}) \otimes t da (x) + \int_{B_{t}} (x - x_{c}) \otimes b ρ dv (x)

(3)

corresponding to the contact and body forces acting on the rigid body $B$ . In the case of Euler angles, the part of the tensor $A$ which comes into play is the axial vector of the skew part of $A$ which is, apart from a constant factor, equal to the moments of the external forces. In the case of quaternions, there is also an additional contribution from the trace of $A$ , and in the case of linear coordinates, the entire astatic load tensor is contributing. These results have not been previously reported in the literature.

A new proof of the well-known relation between angular velocity and unit quaternions is presented and we calculate the mass matrix in terms of quaternions and linear coordinates. In all representations the components of the angular velocity are explicitly presented. This is not usually the case in previous investigations of this subject.

The paper is organized as follows. In section 2, the basic notations and definitions used throughout the paper are summarized. In section 3, the kinematics of the rigid body are reviewed. The definition of a regular coordinate system is given and a proposition giving necessary and sufficient conditions for regularity is presented. In section 4, the equations of motion are stated in the Lagrangian as well as the Eulerian form. The equations are accompanied by appropriate constraint condition corresponding to the specific choice of configuration coordinates. In section 5, the calculation of the mass matrix is discussed and equations of motion, based on the Lagrange-d’Alembert equation and with Lagrangian multipliers eliminated, are accounted for. In sections 6, 7 and 8, a specialization is made to Euler angles, quaternions and linear coordinates, respectively.

To facilitate the reading of this paper some of the more tedious calculations have been filed in an appendix along with basic mathematical results concerning quaternions.

2. Notations

In this paper, $N$ denotes the set of natural numbers (not including 0) and $ℝ$ denotes the set of real numbers. The set of m-dimensional, real column vectors is denoted by $ℝ^{m \times 1}$ and the null vector in $ℝ^{m \times 1}$ is written $0_{m \times 1}$ . $ℝ^{m \times n}$ denotes the set of real matrices of order $m \times n$ with the null matrix written $0_{m \times n}$ . If $A \in ℝ^{m \times n}$ , then $A^{T} \in ℝ^{n \times m}$ is the transpose of $A$ . The rank of a matrix $A \in ℝ^{m \times n}$ is written $rank A$ . If $A$ is a square matrix, then $tr A$ denotes the trace of $A$ and $\det A$ denotes the determinant of $A$ and if $\det A \neq 0$ , then $A^{- 1}$ denotes its inverse. $I_{n \times n}$ denotes the identity matrix in $ℝ^{n \times n}$ . Let $E$ denote a three-dimensional Euclidean point space with the corresponding translation vector space $V$ . Points in $E$ are denoted by $x, y, \dots, X, Y, \dots$ and vectors in $V$ are denoted by $a, b, \dots, u, v, \dots$ . The space of all second order tensors $A$ on $V$ , i.e. linear maps $V \to V$ , is denoted $End (V)$ . We write $1_{V}$ for the unit tensor in $End (V)$ and $tr A, det A, A^{- 1}, A^{T}$ for respectively, the trace, determinant, inverse and transpose of $A \in End (V)$ . The scalar, vector and tensor products of two vectors $a and b$ are denoted $a \cdot b, a \times b$ and $a \otimes b$ , respectively. The norm $| a |$ of a vector $a$ is defined by $| a | = \sqrt{a \cdot a}$ . The set of all rotations on $V$ is denoted $SO (V) =$ ${A \in E n d (V) | A^{T} A = A A^{T} = 1_{V}, d e t A = 1}$ . The space of symmetric tensors is denoted $Sym (V) = {A \in End (V) | A^{T} = A}$ and the space of skew-symmetric tensor is denoted $Skew (V) = {A \in End (V) | A^{T} = - A}$ . If $A \in Skew (V)$ , then $ax (A) \in V$ denotes its axial vector, i.e. $Au = ax (A) \times u, u \in V$ . Given $a \in V$ , then the tensor $a \times 1_{V} \in S k e w (V)$ is defined by $(a \times 1_{V}) u = a \times u, u \in V$ . Let $u \in V$ and $A \in End (V)$ and let $\underline{e} = (\begin{matrix} e_{1} e_{2} e_{3} \end{matrix})$ be an orthonormal basis in $V$ , then $u = \underline{e} {[u]}_{\underline{e}}$ and $A = \underline{e} {[A]}_{\underline{e}} {\underline{e}}^{T}$ where ${[u]}_{\underline{e}}$ and ${[A]}_{\underline{e}}$ are the matrices, in the basis $\underline{e}$ , corresponding to $u$ and $A$ , respectively. See appendix A.4. A right handed orthonormal basis is called a RON-basis.

3. Kinematics

Let $B_{0} \subseteq E$ be a reference placement of the rigid body $B$ . A motion of $B$ is given by a function $q : [0, \infty [\to Ω \subseteq ℝ^{n}, n \geq 6$ , where $Ω$ is an open set, and a mapping $χ : Ω \times B_{0} \to E$ defined by

x = χ (q; X) = x_{c} (q) + R (q) (X - X_{c}), X \in B_{0}

(4)

where $c$ denotes the centre of mass of $B$ with reference place $X_{c}$ and present place $x_{c} = x_{c} (q) = χ (q; X_{c})$ . The reference place of a material particle in Euclidean space is denoted by $X$ and its present place is denoted by $x$ . $R = R (q) \in SO (V)$ denotes the rotation tensor for $B$ and $q = (q^{1}, \dots, q^{n}) \in Ω$ are local configuration coordinates for the body.

The system of configuration coordinates is said to be regular for the body $B$ if, for all fixed $q \in Ω$

\sum_{k = 1}^{n} \frac{\partial χ}{\partial q^{k}} w^{k} = 0, w^{k} \in ℝ, \forall X \in B_{0} \Rightarrow w^{k} = 0, k = 1, \dots, n

This is equivalent to the statement that, for all fixed $q \in Ω$ , the functions

\frac{\partial χ}{\partial q^{k}} = \frac{\partial χ}{\partial q^{k}} (q; X) = \frac{\partial x_{c}}{\partial q^{k}} (q) + \frac{\partial R}{\partial q^{k}} (q) (X - X_{c}), X \in B_{0}, k = 1, \dots, n

are linearly independent vector fields on $B_{0}$ .

The spin tensor of the rigid body is defined by

W = \overset{\cdot}{R} R^{T} = \sum_{k = 1}^{n} \frac{\partial R}{\partial q^{k}} R^{T} {\overset{\cdot}{q}}^{k} = \sum_{k = 1}^{n} W_{, k} {\overset{\cdot}{q}}^{k}

(5)

where

W_{, k} = \frac{\partial R}{\partial q^{k}} R^{T}

The angular velocity of the rigid body is defined by $ω = ax (W)$ and then

ω = \sum_{k = 1}^{n} ax (W_{, k}) {\overset{\cdot}{q}}^{k}

(6)

Furthermore

{\overset{\cdot}{x}}_{c} = \frac{d x_{c} (q (t))}{dt} = \sum_{k = 1}^{n} \frac{\partial x_{c}}{\partial q^{k}} {\overset{\cdot}{q}}^{k}

(7)

Proposition 3.1. The velocity field of the rigid part $B$ may be written

\overset{\cdot}{x} = \overset{\cdot}{x} (q, \overset{\cdot}{q}; x) = {\overset{\cdot}{x}}_{c} (q, \overset{\cdot}{q}) + ω (q, \overset{\cdot}{q}) \times (x - x_{c}), x \in B_{t}

(8)

where $B_{t} = χ (q (t); B_{0})$ is the present placement of $B$

{\overset{\cdot}{x}}_{c} = \sum_{k = 0}^{n} x_{c, k} {\overset{\cdot}{q}}^{k} and ω = \sum_{k = 0}^{n} ω_{, k} {\overset{\cdot}{q}}^{k}

(9(a, b))

x_{c, k} = x_{c, k} (q) = \frac{\partial x_{c}}{\partial q^{k}} and ω_{, k} = ω_{, k} (q) = ax (W_{, k}), k = 0, 1, \dots, n

(10(a, b))

Proof: The proof follows directly from equations (7) and (6), respectively. □

Remark 3.1: As a direct consequence of equation (9) we have

x_{, k} = \frac{\partial {\overset{\cdot}{x}}_{c}}{\partial {\overset{\cdot}{q}}^{k}} and ω_{, k} = \frac{\partial ω}{\partial {\overset{\cdot}{q}}^{k}}, k = 0, 1, \dots, n

■

Remark 3.2: Note that relations in equation (9(b)) and (10(b)) are based on the assumption that $R = R (q) \in SO (V)$ for all $q \in Ω$ .▪

Let ${\underline{e}}^{o} = (e_{1}^{o} e_{2}^{o} e_{3}^{o})$ be a RON-basis fixed in an inertial frame. The rotation tensor may be written as (see appendix A.4 for the notation)

R = {\underline{e}}^{o} {[R]}_{{\underline{e}}^{o}} {\underline{e}}^{o}^{T}

(11)

where ${[R]}_{{\underline{e}}^{o}}$ is the rotation matrix relative to the basis ${\underline{e}}^{o}$ .

Proposition 3.2. We have

ω_{, k} = ax (\underline{e} {[R]}_{{\underline{e}}^{o}}^{T} \frac{\partial {[R]}_{{\underline{e}}^{o}}}{\partial q^{k}} {\underline{e}}^{T})

(12)

where $\underline{e} = (e_{1} e_{2} e_{3})$ defined by $\underline{e} = R {\underline{e}}^{o} = {\underline{e}}^{o} {[R]}_{{\underline{e}}^{o}}$ is a basis fixed to the body in its present placement.

Proof: From equation (11) it follows that

\frac{\partial R}{\partial q^{k}} = {\underline{e}}^{o} \frac{\partial {[R]}_{{\underline{e}}^{o}}}{\partial q^{k}} {\underline{e}}^{o}^{T}

and then

\frac{\partial R}{\partial q^{k}} R^{T} = {\underline{e}}^{o} \frac{\partial {[R]}_{{\underline{e}}^{o}}}{\partial q^{k}} {\underline{e}}^{o}^{T} {\underline{e}}^{o} {[R]}_{{\underline{e}}^{o}}^{T} {\underline{e}}^{o}^{T} = {\underline{e}}^{o} \frac{\partial {[R]}_{{\underline{e}}^{o}}}{\partial q^{k}} {[R]}_{{\underline{e}}^{o}}^{T} {\underline{e}}^{o}^{T} = \underline{e} {[R]}_{{\underline{e}}^{o}}^{T} \frac{\partial {[R]}_{{\underline{e}}^{o}}}{\partial q^{k}} {\underline{e}}^{T},

where we have used the relation ${\underline{e}}^{o} = \underline{e} {[R]}_{{\underline{e}}^{o}}^{T}$ and the fact that ${\underline{e}}^{o}^{T} {\underline{e}}^{o} = I_{3 \times 3}$ . The proposition then follows from (10(b)). □

Remark 3.3: Note that

ω_{, k} = ax (\underline{e} {[R]}_{\underline{e}}^{T} \frac{\partial {[R]}_{\underline{e}}}{\partial q^{k}} {\underline{e}}^{T})

which follows from equation (11) and the fact that the rotation matrices in the bases ${\underline{e}}^{o}$ and $\underline{e}$ are equal; that is ${[R]}_{\underline{e}} = {[R]}_{{\underline{e}}^{o}}$ . ▪

The following proposition expresses the classical Euler–Rodrigues representation formula.

Proposition 3.3. To the rotation tensor $R =$ $R (q) \in SO (V)$ , there exists a unit vector $n = n (q)$ and an angle $φ = φ (q) \in ℝ$ so that

R 1_{V} + 2 \sin^{2} \frac{φ}{2} n \times (n \times 1_{V}) + 2 \sin \frac{φ}{2} \cos \frac{φ}{2} n \times 1_{V}

(13)

This representation is mainly unique. If $R = 1_{V}$ , then $n$ may be arbitrarily chosen, while $φ = 0 (modulo 2 π)$ . If $R \neq 1_{V}$ , then $n$ is uniquely determined, apart from a factor $\pm 1$ . For a given $n$ the angle $φ$ is uniquely determined $(modulo 2 π)$ . We write $R ~ n φ$ .

Proof: The proof has been omitted. A demonstration of this classical result may, for instance, be found in Geradin and Cardona [20]. □

A direct consequence of the previous proposition is the following result relating angular velocity to the rotation axis direction $n$ and the rotation angle $φ$ .

Proposition 3.4. If $R ~ n φ$ , then

ω = Φ \overset{\cdot}{n} + n \overset{\cdot}{φ}

(14)

where

Φ = \sin φ 1_{V} + (1 - \cos φ) n \times 1_{V}

(15)

Proof: See Lidström [21]. □

Taking $ν = n - 3$ , we may split the configuration coordinates into two parts $q_{t} = (q^{1}, q^{2}, q^{3}) \in Ω_{t} \subseteq ℝ^{3}$ representing the translational part of the rigid body and $q_{r} = (q^{4}, \dots, q^{n}) \in Ω_{r} \subseteq ℝ^{ν}$ representing the rotational part so that $q = (q_{t}, q_{r}) \in Ω_{t} \times Ω_{r} = Ω$ and equation (3) may then be written

x = χ (q; X) = x_{c} (q_{t}) + R (q_{r}) (X - X_{c}), X \in B_{0}

(16)

The special orthogonal group $SO (V)$ is a three-dimensional manifold naturally embedded in $End (V)$ , a nine-dimensional vector space. Thus, it is appropriate to take $ν = 3$ . It is possible, and sometimes beneficial, to take $ν$ larger than three in order to avoid singularities and make the coordinate representation global. In the case of $ν > 3$ , the coordinates chosen will be related by $m = ν - 3$ constraint conditions

C (q_{r}) = 0_{m \times 1}

(17)

where it is assumed that the mapping $C : Ω_{r} \to ℝ^{m \times 1}$ is continuously differentiable and the matrix $\frac{\partial C}{\partial q_{r}} \in ℝ^{m \times ν}$ has rank $m$ . From equation (17), it follows the kinematical constraints

g (q_{r}) {\overset{\cdot}{q}}_{r} = 0

(18)

where ${\overset{\cdot}{q}}_{r} \in ℝ^{ν \times 1}$ and the constraint matrix is given by

g (q_{r}) = \frac{\partial C (q_{r})}{\partial q_{r}} \in ℝ^{m \times ν}

Using the basis ${\underline{e}}^{o} = (e_{1}^{o} e_{2}^{o} e_{3}^{o})$ we may write

x_{c} (q) = x_{O} + e_{1}^{o} q^{1} + e_{2}^{o} q^{2} + e_{3}^{o} q^{3}, R = R (q^{4}, q^{5}, \dots, q^{n})

(19)

Then

\frac{\partial χ}{\partial q^{1}} = e_{1}^{o}, \frac{\partial χ}{\partial q^{2}} = e_{2}^{o}, \frac{\partial χ}{\partial q^{3}} = e_{3}^{o}, \frac{\partial χ}{\partial q^{k}} = \frac{\partial R}{\partial q^{k}} (q_{r}) (X - X_{c}), k = 4, \dots, n

(20)

Proposition 3.5. Let $B$ be a rigid body and let $n \geq 6$ . The system of q-coordinates, defined by equation (19), is regular if and only if

\frac{\partial R}{\partial q^{4}}, \dots, \frac{\partial R}{\partial q^{n}} \in End (V)

(21)

are linearly independent tensors for all $q_{r} \in Ω_{r}$ .

Proof: For the ‘if-part’, assume that the tensors in equation (21) are linearly independent and consider the linear combination ( $w^{k} \in ℝ$ )

e_{1}^{o} w^{1} + e_{2}^{o} w^{2} + e_{3}^{o} w^{3} + \sum_{k = 4}^{n} w^{k} \frac{\partial R}{\partial q^{k}} (X - X_{c}) = 0, X \in B_{0}

First take $X = X_{c}$ . Then $e_{1}^{o} w^{1} + e_{2}^{o} w^{2} + e_{3}^{o} w^{3} = 0 \Rightarrow w^{1} = w^{2} = w^{3} = 0$ and consequently

\sum_{k = 4}^{n} w^{k} \frac{\partial R}{\partial q^{k}} (X - X_{c}) = (\sum_{k = 4}^{n} w^{k} \frac{\partial R}{\partial q^{k}}) (X - X_{c}) = 0, X \in B_{0} \Rightarrow \sum_{k = 4}^{n} w^{k} \frac{\partial R}{\partial q^{k}} = 0

and from this it follows that $w^{4} = w^{5} = \dots = w^{n} = 0$ . For the ‘only if-part’ assume that the fields

\frac{\partial χ}{\partial q^{k}} = \frac{\partial R}{\partial q^{k}} (X - X_{c}), X \in B_{0}, k = 4, \dots, n

are linearly independent. Then

\sum_{k = 4}^{n} w^{k} \frac{\partial R}{\partial q^{k}} = 0 \Rightarrow (\sum_{k = 4}^{n} w^{k} \frac{\partial R}{\partial q^{k}}) (X - X_{c}) = 0, X \in B_{0} \Rightarrow \sum_{k = 4}^{n} w^{k} \frac{\partial R}{\partial q^{k}} (X - X_{c}) = 0, X \in B_{0} \Rightarrow

w^{4} = w^{5} = \dots = w^{n} = 0

which proves the proposition. □

Corollary 3.1. The coordinate system is regular only if $n \leq 12$ .

Proof: The proof is an obvious consequence of the previous proposition, since $\dim (End (V)) = 9$ . □

Remark 3.4: From this corollary it follows that the options available for the number of coordinates, in order to obtain a regular coordinate system, are $3 \leq ν \leq 9$ . In this paper we will consider the cases $ν = 3, 4$ and $9$ . ▪

4. Rigid body dynamics using Lagrange–d’Alembert equations

Let $B$ be a rigid body with the mass density, $ρ_{0} = ρ_{0} (X), X \in B_{0}$ . The total mass of the body is given by

m = \int_{B_{0}} ρ_{0} (X) dv (X)

(21)

where $v$ is the Euclidean volume measure.

Proposition 4.1. The kinetic energy of the rigid body may be written

T = T (q, \overset{\cdot}{q}) = \frac{1}{2} {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q}

(22)

where the mass matrix $M_{2}$ is given by

M_{2} = M_{2} (q) = (\begin{matrix} I_{3 \times 3} m & 0_{3 \times ν} \\ 0_{ν \times 3} & M \end{matrix}) \in ℝ^{n \times n}, M = M (q) = (\begin{matrix} m_{11} & m_{12} & \dots & m_{1 ν} \\ m_{21} & m_{22} & \dots & m_{2 ν} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ m_{ν 1} & m_{ν 2} & \dots & m_{ν ν} \end{matrix}) \in ℝ^{ν \times ν}

(23)

Let ${\underline{e}}^{o} = (e_{1}^{o} e_{2}^{o} e_{3}^{o})$ be an orthonormal basis. The sub-matrix $M$ has components

m_{kl} = \sum_{h, i, j = 1}^{3} \frac{\partial R_{hi}}{\partial q^{k}} \frac{\partial R_{hj}}{\partial q^{l}} e_{i}^{o} \cdot E_{c}^{0} e_{j}^{o}, k, l = 1, \dots, ν

(24)

where $R_{ij} = e_{i}^{o} \cdot {Re}_{j}^{o}$ and

E_{c}^{0} = \int_{B_{0}} (X - X_{c}) \otimes (X - X_{c}) ρ_{0} (X) dv (X)

(25)

is the Euler tensor in the reference placement. We have $M \in Sym (ℝ^{ν \times ν})$ and $E_{c}^{0} \in Sym (V)$ .

Proof: By definition

T = \frac{1}{2} \int_{B_{0}} \overset{\cdot}{x} \cdot \overset{\cdot}{x} ρ_{0} dv (X) = \frac{1}{2} \int_{B_{0}} (\sum_{k = 1}^{n} \frac{\partial χ}{\partial q^{k}} {\overset{\cdot}{q}}^{k}) \cdot (\sum_{l = 1}^{n} \frac{\partial χ}{\partial q^{l}} {\overset{\cdot}{q}}^{l}) ρ_{0} dv (X) = \frac{1}{2} \sum_{k, l = 1}^{n} m_{kl} {\overset{\cdot}{q}}^{k} q^{l} = \frac{1}{2} {\overset{\cdot}{q}}^{T} M_{2} \overset{\cdot}{q}

where

m_{kl} = \int_{B_{0}} \frac{\partial χ}{\partial q^{k}} \cdot \frac{\partial χ}{\partial q^{l}} ρ_{0} dv (X), k, l = 1, \dots, n

Now using equation (4) and the representation $R = \sum_{g, i = 1}^{3} e_{g}^{o} \otimes e_{i}^{o} R_{gi}, R_{gi} = R_{gi} (q^{4}, q^{5}, \dots, q^{n})$ , one obtains

\frac{\partial χ}{\partial q^{k}} = \frac{\partial x_{c}}{\partial q^{k}} + \sum_{g, i = 1}^{3} e_{g}^{o} e_{i}^{o} \cdot (X - X_{c}) \frac{\partial R_{gi}}{\partial q^{k}}, X \in B_{0}, k = 1, \dots, n

If $k, l = 1, 2, 3$ , then the mass matrix elements are given by

m_{kl} = \int_{B_{0}} \frac{\partial χ}{\partial q^{k}} \cdot \frac{\partial χ}{\partial q^{l}} ρ_{0} dv (X) = \int_{B_{0}} e_{k}^{o} \cdot e_{l}^{o} ρ_{0} dv (X) = δ_{kl} m

Furthermore, if $k, l = 4, \dots, n$ , then the mass matrix elements are given by

m_{kl} = \int_{B_{0}} \frac{\partial χ}{\partial q^{k}} \cdot \frac{\partial χ}{\partial q^{l}} ρ_{0} dv (X) = \int_{B_{0}} (\sum_{g, i = 1}^{3} e_{g}^{o} e_{i}^{o} \cdot (X - X_{c}) \frac{\partial R_{gi}}{\partial q^{k}}) \cdot (\sum_{h, j = 1}^{3} e_{h}^{o} e_{j}^{o} \cdot (X - X_{c}) \frac{\partial R_{hj}}{\partial q^{l}}) ρ_{0} dv (X) =

\int_{B_{0}} (\sum_{g, i, h, j = 1}^{3} e_{g}^{o} \cdot e_{h}^{o} e_{i}^{o} \cdot (X - X_{c}) \frac{\partial R_{gi}}{\partial q^{k}} e_{j}^{o} \cdot (X - X_{c}) \frac{\partial R_{hj}}{\partial q^{l}}) ρ_{0} dv (X) =

\int_{B_{0}} (\sum_{h, i, j = 1}^{3} e_{i}^{o} \cdot (X - X_{c}) \frac{\partial R_{hi}}{\partial q^{k}} e_{j}^{o} \cdot (X - X_{c}) \frac{\partial R_{hj}}{\partial q^{l}}) ρ_{0} dv (X) =

\sum_{h, i, j = 1}^{3} \frac{\partial R_{hi}}{\partial q^{k}} \frac{\partial R_{hj}}{\partial q^{l}} e_{i}^{o} \cdot \int_{B_{0}} (X - X_{c}) \otimes (X - X_{c}) ρ_{0} dv (X) e_{j}^{o} = \sum_{h, i, j = 1}^{3} \frac{\partial R_{hi}}{\partial q^{k}} \frac{\partial R_{hj}}{\partial q^{l}} e_{i}^{o} \cdot E_{c}^{0} e_{j}^{o}

which proves the proposition. □

Since $E_{c}^{0} \in Sym (V)$ there is a RON-basis ${\underline{i}}^{o} = (i_{1}^{o} i_{2}^{o} i_{3}^{o})$ , the principal basis, fixed to the body in its reference placement such that $E_{c}^{0} i_{k}^{o} = i_{k}^{o} E_{k}, k = 1, 2, 3$ where $E_{1}, E_{2}, E_{3}$ are constant numbers.

Corollary 4.1. Let ${\underline{i}}^{o}$ be a principal basis for $E_{c}^{0}$ , then

m_{kl} = \sum_{j = 1}^{3} E_{j} (\sum_{i = 1}^{3} \frac{\partial R_{ij}}{\partial q^{k}} \frac{\partial R_{ij}}{\partial q^{l}}), k, l = 4, \dots, n

(26)

where

E_{j} = \int_{B_{0}} {(i_{j}^{o} \cdot (X - X_{c}))}^{2} ρ_{0} (X) dv (X), j = 1, 2, 3

are positive real numbers.

Proof: We obtain, by taking ${\underline{e}}^{o} = {\underline{i}}^{o}$ in equation (24)

\begin{array}{l} m_{k l} = \sum_{h, i, j = 1}^{3} \frac{\partial R_{h i}}{\partial q^{k}} \frac{\partial R_{h j}}{\partial q^{l}} i_{i}^{o} \cdot E_{c}^{0} i_{j}^{o} = \sum_{h, i, j = 1}^{3} \frac{\partial R_{h i}}{\partial q^{k}} \frac{\partial R_{h j}}{\partial q^{l}} i_{i}^{o} \cdot i_{j}^{o} E_{j} = \sum_{h, i, j = 1}^{3} \frac{\partial R_{h i}}{\partial q^{k}} \frac{\partial R_{h j}}{\partial q^{l}} δ_{i j} E_{j} = \sum_{h, j = 1}^{3} \frac{\partial R_{h j}}{\partial q^{k}} \frac{\partial R_{h j}}{\partial q^{l}} E_{j} \\ = \sum_{j = 1}^{3} E_{j} (\sum_{i = 1}^{3} \frac{\partial R_{i j}}{\partial q^{k}} \frac{\partial R_{i j}}{\partial q^{l}}) \end{array}

where

E_{j} = i_{j}^{o} \cdot E_{c}^{0} i_{j}^{o} = \int_{B_{0}} {(i_{j}^{o} \cdot (X - X_{c}))}^{2} ρ_{0} (X) dv (X) > 0

This proves the corollary. □

Proposition 4.2. The regularity of the configuration coordinates is equivalent to the positive definiteness of the mass matrix.

Proof: See Lidström [22].

Proposition 4.3. If the system of configuration coordinates is regular, then the matrix $M$ is positive definite.

Proof: By taking ${\overset{\cdot}{q}}_{t} = 0_{3 \times 1}$ in equation (22) one obtains

T = \frac{1}{2} {\overset{\cdot}{q}}_{r}^{T} M {\overset{\cdot}{q}}_{r} \geq 0, \forall {\overset{\cdot}{q}}_{r} \in ℝ^{ν \times 1} .

Furthermore, since $M_{2}$ is positive definite

det M_{2} = det (I_{3 \times 3} m) det M = m^{3} det M \neq 0 \Rightarrow det M \neq 0

and this proves the proposition. □

4.1. The equations of motion

The Lagrange–d’Alembert equations of motion in multibody dynamics are presented in Lidström [10]. When applied to rigid body dynamics with a coordinate representation according to equation (19) the translational and rotational degrees of freedom $(q_{t}, q_{r})$ may be separated according to the following theorem.

Theorem 4.1. We have

{\begin{matrix} m {\ddot{q}}_{t} = Q^{t} \\ M {\ddot{q}}_{r} - (I_{ν \times ν} - g^{T} Γ^{- 1} g M^{- 1}) Q^{r, res} + g_{r}^{T} Γ^{- 1} (\frac{\partial g}{\partial q_{r}} {\overset{\cdot}{q}}_{r}) {\overset{\cdot}{q}}_{r} = 0_{ν \times 1} \\ q (0) = q_{0}, \overset{\cdot}{q} (0) = {\overset{\cdot}{q}}_{0} \\ g (q_{r, 0}) {\overset{\cdot}{q}}_{r, 0} = 0_{m \times 1} \end{matrix}

(27)

where $m$ is the total mass of the rigid body and $M$ is the mass matrix, according to equations (23) and (24), corresponding to the rotational degrees of freedom. Furthermore $Γ = g M^{- 1} g^{T} \in ℝ^{m \times m}$ and

Q^{r, res} = - {\overset{\cdot}{q}}_{r}^{T} (\frac{\partial M^{- 1}}{\partial q_{r}} - \frac{1}{2} {(\frac{\partial M^{- 1}}{\partial q_{r}})}^{t}) {\overset{\cdot}{q}}_{r} + Q^{r}

(28)

is the residual generalized force and $Q = (\begin{matrix} Q^{t} & Q^{r} \end{matrix}) \in ℝ^{1 \times n}$ is the generalized external force on the body. The mass matrix and the external force will be discussed in the following sections.

Proof: This is a direct consequence of theorem 6.1 in Lidström [10] and proposition 4.1. □

Remark 4.1: The Lagrangian multipliers, corresponding to constraint forces, are given by

λ = - Γ^{- 1} (g M^{- 1} Q^{r, resT} + (\frac{\partial g}{\partial q_{r}} {\overset{\cdot}{q}}_{r}) {\overset{\cdot}{q}}_{r})

■

4.2. The generalized external force

The force and moment sums of the external forces on the body may be written

F = \int_{\partial B_{t}} t da (x) + \int_{B_{t}} b ρ dv (x)

(29)

M_{c} = \int_{\partial B_{t}} (x - x_{c}) \times t da (x) + \int_{B_{t}} (x - x_{c}) \times b ρ dv (x)

where $t = t (x, t)$ is the external traction vector, $b = b (x, t)$ is the body force vector, $ρ = ρ (x, t)$ is the mass density and $a$ is the Euclidean area measure.

Proposition 4.4. The components of the external generalized force are given by

Q_{k}^{t} = e_{k}^{o} \cdot F, k = 1, 2, 3, Q_{k}^{r} = W_{, k}^{T} \cdot A, k = 4, \dots, n

(30)

where

A = \int_{\partial B_{t}} (x - x_{c}) \otimes t da (x) + \int_{B_{t}} (x - x_{c}) \otimes b ρ dv (x)

(31)

is the astatic load tensor corresponding to the external forces acting on part $B$ .

Proof: By using equation (20) one obtains for $k = 1, 2, 3$

Q_{k} = \int_{\partial B_{0}} \frac{\partial χ}{\partial q^{k}} \cdot t_{0} da (X) + \int_{B_{0}} \frac{\partial χ}{\partial q^{k}} \cdot b ρ_{0} dv (X) = \int_{\partial B_{0}} e_{k}^{o} \cdot t_{0} da (X) + \int_{B_{0}} e_{k}^{o} \cdot b ρ_{0} dv (X) = e_{k}^{o} \cdot F

Furthermore, for $k = 4, \dots, n$

Q_{k} = \int_{\partial B_{0}} \frac{\partial χ}{\partial q^{k}} \cdot t_{0} da (X) + \int_{B_{0}} \frac{\partial χ}{\partial q^{k}} \cdot b ρ_{0} dv (X) = \int_{\partial B_{0}} \frac{\partial R (q)}{\partial q^{k}} (X - X_{c}) \cdot t_{0} da (X) +

\begin{matrix} \int_{B_{0}} \frac{\partial R (q)}{\partial q^{k}} (X - X_{c}) \cdot b ρ_{0} dv (X) = {(\frac{\partial R (q)}{\partial q^{k}})}^{T} \cdot (\int_{\partial B_{0}} (X - X_{c}) \otimes t_{0} da (X) \\ + \int_{B_{0}} (X - X_{c}) \otimes b ρ_{0} dv (X)) = R (q) {(\frac{\partial R (q)}{\partial q^{k}})}^{T} \cdot (\int_{\partial B_{0}} R (q) (X - X_{c}) \otimes t_{0} da (X) + \int_{B_{0}} R (q) (X - X_{c}) \otimes b ρ_{0} dv (X)) = \end{matrix}

W_{, k}^{T} \cdot (\int_{\partial B_{t}} (x - x_{c}) \otimes t da (x) + \int_{B_{0}} (x - x_{c}) \otimes b ρ dv (x)) = W_{, k}^{T} \cdot A

where $A$ is given by equation (31). □

Remark 4.2: With

A_{0} = \int_{\partial B_{0}} (X - X_{c}) \otimes t_{0} d a (X) + \int_{B_{0}} (X - X_{c}) \otimes b ρ_{0} d v (X)

representing the astatic load tensor in the reference placement, we have $A = R A_{0}$ and

Q_{k}^{r} = {(\frac{\partial R (q)}{\partial q^{k}})}^{T} \cdot A_{0}}, k = 4, \dots, n

▪

5. Rigid body dynamics using Euler equations

The Euler tensor in the present placement of the body, $E_{c}$ , is defined by

E_{c} = \int_{B_{t}} (x - x_{c}) \otimes (x - x_{c}) ρ (x, t) dv (x) \in Sym (V)

(32)

where $ρ (x, t) = ρ_{0} (X), X \in B_{0},$ $x = χ (q; X) \in B_{t}$ is the mass density in the present placement. We have the relation $E_{c} = {RE}_{c}^{0} R^{T}$ , where $E_{c}^{0}$ is defined in equation (25). Let $\underline{i} = (i_{1} i_{2} i_{3})$ denote the principal directions in the present placement, i.e. $\underline{i} = R {\underline{i}}^{0}$ then $E_{c} i_{k} = i_{k} E_{k}, k = 1, 2, 3$ . The inertia tensor, $I_{c}$ , of the body, at the centre of mass, is defined by $I_{c} = (t r E_{c}) 1_{V} - E_{c}$ . It then follows that $I_{c} i_{k} = i_{k} I_{k}, k = 1, 2, 3$ where the principal moments of inertia $I_{1}, I_{2}, I_{3}$ are constant, positive real numbers given by

I_{1} = E_{2} + E_{3}, I_{2} = E_{3} + E_{1}, I_{3} = E_{1} + E_{3}

(33)

The Euler equations for the motion of the rigid body are as follows (see Arnold [23])

F = {\ddot{x}}_{c} m, M_{c} = I_{c} \overset{\cdot}{ω} + ω \times I_{c} ω

(34(a, b))

where $F$ is the force sum and $M_{c}$ the moment sum of the external forces acting on the body according to equation (29). Using the basis $\underline{i}$ we may write

M_{c} = i_{1} M_{1} + i_{2} M_{2} + i_{3} M_{3}, I_{c} = i_{1} \otimes i_{1} I_{1} + i_{2} \otimes i_{2} I_{2} + i_{3} \otimes i_{3} I_{3}, ω = i_{1} ω_{1} + i_{2} ω_{2} + i_{3} ω_{3}

(35)

Equation (34(b)) is then equivalent to

{\begin{matrix} M_{1} = I_{1} {\overset{\cdot}{ω}}_{1} + (I_{3} - I_{2}) ω_{2} ω_{3} \\ M_{2} = I_{2} {\overset{\cdot}{ω}}_{2} + (I_{1} - I_{3}) ω_{3} ω_{1} \\ M_{3} = I_{3} {\overset{\cdot}{ω}}_{3} + (I_{2} - I_{1}) ω_{1} ω_{2} \end{matrix}

(36)

If a coordinate system with $ν > 3$ is used, then equation (36) should be complemented by a kinematical constraint condition according to equation (18).

6. Euler angles

Coordinate representations using Euler angles are found in many textbook on rigid body dynamics. See, for instance, Goldstein et. al. [2] and Arnold [23]. We give a presentation here which is adapted to the general approach of this paper. In this case, $ν = 3$ . If $x_{c} (q) = x_{c} (q^{1}, q^{2}, q^{1}) =$ $x_{O} + i_{1}^{o} q^{1} + i_{2}^{o} q^{2} + i_{3}^{o} q^{3}$ , $q_{r}^{1} = ψ$ , $q_{r}^{2} = θ$ , $q_{r}^{3} = ϕ$ ,

x = χ (q; X) = x_{c} (q^{1}, q^{2}, q^{3}) + R (ψ, θ, ϕ) (X - X_{c}), X \in B_{0}, (ψ, θ, ϕ) \in Ω_{r}

(37)

where $Ω_{r} = {(ψ, θ, ϕ) \in ℝ^{3} | 0 \leq ψ < 2 π, 0 < θ < π, 0 \leq ϕ < 2 π}$ and

{[R]}_{{\underline{i}}^{o}} = (\begin{matrix} \cos ψ \cos ϕ - \sin ψ \cos θ \sin ϕ & - \cos ψ \sin ϕ - \sin ψ \cos θ \cos ϕ & \sin ψ \sin θ \\ \sin ψ \cos ϕ + \cos ψ \cos θ \sin ϕ & - \sin ψ \sin ϕ + \cos ψ \cos θ \cos ϕ & - \cos ψ \sin θ \\ \sin θ \sin ϕ & \sin θ \cos ϕ & \cos θ \end{matrix})

(38)

then the coordinates, $(ψ, θ, ϕ)$ , are called the Euler angles.

Proposition 6.1. The Euler angle coordinate system is regular.

Proof: By a direct computation one finds

\frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ψ} = (\begin{matrix} - \sin ψ \cos ϕ - \cos ψ \cos θ \sin ϕ & \sin ψ \sin ϕ - \cos ψ \cos θ \cos ϕ & \cos ψ \sin θ \\ \cos ψ \cos ϕ - \sin ψ \cos θ \sin ϕ & - \cos ψ \sin ϕ - \sin ψ \cos θ \cos ϕ & \sin ψ \sin θ \\ 0 & 0 & 0 \end{matrix})

\frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial θ} = (\begin{matrix} \sin ψ \sin θ \sin ϕ & \sin ψ \sin θ \cos ϕ & \sin ψ \cos θ \\ - \cos ψ \sin θ \sin ϕ & - \cos ψ \sin θ \cos ϕ & - \cos ψ \cos θ \\ \sin ϕ \cos θ & \cos ϕ \cos θ & - \sin θ \end{matrix})

(39 (a, b, c))

\frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ϕ} = (\begin{matrix} - \cos ψ \sin ϕ - \sin ψ \cos θ \cos ϕ & - \cos ψ \cos ϕ + \sin ψ \cos θ \sin ϕ & 0 \\ - \sin ψ \sin ϕ + \cos ψ \cos θ \cos ϕ & - \sin ψ \cos ϕ - \cos ψ \cos θ \sin ϕ & 0 \\ \cos ϕ \sin θ & - \sin ϕ \sin θ & 0 \end{matrix}),

and then one obtains

w^{4} \frac{\partial R}{\partial ψ} + w^{5} \frac{\partial R}{\partial θ} + w^{6} \frac{\partial R}{\partial ϕ} = 0 \Rightarrow {\begin{matrix} w^{5} \sin ϕ \cos θ + w^{6} \cos ϕ \sin θ = 0 \\ w^{4} \cos ψ \sin θ + w^{5} \sin ψ \cos θ = 0 \\ \forall (ψ, θ, ϕ) \in Ω_{r} \end{matrix}

which implies $w^{4} = w^{5} = w^{6} = 0$ and then, according to proposition 3.5, the Euler angle coordinate system is regular. □

Remark 6.1: Note that if ${[R]}_{{\underline{i}}^{o}} = (R_{ij})$ then

\begin{matrix} \cos θ = R_{33}, \sin θ = \sqrt{1 - \cos^{2} θ}, \cos ψ = - \frac{R_{23}}{\sin θ}, \sin ψ = \frac{R_{13}}{\sin θ} \\ \cos ϕ = \frac{R_{32}}{\sin θ}, \sin ϕ = \frac{R_{31}}{\sin θ} \end{matrix}

These formulas show that singularities are present for values of $θ$ equal to $n π, n = 0, 1, \dots$ ▪

Proposition 6.2. The mass matrix, related to the Euler angle coordinate system, is given by

M = (\begin{matrix} I_{1} \sin^{2} θ \sin^{2} ϕ + I_{2} \sin^{2} θ \cos^{2} ϕ + I_{3} \cos^{2} θ & (I_{1} - I_{2}) \sin θ \sin ϕ \cos ϕ & I_{3} \cos θ \\ (I_{1} - I_{2}) \sin θ \sin ϕ \cos ϕ & I_{1} \cos^{2} ϕ + I_{2} \sin^{2} ϕ & 0 \\ I_{3} \cos θ & 0 & I_{3} \end{matrix})

(40)

where $M$ is positive definite.

Proof: See appendix A.2.1 for the computation of $M$ . Since the Euler angle coordinate system is regular, $M$ is positive definite. □

In the case of Euler angles there are no constraints present which means that $g = 0_{1 \times 3}$ and the Lagrange–d’Alembert equations in equation (27), corresponding to the rotational motion, may be written according to the following proposition.

Proposition 6.3. The rotational motion of the rigid body satisfies

{\begin{matrix} M \ddot{q} = Q^{resT} \\ q (0) = q_{0}, \overset{\cdot}{q} (0) = {\overset{\cdot}{q}}_{0} \end{matrix}

(41)

where $q = {(\begin{matrix} ψ & θ & ϕ \end{matrix})}^{T}$ ; $M$ is given by equation (40) and the components of $Q^{res}$ are

\begin{matrix} Q_{ψ}^{res} = \sin 2 θ (I_{1} \sin^{2} ϕ + I_{2} \cos^{2} ϕ - I_{3}) \overset{\cdot}{θ} \overset{\cdot}{ψ} - \frac{1}{2} (I_{1} - I_{2}) \cos θ \sin 2 ϕ {\overset{\cdot}{θ}}^{2} - (I_{1} - I_{2}) \sin 2 ϕ \sin^{2} θ \overset{\cdot}{ϕ} \overset{\cdot}{ψ} \\ + \sin θ (I_{3} - (I_{1} - I_{2}) \cos 2 ϕ) \overset{\cdot}{ϕ} \overset{\cdot}{θ} + M_{1} \sin θ \sin ϕ + M_{2} \sin θ \cos ϕ + M_{3} \cos θ \end{matrix}

(42)

\begin{array}{l} Q_{θ}^{r e s} = \frac{1}{2} \sin 2 θ (I_{1} \sin^{2} ϕ + I_{2} \cos^{2} ϕ - I_{3}) {\dot{ψ}}^{2} - (I_{2} - I_{1}) \sin 2 ϕ \dot{ϕ} \dot{θ} - ((I_{1} - I_{2}) \sin θ \cos 2 ϕ \\ + \frac{3}{2} I_{3} \sin θ) \dot{ϕ} \dot{ψ} + M_{1} \cos ϕ - M_{2} \sin ϕ \end{array}

Q_{ϕ}^{res} = \frac{1}{2} (I_{1} - I_{2}) \sin 2 ϕ \sin^{2} θ {\overset{\cdot}{ψ}}^{2} + (I_{3} \sin θ + (I_{1} - I_{2}) \sin θ \cos 2 ϕ) \overset{\cdot}{θ} \overset{\cdot}{ψ} + \frac{1}{2} (I_{2} - I_{1}) \sin 2 ϕ {\overset{\cdot}{θ}}^{2} + M_{3}

Proof: In this case, where no constraints are present, we may take $g = 0_{1 \times 3}$ , and then, according to Proposition 5.4

{\begin{matrix} M \ddot{q} - Q^{resT} = 0_{n \times 1} \\ q (0) = q_{0}, \overset{\cdot}{q} (0) = {\overset{\cdot}{q}}_{0} \end{matrix}

For the calculation of components of $Q^{res}$ , see Appendix A.3.1. □

In order to formulate the Euler equations we need the following proposition expressing the angular velocity components in terms of the Euler angles.

Proposition 6.4. The angular velocity of a rigid part is given by $ω = i_{1} ω_{1} + i_{2} ω_{2} + i_{3} ω_{3}$ where

{\begin{matrix} ω_{1} = \overset{\cdot}{ψ} \sin θ \sin ϕ + \overset{\cdot}{θ} \cos ϕ \\ ω_{2} = \overset{\cdot}{ψ} \sin θ \cos ϕ - \overset{\cdot}{θ} \sin ϕ \\ ω_{3} = \overset{\cdot}{ψ} \cos θ + \overset{\cdot}{ϕ} \end{matrix}

(43)

For the proof of this we need the following lemma.

Lemma 6.1. We have

ω_{, ψ} = i_{1} \sin θ \sin ϕ + i_{2} \sin θ \cos ϕ + i_{3} \cos θ, ω_{, θ} = i_{1} \cos ϕ - i_{2} \sin ϕ, ω_{, ϕ} = i_{3}

Proof: From proposition 3.2 one obtains

ω_{, ψ} = ax (\underline{i} {[R]}_{{\underline{i}}^{o}}^{T} \frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ψ} {\underline{i}}^{T})

where, by using equation (38) and (39(a)), one obtains

{[R]}_{{\underline{i}}^{o}}^{T} \frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ψ} = (\begin{matrix} 0 & - \cos θ & \sin θ \cos ϕ \\ \cos θ & 0 & - \sin θ \sin ϕ \\ - \sin θ \cos ϕ & \sin θ \sin ϕ & 0 \end{matrix})

and consequently $ω_{, ψ} = i_{1} \sin θ \sin ϕ + i_{2} \sin θ \cos ϕ + i_{3} \cos θ$ . Expressions for $ω_{, θ}$ and $ω_{, ϕ}$ are obtained by exploiting the following expressions

{[R]}_{{\underline{i}}^{o}}^{T} \frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial θ} = (\begin{matrix} 0 & 0 & - \sin ϕ \\ \cos θ & 0 & - \cos ϕ \\ \sin ϕ & \cos ϕ & 0 \end{matrix}), {[R]}_{{\underline{i}}^{o}}^{T} \frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ϕ} = (\begin{matrix} 0 & - \overset{\cdot}{ϕ} & 0 \\ \overset{\cdot}{ϕ} & 0 & 0 \\ 0 & 0 & 0 \end{matrix})

where use has been made of equations (39(b)) and (39(c)). This proves the lemma. □

Proof (of Proposition 6.4): We have, according to equation (9(b)), and the previous lemma

ω = ω_{, ψ} \overset{\cdot}{ψ} + ω_{, θ} \overset{\cdot}{θ} + ω_{, ϕ} \overset{\cdot}{ϕ} = i_{1} (\overset{\cdot}{ψ} \sin θ \sin ϕ + \overset{\cdot}{θ} \cos ϕ) + i_{2} (\overset{\cdot}{ψ} \sin θ \cos ϕ - \overset{\cdot}{θ} \sin ϕ) + i_{3} (\overset{\cdot}{ψ} \cos θ + \overset{\cdot}{ϕ})

□

The Euler equations for the rigid part are then given by equations (36) and (43). One then obtains the following result.

Proposition 6.5. The Euler equations for the rigid part may be written $A \ddot{q} = F (q, \overset{\cdot}{q})$ where $q = {(\begin{matrix} ψ & θ & ϕ \end{matrix})}^{T}$ and

\begin{array}{l} A = (\begin{matrix} I_{1} \sin θ \sin ϕ & I_{1} \cos ϕ & 0 \\ I_{2} \sin θ \cos ϕ & - I_{2} \sin ϕ & 0 \\ I_{3} \cos θ & 0 & I_{3} \end{matrix}) \\ F (q, \dot{q}) = (\begin{matrix} M_{1} - I_{1} (\dot{ψ} \dot{θ} \cos θ \sin ϕ + \dot{ψ} \dot{ϕ} \sin θ \cos ϕ - \dot{θ} \dot{ϕ} \sin ϕ) - (I_{3} - I_{2}) (\dot{ψ} \sin θ \cos ϕ - \dot{θ} \sin ϕ) (\dot{ψ} \cos θ + \dot{ϕ}) \\ M_{2} - I_{2} (\dot{ψ} \dot{θ} \cos θ \cos ϕ - \dot{ψ} \dot{ϕ} \sin θ \sin ϕ - \dot{θ} \dot{ϕ} \cos ϕ) - (I_{1} - I_{3}) (\dot{ψ} \cos θ + \dot{ϕ}) (\dot{ψ} \sin θ \sin ϕ + \dot{θ} \cos ϕ) \\ M_{3} + I_{3} \dot{ψ} \dot{θ} \sin θ - (I_{2} - I_{1}) (\dot{ψ} \sin θ \sin ϕ + \dot{θ} \cos ϕ) (\dot{ψ} \sin θ \cos ϕ - \dot{θ} \sin ϕ) \end{matrix}) \end{array}

Remark 6.2: Since

det A = - I_{1} I_{2} I_{3} \sin θ \neq 0, (ψ, θ, ϕ) \in Ω_{r}

the matrix $A$ is non-singular. It is however not symmetric. ▪

Remark 6.3: Obviously, both $A = A (ψ, θ, ϕ)$ and $F = F (ψ, θ, ϕ)$ are highly non-linear functions due to the presence of trigonometric functions. ▪

Remark 6.4: The rotational part of the kinetic energy is, according to equation (43), given by

\begin{array}{l} T = \frac{1}{2} ω \cdot I_{c} ω = \frac{1}{2} (ω_{1}^{2} I_{1} + ω_{2}^{2} I_{2} + ω_{2}^{2} I_{2}) \\ = \frac{1}{2} ({(\dot{ψ} \sin θ \sin ϕ + \dot{θ} \cos ϕ)}^{2} I_{1} + {(\dot{ψ} \sin θ \cos ϕ - \dot{θ} \sin ϕ)}^{2} I_{2} + {(\dot{ψ} \cos θ + \dot{ϕ})}^{2} I_{3}) \end{array}

and the generalized forces are given by

Q_{ψ} = ω_{, ψ} \cdot M_{c} = M_{1} \sin θ \sin ϕ + M_{2} \sin θ \cos ϕ + M_{3} \cos θ

Q_{θ} = ω_{, θ} \cdot M_{c} = M_{1} \cos ϕ - M_{2} \sin ϕ, Q_{ϕ} = ω_{, ϕ} \cdot M_{3} = M_{3}

Lagrange’s equations may then be written

{\begin{matrix} \frac{d}{dt} (\frac{\partial T}{\partial \overset{\cdot}{ψ}}) - \frac{\partial T}{\partial ψ} = M_{1} \sin θ \sin ϕ + M_{2} \sin θ \cos ϕ + M_{3} \cos θ \\ \frac{d}{dt} (\frac{\partial T}{\partial \overset{\cdot}{θ}}) - \frac{\partial T}{\partial θ} = M_{1} \cos ϕ - M_{2} \sin ϕ \\ \frac{d}{dt} (\frac{\partial T}{\partial \overset{\cdot}{ϕ}}) - \frac{\partial T}{\partial ϕ} = M_{3} \end{matrix}

which is, of course, equivalent to equation (41). ▪

7. Quaternions

In this case $ν = 4$ and $q_{r} \in Ω_{r} = ℝ^{4}$ . If $q_{r}^{1} = ε_{0}$ , $q_{r}^{2} = ε_{1}$ , $q_{r}^{3} = ε_{2}$ , $q_{r}^{4} = ε_{3}$ and the transplacement of the body is given by

x = χ (q; X) = x_{c} (q^{1}, q^{2}, q^{3}) + R (ε_{0}, ε_{1}, ε_{2}, ε_{3}) (X - X_{c}), X \in B_{0},

(44)

where

{[R]}_{{\underline{i}}^{o}} = (\begin{matrix} 1 - 2 (ε_{2}^{2} + ε_{3}^{2}) & 2 (ε_{1} ε_{2} - ε_{0} ε_{3}) & 2 (ε_{1} ε_{3} + ε_{0} ε_{2}) \\ 2 (ε_{1} ε_{2} + ε_{0} ε_{3}) & 1 - 2 (ε_{1}^{2} + ε_{3}^{2}) & 2 (ε_{2} ε_{3} - ε_{0} ε_{1}) \\ 2 (ε_{1} ε_{3} - ε_{0} ε_{2}) & 2 (ε_{2} ε_{3} + ε_{0} ε_{1}) & 1 - 2 (ε_{1}^{2} + ε_{2}^{2}) \end{matrix})

(45)

then the coordinates, $(ε_{0}, ε_{1}, ε_{2}, ε_{3})$ are called quaternions. The matrix in equation (45) is introduced in Appendix A.1 where it is shown that ${[R]}_{{\underline{i}}^{o}}$ is orthonormal if and only if the constraint condition

C (ε_{0}, ε_{1}, ε_{2}, ε_{3}) = ε_{0}^{2} + ε_{1}^{2} + ε_{2}^{2} + ε_{3}^{2} - 1 = 0

(46)

is satisfied. Thus, in order for equations (44) and (45) to represent a rigid transplacement condition, equation (46) has to be satisfied. Note that $det {[R]}_{{\underline{i}}^{o}} = {(ε_{0}^{2} + ε_{1}^{2} + ε_{2}^{2} + ε_{3}^{2})}^{3}$ .

Remark 7.1: If ${[R]}_{{\underline{i}}^{o}} = (R_{ij})$ , then

ε_{0}^{2} = \frac{R_{11} + R_{22} + R_{33} + 1}{4}, ε_{i}^{2} = \frac{1 + 2 R_{ii} - R_{11} - R_{22} - R_{33}}{4}, i = 1, 2, 3

and the coordinate representation using quaternions is thus free of singularities. This is obtained at the expense of a double cover of $SO (V)$ ; see remark A.3. ▪

Proposition 7.1. The quaternion coordinate system is regular.

Proof: By a direct computation one finds

\frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ε_{0}} = (\begin{matrix} 0 & - 2 ε_{3} & 2 ε_{2} \\ 2 ε_{3} & 0 & - 2 ε_{1} \\ - 2 ε_{2} & 2 ε_{1} & 0 \end{matrix}), \frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ε_{1}} = (\begin{matrix} 0 & 2 ε_{2} & 2 ε_{3} \\ 2 ε_{2} & - 4 ε_{1} & - 2 ε_{0} \\ 2 ε_{3} & 2 ε_{0} & - 4 ε_{1} \end{matrix})

\frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ε_{2}} = (\begin{matrix} - 4 ε_{2} & 2 ε_{1} & 2 ε_{0} \\ 2 ε_{1} & 0 & 2 ε_{3} \\ - 2 ε_{0} & 2 ε_{3} & - 4 ε_{2} \end{matrix}), \frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ε_{3}} = (\begin{matrix} - 4 ε_{3} & - 2 ε_{0} & 2 ε_{1} \\ 2 ε_{0} & - 4 ε_{3} & 2 ε_{2} \\ 2 ε_{1} & 2 ε_{2} & 0 \end{matrix})

and it follows that if

\begin{matrix} w_{0} \frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ε_{0}} + w_{1} \frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ε_{1}} + w_{2} \frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ε_{2}} + w_{3} \frac{\partial {[R]}_{{\underline{i}}^{o}}}{\partial ε_{3}} = w_{0} (\begin{matrix} 0 & - 2 ε_{3} & 2 ε_{2} \\ 2 ε_{3} & 0 & - 2 ε_{1} \\ - 2 ε_{2} & 2 ε_{1} & 0 \end{matrix}) \\ + w_{1} (\begin{matrix} 0 & 2 ε_{2} & 2 ε_{3} \\ 2 ε_{2} & - 4 ε_{1} & - 2 ε_{0} \\ 2 ε_{3} & 2 ε_{0} & - 4 ε_{1} \end{matrix}) + w_{2} (\begin{matrix} - 4 ε_{2} & 2 ε_{1} & 2 ε_{0} \\ 2 ε_{1} & 0 & 2 ε_{3} \\ - 2 ε_{0} & 2 ε_{3} & - 4 ε_{2} \end{matrix}) + w_{3} (\begin{matrix} - 4 ε_{3} & - 2 ε_{0} & 2 ε_{1} \\ 2 ε_{0} & - 4 ε_{3} & 2 ε_{2} \\ 2 ε_{1} & 2 ε_{2} & 0 \end{matrix}) \\ = 2 (\begin{matrix} - 2 w_{2} ε_{2} - 2 w_{3} ε_{3} & - w_{0} ε_{3} + w_{1} ε_{2} + w_{2} ε_{1} - w_{3} ε_{0} & w_{0} ε_{2} + w_{1} ε_{3} + w_{2} ε_{0} + w_{3} ε_{1} \\ w_{0} ε_{3} + w_{1} ε_{2} + w_{2} ε_{1} + w_{3} ε_{0} & - 2 w_{1} ε_{1} - 2 w_{3} ε_{3} & - w_{0} ε_{1} - w_{1} ε_{0} + w_{2} ε_{3} + w_{3} ε_{2} \\ w_{0} ε_{2} + w_{1} ε_{3} - w_{2} ε_{0} + w_{3} ε_{1} & w_{0} ε_{1} + w_{1} ε_{0} + w_{2} ε_{3} + w_{3} ε_{2} & - 2 w_{1} ε_{1} - 2 w_{2} ε_{2} \end{matrix}) = 0_{3 \times 3} \end{matrix}

for all $\forall ε_{i} \in ℝ$ , then, for instance, $w_{0} ε_{2} + w_{1} ε_{3} + w_{2} ε_{0} + w_{3} ε_{1} = 0, \forall ε_{i} \in ℝ \Rightarrow$ $w_{i} = 0, i = 0, 1, 2, 3$ and thus the quaternion coordinate system is regular. □

Proposition 7.2. The mass matrix, related to the quaternion coordinate system, is given by

M = 4 (\begin{matrix} I_{1} ε_{1}^{2} + I_{2} ε_{2}^{2} + I_{3} ε_{3}^{2} & I_{1} ε_{0} ε_{1} + (I_{2} - I_{3}) ε_{2} ε_{3} & I_{2} ε_{0} ε_{2} + (I_{3} - I_{1}) ε_{1} ε_{3} & I_{3} ε_{0} ε_{3} + (I_{1} - I_{2}) ε_{1} ε_{2} \\ I_{1} ε_{0} ε_{1} + (I_{2} - I_{3}) ε_{2} ε_{3} & I_{1} ε_{0}^{2} + I_{3} ε_{2}^{2} + 4 I_{1} ε_{1}^{2} + I_{2} ε_{3}^{2} & (2 (I_{1} + I_{2}) - I_{3}) ε_{1} ε_{2} + (I_{1} - I_{2}) ε_{0} ε_{3} & (2 (I_{1} + I_{3}) - I_{2}) ε_{1} ε_{3} + (I_{3} - I_{1}) ε_{0} ε_{2} \\ I_{2} ε_{0} ε_{2} + (I_{3} - I_{1}) ε_{1} ε_{3} & (2 (I_{1} + I_{2}) - I_{3}) ε_{1} ε_{2} + (I_{1} - I_{2}) ε_{0} ε_{3} & I_{2} ε_{0}^{2} + I_{3} ε_{1}^{2} + 4 I_{2} ε_{2}^{2} + I_{1} ε_{3}^{2} & (2 (I_{2} + I_{3}) - I_{1}) ε_{2} ε_{3} + (I_{2} - I_{3}) ε_{0} ε_{1} \\ I_{3} ε_{0} ε_{3} + (I_{1} - I_{2}) ε_{1} ε_{2} & (2 (I_{1} + I_{3}) - I_{2}) ε_{1} ε_{3} + (I_{3} - I_{1}) ε_{0} ε_{2} & (2 (I_{2} + I_{3}) - I_{1}) ε_{2} ε_{3} + (I_{2} - I_{3}) ε_{0} ε_{1} & I_{3} ε_{0}^{2} + I_{2} ε_{1}^{2} + I_{1} ε_{2}^{2} + 4 I_{3} ε_{3}^{2} \end{matrix})

(47)

where $M$ is positive definite.

Proof: See Appendix A.2.2 for the computation of $M$ . Since the quaternion coordinate system is regular, $M$ is positive definite. □

Theorem 7.1. The rotational motion of the rigid body satisfies

{\begin{matrix} M \underset{~}{\ddot{ε}} = (I_{4 \times 4} - \frac{1}{{\underset{~}{ε}}^{T} M^{- 1} \underset{~}{ε}} \underset{~}{ε} {\underset{~}{ε}}^{T} M^{- 1}) Q^{resT} - \frac{{ε_{~}^{.}}^{T} ε_{~}^{.}}{{\underset{~}{ε}}^{T} M^{- 1} \underset{~}{ε}} \underset{~}{ε} = 0_{ν \times 1} \\ \underset{~}{ε} (0) = {\underset{~}{ε}}_{0}, \underset{~}{\overset{\cdot}{ε}} (0) = {\underset{~}{\overset{\cdot}{ε}}}_{0} \\ {\underset{~}{ε}}_{0} {\underset{~}{\overset{\cdot}{ε}}}_{0} = 0 \end{matrix},

where $\underset{~}{ε} = {(\begin{matrix} ε_{0} & ε_{1} & ε_{2} & ε_{3} \end{matrix})}^{T} \in ℝ^{4 \times 1}$ and $M$ is given by equation (47) and the components of $Q^{res}$ are given by

Q_{k}^{r e s} = - \sum_{i, j = 0}^{3} (\frac{\partial m_{j k}}{\partial ε_{i}} - \frac{1}{2} \frac{\partial m_{i j}}{\partial ε_{k}}) {\dot{ε}}_{i} {\dot{ε}}_{j} + 2 ε_{k} t r A + 2 ε_{k} \cdot M_{c}, k = 0, \dots, 3,

where

ε_{0} = i^{o} (- ε_{1}) + i_{2}^{o} (- ε_{2}) + i_{3}^{o} (- ε_{3}), ε_{1} = i_{1}^{o} ε_{0} + i_{2}^{o} ε_{3} + i_{3}^{o} (- ε_{2})

(48)

ε_{2} = i_{1}^{o} (- ε_{3}) + i_{2}^{o} ε_{0} + i_{3}^{o} ε_{1}, ε_{3} = i_{1}^{o} ε_{2} + i_{2}^{o} (- ε_{1}) + i_{3}^{o} ε_{0}

Proof: In this case, $q_{r} = \underset{~}{ε}$ , $g = \overset{T}{\underset{~}{ε}}$ and $Γ = \overset{T}{\underset{~}{ε}} M^{- 1} \underset{~}{ε} \in ℝ^{+}$ . Then

g^{T} Γ^{- 1} g M^{- 1} = \frac{\underset{~}{ε} {\underset{~}{ε}}^{T}}{{\underset{~}{ε}}^{T} M^{- 1} \underset{~}{ε}} M^{- 1}, \frac{\partial g}{\partial q_{r}} = I_{4 \times 4}

and

g^{T} Γ^{- 1} (\frac{\partial g}{\partial q_{r}} {\dot{q}}_{r}) {\dot{q}}_{r} = \frac{{\underset{~}{\dot{ε}}}^{T} \underset{~}{\dot{ε}}}{\underset{~}{ε} M^{- 1} \underset{~}{ε}} \underset{~}{ε}

From Proposition 4.1 we then obtain

M \underset{~}{\ddot{ε}} - (I_{4 \times 4} - \frac{1}{{\underset{~}{ε}}^{T} M^{- 1} \underset{~}{ε}} \underset{~}{ε} {\underset{~}{ε}}^{T} M^{- 1}) Q^{r e s T} + \frac{{\underset{~}{\dot{ε}}}^{T} \underset{~}{\dot{ε}}}{{\underset{~}{ε}}^{T} M^{- 1} \underset{~}{ε}} \underset{~}{ε} = 0_{ν \times 1}

where

Q^{res} = - \underset{~}{{\overset{\cdot}{ε}}^{T}} (\frac{\partial M}{\partial \underset{~}{ε}} - \frac{1}{2} {(\frac{\partial M}{\partial \underset{~}{ε}})}^{t}) \underset{~}{\overset{\cdot}{ε}} + Q

is the residual generalized force and the components of the external generalized force $Q$ are given by $Q_{k} = 2 ε_{k} t r A + 2 ε_{k} \cdot M_{c}, k = 0, 1, 2, 3$ . See Appendix A.3.2. □

Remark 7.2: The explicit calculation of the components

{(- {\underset{~}{\overset{\cdot}{ε}}}^{T} (\frac{\partial M}{\partial \underset{~}{ε}} - \frac{1}{2} {(\frac{\partial M}{\partial \underset{~}{ε}})}^{t}) \underset{~}{\overset{\cdot}{ε}})}_{k} = - \sum_{i, j = 0}^{3} (\frac{\partial m_{jk}}{\partial ε_{i}} - \frac{1}{2} \frac{\partial m_{ij}}{\partial ε_{k}}) {\overset{\cdot}{ε}}_{i} {\overset{\cdot}{ε}}_{j}, k = 0, \dots, 3

is straightforward, but somewhat lengthy, and is not undertaken here. ▪

Remark 7.3: Note that

\underset{~}{ε} {\underset{~}{ε}}^{T} = (\begin{matrix} ε_{0}^{2} & ε_{0} ε_{1} & ε_{0} ε_{2} & ε_{0} ε_{3} \\ ε_{1} ε_{0} & ε_{1}^{2} & ε_{1} ε_{2} & ε_{1} ε_{3} \\ ε_{2} ε_{0} & ε_{2} ε_{1} & ε_{2}^{2} & ε_{2} ε_{3} \\ ε_{3} ε_{0} & ε_{3} ε_{1} & ε_{3} ε_{2} & ε_{3}^{2} \end{matrix}) .

▪

Remark 7.4: An alternative approach to the presentation above is given in Udwadia and Schutte [15] where quaternion algebra is further exploited by introducing a ‘four-vector’ representing the torque and a ‘ $4 \times 4$ -matrix’ representing the inertia tensor. This seems, however to involve some arbitrariness and leads to unnecessary complications. ▪

Remark 7.5: The Lagrangian multiplier reads

λ = - \frac{1}{{\underset{~}{ε}}^{T} M^{- 1} \underset{~}{ε}} ({\underset{~}{ε}}^{T} M^{- 1} Q^{resT} + \underset{~}{{\overset{\cdot}{ε}}^{T}} \underset{~}{\overset{\cdot}{ε}})

There does not seem to be any natural physical interpretation of this multiplier in terms of constraint forces. ▪

Let $R ~ n φ$ and associate to this rotation tensor the unit quaternion $q_{n} = \cos \frac{φ}{2} + n \sin \frac{φ}{2}$ .

Then, according to Proposition A.1

R u = q_{n} u {q_{n}}^{- 1}, u \in V

(49)

The definitions and properties of quaternions are summarized in Appendix A.1. We have the following well-known result, where the proof given here seems to be new.

Proposition 7.3. If $q_{n} = \cos \frac{φ}{2} + n \sin \frac{φ}{2}$ then

{\overset{\cdot}{q}}_{n} = \frac{1}{2} ω q_{n}

(50)

where the product on the right-hand side is quaternionic.

Proof: If $q_{n} = \cos \frac{φ}{2} + n \sin \frac{φ}{2}$ , then

{\overset{\cdot}{q}}_{n} = - \frac{\overset{\cdot}{φ}}{2} \sin \frac{φ}{2} + \overset{\cdot}{n} \sin \frac{φ}{2} + \frac{\overset{\cdot}{φ}}{2} n \cos \frac{φ}{2}

(51)

and

ω q_{n} = - ω \cdot n \sin \frac{φ}{2} + ω \cos \frac{φ}{2} + ω \times n \sin \frac{φ}{2}

(52)

From Proposition 3.3 we have $ω = Φ \overset{\cdot}{n} + n \overset{\cdot}{φ}$ where $n \cdot Φ \overset{\cdot}{n} = 0$ and then

ω \cdot n = (Φ \overset{\cdot}{n} + n \overset{\cdot}{φ}) \cdot n = \overset{\cdot}{φ}

(53)

\begin{matrix} ω \cos \frac{φ}{2} = (Φ \overset{\cdot}{n} + n \overset{\cdot}{φ}) \cos \frac{φ}{2} = \cos \frac{φ}{2} \sin φ \overset{\cdot}{n} + \cos \frac{φ}{2} (1 - \cos φ) n \times \overset{\cdot}{n} + \overset{\cdot}{φ} \cos \frac{φ}{2} n \\ = 2 \cos^{2} \frac{φ}{2} \sin \frac{φ}{2} \overset{\cdot}{n} + 2 \cos \frac{φ}{2} \sin \frac{φ}{2} n \times \overset{\cdot}{n} + \overset{\cdot}{φ} \cos \frac{φ}{2} n \end{matrix}

(54)

Furthermore

\begin{array}{l} ω \times n = (Φ \dot{n} + n \dot{φ}) \times n = - n \times Φ \dot{n} = - n \times (\sin φ 1 + (1 - \cos φ) n \times 1) \dot{n} = \\ - (\sin φ n \times \dot{n} + (1 - \cos φ) n \times (n \times \dot{n})) = - (\sin φ n \times \dot{n} - (1 - \cos φ) \dot{n}) \end{array}

and then

ω \times n \sin \frac{φ}{2} = - \sin \frac{φ}{2} \sin φ n \times \overset{\cdot}{n} + \sin \frac{φ}{2} (1 - \cos φ) \overset{\cdot}{n} = - 2 \sin^{2} \frac{φ}{2} \cos \frac{φ}{2} n \times \overset{\cdot}{n} + 2 \sin^{3} \frac{φ}{2} \overset{\cdot}{n}

(55)

By combining equations (51)–(55) one obtains

\begin{matrix} ω q_{n} = - \overset{\cdot}{φ} \sin \frac{φ}{2} + \cos \frac{φ}{2} \sin φ \overset{\cdot}{n} + \cos \frac{φ}{2} (1 - \cos φ) n \times \overset{\cdot}{n} + \cos \frac{φ}{2} n \overset{\cdot}{φ} \\ - \sin \frac{φ}{2} \sin φ n \times \overset{\cdot}{n} + \sin \frac{φ}{2} (1 - \cos φ) \overset{\cdot}{n} = - \overset{\cdot}{φ} \sin \frac{φ}{2} + 2 \cos^{2} \frac{φ}{2} \sin \frac{φ}{2} \overset{\cdot}{n} + \overset{\cdot}{φ} \cos \frac{φ}{2} n + 2 \sin^{3} \frac{φ}{2} \overset{\cdot}{n} \\ = - \overset{\cdot}{φ} \sin \frac{φ}{2} + 2 \sin \frac{φ}{2} \overset{\cdot}{n} + \cos \frac{φ}{2} n \overset{\cdot}{φ} \end{matrix}

Thus

\frac{1}{2} ω q_{n} = - \frac{1}{2} \overset{\cdot}{φ} \sin \frac{φ}{2} + \sin \frac{φ}{2} \overset{\cdot}{n} + \frac{1}{2} \cos \frac{φ}{2} n \overset{\cdot}{φ}

and by comparing this with equation (51) the proposition is proved. □

A pull-back of the angular velocity to the reference placement gives

Ω = R^{T} ω = i_{1}^{o} ω_{1} + i_{2}^{o} ω_{2} + i_{3}^{o} ω_{3}

where $\underline{i}^{0} = (i_{1}^{0} i_{2}^{0} i_{3}^{0}) = R^{T} \underline{i}$ is a fixed (time-independent) basis in the reference. Using quaternions we have

Ω = R^{T} ω = q_{n}^{- 1} ω q_{n} = 2 q_{n}^{- 1} {\overset{\cdot}{q}}_{n}

With $q_{n} = ε_{0} + i_{1}^{o} ε_{1} + i_{2}^{o} ε_{2} + i_{3}^{o} ε_{3}$ then

q_{n}^{- 1} = ε_{0} - i_{1}^{o} ε_{1} - i_{2}^{o} ε_{2} - i_{3}^{o} ε_{3}, {\overset{\cdot}{q}}_{n} = {\overset{\cdot}{ε}}_{0} + i_{1}^{o} {\overset{\cdot}{ε}}_{1} + i_{2}^{o} {\overset{\cdot}{ε}}_{2} + i_{3}^{o} {\overset{\cdot}{ε}}_{3}

and consequently (see appendix A.1)

\begin{matrix} Ω = 2 q_{n}^{- 1} {\overset{\cdot}{q}}_{n} = 2 (1 i_{1}^{o} i_{2}^{o} i_{3}^{o} \end{matrix}) (\begin{matrix} ε_{0} & ε_{1} & ε_{2} & ε_{3} \\ - ε_{1} & ε_{0} & ε_{3} & - ε_{2} \\ - ε_{2} & - ε_{3} & ε_{0} & ε_{1} \\ - ε_{3} & ε_{2} & - ε_{1} & ε_{0} \end{matrix}) (\begin{matrix} {\overset{\cdot}{ε}}_{0} \\ {\overset{\cdot}{ε}}_{1} \\ {\overset{\cdot}{ε}}_{2} \\ {\overset{\cdot}{ε}}_{3} \end{matrix}) = i_{1}^{o} 2 (- ε_{1} {\overset{\cdot}{ε}}_{0} + ε_{0} {\overset{\cdot}{ε}}_{1} + ε_{3} {\overset{\cdot}{ε}}_{2} - ε_{2} {\overset{\cdot}{ε}}_{3}) + i_{2}^{o} 2 (- ε_{2} {\overset{\cdot}{ε}}_{0} - ε_{3} {\overset{\cdot}{ε}}_{1} + ε_{0} {\overset{\cdot}{ε}}_{2} + ε_{1} {\overset{\cdot}{ε}}_{3}) + i_{3}^{o} 2 (- ε_{3} {\overset{\cdot}{ε}}_{0} + ε_{2} {\overset{\cdot}{ε}}_{1} - ε_{1} {\overset{\cdot}{ε}}_{2} + ε_{0} {\overset{\cdot}{ε}}_{3})

From this we may draw the following conclusion.

Proposition 7.4. The components of the angular velocity relative to the basis $\underline{i}$ are given by

{\begin{matrix} ω_{1} = 2 (- ε_{1} {\overset{\cdot}{ε}}_{0} + ε_{0} {\overset{\cdot}{ε}}_{1} + ε_{3} {\overset{\cdot}{ε}}_{2} - ε_{2} {\overset{\cdot}{ε}}_{3}) \\ ω_{2} = 2 (- ε_{2} {\overset{\cdot}{ε}}_{0} - ε_{3} {\overset{\cdot}{ε}}_{1} + ε_{0} {\overset{\cdot}{ε}}_{2} + ε_{1} {\overset{\cdot}{ε}}_{3}) \\ ω_{3} = 2 (- ε_{3} {\overset{\cdot}{ε}}_{0} + ε_{2} {\overset{\cdot}{ε}}_{1} - ε_{1} {\overset{\cdot}{ε}}_{2} + ε_{0} {\overset{\cdot}{ε}}_{3}) \end{matrix} .

Furthermore, we have the following proposition.

Proposition 7.5. The components of the angular velocity relative to the basis $\underline{i}^{0}$ are given by

{\begin{matrix} ω_{1}^{0} = 2 (- ε_{1} {\overset{\cdot}{ε}}_{0} + ε_{0} {\overset{\cdot}{ε}}_{1} - ε_{3} {\overset{\cdot}{ε}}_{2} + ε_{2} {\overset{\cdot}{ε}}_{3}) \\ ω_{2}^{0} = 2 (- ε_{2} {\overset{\cdot}{ε}}_{0} + ε_{3} {\overset{\cdot}{ε}}_{1} + ε_{0} {\overset{\cdot}{ε}}_{2} - ε_{1} {\overset{\cdot}{ε}}_{3}) \\ ω_{3}^{0} = 2 (- ε_{2} {\overset{\cdot}{ε}}_{0} + ε_{3} {\overset{\cdot}{ε}}_{1} + ε_{0} {\overset{\cdot}{ε}}_{2} - ε_{1} {\overset{\cdot}{ε}}_{3}) \end{matrix}

and then

ω_{, 0} = - 2 (i_{1}^{o} ε_{1} + i_{2}^{o} ε_{2} + i_{3}^{o} ε_{3}), ω_{, 1} = 2 (i_{1}^{o} ε_{0} + i_{2}^{o} ε_{3} + i_{3}^{o} (- ε_{2})),

(56)

ω_{, 2} = 2 (i_{1}^{o} (- ε_{3}) + i_{2}^{o} ε_{0} + i_{3}^{o} ε_{1}), ω_{, 3} = 2 (i_{1}^{o} ε_{2} + i_{2}^{o} (- ε_{1}) + i_{3}^{o} ε_{0})

Proof: We have $ω = \underline{i} {[ω]}_{\underline{i}} = R {\underline{i}}^{o} {[ω]}_{\underline{i}} = {\underline{i}}^{o} {[R]}_{{\underline{i}}^{o}} {[ω]}_{\underline{i}} \Rightarrow {[ω]}_{{\underline{i}}^{o}} = {[R]}_{{\underline{i}}^{o}} {[ω]}_{\underline{i}} =$

\begin{array}{l} (\begin{matrix} 1 - 2 (ε_{2}^{2} + ε_{3}^{2}) & 2 (ε_{1} ε_{2} - ε_{0} ε_{3}) & 2 (ε_{1} ε_{3} + ε_{0} ε_{2}) \\ 2 (ε_{1} ε_{2} + ε_{0} ε_{3}) & 1 - 2 (ε_{1}^{2} + ε_{3}^{2}) & 2 (ε_{2} ε_{3} - ε_{0} ε_{1}) \\ 2 (ε_{1} ε_{3} - ε_{0} ε_{2}) & 2 (ε_{2} ε_{3} + ε_{0} ε_{1}) & 1 - 2 (ε_{1}^{2} + ε_{2}^{2}) \end{matrix}) (\begin{matrix} 2 (- ε_{1} {\dot{ε}}_{0} + ε_{0} {\dot{ε}}_{1} + ε_{3} {\dot{ε}}_{2} - ε_{2} {\dot{ε}}_{3}) \\ 2 (- ε_{2} {\dot{ε}}_{0} - ε_{3} {\dot{ε}}_{1} + ε_{0} {\dot{ε}}_{2} + ε_{1} {\dot{ε}}_{3}) \\ 2 (- ε_{3} {\dot{ε}}_{0} + ε_{2} {\dot{ε}}_{1} - ε_{1} {\dot{ε}}_{2} + ε_{0} {\dot{ε}}_{3}) \end{matrix}) = \\ (\begin{matrix} 2 (- ε_{1} {\dot{ε}}_{0} + ε_{0} {\dot{ε}}_{1} - ε_{3} {\dot{ε}}_{2} + ε_{2} {\dot{ε}}_{3}) \\ 2 (- ε_{2} {\dot{ε}}_{0} + ε_{3} {\dot{ε}}_{1} + ε_{0} {\dot{ε}}_{2} - ε_{1} {\dot{ε}}_{3}) \\ 2 (- ε_{2} {\dot{ε}}_{0} + ε_{3} {\dot{ε}}_{1} + ε_{0} {\dot{ε}}_{2} - ε_{1} {\dot{ε}}_{3}) \end{matrix}) \end{array}

Thus

\begin{matrix} ω & = i_{1}^{o} 2 (- ε_{1} {\overset{\cdot}{ε}}_{0} + ε_{0} {\overset{\cdot}{ε}}_{1} - ε_{3} {\overset{\cdot}{ε}}_{2} + ε_{2} {\overset{\cdot}{ε}}_{3}) + i_{2}^{o} 2 (- ε_{2} {\overset{\cdot}{ε}}_{0} + ε_{3} {\overset{\cdot}{ε}}_{1} + ε_{0} {\overset{\cdot}{ε}}_{2} - ε_{1} {\overset{\cdot}{ε}}_{3}) + i_{3}^{o} 2 (- ε_{3} {\overset{\cdot}{ε}}_{0} - ε_{2} {\overset{\cdot}{ε}}_{1} + ε_{1} {\overset{\cdot}{ε}}_{2} + ε_{0} {\overset{\cdot}{ε}}_{3}) \\ = - 2 (i_{1}^{o} ε_{1} + i_{2}^{o} ε_{2} + i_{3}^{o} ε_{3}) {\overset{\cdot}{ε}}_{0} + 2 (i_{1}^{o} ε_{0} + i_{2}^{o} ε_{3} - i_{3}^{o} ε_{2}) {\overset{\cdot}{ε}}_{1} + 2 (- i_{1}^{o} ε_{3} + i_{2}^{o} ε_{0} + i_{3}^{o} ε_{1}) {\overset{\cdot}{ε}}_{2} + 2 (i_{1}^{o} ε_{2} - i_{2}^{o} ε_{1} + i_{3}^{o} ε_{0}) {\overset{\cdot}{ε}}_{3} \end{matrix}

Then, according to Corollary 3.1

ω_{, 0} = - 2 (i_{1}^{o} ε_{1} + i_{2}^{o} ε_{2} + i_{3}^{o} ε_{3}), ω_{, 1} = 2 (i_{1}^{o} ε_{0} + i_{2}^{o} ε_{3} - i_{3}^{o} ε_{2})

ω_{, 2} = 2 (- i_{1}^{o} ε_{3} + i_{2}^{o} ε_{0} + i_{3}^{o} ε_{1}), ω_{, 3} = 2 (i_{1}^{o} ε_{2} - i_{2}^{o} ε_{1} + i_{3}^{o} ε_{0})

and this proves the proposition. □

Corollary 7.1. We have $ω_{, k} = 2 ε_{k}, k = 0, 1, 2, 3$ .

Proof: The proof follows directly by comparing equations (56) and (48). □

Remark 7.6: The components of the external generalized force in theorem 7.1 may now be written

Q_{k} = 2 ε_{k} t r A + ω_{, k} \cdot M_{c}, k = 0, 1, 2, 3

to be compared with a similar expression, in the case of Euler angles, given in proposition A.2. ▪

The Euler equations, in terms of quaternions, are then given by equations (36) and (56) along with the kinematical constraint condition

ε_{0} {\overset{\cdot}{ε}}_{0} + ε_{1} {\overset{\cdot}{ε}}_{1} + ε_{2} {\overset{\cdot}{ε}}_{2} + ε_{3} {\overset{\cdot}{ε}}_{3} = 0

(57)

Introducing ${\underset{~}{i}}^{o} = (\begin{matrix} 1 & i_{1}^{o} & i_{2}^{o} & i_{3}^{o} \end{matrix})$ and $\underset{~}{ε} = {(\begin{matrix} ε_{0} & ε_{1} & ε_{2} & ε_{3} \end{matrix})}^{T}$ and the matrices $A$ and $B$ defined by

A = A (\underset{~}{ε}) = (\begin{matrix} ε_{0} & ε_{1} & ε_{2} & ε_{3} \\ - ε_{1} & ε_{0} & - ε_{3} & ε_{2} \\ - ε_{2} & ε_{3} & ε_{0} & - ε_{1} \\ - ε_{3} & - ε_{2} & ε_{1} & ε_{0} \end{matrix}), B = B (\underset{~}{ε}) = (\begin{matrix} - ε_{1} & ε_{0} & ε_{3} & - ε_{2} \\ - ε_{2} & - ε_{3} & ε_{0} & ε_{1} \\ - ε_{3} & ε_{2} & - ε_{1} & ε_{0} \end{matrix})

we have

Ω = {\underset{~}{i}}^{o} 2 A \underset{~}{\overset{\cdot}{ε}} = {\underline{i}}^{o} 2 B \underset{~}{\overset{\cdot}{ε}}

Thus ${[ω]}_{\underline{i}} = 2 B \underset{~}{\overset{\cdot}{ε}}$ and the Euler equations may be written

{[I_{c}]}_{\underline{i}} \overset{\cdot}{{[ω]}_{\underline{i}}} + {[ω]}_{\underline{i}} \times {[I_{c}]}_{\underline{i}} {[ω]}_{\underline{i}} = I (2 \dot{B} \underset{~}{\dot{ε}} + 2 B \underset{~}{\ddot{ε}}) + 2 B \underset{~}{\dot{ε}} \times I 2 B \underset{~}{\dot{ε}} = M_{c}

where $I = {[I_{c}]}_{\underline{i}} = d i a g (\begin{matrix} I_{1} & I_{2} & I_{3} \end{matrix}) \in ℝ^{3 \times 3}$ and $M_{c} = {[M_{c}]}_{\underline{i}} = {(\begin{matrix} M_{1} & M_{2} & M_{3} \end{matrix})}^{T}$ . The Euler equations may then be written

B \underset{~}{\ddot{ε}} = \frac{1}{2} I^{- 1} M_{c} - \dot{B} \underset{~}{\dot{ε}} - 2 I^{- 1} B \underset{~}{\dot{ε}} \times I B \underset{~}{\dot{ε}}

(59)

To this we add the constraint condition $\overset{T}{\underset{~}{ε}} \underset{~}{\overset{\cdot}{ε}} = 0$ .

Theorem 7.2. The Euler equations and the constraint condition may be summarized in the following equation

A (\underset{~}{ε}) \underset{~}{\ddot{ε}} = F (\underset{~}{ε}, \underset{~}{\dot{ε}}), F (\underset{~}{ε}, \underset{~}{\dot{ε}}) = (\begin{matrix} - {| \underset{~}{\dot{ε}} |}^{2} \\ \frac{1}{2} I^{- 1} M_{c} - \dot{B} \underset{~}{\dot{ε}} - 2 I^{- 1} B \underset{~}{\dot{ε}} \times I B \underset{~}{\dot{ε}} \end{matrix})

(60)

Proof: The constraint ${\underset{~}{ε}}^{T} \underset{~}{\overset{\cdot}{ε}} = 0$ implies ${\underset{~}{\overset{\cdot}{ε}}}^{T} \underset{~}{\overset{\cdot}{ε}} + {\underset{~}{ε}}^{T} \underset{~}{\ddot{ε}} = 0$ and then ${\underset{~}{ε}}^{T} \underset{~}{\ddot{ε}} = - | \underset{~}{\overset{\cdot}{ε}} |^{2}$ . Thus, by combining this and the Euler equation in equation (59), we have, since $A = {(\begin{matrix} \underset{~}{ε} & B^{T} \end{matrix})}^{T}$ , the equation

(\begin{matrix} {\underset{~}{ε}}^{T} \underset{~}{\ddot{ε}} \\ B \underset{~}{\ddot{ε}} \end{matrix}) = (\begin{matrix} {\underset{~}{ε}}^{T} \\ B \end{matrix}) \underset{~}{\ddot{ε}} = A \underset{~}{\ddot{ε}} = (\begin{matrix} - {| ε_{~}^{\cdot} |}^{2} \\ \frac{1}{2} I^{- 1} M_{c} - \overset{\cdot}{B} \underset{~}{\overset{\cdot}{ε}} - 2 I^{- 1} B \underset{~}{\overset{\cdot}{ε}} \times I B \underset{~}{\overset{\cdot}{ε}} \end{matrix})

and this proves the theorem. □

Remark 7.7: Note that since $A$ , according to appendix A.1, is orthonormal, equation (60) may be written

\underset{~}{\ddot{ε}} = A^{T} (\begin{matrix} - {| \underset{~}{\overset{\cdot}{ε}} |}^{2} \\ \frac{1}{2} I^{- 1} M_{c} - \overset{\cdot}{B} \underset{~}{\overset{\cdot}{ε}} - 2 I^{- 1} B \underset{~}{\overset{\cdot}{ε}} \times I B \underset{~}{\overset{\cdot}{ε}} \end{matrix})

(61)

and that the right-hand side of equation (61) contains only polynomial expressions (of the fifth degree) in the configuration coordinates and their time derivatives. ▪

Remark 7.8: Approaches similar to the one presented above are given in Chou [12] and Koshlyakov [13]. ▪

Remark 7.9: We note that the rotational part of the kinetic energy is given by

T = 2 ({(- ε_{1} {\overset{\cdot}{ε}}_{0} + ε_{0} {\overset{\cdot}{ε}}_{1} + ε_{3} {\overset{\cdot}{ε}}_{2} - ε_{2} {\overset{\cdot}{ε}}_{3})}^{2} I_{1} + {(- ε_{2} {\overset{\cdot}{ε}}_{0} - ε_{3} {\overset{\cdot}{ε}}_{1} + ε_{0} {\overset{\cdot}{ε}}_{2} + ε_{1} {\overset{\cdot}{ε}}_{3})}^{2} I_{2} + {(- ε_{3} {\overset{\cdot}{ε}}_{0} + ε_{2} {\overset{\cdot}{ε}}_{1} - ε_{1} {\overset{\cdot}{ε}}_{2} + ε_{0} {\overset{\cdot}{ε}}_{3})}^{2} I_{3})

(62)

This is, if $\overset{T}{\underset{~}{ε}} \underset{~}{ε} = 1$ , identical to

T = \frac{1}{2} \underset{~}{{\overset{\cdot}{ε}}^{T}} M \underset{~}{\overset{\cdot}{ε}}

where $M$ is given by equation (47). ▪

Remark 7.10: Note that the kinetic energy may be written in the form $T = \frac{1}{2} {[ω]}_{\underline{i}}^{T} I {[ω]}_{\underline{i}} = \frac{1}{2} {(2 B \underset{~}{\overset{\cdot}{ε}})}^{T} I 2 B \underset{~}{\overset{\cdot}{ε}} = 2 \underset{~}{{\overset{\cdot}{ε}}^{T}} J \underset{~}{\overset{\cdot}{ε}}$ where $J = B^{T} I B \in ℝ^{4 \times 4}$ is a symmetric, but singular matrix. ▪

8. Linear coordinates

Linear coordinates are obtained by using the matrix elements $R_{ij}$ of the matrix ${[R]}_{{\underline{i}}^{o}}$ as generalized coordinates ( $q_{r}^{3 (i - 1) + j} = R_{ij}, i, j = 1, 2, 3$ ). In this case $ν = 9$ and $q_{r} \in Ω_{r} = ℝ^{9}$ and the coordinates $R_{ij}$ are called linear coordinates. The transplacement of the rigid part may then be written

x = χ (q; X) = x_{O} + i_{1}^{o} q^{1} + i_{2}^{o} q^{2} + i_{3}^{o} q^{3} + \sum_{i, j = 1}^{3} i_{i}^{o} i_{j}^{o} \cdot (X - X_{c}) R_{ij}, X \in B_{0}

(63)

Remark 8.1: The designation linear coordinates is referring to configuration coordinates where the transplacement may be written [18]

x = χ (q; X) = x_{O} + \sum_{k = 1}^{n} a_{k} (X) q^{k}, X \in B_{0}

and where the shape functions, $a_{k} = a_{k} (X), X \in B_{0}, k = 1, \dots, n$ , are differentiable. In this case the transplacement $χ$ is, for fixed $X \in B_{0}$ , a linear function of the configuration coordinates. ▪

Proposition 8.1. The system of coordinates, defined by equation (63) is regular.

Proof: The tensors

\frac{\partial R}{\partial R_{ij}} = i_{i}^{o} \otimes i_{j}^{o}, i, j = 1, 2, 3

are linearly independent. Thus, according to proposition 3.5, the coordinate system is regular. □

The rotation matrix, satisfies ${[R]}_{{\underline{i}}^{o}} {[R]}_{{\underline{i}}^{o}}^{T} = I_{3 \times 3}$ . This condition is equivalent to the constraint condition $C (q_{r}) = 0_{6 \times 1}$ where $C : Ω_{r} \to ℝ^{6 \times 1}$ and

C_{1} (q_{r}) = \sum_{k = 1}^{3} R_{1 k}^{2} - 1, C_{2} (q_{r}) = \sum_{k = 1}^{3} R_{2 k}^{2} - 1, C_{3} (q_{r}) = \sum_{k = 1}^{3} R_{3 k}^{2} - 1

C_{4} (q_{r}) = \sum_{k = 1}^{3} R_{1 k} R_{2 k}, C_{5} (q_{r}) = \sum_{k = 1}^{3} R_{1 k} R_{3 k}, C_{6} (q_{r}) = \sum_{k = 1}^{3} R_{2 k} R_{3 k}

The corresponding kinematical constraint reads

g {\overset{\cdot}{q}}_{r} = 0_{6 \times 1}

(64)

where the constraint matrix $g = (g_{ν k}) \in ℝ^{6 \times 9}$ is given by

g = \frac{\partial C}{\partial q_{r}} = (\begin{matrix} R_{11} & R_{12} & R_{13} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & R_{21} & R_{22} & R_{23} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & R_{31} & R_{32} & R_{33} \\ R_{21} & R_{22} & R_{23} & R_{11} & R_{12} & R_{13} & 0 & 0 & 0 \\ R_{31} & R_{32} & R_{33} & 0 & 0 & 0 & R_{11} & R_{12} & R_{13} \\ 0 & 0 & 0 & R_{31} & R_{32} & R_{33} & R_{21} & R_{22} & R_{23} \end{matrix})

(65)

The constraint matrix in equation (65) has full rank, i.e. $rank (g) = 6$ . See Lidström [18].

Proposition 8.2. The mass matrix is given by

M = (\begin{matrix} M_{3 \times 3} & 0_{3 \times 3} & 0_{3 \times 3} \\ 0_{3 \times 3} & M_{3 \times 3} & 0_{3 \times 3} \\ 0_{3 \times 3} & 0_{3 \times 3} & M_{3 \times 3} \end{matrix}) \in ℝ^{9 \times 9}

(66)

where

M_{3 \times 3} = \frac{1}{2} (\begin{matrix} I_{2} + I_{3} - I_{1} & 0 & 0 \\ 0 & I_{3} + I_{1} - I_{2} & 0 \\ 0 & 0 & I_{3} + I_{1} - I_{2} \end{matrix})

(67)

$M$ is positive definite.

Proof: See Appendix A.2.3 for the computation of $M$ . Since the quaternion coordinate system is regular, $M$ is positive definite. □

Remark 8.2: Note that the mass matrix is constant. The positive definiteness of $M$ is related to the inequalities

I_{1} + I_{2} > I_{3}, I_{3} + I_{1} > I_{2}, I_{2} + I_{3} > I_{1}

We have

det M = {(det M_{3 \times 3})}^{3} = \frac{1}{512} {(I_{2} + I_{3} - I_{1})}^{3} {(I_{3} + I_{1} - I_{2})}^{3} {(I_{3} + I_{1} - I_{2})}^{3} > 0

▪

Remark 8.3: The rotational part of the kinetic energy may be written

T = \frac{1}{2} {\overset{\cdot}{q}}_{r} M {\overset{\cdot}{q}}_{r}^{T} = \frac{1}{2} \sum_{i, j = 1}^{3} {\overset{\cdot}{R}}_{ij}^{2} E_{j}

▪

Theorem 8.1. The rotational motion of the rigid body satisfies

{\begin{matrix} M {\ddot{q}}_{r} = Q^{T} - g^{T} Γ^{- 1} (\frac{\partial g}{\partial q_{r}} {\overset{\cdot}{q}}_{r}) {\overset{\cdot}{q}}_{r} \\ q_{r} (0) = q_{r, 0}, {\overset{\cdot}{q}}_{r} (0) = {\overset{\cdot}{q}}_{r, 0} \\ g (q_{r, 0}) {\overset{\cdot}{q}}_{r, 0} = 0_{6 \times 1} \end{matrix}

(68)

where $M$ is given by equation (66), $Γ = g M^{- 1} g^{T} \in ℝ^{6 \times 6}$ and $(\frac{\partial g}{\partial q_{r}} {\overset{\cdot}{q}}_{r}) {\overset{\cdot}{q}}_{r} \in ℝ^{6 \times 1}$ with components

{((\frac{\partial g}{\partial q_{r}} {\dot{q}}_{r}) {\dot{q}}_{r})}_{i} = \sum_{j = 1}^{9} \sum_{k = 1}^{9} \frac{\partial g_{i j}}{\partial q_{r}^{k}} {\dot{q}}_{r}^{j} {\dot{q}}_{r}^{k}, i = 1, …, 6

The components of the external generalized force $Q$ are given by $Q_{k} = i_{j} \cdot A i_{i}^{o}$ , $k = 3 (i - 1) + j, i, j = 1, 2, 3$ .

Proof: since the mass matrix is constant it follows that

Q^{res} = - {\overset{\cdot}{q}}_{r}^{T} (\frac{\partial M}{\partial q_{r}} - \frac{1}{2} {(\frac{\partial M}{\partial q_{r}})}^{t}) {\overset{\cdot}{q}}_{r} + Q = Q

For the calculation of components of $Q$ see Appendix A.3.3. The theorem is then a direct consequence of Theorem 4.1. □

Remark 8.4: The Lagrangian multiplier reads

λ = - Γ^{- 1} (g M^{- 1} Q^{T} + (\frac{\partial g}{\partial q_{r}} {\overset{\cdot}{q}}_{r}) {\overset{\cdot}{q}}_{r}) \in ℝ^{6 \times 1}

A physical interpretation of the six multiplier components contained in $λ$ may be accomplished if we consider a second body $B'$ , consisting of four point masses, which is inertia equimoment to the rigid part $B$ (same total mass and same principal moments of inertia with respect to the centre of mass). This body $B'$ always exists but it is not unique. If the point masses are connected with rigid, mass-less rods, then this body becomes rigid and the forces in these rod will be given by the components of $λ$ . ▪

Proposition 8.3. The angular velocity is $ω = i_{1} ω_{1} + i_{2} ω_{2} + i_{3} ω_{3}$ where

{\begin{matrix} ω_{1} = R_{11} (R_{21} {\overset{\cdot}{R}}_{31} + R_{22} {\overset{\cdot}{R}}_{32} + R_{23} {\overset{\cdot}{R}}_{33}) + R_{21} ({\overset{\cdot}{R}}_{11} R_{31} + {\overset{\cdot}{R}}_{12} R_{32} + {\overset{\cdot}{R}}_{13} R_{33}) + \\ R_{31} (R_{11} {\overset{\cdot}{R}}_{21} + R_{12} {\overset{\cdot}{R}}_{22} + R_{13} {\overset{\cdot}{R}}_{23}) \\ ω_{2} = R_{12} (R_{21} {\overset{\cdot}{R}}_{31} + R_{22} {\overset{\cdot}{R}}_{32} + R_{23} {\overset{\cdot}{R}}_{33}) + R_{22} ({\overset{\cdot}{R}}_{11} R_{31} + {\overset{\cdot}{R}}_{12} R_{32} + {\overset{\cdot}{R}}_{13} R_{33}) + \\ R_{32} (R_{11} {\overset{\cdot}{R}}_{21} + R_{12} {\overset{\cdot}{R}}_{22} + R_{13} {\overset{\cdot}{R}}_{23}) \\ ω_{3} = R_{13} (R_{21} {\overset{\cdot}{R}}_{31} + R_{22} {\overset{\cdot}{R}}_{32} + R_{23} {\overset{\cdot}{R}}_{33}) + R_{23} ({\overset{\cdot}{R}}_{11} R_{31} + {\overset{\cdot}{R}}_{12} R_{32} + {\overset{\cdot}{R}}_{13} R_{33}) + \\ R_{33} (R_{11} {\overset{\cdot}{R}}_{21} + R_{12} {\overset{\cdot}{R}}_{22} + R_{13} {\overset{\cdot}{R}}_{23}) \end{matrix}

(69)

Proof: We have $\overset{\cdot}{R} R^{T} = {\underline{i}}^{o} {[\overset{\cdot}{R}]}_{{\underline{i}}^{o}} {[R]}_{{\underline{i}}^{o}}^{T} {\underline{i}}^{o}^{T}$ where

{[\overset{\cdot}{R}]}_{{\underline{i}}^{o}} {[R]}_{{\underline{i}}^{o}}^{T} = (\begin{matrix} {\overset{\cdot}{R}}_{11} & {\overset{\cdot}{R}}_{12} & {\overset{\cdot}{R}}_{13} \\ {\overset{\cdot}{R}}_{21} & {\overset{\cdot}{R}}_{22} & {\overset{\cdot}{R}}_{23} \\ {\overset{\cdot}{R}}_{31} & {\overset{\cdot}{R}}_{32} & {\overset{\cdot}{R}}_{33} \end{matrix}) (\begin{matrix} R_{11} & R_{21} & R_{31} \\ R_{12} & R_{22} & R_{32} \\ R_{13} & R_{23} & R_{33} \end{matrix}) = (\begin{matrix} 0 & - (R_{11} {\overset{\cdot}{R}}_{21} + R_{12} {\overset{\cdot}{R}}_{22} + R_{13} {\overset{\cdot}{R}}_{23}) & {\overset{\cdot}{R}}_{11} R_{31} + {\overset{\cdot}{R}}_{12} R_{32} + {\overset{\cdot}{R}}_{13} R_{33} \\ - ({\overset{\cdot}{R}}_{11} R_{21} + {\overset{\cdot}{R}}_{12} R_{22} + {\overset{\cdot}{R}}_{13} R_{23}) & 0 & - (R_{21} {\overset{\cdot}{R}}_{31} + R_{22} {\overset{\cdot}{R}}_{32} + R_{23} {\overset{\cdot}{R}}_{33}) \\ R_{11} {\overset{\cdot}{R}}_{31} + R_{12} {\overset{\cdot}{R}}_{32} + R_{13} {\overset{\cdot}{R}}_{33} & - ({\overset{\cdot}{R}}_{21} R_{31} + {\overset{\cdot}{R}}_{22} R_{32} + {\overset{\cdot}{R}}_{23} R_{33}) & 0 \end{matrix})

where the constraint condition in equation (64) has been used. Thus the components of the angular velocity relative to the basis ${\underline{i}}^{o}$ are

{[ω]}_{{\underline{i}}^{o}} = (\begin{matrix} ω_{1}^{0} \\ ω_{2}^{0} \\ ω_{3}^{0} \end{matrix}) = (\begin{matrix} R_{21} {\overset{\cdot}{R}}_{31} + R_{22} {\overset{\cdot}{R}}_{32} + R_{23} {\overset{\cdot}{R}}_{33} \\ {\overset{\cdot}{R}}_{11} R_{31} + {\overset{\cdot}{R}}_{12} R_{32} + {\overset{\cdot}{R}}_{13} R_{33} \\ R_{11} {\overset{\cdot}{R}}_{21} + R_{12} {\overset{\cdot}{R}}_{22} + R_{13} {\overset{\cdot}{R}}_{23} \end{matrix}) .

However, $ω = {\underline{i}}^{o} {[ω]}_{{\underline{i}}^{o}} = \underline{i} {[ω]}_{\underline{i}} = R {\underline{i}}^{o} {[ω]}_{\underline{i}} = {\underline{i}}^{o} {[R]}_{{\underline{i}}^{o}} {[ω]}_{\underline{i}}$ . Consequently

\begin{array}{l} {[ω]}_{\underline{i}} = (\begin{matrix} ω_{1} \\ ω_{2} \\ ω_{3} \end{matrix}) = {[R]}_{{\underline{i}}^{o}}^{T} {[ω]}_{{\underline{i}}^{o}} = (\begin{matrix} R_{11} ω_{1}^{0} + R_{21} ω_{2}^{0} + R_{31} ω_{3}^{0} \\ R_{12} ω_{1}^{0} + R_{22} ω_{2}^{0} + R_{32} ω_{3}^{0} \\ R_{13} ω_{1}^{0} + R_{23} ω_{2}^{0} + R_{33} ω_{3}^{0} \end{matrix}) = \\ (\begin{matrix} R_{11} (R_{21} {\dot{R}}_{31} + R_{22} {\dot{R}}_{32} + R_{23} {\dot{R}}_{33}) + R_{21} ({\dot{R}}_{11} R_{31} + {\dot{R}}_{12} R_{32} + {\dot{R}}_{13} R_{33}) + R_{31} (R_{11} {\dot{R}}_{21} + R_{12} {\dot{R}}_{22} + R_{13} {\dot{R}}_{23}) \\ R_{12} (R_{21} {\dot{R}}_{31} + R_{22} {\dot{R}}_{32} + R_{23} {\dot{R}}_{33}) + R_{22} ({\dot{R}}_{11} R_{31} + {\dot{R}}_{12} R_{32} + {\dot{R}}_{13} R_{33}) + R_{32} (R_{11} {\dot{R}}_{21} + R_{12} {\dot{R}}_{22} + R_{13} {\dot{R}}_{23}) \\ R_{13} (R_{21} {\dot{R}}_{31} + R_{22} {\dot{R}}_{32} + R_{23} {\dot{R}}_{33}) + R_{23} ({\dot{R}}_{11} R_{31} + {\dot{R}}_{12} R_{32} + {\dot{R}}_{13} R_{33}) + R_{33} (R_{11} {\dot{R}}_{21} + R_{12} {\dot{R}}_{22} + R_{13} {\dot{R}}_{23}) \end{matrix}) \end{array}

This proves the proposition. □

The Euler equations, in terms of linear coordinates, are then given by equations (36) and (69) along with the kinematical constraint condition in equation (64). By introducing the matrix $B$ defined by

B (q_{r}) = (\begin{matrix} R_{21} R_{31} & R_{21} R_{32} & R_{21} R_{33} & R_{31} R_{11} & R_{31} R_{12} & R_{31} R_{13} & R_{11} R_{21} & R_{11} R_{22} & R_{11} R_{23} \\ R_{22} R_{31} & R_{22} R_{32} & R_{22} R_{33} & R_{32} R_{11} & R_{32} R_{12} & R_{32} R_{13} & R_{12} R_{21} & R_{12} R_{22} & R_{12} R_{23} \\ R_{23} R_{31} & R_{23} R_{32} & R_{23} R_{33} & R_{33} R_{11} & R_{33} R_{12} & R_{33} R_{13} & R_{13} R_{21} & R_{13} R_{22} & R_{13} R_{23} \end{matrix})

we have ${[ω]}_{\underline{i}} = B {\overset{\cdot}{q}}_{r}$ and the Euler equations may be written

{[I_{c}]}_{\underline{i}} \overset{\cdot}{{[ω]}_{\underline{i}}} + {[ω]}_{\underline{i}} \times {[I_{c}]}_{\underline{i}} {[ω]}_{\underline{i}} = I (\dot{B} {\dot{q}}_{r} + B {\ddot{q}}_{r}) + B {\dot{q}}_{r} \times I B {\dot{q}}_{r} = M_{c}

B {\ddot{q}}_{r} = I^{- 1} M_{c} - \dot{B} {\dot{q}}_{r} - I^{- 1} \dot{B} {\dot{q}}_{r} \times I \dot{B} {\dot{q}}_{r}

(70)

To this we add the constraint condition $g {\overset{\cdot}{q}}_{r} = 0_{6 \times 1}$ .

Theorem 8.2. The Euler equations and the constraint condition may be summarized in the equation

A (q_{r}) {\ddot{q}}_{r} = F (q_{r}, {\overset{\cdot}{q}}_{r})

A = (\begin{matrix} g \\ B \end{matrix}), F (q_{r}, {\dot{q}}_{r}) = (\begin{matrix} - {\dot{g}}_{r} {\dot{q}}_{r} \\ I^{- 1} M_{c} - \dot{B} {\dot{q}}_{r} - I^{- 1} \dot{B} {\dot{q}}_{r} \times I \dot{B} {\dot{q}}_{r} \end{matrix})

Proof: The constraint $g {\overset{\cdot}{q}}_{r} = 0_{6 \times 1}$ implies $g_{r} {\ddot{q}}_{r} = - {\overset{\cdot}{g}}_{r} {\overset{\cdot}{q}}_{r}$ . Thus, by combining this and the Euler equation in equation (67), we have

(\begin{matrix} g {\ddot{q}}_{r} \\ B {\ddot{q}}_{r} \end{matrix}) = (\begin{matrix} g \\ B \end{matrix}) {\ddot{q}}_{r} = A {\ddot{q}}_{r} = (\begin{matrix} - \overset{\cdot}{g} {\overset{\cdot}{q}}_{r} \\ I^{- 1} M_{c} - \overset{\cdot}{B} {\overset{\cdot}{q}}_{r} - I^{- 1} \overset{\cdot}{B} {\overset{\cdot}{q}}_{r} \times I \overset{\cdot}{B} {\overset{\cdot}{q}}_{r} \end{matrix})

and this proves the theorem. □

Remark 8.5. The matrix $A$ reads

A = (\begin{matrix} R_{11} & R_{12} & R_{13} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & R_{21} & R_{22} & R_{23} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & R_{31} & R_{32} & R_{33} \\ R_{21} & R_{22} & R_{23} & R_{11} & R_{12} & R_{13} & 0 & 0 & 0 \\ R_{31} & R_{32} & R_{33} & 0 & 0 & 0 & R_{11} & R_{12} & R_{13} \\ 0 & 0 & 0 & R_{31} & R_{32} & R_{33} & R_{21} & R_{22} & R_{23} \\ R_{21} R_{31} & R_{21} R_{32} & R_{21} R_{33} & R_{31} R_{11} & R_{31} R_{12} & R_{31} R_{13} & R_{11} R_{21} & R_{11} R_{22} & R_{11} R_{23} \\ R_{22} R_{31} & R_{22} R_{32} & R_{22} R_{33} & R_{32} R_{11} & R_{32} R_{12} & R_{32} R_{13} & R_{12} R_{21} & R_{12} R_{22} & R_{12} R_{23} \\ R_{23} R_{31} & R_{23} R_{32} & R_{23} R_{33} & R_{33} R_{11} & R_{33} R_{12} & R_{33} R_{13} & R_{13} R_{21} & R_{13} R_{22} & R_{13} R_{23} \end{matrix}),

with $det A = {(det R)}^{4} = 1$ .

9. Concluding remarks

Coordinate representations for rigid parts in multibody dynamics have been discussed in this paper. General configuration coordinates are introduced along with the requirement of regularity. The regularity of the configuration coordinates is equivalent to the positive definiteness of the mass matrix. This observation is often something left unnoticed in the mechanics literature. Focus is then placed on Euler angles, quaternions and linear coordinates. These coordinate systems are all demonstrated to be regular. Mass matrices corresponding to the different coordinate systems are derived and it is found that in terms of the Euler angles the mass matrix components are highly non-linear functions, since they are products of trigonometric functions. In terms of quaternions the matrix component are quadratic polynomials and using linear coordinates the mass matrix obtained is constant. Thus, there is a decreasing degree of complexity concerning mass matrix components when going from Euler angles to linear coordinates. This simplification is accomplished at the expense of an increasing gross number of degrees of freedom and the introduction of necessary constraint conditions. Equations of motion are formulated, using Lagrange’s as well as Euler’s methods. The equations are supplemented by the appropriate constraint conditions in the cases of quaternion and linear coordinates. The equations are compared with respect to their structural complexity. In the case of Euler angles the two methods give equations of similar complexity, but with a small advantage for Lagrange’s method, since it delivers $A$ as the mass matrix, which is symmetric. Using quaternions does favour Euler’s method through a simpler and more direct incorporation of the constraint condition. Furthermore the matrix $A$ is, in this case, orthonormal. Linear coordinates seem to be more convenient when using Lagrange’s equations. This is partly due to the fact that the mass matrix, in this case, is constant.

The paper also gives a new proof of the well-known relation between angular velocity and unit quaternions. In all representations the components of the angular velocity are explicitly presented. This is not always the case in previously published papers on this subject.

Footnotes

Appendix

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

References

Euler

. Theoria motus corporum solidorum seu rigidorum. Rostock und Greifswald: Röse, 1765.

Goldstein

Pool

Safco

. Classical mechanics. 3rd ed. San Francisco, California, USA: Addison Wesley, 2002.

Afzali-Far

Lidström

Nilsson

. Parametric damped vibrations of Gough–Stewart platforms for symmetric configurations. Mech Machine Theory 2014; 80: 52–69.

Hamilton

. On quaternions and the rotation of a solid body. Proceed Royal Irish Acad 1850; 4: 38–56.

Euler

. Problema algebraicum ob affectiones prorsus singulares memorabile. Novi Commentari Academiae Scientiarum Imperalis Petropolitanae 1771; 15: 75–106.

Rodrigues

. Des lois géométriques qui régissant les déplacements d’un systèm solide dans l’espace, et de la variation des coordonnées provenant de ses déplacement considérés indépendamment des causes qui peuvent les produire. J de Mathématique Pures et Appliquées 1840; 5: 330–440.

Selig

. Geometric fundamentals of robotics. New York: Springer, 2005.

Shuster

. A survey of attitude representations. J Astro Sci 1993; 41(4): 439–517.

Altmann

. Rotations, quaternions, and double groups. New York: The Clarendon Press, Oxford University Press, 1986.

10.

Lidström

. On the equations of motion in constrained multibody dynamics. Math Mech Solids 2012; 17(3): 209–242.

11.

Rapaport

. Dynamics simulations using quaternions. J Comput Physics 1985; 60: 306–314.

12.

Chou

JCK

. Quaternion kinematic and dynamic differential equations. IEEE Trans Robotics Automat 1992; 8(1); 1992.

13.

Koshlyakov

. Generalized Euler equations in quaternions. Ukrainian Math J 1994; 46(10): 1561–1564.

14.

Betsch

Siebert

. Rigid body dynamics in terms of quaternions: Hamiltonian formulation and conserving numerical integration. Int J Numer Meth Engr 2009; 79: 444–473.

15.

Udwadia

Schutte

. An alternative derivation of the quaternion equations of motion for rigid-body rotational dynamics. J Appl Mech 2010; 77: 1–4.

16.

Zoric

Lazarevic

Simonovic

. Multi-body kinematics and dynamics in terms of quaternions: Lagrange formulation in covariant form—Rodrigues approach. FME Trans 2010; 38: 19–28.

17.

Choi

Yang

. On singularity of rigid-body dynamics using quaternion-based models. J Appl Mech 2012; 79: 024502-1–024502-7.

18.

Lidström

. Linear coordinates in multibody dynamics. Math Mech Solids 2012; 17(6): 587–623.

19.

Krenk

Nielsen

. Conservative rigid body dynamics by convected base vectors with implicit constraints. Comput Methods Appl Mech Engr 2013; 269: 437–453.

20.

Geradin

Cardona

. Flexible multibody dynamics. New York: J Wiley & Sons, 2001.

21.

Lidström

. On the relative rotation of rigid parts and the visco-elastic torsion bushing element. Math Mech Solids 2013; 18(8): 788–802.

22.

Lidström

. On the equations of motion in multibody dynamics. Math Mech Solids 2012; 17(2): 165–205.

23.

Arnold

. Mathematical methods of classical mechanics. 2nd ed. New York: Springer, 1989.

24.

Jacobson

. Lectures in abstract algebra, volume I—basic concepts. Princeton, New Jersey, USA: D Van Nostrand Company, Inc., 1951.