SLIDE 1 Basilio Bona - Dynamic Modelling 71 If we use the quaternions, the Euler parameters, the Cayley-Klein parameters or the Rodrigues vectors, it in unclear what are the parameters that take the place
Using these alternative representations, we cannot say that the αi are “angles”, but it is necessary to speak in more general terms of “angular parameters”, whose knowledge nevertheless allows to compute the chosen (Euler or RPY or Cardan) angles. For instance, using the quaternions we can define: α1 = u1 sin θ 2, α2 = u2 sin θ 2, α3 = u3 sin θ 2 (2.106) and implicitly assume the unit norm. However it is a very common practice to associate to the αi the Euler or the RPY angles, although in the last few years quaternion representation has gained much attention, mainly in the aerospace and computer graphics applications.
2.13 Point Kinematics
Having introduced in the previous Sections the different representations of a geomet- rical point or of a vector, both polar or axial, in a given reference frame, and having also characterized the various types of rigid displacements in the Euclidean space, we are now ready to describe the motion of geometrical points, assumed massless. This description takes the name of kinematics and is distinct from dynamics since the former studies the motion of points or rigid bodies establishing the relations among positions, velocities and accelerations from a pure geometrical point of view, while the latter studies the influence of the forces and torques on the bodies on these quantities. In few words, kinematics study the motion without considering its causes and its effects, while dynamics studies how the external actions on the given bodies are related to their motion. Suppose to have a geometrical point P moving in the 3D Euclidean space with a time law specified by a function P(t); we can assume a vector representation of P in a given reference frame x P(t) = x1P(t) x2P(t) x2P(t)
- r, to keep notation simple
x(t) = x1(t) x2(t) x2(t)
SLIDE 2 72 Basilio Bona - Dynamic Modelling Recalling the vector derivative in A.5, we write d dtx(t) = d dtx1(t) d dtx2(t) d dtx2(t) ≡ ˙ x(t) = ˙ x1(t) ˙ x2(t) ˙ x3(t) ≡ v(t) Since the velocity is the limit dP(t) dt = lim
∆T→0
− − → ∆P ∆T the signed segment − − → ∆P can be different in different reference frames. If the point P is fixed in a local reference frame Rm, we have [− − → ∆P ]
Rm = 0; but if
Rm moves with respect to a world frame R0, then [− − → ∆P ]
R0 ̸= 0. In the first case
we speak on “local” velocity, in the second case we speak of “global” or absolute velocity or total velocity. Usually when we do not specify otherwise, we will consider always the total velocity of a point. When the vector x(t) is subject to a time-varying rotation R(t)x(t), its derivate is computed according to the normal derivative product rule, i.e., d dt [R(t)x(t)] = [ d dtR(t) ] x(t) + R(t) [ d dtx(t) ] therefore it is necessary to derive the rotation matrix R(t). To understand the consequences of this operation, we have to study its properties.
2.13.1 Rotation Matrix Derivative and Angular Velocity
To start we assume that the generic rotation matrix R is function of a generic variable x, that is R = R(x). We recall that the orthogonality of the matrix implies R(x)RT(x) = I We derive this relation with respect to a generic variable x obtaining: dR(x) dx RT(x) + R(x)dRT(x) dx = O (2.107) It is evident that the first term is the transpose of the second term, therefore, recalling the definition of skew-symmetric matrix in Appendix B.6, we can define the following skew-symmetric matrix S (u(x))
def
= dR(x) dx RT(x), (2.108)
SLIDE 3 Basilio Bona - Dynamic Modelling 73 we have introduced a generic vector u(x) because we know that a skew-symmetric matrix in R3 embeds in its structure the components of a vector. We will later discuss its precise meaning. Taking both terms of (2.107) and post-multiplying them by R(x) we obtain dR (x) dx = S (u (x)) R (x) (2.109) and it results that the derivative of an orthonormal matrix is the matrix itself pre- multiplied by a suitable anti-symmetric matrix. S(u) in (2.108) is function of a generic vector u(x), itself a function of a generic scala variable x. In particular, when the rotation matrix R depends on an angle θ(t) around an axis given by the unit vector u, we can write: dR (u, θ(t)) dt ≡ ˙ R (u, θ(t))
def
= S (ω(t)) R (u, θ(t)) (2.110) where ω(t) is the instantaneous total angular velocity vector of the reference frame represented by the matrix R(t). The angular velocity is represented in the world reference frame, but we omit to write it as ω0 for notational simplicity. Now we take into consideration some simple cases of vector transformation R(t)r(t): we start with r(t) = r constant in time in a “mobile” reference frame Rm a we transform it n the world reference frame R0 r 0(t) = R0
m(t)r m
the time derivative is computed as: ˙ r 0(t) = ˙ R
m(t)r m
(2.111) the vectors in this expression are represented in two different reference frames, so it is necessary to represent both in R0; considering that r m = (R0
m)Tr 0
we have at the end, omitting for notational convenience the dependence of R on t, ˙ r 0(t) = ˙ R
1(R0 1)Tr 0(t)
(2.112) In mechanics textbooks one often encounters the following notation Ω(t) ≡ ˙ R
1(R0 1)T
that coincide with our previous skew-symmetric matrix S(ω); indeed Ω = ˙ R
1(R0 1)T = S(ω)R0 1(R0 1)T = S(ω)
SLIDE 4 74 Basilio Bona - Dynamic Modelling This matrix is often calle the angular velocity matrix; in some textbooks the skew-symmetric matrix Ω(t) = S(ω(t) is denoted with the symbol ω(t). Considering again the identity (2.111), we observe that it expresses the well known formula that gives the linear velocity of a point P fixed in a reference frame that rotates with an angular velocity ω: ˙ r 0(t) = Ω(t)r 0(t) = S(ω(t))r 0(t) = ω(t) × r 0(t) (2.113) Other interesting properties of the skew-symmetric matrices are related to the rigid body representation in 3D space. As shown in (2.40), it is possible to compute the rotation matrix R from the rotation axis represented by the non-unit vector v and the rotation angle θ = ∥v∥, considering the skew-symmetric matrix S(v). This equality is a theoretical consequence of the following property, in general valid for Lie groups and algebras (refer to [34] for details): R(v) ≡ R(u, θ) = eS(v) ≡ eS(θ, u) =
∞
∑
k=0
1 k!S k(θu) (2.114) that relates the rotation matrix R to the exponential of the skew-symmetric matrix S(v) = S(u, θ), function of the unit axis u and angle θ. Taking the Taylor series of the matrix exponential and considering that the rotations are elements of a cyclic group, and this reflects on the strusture of S i.e., S 3 = −S, S 4 = −S 2 we obtain: R = I + sin θ θ S(v) + 1 − cos θ θ2 S 2(v) (2.115) from which one gets the identities (2.40) and (2.41). Without entering into details, that can be found in the cited textbook, we write here the fundamental identity that relates the angular velocity vector ω(t), the unit vector u(t) and the angle θ(t): ω(t) = ˙ θ(t)u(t) + sin θ(t) ˙ u(t) + (1 − cos θ(t)) S (u(t)) ˙ u(t) (2.116) From (2.116) we see that the vector ω(t) is not the formal derivative of another vector, except from the simple case that u is constant, i.e., u(t) = c. In such a case we have: ω(t) = ˙ θ(t)c = ˙ r(t) (2.117) that represents the total time derivative of the vector r(t) = θ(t)c. Another useful formula to get ω(t) from a rotation matrix R = [ r g b ] is the following (see [17]): ω(t) = 1 2(r × ˙ r + g × ˙ g + b × ˙ b) (2.118) Notice the similarity of this equality with (2.105).
SLIDE 5
Basilio Bona - Dynamic Modelling 75
2.13.2 Infinitesimal Rotations
We have seen that, given two rotations respectively represented by two angular pa- rameters α1 = [ α11 α12 α13 ]T and α2 = [ α21 α22 α23 ]T, the rotation resulting from their composition cannot be obtained (except particular cases) from the simple sum of the two angles: α ̸= α1 + α2 Now let us consider an infinitesimal rotation defined by the infinitesimal angle vector dα = [ dα1 dα2 dα3 ]T the associated rotation matrix can be written as R(dα) = R(u(t), dθ(t)) = Ru(t),dθ(t) (2.119) From (2.115), and taking into account that for dθ → 0 we have 1 − cos (dθ) (dθ)2 → 0 sin (dθ) (dθ) → 0 the following approximate relation holds R (dα) ≃ I + S (dα) . (2.120) The inverse of R(dα) is computed as R(dα)−1 = R(dα)T ≃ I − S(dα) = I + S(−dα) ≃ R(−dα) This relation can be verified writing R(dα)TR(dα) = (I + S)(I − S) = I + S − S + S 2 ≃ I where S 2(dα) ≃ O since it contains second order infinitesimals, we conclude that the inverse of R(dα) is R(−dα). Using (2.120) it is possible to show that dα behaves like a proper vector. For instance, the commutative property holds; indeed the following identity is true: R (dα1) R (dα2) = (I + S (dα1)) (I + S (dα2)) = I + S (dα1) + S (dα2) + S (dα1) S (dα2) Since the last term goes to zero, and recalling (B.3), we have R (dα1) R (dα2) = I + S (dα1) + S (dα2) = I + S (dα1 + dα2) = R (dα1 + dα2) = R (dα2 + dα1) = R (dα2) R (dα1) (2.121)
SLIDE 6 76 Basilio Bona - Dynamic Modelling We know that from (2.110)) it results ˙ R(t) = S (ω(t)) R (u, θ(t)) (2.122) Since R does not depend on the choice of the representation of the parameters α used to characterize the rotation, the formula (2.122) can be re-written as: ˙ R (α(t)) = S (ω(t)) R (α(t)) (2.123) Considering the incremental rate we have ∆R = R (α + dα) − R (α) ≃ ˙ R(α)dt = S (ω(t)) R (α) dt = S (ω(t)dt) R (α) = S (dα) R (α) from which it follows R (α + dα) ≃ R (α) + S (dα) R (α) = (I + S (dα)) R (α) = R (dα) R (α) (2.124) We recall that the last term in the above identity does not commute. We conclude this Section with an observation: the role played by the skew-symmetric matrix S(ω) is fundamental in the definition of the angular velocity of a reference frame, and the two relations (2.114) and (2.122) are an example of this role. Now, if we introduce the following three matrices M i: M 1 = −1 1 M 2 = 1 −1 M 3 = −1 1 (2.125) thta are commonly called infinitesimal rotation generators), we are able to build S(ω) as a weighted sum S(ω) = ∑
i=1,3
ωiM i and observe that the M i form a basis for S(ω). Moreover the M i have the property that the differences between their mutual products obey to M iM j − M jM i ≡ [M i, M j] = ϵijkM k (2.126) where the ϵijk are called permutation symbols or Levi-Civita symbols defined as ϵijk = i = j, j = k, k = i; two symbols out of three are equal +1 (i, j, k) ∈ {(1, 2, 3), (2, 3, 1), (3, 1, 2)} −1 (i, j, k) ∈ {(1, 3, 2), (2, 1, 3), (3, 2, 1)} (2.127) The matrix operator in (2.126) is referred as a commutator or Lie bracket and it defines the Lie algebra on the orthonormal matrix group. The properties and the importance of such an algebra will not be detailed here; the interested reader can find more details in [34].
SLIDE 7
Basilio Bona - Dynamic Modelling 77
2.13.3 Infinitesimal Rotations and Quaternions
If the rotation angles are small it is possible to directly compute the correspond- ing quaternions and vice-versa. Let us assume to have a small rotation ∆θ rep- resented by the (small) angle parameters ∆α, and the corresponding quaternion h = [ h0 h1 h2 h3 ] ; from (2.120) we have R(∆α) = I + S(∆α) = 1 −∆α3 ∆α2 ∆α3 1 −∆α1 −∆α2 ∆α1 1 = r11 r12 r13 r21 r22 r23 r31 r32 r33 (2.128) comparing (2.128) with (2.89), we obtain: ∆α1 = 1 2(r32 − r23) = 2h1h0 ≃ 2h1 ∆α2 = 1 2(r13 − r31) = 2h2h0 ≃ 2h2 ∆α3 = 1 2(r21 − r12) = 2h3h0 ≃ 2h3 Taking into consideration definition (2.119), the approximation is reasonable; in- deed, applying (2.114), we obtain h0 = cos ∆θ 2 ≃ 1 since ∆θ is a small angle. Hence, if we know the value of an infinitesimal rotation dα = [ dα1 dα2 dα3 ]T, we can easily compute the corresponding quaternions as h ≈ [ 1 dα1 2 dα2 2 dα3 2 ]
2.14 Total Velocity and Acceleration of Points
When we speak of linear or angular velocites of a geometrical point or a rigid body, we mean the relative velocity of the geometrical point with respect to a well defined reference frame, or the velocity of the reference frame describing the rigid body with respect to another well defined reference frame. To be precise, when we say that a rigid body has at the time t an angular velocity ω(t), we mean that at t the relation (2.116), that we report below for the sake of clarity, is true: ω(t) = ˙ θ(t)u(t) + sin θ(t) ˙ u(t) + (1 − cos θ(t)) S (u(t)) ˙ u(t) This means that we have to give non only ˙ θ(t), but also the instantaneous rotation axis u(t) at time t. If the rotation axis remains constant, then (2.116) reduces to (2.117).
SLIDE 8 78 Basilio Bona - Dynamic Modelling It is important to note that the vector ω(t) cannot be translated parallel to itself without modifying the physical significance of the described phenomenon, since the application point of ω(t) gives the position of the axis around which the rotation takes place; this is not valid for the linear velocity given by the vector v(t): this one is a free vector, while ω(t) is an applied vector. We recall that an applied vector is defined by the couple (P, w), where P is the geometrical point of application and w is defined such that − → PQ = w To be more precise, when we speak of angular (or linear) velocities we should use the following (cumbersome) notation: ωk
ij(P)
(2.129) where we denote the relative angular velocity of the reference frame (or of the associated rigid body) Rj with respect to the reference frame Ri, represented in the reference frame Rk, and where P is the application point of ω; when P is known
- r superfluous, we only write ωk
- ij. Usually we implicitly assume to represent the
angular velocity as ωi
ij(P).
Whit this notation it is always true that ωk
jj = 0
∀j, k (2.130) since no relative motion will take place, while, in general, it results ωj
ij ̸= ωi ij
∀i, j (2.131) provided that Ri and Rj are non-coincident. Moreover, given the vectorial nature
- f the velocity, the additive property holds
ωk
im + ωk mj = ωk ij
from which we have the following property ωk
ij + ωk ji = ωk ii = 0
⇒ ωk
ij = −ωk ji
(2.132) that says that, if a reference frame rotates with an angular velocity with respect to another reference frame, this last rotates with an angular velocity that is the same in module, but opposite in direction to the previous one. Now assume that x 1 is the representation of the geometrical point A in R1 and that this frame roto-translates with respect to R0; assume also that R0
1(t) and t0 1(t) are
respectively the rotation matrix and the translation vector between the two frames, and that ω0
01(P) is the relative angular velocity of R1 with respect to R0, represented
in R0. We have still to define the application point P; it varies according to the fact that the rotation takes place with respect to the fixed frame R0 axes, or with respect to the translated local frame R1 axes. In the following we will consider the two cases.
SLIDE 9 Basilio Bona - Dynamic Modelling 79 Case A We consider here a roto-translation described by T 0
1 = T(d 0 1)T(R0 1)
(2.133) that represents a translation d 0
1, followed by a rotation R0 1, with instantaneous
rotational velocity ω0
01, performed with respect to the mobile frame (post-mobile
rule). The motion defined in (2.133) allows the following alternative description: “a rota- tion followed by a translation performed with respect to the fixed frame (pre-fixed rule)”. In this case we can state that the application point P of the applied vector(P, ω) is the origin of R1. In this case the vector r 1(t) represented in R1 is described in R0 by: r 0(t) = R0
1(t)r 1(t) + d 0 1(t);
(2.134) Computing the time derivative of this identity, recalling (2.113) and omitting for brevity to write the time dependence in the various terms, we have: ˙ r 0 = ˙ R
1r 1 + R0 1 ˙
r 1 + ˙ d
1
= S (ω0
01) R0 1r 1 + R0 1 ˙
r 1 + ˙ d
1
= ω0
01 ×
( R0
1r 1)
+ R0
1 ˙
r 1 + ˙ d
1
(2.135) If now we define ρ0(t) = R0
1r 1(t)
(2.136) we have ˙ r 0(t) = ω0
01(t) × ρ0(t)
+ R0
1 ˙
r 1(t)
+ ˙ d
1(t) C
(2.137) This is the classical relation that gives the total linear velocity in R0 of a point P that moves with a linear velocity ˙ r 1(t) in the reference system R1, itself moving at a linear velocity ˙ d
1(t) and angular velocity ω0 01(t) applied at the origin of R1.
The three terms are: A) the tangential velocity due to the rotational velocity of the frame; B) the proper velocity of the point, represented in R0; C) the translational velocity of the frame. The acceleration ¨ r 0(t) of the point P expressed in R0 can be computed considering the identity: d (a × b) dt = (da dt × b ) + ( a × db dt )
SLIDE 10 80 Basilio Bona - Dynamic Modelling Taking the time derivative of ˙ r 0(t) in (2.137), using (2.110) and (2.136), and omitting again to write the time dependence, we obtain: ¨ r 0(t) = d dt (ω0
01(t) × ρ0(t)) + d
dt ( R0
1 ˙
r 1(t) ) + d dt ˙ d
1(t)
= ˙ ω0
01 × ρ0 + ω0 01 × ˙
r 1 + ˙ R
1 ˙
r 1 + R0
1¨
r 1 + ¨ d
1
= ˙ ω0
1 × ρ0 + ω0 01 × d
dt ( R0
1r 1)
+ ˙ R
1 ˙
r 1 + R0
1¨
r 1 + ¨ d
1
= ˙ ω0
01 × ρ0 + ω0 01 ×
( ˙ R
1r 1 + R0 1 ˙
r 1) + ˙ R
1 ˙
r 1 + R0
1¨
r 1 + ¨ d
1
= ˙ ω0
01 × ρ0 + ω0 01 × (ω0 01 × ρ0) + ω0 01 ×
( R0
1 ˙
r 1) + ˙ R
1 ˙
r 1 + R0
1¨
r 1 + ¨ d
1
= ˙ ω0
01 × ρ0 + ω0 01 × (ω0 01 × ρ0) + ω0 01 ×
( R0
1 ˙
r 1) + ω0
01 ×
( R0
1 ˙
r 1) +R0
1¨
r 1 + ¨ d
1
= ˙ ω0
01 × ρ0
+ ω0
01 ×
( ω0
01 × ρ0)
+ 2ω0
01 ×
( R0
1 ˙
r 1)
+ R0
1¨
r 1 + ¨ d
1
(2.138) where A) is the tangential acceleration; B) is the centripetal acceleration; C) is the Coriolis acceleration; D) is the sum of the linear acceleration of the point represented in the fixed frame and the acceleration of the frame (already represented in R0). Case B We consider here a roto-translation described by
1(t) = T
( R0
1
) T ( d 0
1
) (2.139) that represents a rotation R0
1 with instantaneous rotational velocity ω0 01, followed
by a translation d 0
1 performed with respect to the mobile frame (post-mobile rule).
The motion defined in (2.139) allows the following alternative description: “a trans- lation followed by a rotation performed with respect to the fixed frame (pre-fixed rule)”. In this case we can state that the application point P of the applied vector(P, ω) is the origin of R0. The point r 1 is described in R0 by: r 0(t) = R0
1(t)r 1(t) + R0 1(t)d 0 1(t);
(2.140) Taking the time derivative and omitting again to write the time dependence, we
˙ r 0 = ˙ R
1r 1 + R0 1 ˙
r 1 + ˙ R
1d 0 1 + R0 1 ˙
d
1
= S(ω0
01)R0 1
( r 1 + d 0
1
) + R0
1 ˙
r 1 + R0
1 ˙
d
1
= ω0
01 × (R0 1r 1 + R0 1d 0 1) + R0 1 ˙
r 1 + R0
1 ˙
d
1
(2.141) If now we define ρ0 = R0
1r 1 + R0 1d 0 1
(2.142)
SLIDE 11 Basilio Bona - Dynamic Modelling 81 we have ˙ r 0 = ω0
01 × ρ0 + R0 1 ˙
r 1 + R0
1 ˙
d
1
(2.143) that is identical to (2.137), with R0
1 ˙
d
1 instead of ˙
d
- 1. Taking the time derivative of
(2.143) we obtain the total acceleration: ¨ r 0(t) = ˙ ω0
1 × ρ0 A
+ ω0
01 ×
( ω0
01 × ρ0)
+ 2ω0
01 ×
( R0
1 ˙
r 1 + R0
1 ˙
d
1
)
+ R0
1¨
r 1 + R0
1¨
d
1
(2.144) The formula has the A), B), C) and D) terms identical with those in (2.138), the only difference being the translational velocity and acceleration of the frame R1 expressed in the fixed reference R0. In conclusions, the total acceleration of a point that moves in a local reference frame, that itself moves with respect to a fixed frame, is the sum of the following terms (time dependency omitted for brevity): Tangential acceleration Due to rotational acceleration of the local frame with respect to the fixed frame: A = ˙ ω0
01 × ρ0
where ρ0 is the distance of the point P from the origin of the fixed frame. Centripetal acceleration Due to the variation of the tangential acceleration, that produce an acceleration directed toward the instantaneous rotation center: B = ω0
01 ×
( ω0
01 × ρ0)
Coriolis acceleration Due to the linear velocity in a rotating frame: C = 2ω0
01 × v 0
where v 0 = R0
1 ˙
r 1
v 0 = R0
1( ˙
r 1 + ˙ d
1)
is the total linear velocity; Linear acceleration Due to the linear acceleration of the point P and of the reference frame: a0 = R0
1¨
r 1 + ¨ d
1
a0 = R0
1(¨
r 1 + ¨ d
1)
Notice that the symbol ρ0 has a different meaning in the two cases, as put in evidence by (2.136) and (2.142).
