Rigid Body Dynamics If the aim of kinematics is to describe the body - - PDF document

SLIDE 1

Rigid Body Dynamics

If the aim of kinematics is to “describe” the body motion, the aim of dynamics is to “explain” it; the history of mechanics shows that the passage from description to explanation requires the introduction of a new physical entity, that of mass or, in alternative, that of force. The fundamental law of mechanics, due to Newton, is analytically expressed as a vectorial equation between force, mass and acceleration of a translating rigid mass f = ma

• r else

  fx fy fz   = m   ax ay az   where the linear acceleration a = [ ax ay az ]T is a kinematic quantity, obtained ad the time derivative of the body linear velocity. This law is true in principle only for a single point-mass particle with mass m, where the applied force f and the acceleration a can exchange the role of cause and effect: if a force f is applied to the particle, the particle accelerates with a linear acceleration equal to a, and conversely, if a particle has a linear acceleration a, the particle is subject to a force f proportional to its mass. If we jointly know the two vectors f = [ fx fy fz ]T and a = [ ax ay az ]T, we can compute the mass from one or any of the following scalar relations as m = fx ax = fy ay = fz az Though, in this last case, we are applying a circular argument, since no definitions

• f “acceleration” or “force” exist that are independent from the measurement of a

mass; we are therefore compelled to use some trick, as clearly illustrated in [?]. If a single point mass is connected to others to form a rigid body, the Newton law is still valid, given that we observe some precautions. Every point-mass shall be isolated and we must consider and deduce the forces applied on it by the other masses; that is, we must introduce the constraint forces in addition to the external forces applied on the body. It is important to notice that the vector equations are dependent on the representa- tion used to characterize its components, and change changing the reference frame 1

SLIDE 2

2 used to represent the force and acceleration vectors, although the masses, at least in a non-relativistic motion state, do not vary changing the reference frame. In this textbook we will use an analytical approach, as defined in [?], instead of the vectorial one due to Newton. We consider the multibody system as a system in which the dynamical equations derive from a unifying principle that implicitly includes and generates these equations. This principle is based on the fact that, in order to describe the motion of a multi- body system is sufficient to consider and use in a proper way some suitable scalar quantities; these were in origin called by Leibnitz vis viva and work function, nowa- days take the name of kinetic energy and potential energy. They are an example of the so-called state functions, since to each value of the state vector1 This general principle takes the name of principle of least action; it can be roughly described in the following way. Let us consider the space Q of the generalized coordinates q ∈ Q, as sketched in Figure 8 for a two-dimensional space Q; let us assume that a particle starts its motion at time t1 in the state Q1 = q(t1) and ends it motion at time t2 having reached the state Q2 = q(t2). Let us further assume that the motion keeps constant the sum E = C + P of the kinetic energy C and the potential energy P that the particle has at time t1. Given the continuity of motion, the two points Q1 and Q2 are connected by a continuous path (or trajectory), as the one represented by a continuous line in Figure 8; this trajectory is called the true trajectory, and at least in principle, it is unknown, since it is exactly what we want to compute as the result of the dynamical equation analysis. If we had chosen at random a different trajectory, with the only condition that the two boundary point remain fixed, and called this trajectory a perturbed trajectory, the chance to obtain exactly the true trajectory would have been minimal. Now we ask what is that characterize the true trajectory with respect to all possible

• ther perturbed trajectories. Euler was the first mathematician to contribute to

the solution of this problem, but it was Lagrange that developed a complete theory, that was later extended by Hamilton: the true trajectory is the one that minimizes the integral of the so-called vis-viva (i.e., two times the kinetic energy) of the entire motion between Q1 and Q2. This integral is called action and has a constant and well defined value for each perturbed trajectory at constant E. The least action principle states that the nature “chooses”, among the infinite number of trajectories starting in q(t1) and ending in q(t2), the trajectory that

1The concept of state will be defined in Section XXX; for the moment we simply consider that

the state corresponds to the two vectors q(t) and ˙ q(t).

SLIDE 3

3 Figure 1: True and perturbed trajectories in the configuration space Q. minimizes the definite integral S of a particular state function C∗(q(t), ˙ q(t))2 that depends on the generalized coordinates q(t) and the generalized velocities, where ˙ q(t) S = ∫ t2

t1

C∗(q(t), ˙ q(t))dt (1) In order to apply this principle it is therefore necessary to compute the trajectory in the space Q that minimizes S. The integral in (1) is computed between the initial time t1 and the final time t2 and must obey to the boundary constraints active at these two times. The minimization of a functional3 is based on a particular mathematical technique, called calculus of variations, whose illustration goes beyond the scope of the present text. The interested reader can find additional treatment in [?], Appendix D. FARE APPENDICE? The conditions that assure the attainment of the minimum of S provide a set of differential equations that contain the first and second time derivatives of the qi(t); this set completely describes the system motion dynamics. As a final comment, we can say that the least action principle makes use of a concept

• the action - that overturns the usual representation of the physical system dynamics

2We will see later what this state function C∗ is. 3A functional is a mapping between a function and a real number; the function shall be con-

sidered as a whole, i.e., not a single particular value; in this sense a functional is often the integral

• f the function, as in the above equation for S.

SLIDE 4

4 by means of a set of differential equations. With some approximation we can say that the differential equations specify the evolution of a physical quantity as the result of infinitesimal increments of time or position; summing up this infinitesimal variations we obtain the physical variables at every instant, knowing only their initial value and possibly some initial derivative: we can say that the motion has a local representation. On the contrary, the action characterizes the motion dynamics requiring only the knowledge of the states at the initial and final times; every intermediate value of the variables can be determined by the minimization of the action, that is a global, rather than local, measure. In any case the implementation of the principle of least action produces as a result a set of differential equations, different from these obtained with the local representation. In the following Sections we will illustrate how to obtain these equations.

0.1 Point Mass System

Yo introduce the Lagrangian approach and compare it with the Newtonian approach we must define some general physical entities Let us start to consider a sate of N point masses mi, with i = 1, · · · , N), as schemat- ically represented inn Figure 2. These systems are called multi-point systems. Figure 2: A system of N point masses mi. The position r i and the velocity v i are represented in a generic reference frame Rb. Now we consider a pure rotation motion around a point O, that, for simplicity, is the origin of the reference frame Rb; all vectors will be represented in this frame. If necessary we indicate the vectors as r b

i, otherwise with r i.

The position of the generic mass mi is defined by the vector r i = [ xi yi zi ]T. If the rotation velocity of the system is given by the vector ω, every mass will acquire

SLIDE 5

0.1. POINT MASS SYSTEM 5 a linear velocity v i, as the result of the rotation, whose value is v i(t) = ω(t) × r i(t) = S(ω(t))r i(t) (2) We define the linear momentum4 pi(t) as the product between the mass mi and the velocity v i(t) pi(t) = miv i(t) In this Chapter, the symbol p indicates the linear momentum an not the cartesian pose of a body or the position of a point, as in the preceding Chapters.

0.1.1 Moment of a Force

We briefly recall that, given a point mass located in a point P represented by the vector r in the reference frame R, and a force f applied to it, the moment of the force with respect to a point O is given by the cross product r × f = r × mdv dt = r × dp dt (3) where the terms on the right side arise from the Newton equation. It is important to recall that the term “moment” indicates the cross product, while the term “momentum” indicates the vector p = mv. Applying the derivative rules in (??) to the right end side of (3), we obtain r × f = r × dp dt = d dt(r × p) − dr dt × p = d dt(r × p) − v × p = dh dt where the last term is always zero, since the moment p and the linear velocity of the point mass v are always collinear, and where the vector h = r × p is called angular moment or moment of the momentum. Since the moment r × f is often referred as the torque applied by the force with respect to the point O, we can write τ = d dth Note the analogy f = d dtp = ˙ p τ = d dth = ˙ h

4In Italian quantit`

a di moto (lineare).

SLIDE 6

6 that shows the relations between the generalized force/torque and the generalized momentum. Figure 3: A point mass m with the applied force f and it linear moment p.

0.1.2 Angular Momentum and Inertia Matrix

Given a point O in space and a point mass mi with position r i, belonging to a set of point masses, we have define the angular momentum, or moment of momentum5 hi, as hi(t) = r i(t) × pi(t) = r i(t) × (miv i(t)) Replacing (2) in the above relation and omitting for simplicity the time dependency, we obtain hi = mi(r i × (ω × r i)) Summing up all the N point masses contributions, we have the total angular mo- mentum h =

N

∑

i=1

hi = ∑

i

mi(r i × (ω × r i)) (4) Recalling from (??) that the triple cross product enjoys of the following property a × (b × c) = (aTc)b − (aTb)c, where aTc and aTb are scalar products, it follows that the momentum can be written as h = ∑

i

mi ( (r T

i r i)ω − (r T i ω)r i

)

5In Italian momento della quantit`

a di moto or quantit` a di moto angolare.

SLIDE 7

0.1. POINT MASS SYSTEM 7 Replacing with r2

i the square norm ∥r i∥2 = r T i r i = (x2 i + y2 i + z2 i ) and writing

explicitely r T

i ω = (xiωx + yiωy + ziωz) we obtain

h = ∑

i

mi  r2

i

  ωx ωy ωz   − (xiωx + yiωy + ziωz)   xi yi zi    

• r

h =   hx hy hz   = ∑

i

   mi(r2

i − x2 i )ωx − mixiyiωy − mixiziωz

−mixiyiωx + mi(r2

i − y2 i )ωy − miyiziωz

−mixiziωx − miyiziωy + mi(r2

i − z2 i )ωz

   (5) This relation can be written in matrix form as h = ∑

i

  Γxx,i Γxy,i Γxz,i Γyx,i Γyy,i Γyz,i Γzx,i Γzy,i Γzz,i     ωx ωy ωz   = ∑

i

Γ iω (6) where Γxx,i = mi(r2

i − x2 i ) = mi(y2 i + z2 i )

Γyy,i = mi(r2

i − y2 i ) = mi(x2 i + z2 i )

Γzz,i = mi(r2

i − z2 i ) = mi(x2 i + y2 i )

and Γxy,i = Γyx,i = −mixiyi Γxz,i = Γzx,i = −mixizi Γyz,i = Γzy,i = −miyizi If we now introduce a new matrix defined as Γ =    Γxx Γxy Γxz Γyx Γyy Γyz Γzx Γzy Γzz    = ∑

i

   Γxx,i Γxy,i Γxz,i Γyx,i Γyy,i Γyz,i Γzx,i Γzy,i Γzz,i    = ∑

i

Γ i from (6) we have h =   Γxx Γxy Γxz Γyx Γyy Γyz Γzx Γzy Γzz     ωx ωy ωz   = Γω (7)

• r better, restating the time dependency

h(t) = Γ(t)ω(t)

SLIDE 8

8 The matrix Γ is called inertia matrix or inertia tensor and has on the main diagonal the inertia moments Γxx = ∑

i

mi(r2

i − x2 i ) =

∑

i

mi(y2

i + z2 i ) =

∑

i

mid2

x

Γyy = ∑

i

mi(r2

i − y2 i ) =

∑

i

mi(x2

i + z2 i ) =

∑

i

mid2

y

(8) Γzz = ∑

i

mi(r2

i − z2 i ) =

∑

i

mi(x2

i + y2 i ) =

∑

i

mid2

z

and outside the diagonal the inertia products Γxy = Γyx = − ∑

i

mixiyi Γxz = Γzx = − ∑

i

mixizi (9) Γyz = Γzy = − ∑

i

miyizi Inertia moments As one can notice, the inertia moments are obtained summing the contribu- tions of the products of each mass by the square euclidean distance d2

x, d2 y, d2 z

from the three axes x, y and z. Prodotti d’inerzia When the inertia products are zero, the inertia matrix becomes diagonal Γ =   Γxx Γyy Γzz   In this case the reference frame axes are aligned along what are called the principal inertia axes of the body, and the matrix is the so called principal inertia matrix. Another possible representation of h and Γ can be obtained using (4) and the properties of the skew-symmetric matrices: h = ∑

i

mi(r i × (ω × r i)) = ∑

i

miS(r i)S(ω)r i (10) Recalling that S(ω)r i = −S(r i)ω, one obtains: h = ∑

i

miS(r i)S(ω)r i = ∑

i

−miS(r i)S(r i)ω (11)

SLIDE 9

0.1. POINT MASS SYSTEM 9 and being ω a common term, at the end one has h = ( − ∑

i

miS(r i)S(r i) ) ω = Γω (12) with Γ = − ∑

i

miS(r i)S(r i) (13) The inertia matrix Γ depends on the r i vectors that characterize the positions of the masses mi in a given reference frame, with respect to a given point; this point can be fixed, as the origin O of the frame, or mobile (i.e., local) as the center of mass C of the body; in the first case the inertia matrix will be indicated as Γ o, in the second case as Γ c. Considering two reference frames with a common origin, but one rotated with respect to the other, the inertia matrices will be different; the rule to transform one into the

• ther will be detailed later on.

The inertia matrix Γ represents the inertial properties of a rigid body subject to a rotation, in the same way as the mass represents the inertia properties of a body subject to a translation. Equation (7) is qualitatively similar to the expression of the total linear momentum p(t) = mv(t). Since we can write it as p(t) = M v(t), where M = mI , the two momentum are p(t) = M v(t) h(t) = Γ(t)ω(t) with the difference that usually the mass matrix M does not depend on time, while Γ(t) does. Example 0.1.1 Let us consider a simple body system composed by two equal masses as in Figure 4; we want to compute the angular momenta and the inertia matrices with respect to the points O and C. Assume also that Ra is fixed, while Rb is on the barycenter (or center of mass) of the body having the axis i aligned with the segment − − → CB. Assume r 1 =   1 2   ; r c =   2 2   ; r 2 =   3 2   ; c1 =   −1   ; c2 =   1  

SLIDE 10

10 Figure 4: The two equal masses system considered in Example 0.1.1. and m1 = m2 = 1. Since Ra e Rb are parallel, the rotation matrix Ra

b = I and we can neglect it.

We start computing Γ o as Γ o = −m1S(r 1)S(r 1) − m2S(r 2)S(r 2) =   8 −8 −8 10 18   and then Γ c as Γ c = −m1S(c1)S(c1) − m2S(c2)S(c2) =   2 2   Recalling the parallel axes theorem, that we will introduce in Paragraph 0.3.3, and assuming also (23) and (24), we obtain Γ o − Γ c =   8 −8 −8 8 16   = −(m1 + m2)S(r c)S(r c) Now assuming that the angular velocity of the system is ωa

ab = ωb ab =

  1   we can compute the momenta as follows ho =   8 −8 −8 10 18     1   =   18   ; hc =   2 2     1   =   2  

SLIDE 11

0.2. DISTRIBUTED MASS SYSTEMS 11

0.2 Distributed Mass Systems

We make a further step assuming that the body whose angular moment we want to compute is composed by a continuous distribution of infinitesimal point masses dm, having point density ρ(r) = ρ(x, y, z), depending on the position (x, y, z), and where the infinitesimal mass is dm = ρ(x, y, z)dV . The coordinate r takes its values in V, where V is a finite space region that in- cludes the body, with a total volume V = ∫

V dV and total mass mtot =

∫

V dm =

∫

V ρ(x, y, z)dV .

We can apply the same arguments used above for N point masses to a system of distributed masses, replacing in (8) and (9) the summation operator with the integral

• perator in the volume V. We can now write in compact form the inertia matrix as:

Γ = ∫

V

ρ(x, y, z)   y2 + z2 −xy −xz −xy x2 + z2 −yz −xz −yz x2 + y2   dV (14)

• btaining the inertia moments

Γxx = ∫

V

ρ(r)(y2 + z2) dV Γyy = ∫

V

ρ(r)(x2 + z2) dV (15) Γzz = ∫

V

ρ(r)(x2 + y2) dV and the inertia products Γxy = Γyx = − ∫

V

ρ(r)xy dV Γxz = Γzx = − ∫

V

ρ(r)xz dV (16) Γyz = Γzy = − ∫

V

ρ(r)yz dV Notice that ρ(r)dV is equivalent to an infinitesimal mass dm(r), therefore we can also write relation (14) as Γ = ∫

V

  y2 + z2 −xy −xz −xy x2 + z2 −yz −xz −yz x2 + y2   dm(r) (17)

0.2.1 Some Observations

The inertia matrix is implicitly defined when we define the angular momentum h as in (10); this one in turn is defined introducing a rotation around a point O,

SLIDE 12

12 characterized by the angular velocity ω. The cartesian components of h are represented in the chosen reference frame, that

• ften, but not always, coincides with a reference frame with its origin in O; the

inertia matrix (or inertia tensor) describes the way the mass is distributed with respect to the axes of R0(O; x, y, z). If the frame is rigidly attached to the body, the inertia matrix is time independent, otherwise it changes with time. This is the reason why it is common to read in textbooks “inertia moment with respect to a point” and/or “inertia matrix with respect to the axes”. This difference shall not generate confusion since it is an example of a casual use of mathematical concepts, that, once understood, does not represents a problem anymore. A reference frame attached to the body is called a body-frame and any inertia matrix with respect to this frame is always constant. If one chooses a different reference frame attached to the body, with the same ori- gin O, but different axes, the components of h do change and therefore also the components of Γ, but they remain time-invariant. If the frame is translated, Γ changes according to the parallel axes theorem (see Section 0.3.3) If the body rotates with respect to the reference frame (and therefore this one is not a body-frame), the components of Γ vary with time. To conclude, we recall that what has been presented here is valid only for rigid bodies.

0.3 Mathematical Formulation

We have seen that a inertia matrix or tensor6 of a rigid body specifies and summarize the inertial characteristics of the body with respect to the variations of the angular velocity associated to the body rotation. The matrix is always defined specifying explicitly or implicitly a point with respect to which it is computed; usually this point is the center of mass of the body, but this choice it is not strictly necessary. In order to define correctly the inertia matrix of a rigid body it is proper to introduce a number of subscripts and superscripts that make clear the meaning with respect to the reference frames defined by the user. Given a rigid body B, identified by the subscript b, and a reference system Rb, called

6The term tensor is used to indicate an abstract mathematical entity that generalizes the

vectors, while we use the term matrix to indicate its representation in some defined reference

• frame. We prefer to use “matrix” instead of “tensor” because it is more common in modelling and

control literature

SLIDE 13

0.3. MATHEMATICAL FORMULATION 13 local frame, the position of the center of mass C of the body is defined by the vector r b

c, computed as follows:

r b

c

∫

B

dm = r b

c mb tot =

∫

B

r bdm (18) where r b represents the position in Rb of the generic point mass dm belonging to the body B; moreover mb

tot =

∫

B dm is the total mass of the body B.

The inertia matrix Γ b

c ∈ R3×3 around the center of mass C is defined, as seen in

Section 0.1.2, resorting to the angular momentum definition, that we write again here, with the addition of the introduced subscripts and superscripts: hb

c = Γ b cωb

(19) where hb

c is the momentum of the rigid body with respect to C and ωb 0 ≡ ωb 0b is the

total angular velocity of the rigid body; both are 3D vectors represented in Rb. Take care not to confuse the total angular velocity ωb

0b represented in Rb and the

same velocity represented in R0, that we shall indicate with the symbol ω0

0b.

The use of various subscripts and superscript in Γ b

c is non a useless pedantry; on

the contrary it highlights the fundamental points of its definition: Γ b

c is the inertia

matrix of the rigid body B with respect to its center of mass C, represented in the reference frame Rb, that in general can be placed everywhere one likes, but usually attached to the body. The inertia matrix changes changing the reference frame, as we will see later. If now we attach Rb to the rigid body with its origin in C we have from (18) r b

c = 0

→ ∫

B

r bdm = 0 (20) and the inertia matrix Γ b

c will be defined as:

Γ b

c = −

∫

B

S(r b)S(r b)dm = ∫

B

[ r b 2 I − r b(r b)T] dm =   Γxx Γxy Γxz Γyx Γyy Γyz Γzx Γzy Γzz   (21) where we notice that r b(r b)T is a rank one dyadic product (see Appendix ??), and that S(x) is a skew-symmetric matrix (see Appendix ??). Observe the analogy of the second term in (21) with the second term in (13).

0.3.1 Rotation of the Reference Frame

If now we consider another reference frame Rk, rotated around the common origin (in this case C) with respect to Rb, and the rotation matrix given by Rk

b, then the

SLIDE 14

14 inertia matrix is computed according the following relation Γ k

c = Rk bΓ b cRb k = Rk bΓ b c(Rk b)T

(22)

• r

Γ k

cRk b = Rk bΓ b c.

using the usual notation Ra

b = (Rb a)T.

When the generic reference frame Rb is attached at the body B and moving with it, i.e., it is a body-frame, the inertia matrix relative to Rb is time invariant; in the

• ther case it is time variant.

0.3.2 Principal Inertia Matrix

Often it is appropriate to refer to a “privileged” reference frame R∗, implicitly defined as that frame having its origin in the body center of mass and with respect to which the inertia matrix is diagonal, i.e., Γ ∗

c =

  Γx Γy Γz   ; This frame has the three unit vectors i ∗, j ∗ andk ∗, aligned along the so-called allineati principal inertia axes; the matrix Γ ∗

c is called prende il nome di prin-

cipal inertia matrix.

0.3.3 Parallel Axes Theorem

Assume that we know the inertia matrix with respect to a reference frame with

• rigin in the center of mass C, and we want to compute it with respect to another

frame, with different origin O, but parallel axes. To do so we must express the relation between the infinitesimal mass ρ dV = dm with respect to the two points O e C is (vedi Fig. 6) It is straightforward to see that: r b

i = r b c + cb i,

∀i = 1, . . . , N where r b

c =

[ tx ty tz ]T represents the translation of the reference frame from O to C. As we already know, the inverse translation, from C to O, is given by −r b

c.

This relation is valid for any vector represented in a couple of translated reference frames.

SLIDE 15

0.3. MATHEMATICAL FORMULATION 15 Figure 5: An example of parallel axes inertia computation. The point mass element dmi that was in the position cb

i, is now in position r b i.

It is now possible to compute the new inertia matrix Γ b

• , taking into account (21),

as Γ b

• = Γ b

c − mtotS(r b c)S(r b c) = Γ b c + mtot

[ r b

c

• 2 I − r b

c(r b c)T]

(23) Simplifying the notation as follows Γ b

c ≡ Γ =

  Γxx Γxy Γxz Γyx Γyy Γyz Γzx Γzy Γzz   e Γ b

• ≡ Γ ′ =

  Γ ′

xx

Γ ′

xy

Γ ′

xz

Γ ′

yx

Γ ′

yy

Γ ′

yz

Γ ′

zx

Γ ′

zy

Γ ′

zz

  with t ≡ r b

c and with m ≡ mtot, we will obtain

Γ ′ = Γ − mS(t)S(t) (24) = Γ + m [ ∥t∥2 I − ttT] = Γ + m   (t2

y + t2 z)

−txty −txtz −txty (t2

x + t2 z)

−tytz −txtz −tytz (t2

x + t2 y)

  that gives Γ ′

xx = Γxx + m

( t2

y + t2 z

) Γ ′

yy = Γyy + m (t2 x + t2 z)

Γ ′

zz = Γzz + m

( t2

x + t2 y

) . (25)

SLIDE 16

16 Figure 6: Parallel axes theorem: relation between the infinitesimal mass position in the two cases. The inertia products can be written as Γ ′

xy = Γxy − m txty

Γ ′

xz = Γxz − m txtz

Γ ′

yz = Γyz − m tytz

(26) The two relations (22) and (23) allow to compute the inertia matrix with respect to any other point O of choice, knowing the inertia matrix with respect to the point C. The relation introduced above are valid both for discrete point-mass systems and for distribute mass systems.

0.4 Angular Momentum and Euler Equation

The expression of the angular momentum has been introduced in (19) without any additional explanation; here we justify it in a more formal way. We can assume that a rigid body B is composed by a continuum of point-mass particles with mass dm located at r k, with velocity v k = ˙ r k and acceleration ak = ˙ v k; the superscript k identifies the generic reference frame Rk with respect to which the vectors are represented. Since B is a rigid body, there will be a number of constraint relations linking positions, velocities and accelerations, as detailed in Section ??. To simplify the mathematical treatment of the rigid body dynamics, we can reduce the action of the forces and torques acting on the body B to two equivalent effects:

SLIDE 17

0.4. ANGULAR MOMENTUM AND EULER EQUATION 17

• an equivalent linear force f k

c, with its line of action across the body center

• f mass C, and with a module equal to the module of the vector sum of the

applied forces;

• an equivalent torque τ k

c, equal to the total moment of produced by the forces

acting on B, relative to C. The second Newton Law states that the resulting inertial force is equal to the time derivative of the linear momentum f k

c =

( d dtmv k

c

) = m ( d dtv k

c

) = mak

c

where v k

c is the instantaneous velocity of C, m is the total mass of the body and

ak

c = ˙

v k

c is the total acceleration of the center of mass; the relation is valid for any

choice of the reference frame Rk used to describe the vectors. The Euler equation states that the inertial torque is equal to the time derivative

• f the angular momentum

τ k

c =

( d dthk

c

) where in both cases we have pointed out that the total derivative must be performed with respect to an inertial7 (or pseudo-inertial) reference frame R0; It is not neces- sary that Rk is an inertial frame, but only that the derivative, seen as the limit of the difference quotient between vectors, is computed in the inertial frame (see also Section ??). The vector hk

c is the angular momentum with respect to the mass center C and is

defined, in analogy to what has been introduce in Section 0.1.2, as follows: hk

c =

∫

B

(ck × v k)dm (27) where ck indicates the position of an element of infinitesimal mass dm with respect to the center of mass C of the body (see also Figure 7). Since the position of an element dm with respect to the origin O of the reference frame Rk can be decomposed into the sum r k = ck + r k

c

it results ck = r k − r k

c

where we have indicated with r k

c the position of the center of mass C and with r k

the position of dm with respect to O.

7We recall that an inertial reference frame must be fixed or in linear uniform motion, i.e.,

SLIDE 18

18 Figure 7: The position r k of on infinitesimal mass element dm is decomposed into the sum r k

c + ck. The reference frame Rk is assumed to be inertial and all vectors

are represented in it. Corretto fin qua Therefore, the velocity v k = ˙ r k of the generic elementary mass can be expressed in relation to the velocity of the center of mass v k

c as follows

v k ≡ ˙ r k = ˙ ck + ˙ r k

c = v k c + ωk kb × ck

(28) The second term of the right hand side is obtained considering the relation (??), rewritten here for convenience comodit` a ˙ r 0(t) = ω0

01(t) × ρ0(t) + R0 1 ˙

r 1(t) + ˙ d

1(t)

(29) where ˙ r 1(t) = ˙ d

1(t) = 0, ρ0(t) ≡ ck and ωk kb ≡ ω0 01(t)

As usual, the angular velocity ωk

kb, is the total angular velocity of B with respect to

Rk and is represented in Rk itself. Replacing (28) in the (27), we obtain hk

c =

∫

B

ckdm × v k

c +

∫

B

ck × (ωk

kb × ck)dm =

the term ∫

B ckdm goes to zero since, in accordance with the definition of the mass

center in (18), one can write: ∫

B

ckdm = ∫

B

(r k − r k

c)dm =

∫

B

r kdm − mr k

c = 0

• ne in which the accelerations acting on it are zero or negligible; in this last case the reference

frame is pseudo-inertial. For example, a reference frame on your study table experiences only the accelerations due to the Earth rotation around its axis and around the Sun, to neglect other more complex motions of the Sun in the local galaxy, etc. In this textbook we call “inertial” any pseudo-inertial reference frame.

SLIDE 19

0.4. ANGULAR MOMENTUM AND EULER EQUATION 19 being ∫ r kdm = mr k

c

In conclusion, we have hk

c =

∫

B

ck × (ωk

kb × ck)dm = −

∫

B

ck × (ck × ωk

kb)dm =

= − ∫

B

S(ck)S(ck)ωk

kbdm = −

∫

B

S(ck)S(ck)dm ωk

kb =

= definition (21) = = Γ k

cωk kb

(30) If we represent the vectors and the inertia matrix in a different reference frame Rℓ we write hℓ

c = Γ ℓ cωℓ kb.

with the relation between the inertia matrices Γ ℓ

c = Rℓ kΓ k cRk ℓ

(31) and the angular velocities ωℓ

kb = Rℓ kωk kb

therefore we can observe that hℓ

c = Γ ℓ cωℓ kb = Rℓ kΓ k c Rk ℓRℓ k I

ωk

kb = Rℓ kΓ k cωk kb = Rℓ khk c

To complete the Euler equation, we should define the time derivative of the angular

• momentum. Recalling the expression (??), we can write

˙ h

k c =

( d dt ) (Γ k

cωk kb) =

( d dt )

k

(Γ k

cωk kb) + ωk cb × (Γ k cωk cb)

Since Γ k

c is a constant in Rk, it follows that

˙ h

k b = Γ k bαk 0b + ωk 0b × (Γ k bωk 0b)

(32) where αk

0b ≡ ˙

ωk

0b is the total angular acceleration of the body b represented in Rk.

In conclusion, the time derivative of the angular momentum, i.e., the inertial torques that appears in the Euler equation, can be simply written as: ˙ h = Γ ˙ ω + ω × (Γω) (33) where the simplified symbols are to be correctly interpreted.

SLIDE 20

20 Alternative Formulation Another way to obtain (33) is the following: we assume that the body b rotates with a total angular velocity ω0 with respect to an inertial reference frame R0, whose

• rigin is now located at the center-of-mass of the body.

The Euler equation, that describes the dynamics of a rotating body, states that the variation of the angular momentum h0 is equal to the vector τ 0, that is the sum of all torques applied to the body: ˙ h0 ≡ d dt(Γ 0ω0) = τ 0 (34) Both Γ 0 and ω0 are time functions, so we have ˙ h0 = ˙ Γ 0ω0 + Γ 0 ˙ ω0 (35) In order to easily compute ˙ Γ 0 it is convenient to change the reference frame; if Rℓ is the local frame attached to the body, i.e., the body-frame, having its origin in the center-of-mass, if R ≡ R0

ℓ is the rotation matrix representing Rℓ in R0, and if Γ ℓ

is the inertia tensor expressed in Rℓ, then we have seen from (22) that Γ 0 = RΓ ℓRT (36) Calling ωℓ the total angular velocity of the body represented in the frame Rℓ, i.e., ωℓ ≡ ωℓ

0b, and calling ω0 the same velocity represented in R0, i.e., ω0 0b ≡ ω0, we

can write ω0 = Rωℓ and consequently h0 = Γ 0ω0 = RΓ ℓRT Rωℓ = RΓ ℓωℓ = Rhℓ (37) where we have introduced the new vector hℓ = Γ ℓωℓ that is the angular momen- tum represented in Rℓ. We can see that the angular momentum is a vector that transforms from one reference frame to another with the same rule of other vectors. Replacing h0 in (37) and computing the derivative, we obtain: ˙ h0 = d dt(RΓ ℓωℓ) = ˙ RΓ ℓωℓ + R ˙ Γ ℓωℓ + RΓ ℓ ˙ ωℓ (38) The inertia tensor Γ ℓ is time invariant, since it is relative to a body-frame and the body does not change its shape, therefore ˙ Γ ℓ = O, and (38) simplifies as follwos ˙ h0 = S(ω0)RΓ ℓωℓ + RΓ ℓ ˙ ωℓ = ω0 × (RΓ ℓωℓ) + RΓ ℓ ˙ ωℓ (39) This vector is represented in R0 but includes both ωℓ and ω0; if we want to represent everything in Rℓ we should pre-multiply by RT. Recalling the properties of the skew-symmetric matrices, we obtain ˙ hℓ = RT ˙ h0 = RTS(ω0)RΓ ℓωℓ + RTRΓ ℓ ˙ ωℓ = S(RTω0)Γ ℓωℓ + Γ ℓ ˙ ωℓ = S(ωℓ)Γ ℓωℓ + Γ ℓ ˙ ωℓ = ωℓ × Γ ℓωℓ + Γ ℓ ˙ ωℓ (40)

SLIDE 21

0.5. NEWTON-EULER EQUATIONS 21 Poich´ e ˙ hℓ = τ ℓ, segue: τ ℓ = ωℓ × Γ ℓωℓ + Γ ℓ ˙ ωℓ (41) where τ ℓ = RTτ 0 is the sum of the applied torques, represented in Rℓ. If we prefer to represent the torque in R0, we use (39), replacing Γ 0R instead of RΓ ℓ, obtained from (36); in this way we obtain τ 0 = ω0 × Γ 0ω0 + Γ 0 ˙ ω0 (42) The equations (33), (41) and (42) are identical, provided that all vectors and tensor matrices are expressed in the same reference frame.

0.5 Newton-Euler Equations

The Newton-Euler equations are vector equilibrium equations that state that all forces and torques acting on the body B, including the inertial ones due to the motion, the external ones and those due to the constraint, are in (dynamical) equi- librium. There are two vector equations, the first, due to Newton, describes the equilibrium

• f the linear forces, while the second, due mainly to Euler, describes the equilibrium
• f the linear momenta.

Assuming that the rigid body is connected with other rigid bodies, we must consider also the reaction forces and torques that the other bodies transmit to the considered body. For simplicity, from now on, all vectors will be represented in a common reference frame Rk, but omitting the subscript k. The motion of the rigid body B is decomposed into two distinct motions

• 1. A translation motion of the center-of-mass C, described by the Newton

equilibrium vector equation f aC + f vC + mg C − maC = 0 (43) where aC = ˙ v C = ¨ r C is the total acceleration, with respect to an inertial reference frame R0, of the center-of-mass C, g C is the acceleration due to the gravitational field, applied to C, f aC is the resultant of all the external forces acting on the body applied in C, and f vC is the resultant of the constraint forces, with the application point brought in C.

• 2. A rotation motion around a generic point P, described by the Euler equi-

librium vector equation τ aP + τ vP + τ f

aP + τ f vP + τ gP + τ iP − Γ b/Pα0 − ω0 × Γ b/Pω0 = 0

(44)

SLIDE 22

22 where τ aP is the resultant of the active applied external torques, τ vP is the resultant of the constraint torques, τ f

aP is is the resultant of the torques due to

the applied external forces, τ f

vP is the resultant of torques due to the constraint

forces, τ gP is the torque due to the gravitational field, and τ iP is the torque due to the inertial force −maP, ω0 is the total angular velocity ω0 ≡ ω0

0b, and

α0 is the total angular acceleration α0 = ( d dt ) ω0

0b.

Equation (44) is simplified if the point P coincides with the center-of-mass C; indeed in this case we have τ gP = τ iP = 0 and (44)reduces to τ aC + τ vC + τ f

aC + τ f vC − Γ bα0 − ω0 × Γ bω0 = 0

(45) Notice that equations (43) and (45) are functions of aC, the linear acceleration of the center-of-mass, of the total angular acceleration α0 and of the total angular velocity ω0. If the acceleration due to the gravitational field is negligible, then g ≈ 0. The Newton-Euler equations can be also written in a more compact form, introduc- ing suitable matrices, as follows [mI Γ b ] [aC α0 ] + [ S(ω0)Γ b ] = [f totC τ totC ] (46) where f totC = f aC + f vC + mg C τ totC = τ aC + τ vC + τ f

aC + τ f vC

Recalling that aC = ˙ v C + ω0 × v C, where v C is the total linear velocity of the center-of-mass, a second form of the equation is possible [mI Γ b ] [ ˙ v C α0 ] + [mS(ω0) S(ω0)Γ b ] [v C ω0 ] = [f totC τ totC ] (47) If instead of the center-of-mass C we want to express the equation with respect to another point P, the resulting equation become [ mI mS(r pc)T mS(r pc) Γ bω0 ] [ap α0 ] + [mS(ω0)S(ω0)r pc S(ω0)Γ b/pω0 ] = [f totP τ totP ] (48) where r pc is the oriented segment from P to C, and ap = ¨ r p is the total acceleration

• f P.

In conclusion, we want again point out that the difficulties arising with the use of these equations is due to the necessity to include the constraint forces and torques; these do not contribute to characterize the body motion, but have the only function to replace the geometrical constraints with vector quantities to be included in the equations However, if one is interested in numerical computation of the body dynamics the Newton-Euler equations, suitably organized in a recursive form, are an efficient way to do so, as we will see in Section 0.8.

SLIDE 23

0.6. VIRTUAL WORK PRINCIPLE 23

0.6 Virtual Work Principle

Virtual displacements, as seen in Section ??, give the degrees-of-freedom of a rigid body and are used to define the virtual work. Given a system composed by N point masses, each located in a position defined by the vector r i, and N forces (possibly having zero magnitude) acting on them from the external world, the virtual work δW is defined as in relation (??), reported here for completeness δW =

N

∑

i=1

f i · δr i If the system is a continuous rigid body, on which both forces and torques act, we can write a more complete expression, as in (??) δW =

Nf

∑

i=1

f i · δr i +

Nτ

∑

i=1

τ i · δαi = 0 (49) where τ i are the applied torques and δαi the respective angular virtual displace- ments. The virtual work principle can be applied also to dynamic equilibrium; indeed, if we introduce the inertia forces −miai, and if we define with the name of effective forces the sum of the applied and inertia forces F i = f i − miai acting on each mass particle, the D’Alembert principle [?] states that the total virtual work of the effective forces is zero for all invertible variations that satisfy the kinematic conditions given by : ∑

i

F iδ · r i = ∑

i

(f i − miai) · δr i = 0 (50)

0.7 Lagrange Equations in Mechanical Systems

We have already seen in Section 0.5 how the Newton-Euler approach determines the dynamical equations of the system. We will see now how the Lagrange approach can produce an equivalent set of equations, starting from a very different viewpoint. The Lagrange approach is based on the definition of two scalar quantities, namely the total kinetic co-energy 8 and the total potential energy associated with the body. In general the Lagrange method allows to define a set of Lagrange equations, that have some advantages with respect to the vector equations provided by the Newton-Euler approach.

8The reason of using the term co-energy instead of the term energy, will be clarified later.

SLIDE 24

24

• The approach provides n second order scalar differential equations, directly ex-

press in the generalized coordinates ˙ qi(t) e qi(t), that provide an easier starting point for writing the state variable equations than the Newton-Euler equations.

• If holonomic constraint are present, in the equations the constraint force do

not appear: This is very convenient when the equations are used for control design or simulation problems, but may be annoying when we want to explicitly compute constraint forces. We will discuss this topic in Section 0.7.5.

• The kinetic co-energies and the potential energies are independent of the ref-

erence frame used to represent the body motion.

• The kinetic co-energies and the potential energies are additive scalars: in a

multi-body system the total energies/co-energies is the sum of each component energy/co-energy. We will now formally define these two quantities.

0.7.1 Kinetic Energy and Co-energy

First we will define these energies for a point mass and then for a rigid extended body. Point-mass The kinetic energy associated to a point-mass m is defined as the work necessary to increase the linear or angular momentum from 0 to h, i.e., C(h) = ∫ h dW (51) The infinitesimal work associated to the mass is given by dW = f · dr where f is the resultant of the applied forces on the mass and dr is the infinitesimal position increment. The infinitesimal increment dr differs from the virtual displace- ment δr since the former is consistent with the laws of dynamics, i.e., dr = vdt, while the latter if free or restricted only by the geometry of the constraints. In order to have as integrand in (51) only a quantity that depends on the integration limits, we recall that f = dh dt

SLIDE 25

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 25 that relates the force to the variation of the linear momentum; the resulting in- finitesimal work is therefore dW = f · dr = dh dt · dr = dh dt · vdt = v · dh and so we can write C(h) = ∫ h v · dh (52) The kinetic energy is a scalar state function associated to the particle, that de- pends only of the state variables v and h associated to the point mass. It is possible to define another state function associated to the point-mass, that we will call kinetic co-energy, defined as C∗(v) = ∫ v h · dv (53) As shown in Figure 8, between the energy and the co-energy, a relation exists, namely C∗(v) = h · v − C(h) (54) This relation is an example of the Legendre transform. In Section ?? we will present a detailed discussion of this topic. Figure 8: The relation existing between the linear momentum h and the linear velocity v; in this particular case the two quantities are scalars and the relation assumes a constant mass. In particular, if the mass particle is moving at a velocity significantly smaller that the light speed c, i.e., it is not a relativistic mass, we can write h = mv and the two

SLIDE 26

26 “energies” become C(h) = ∫ h 1 mh · dh = 1 2mh · h = 1 2m ∥h∥2 C∗(v) = ∫ v mv · dv = 1 2mv · v = 1 2m ∥v∥2 (55) As one can see, in this case the kinetic energy and co-energy are the same; this does not happens for relativistic particles. We conclude this Section pointing out that in many textbooks the kinetic energy is indicated by the symbol T, but we prefer to use this not-so-common symbol C or C∗, and since these are functions of the generalized coordinates and velocities, often we write C(q, ˙ q) or C∗(q, ˙ q) Multiple Point-masses Bodies Now we take into consideration an extended body composed by N masses mi. The kinetic co-energy will be the sum of the kinetic co-energy of each mass, C∗(v) = 1 2

N

∑

i=1

miv i · v i (56) In Figure 9 we have an inertial reference frame R0, with origin in a fixed point O and a body-frame Rb; if we consider the velocity v i = ˙ r i with respect to R0, we consider first the linear velocity and then the angular velocity. The velocities in R0 can be computed from the general relation (??), rewritten here for convenience ˙ r 0(t) = ω0

01(t) × ρ0(t) + R0 1 ˙

r 1(t) + ˙ d

1(t)

where the second term R0

b ˙

r i(t) at the right hand side is zero, since the point-masses are fixed with respect to the body-frame. Now we consider a a purely translatory motion and then a purely rotational motion. Translational Motion If the motion is purely translational, (??) gives ˙ r 0

i (t) = ˙

d

b(t) ≡ v 0(t)

where v 0 is the total linear velocity with respect to R0; therefore each mass mi will have the same velocity v 0 and the relation (56) will reduce to C∗ = 1 2v 0 · v 0

N

∑

i=1

mi = 1 2mtot v 0 · v 0 = 1 2mtot ∥v 0∥2 = 1 2v T

0 (mtotI )v 0

(57)

SLIDE 27

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 27 Figure 9: A rigid system composed by N point masses, and the various position

• vectors. C is the body center-of-mass.

where the mass mtot is the total body mass. So this equation states that the kinetic (co)energy of a purely translational motion is equivalent to that of one single particle with mass mtot having the total translational velocity v 0. We recall that the total mass mtot can be ideally concentrated in the body center-

• f-mass C del corpo, that obeys to the relation

r cmtot = ∑

i

r imi → r c = ∑

i

mi mtot r i Since the velocity is equal for all points of the body, v 0 is also the velocity of the center-of-mass C; if we use the symbol v 0c ≡ v 0 for this velocity, we can write C∗ = 1 2mtot v 0c · v 0c = 1 2mtot ∥v 0c∥2 = 1 2v T

0c(mtotI )v 0c

that gives the usual rule: “the kinetic energy is half the product of the total mass for total velocity squared”. Rotational Motion If the motion is purely rotational, (??) reduces to ˙ r 0

i (t) = ω0 01(t) × r b i

where ω0

01 ≡ ω0 is the total angular velocity; therefore, using a simpler notation,

the velocity of each mass follows the relation v i = ω0 × r i

SLIDE 28

28 where r i is the position of the i-th mass in R0. Replacing this expression in (56) we obtain C∗ = 1 2

N

∑

i=1

mi(ω0 × r i) · (ω0 × r i) (58) Recalling that, given three generic vectors a, b, c, we can write a · (b × c) = b · (c × a) and calling a ≡ (ω0 × r i), b ≡ ω0, c ≡ r i we obtain C∗ = 1 2

N

∑

i=1

miω0 · r i × (ω0 × r i) = 1 2ω0 · ( N ∑

i=1

mir i × (ω0 × r i) ) Now, recalling (4), we obtain C∗ = 1 2ω0 · ( N ∑

i=1

hi ) = 1 2ω0 · h0 where h0 is the total angular momentum with respect to O. Since h0 = Γ 0ω0, we can conclude writing this expression C∗ = 1 2ωT

0 Γ 0ω0

(59) So we interpret this equation stating that the kinetic (co)energy of purely rotating body is dependent on the total angular velocity and the inertia matrix with respect to the point O. Notice the similarity between (57) and (59): C∗ = 1 2v T

0c(mtotI )v 0c

C∗ = 1 2ωT

0 Γ 0ω0

In both cases the energy is proportional to the “weighted” scalar product of the velocities, where the weighing matrix is the fixed mass diagonal matrix, or the variable inertia symmetric matrix. Alternative Formulation An alternative formulation is obtained considering the center-of-mass C, whose po- sition in R0 is given by r c, and representing the position of each point-mass as r i = r c + ci, where ci is the oriented segment from C to the i-th mass, as in Figure 9.

SLIDE 29

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 29 We compute the velocity v i = ω0 × r i as v i = ω0 × (r c + ci) = ω0 × r c + ω0 × ci = v c + ω0 × ci and, from (56), we have C∗ = 1 2 ∑

i

miv i · v i = 1 2 ∑

i

mi(v c + ω0 × ci) · (v c + ω0 × ci) = 1 2 {∑

i

miv c · v c + ∑

i

mi(ω0 × ci) · (ω0 × ci) + . . . } (60) In this case v c is the (tangential) velocity of the center-of-mass due to the angular motion; considering (58), where now in place of r i we have ci, we can rewrite the above equation as C∗ = 1 2 { v T

c (mtotI )v c + ωT 0 Γ cω0

} (61) where Γ c is the inertia matrix with respect to the body center-of-mass. Natural Systems In general, one can extend the arguments used above to the case in which each mass has a position that depends, in addition to the coordinates qi, also explicitly from the time t; therefore a system consisting of N point-mass particles defined by a position r i(q(t), t) and by the velocity v i(q(t), t) ≡ ˙ r i(q(t), t) =

n

∑

k=1

∂r i ∂qk ˙ qk + ∂r i ∂t , have a kinetic (co)energy that is defined as C∗ = 1 2

N

∑

i=1

mi ˙ r i · ˙ r i = 1 2

N

∑

i=1

[ n ∑

j=i n

∑

k=1

∂r i ∂qj · ∂r i ∂qk ˙ qj ˙ qk + 2∂r i ∂t ·

n

∑

j=i

∂r i ∂qj ˙ qj + ∂r i ∂t · ∂r i ∂t ] (62) If we introduce the following coefficients αjk =

N

∑

i=1

mi ∂r i ∂qj · ∂r i ∂qk βj =

N

∑

i=1

mi ∂r i ∂t · ∂r i ∂qj γ = 1 2

N

∑

i=1

mi ∂r i ∂t · ∂r i ∂t

SLIDE 30

30 we can write (62) as C∗ = C∗

2 + C∗ 1 + C∗

(63) where C∗

2

= 1 2

n

∑

j=1 n

∑

k=1

αjk ˙ qj ˙ qk C∗

1

=

n

∑

j=1

βj ˙ qj C∗ = γ Therefore, in general, the kinetic (co)energy is the sum of three terms: a quadratic, a linear and a constant one. If there is no direct time dependence the linear and the constant terms disappear, and the energy is a purely quadratic function. A system where

• the vectors r i are not explicitly depending on time, but depend on time only

through the generalized coordinates q(t),

• the field of external forces is conservative,
• no constraints are present,

has a kinetic (co)energy that is quadratic, as C = C2 = 1 2

n

∑

j=1 n

∑

k=1

αjk ˙ qj ˙ qk These systems are a called natural systems; the others are called non natural systems. In a non natural system we observe that the term C∗

0 behaves as a sort of potential

energy and is naturally connected to the centrifugal forces9, that are proportional to the square of the angular velocity ∂r i ∂t · ∂r i ∂t , while the term C1 produces the so called Coriolis forces, that are proportional to the product of two velocities Extended Rigid Body In analogy with the derivation used to compute the kinetic (co)energy of a multiple point-masses body, it is possible to compute the kinetic (co)energy of an extended

9Consider, as an example, a satellite orbiting the Earth, where the gravitational (conservative)

field forces are kept in equilibrium by the centrifugal force due to the angular velocity of the satellite around the Earth.

SLIDE 31

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 31 rigid body, as that represented in Figure 10. The relation is C∗(v) = 1 2 ∫

B

v · v dm = 1 2 ∫

B

ρ(r) v · v dV where v is the total velocity of every element with mass dm; the integral is com- puted considering the entire body volume V ; ρ(r) is the local density, function of the position r of the elements; the density is assumed to be uniform in time, but generally not in position. Figure 10: da definire Again we have an inertial reference frame R0, with origin in a fixed point O and a body-frame Rb; with arguments that are similar to those used to derive relation (61), we have now C∗(v, ω) = 1 2mv c · v c + 1 2ωT

0 Γ cω0 = 1

2 [ v T

c (mI )v c + ωT 0 Γ cω0

] (64) where m = ∫

B dm =

∫

B ρ(r)dV is the total body mass, Γ c is the inertia moment

with respect to the center-of-mass C, v c is the total velocity of the center of mass and ω0 is the total angular velocity of the body. From (64) we can state the well known law for computing the kinetic (co)energy of a rigid body The kinetic (co)energy of an extended rigid body is the sum of two terms: a transla- tional part, equal to the kinetic energy equivalent to the entire mass concentrated in the body center-of-mass, and an rotational part, equal to the kinetic energy due to the rotational inertia of the body with respect to its center-of-mass. The reference system with respect to which the inertia matrix Γ c is computed must be considered as a body-frame, with its origin in C; in this way the inertia matrix

SLIDE 32

32 is constant in time. However this frame may be rotate with respect to the frame where ω0 is represented; in this case, as given by (31) and (36), the Γ c value must be expressed in the common frame. Example 0.7.1 Let us consider a prismatic rigid body B (as a rod); we want to compute the kinetic

• energy. As shown in Figure 11 we have two possible reference frames, Ra e Rb.

Figure 11: Computation of the inertia matrix for a rigid body and two reference frames. The first one is inertial and the total angular velocity of the body is ωa = [ ˙ θ ] ; the second is a body-frame with origin in the center-of-mass C. At a given time t the two systems are related by the following rotation matrix Ra

b = Rot(i, 90◦) =

  1 −1 1   The axes of Rb are aligned with the principal inertia axes, and the inertia matrix is Γ b =   Γx Γy Γz   Since the kinetic (co)energy is independent of the reference frame, provided that we use the total angular velocity, we can compute it in two different ways.

• 1. We can represent Γ in Ra and use the total angular velocity ωa represented

in that frame C∗ = 1 2(ωa)TΓ aωa

SLIDE 33

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 33

• 2. We can represent Γ in Rb and use the total angular velocity ωb represented

in that frame C∗ = 1 2(ωb)TΓ bωb In the first case we have Γ a = Ra

bΓ bRb a =

  1 −1 1     Γx Γy Γz     1 1 −1   =   Γx Γz Γy   and consequently C∗ = 1 2(ωa)TΓ aωa = 1 2 [ ˙ θ ]   Γx Γz Γy     ˙ θ   = 1 2Γy ˙ θ2 In the second case we have ωb = Rb

aωa =

  1 1 −1     ˙ θ   =   ˙ θ   and consequently C∗ = 1 2(ωb)TΓ bωb = 1 2 [ ˙ θ ]   Γx Γy Γz     ˙ θ   = 1 2Γy ˙ θ2 In both cases, as expected, the kinetic (co)energy is the same.

• Generalized Coordinates

If we want to write the kinetic (co)energy as a function of the generalized velocities ˙ qi(t), we shall recall the direct kinematic velocity function, that requires to know the linear Jacobian and the geometrical angular Jacobian: v c = J ℓ(q) ˙ q

• r

ω = J ω(q) ˙ q (65) Here we will use the geometrical Jacobian, since we are interested in ω0 and not in ˙ α0 to express C∗(v, ω) in (64) as a function of q and ˙

• q. Using the linear (??) and

the angular geometrical jacobians we obtain C∗( ˙ q, q) = 1 2 [ ˙ q TJ T

ℓ (mI )J ℓ ˙

q + ˙ q TJ T

ωΓ cJ ω ˙

q ] (66)

SLIDE 34

34 and finally C∗( ˙ q, q) = 1 2 ˙ q T [ J T

ℓ (mI )J ℓ + J T ωΓ cJ ω

] ˙ q = 1 2 ˙ q TΓ tot ˙ q (67) This relation shows that the kinetic (co)energy is a quadratic function of the overall inertial mass of the body; since the kinetic energy of a body with a non-zero mass is always positive, the matrix Γ tot is always positive definite and hence invertible.

0.7.2 Potential Energy

Gravitational Energy A force f is said to be conservative or potential when it the negative gradient of a potential function P(r); the potential function is a scalar function and depends

• n the position r =

[ x y z ]T, defined as: f (r) = −∇P(r) = − [∂P(r) ∂x ∂P(r) ∂y ∂P(r) ∂z ]T (68) The symbol ∇ indicates the gradient operator that transforms a scalar function into a vector10. If a potential function exists, it is conventionally called potential energy of the system; it must be observed that, if it exists, it is unique apart from an additive

• constant. This implies, as we will see later, that the effects on the body dynamics

depend only from the potential energy variation, and not on its absolute value. The dynamical systems in which all the forces doing work are due to a potential function, are called conservative systems. In a conservative force field the potential energy P is obtained integrating the con- servative forces, i.e., P(r) = − ∫ r

r 0

f (r) · dr = ∫ r

r0

∇P(r) · dr (69) where r 0 represents a generic initial position of the mass, and the minus sign estab- lish that the potential energy decreases if the force field perform the work on the

• body. A typical example of conservative force filed is the gravitational field around

Earth, or in general any mass. The potential field produces the so-called weight forces. The potential energy due to a gravitational field and associated to a generic mass m is given by the following relation: Pg = −mg · r c0 (70)

10Often the gradient is interpreted as a row vector, but we prefer to follow the usual column

representation rule, so we introduce the transpose symbol T.

SLIDE 35

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 35 where g is the local gravitational acceleration vector and r c0 is the body center-of- mass position vector with respect to an orthogonal plane to g that has the conven- tional zero value of potential energy (zero energy plane); it is the variation of energy with respect to this plane that is considered. Elastic Energy Moreover there is another force field related to potential energy; it is that due to elastic elements, as the ideal springs: these elements are the abstract model of a proportional relation between elongation and force, so that force is proportional to the deformation of the ideal spring element. If we assume a one-dimensional linear spring we can write f = kee (71) while if we assume a one-dimensional torsional or torsion spring we write τ = k′

eδ

(72) where e and δ are, respectively, the linear elongation, and the angular defor- mation from the rest position of the spring; ke and k′

e are the so-called elastic

constants of the springs. The elongation and the deformation physically represent a difference in position or angle between a rest state of the elastic element and the deformed state of the same element. In one dimensional linear springs the force and the deformation are co-linear, while in one dimensional torsion spring the torque and the rotation axis are co-linear as well. Figure 13 depicts a one dimensional linear spring, where the force and the deforma- tion are co-linear. It is also possible to find elastic elements where the elastic constant depends on the deformation itself, as in f = ke(e)e τ = k′

e(δ)δ

(73) Figure 12: The scheme of a one dimensional linear spring. The potential energy is the integral of the virtual work performed by the spring deformation P(e) = ∫ e f · de (74)

SLIDE 36

36 Figure 13: The scheme of a one dimensional torsion spring.

• r

P(δ) = ∫ δ τ · dδ (75) We can also in this case define the potential co-energies, that are P∗(f ) = ∫ f e · df

• r

P∗(τ) = ∫ τ δ · dτ The relation between P∗(f ) and P(e) is similar to the co-energy Legendre transfor- mation given in (54) P∗(f ) = f · e − P(e) (76) as illustrated in Figure 14 Figure 14: da definire When the elastic elements are linear with constant ke, the potential energy and co-energy are as follows P(e) = 1 2eT(keI )e = 1 2ke ∥e∥2 (77)

SLIDE 37

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 37 and P∗(f ) = 1 2f T(keI )−1f = 1 2ke ∥f ∥2 (78) If When torsion springs are considered with constant k′

e, the potential energy and

co-energy are as P(δ) = 1 2δT(k′

eI )δ = 1

2ke ∥δ∥2 (79) and P∗(τ) = 1 2τ T(k′

eI )−1τ =

1 2k′

e

∥τ∥2 (80) When the linear relations (71) and (72) are valid, the energies and co-energies are equal P(e) = P∗(f ) and P(δ) = P∗(τ) If the elastic constants are different along the three directions, as, for example, in the case of the rubber-band ball of Figure 15, we have the following, more general, relations: P(e) = 1 2eTK ee; P∗(f ) = 1 2f TK −1

e f

and P(δ) = 1 2δTK ′

eδ;

P∗(τ) = 1 2τ T(K ′

e)−1τ

where K e = diag(ke1, ke2, ke3) e K ′

e = diag(k′ e1, k′ e2, k′ e3) are the elastic constant

matrices along the three dimensional components. Figure 15: A rubber-band sphere, an example of a three dimensional elastic element. We observe that, in general, the potential energy is function of the position coordi- nates alone, so we will write P(q)11. Nevertheless, there are some cases where the potential energy is also dependent

• n the velocity coordinates ˙

q; in this case it will be designated as generalized potential energy: the most notable example is that involving the motion of a charged particle in an electromagnetic field, where we can write P(q, ˙ q) = e ( ϕ(r) − 1 cv(r, ˙ r) · A ) where e is the electric charge of the particle, ϕ is the scalar electric potential, A is the vector potential of the electromagnetic field, v is the particle speed and c is the light velocity.

11Many textbooks write the symbol V (q) for potential energy, but we prefer to use the symbol

P.

SLIDE 38

38

0.7.3 Generalized forces in Holonomic Systems

Let us assume that we have a system composed by N point-masses, subject only to holonomic, ideal, frictionless constraints. The forces f i acting on the i-th mass can be classified according to three groups:

• N ′′′ forces f c

i due to constraint reactions.

• N ′′ forces f c

i due to a conservative field.

• N ′ non conservative forces f nc

i .

i.e., f i =

N′

∑

i=1

f nc

i + N′′

∑

i=1

f c

i + N′′′

∑

i=1

f c

i

Constraint Forces The allowed virtual displacements δr i are always tangent to the constraints, while the constraint forces f v

i are always orthogonal to te constraints, as illustrated in

Figure 16; from this assumption it follows that f c

i · dr i = 0

and f c

i · δr i = 0

Therefore the work done by the constraint forces is zero (the forces “do not work”) δWv =

N′′′

∑

i=1

f c

i · δr i = 0.

Figure 16: Vincoli geometrici, spostamenti virtuali e forze vincolari.

SLIDE 39

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 39 Conservative Forces The N ′′ conservative forces do a work that, considering (68), results to be δWc =

N′′

∑

i=1

f c

i · δr i = − N′′

∑

i=1

∇Pi · δr i = −

N′′

∑

i=1

∇Pi · [ n ∑

j=1

∂r i ∂qj δqj ] =

n

∑

j=1

[ N′′ ∑

i=1

−∇Pi · ∂r i ∂qj ]

• Fc

j

δqj (81) This last expression highlights the so called generalized conservative forces Fc

j and

allows to transform the virtual work from a function of the positions r to a function

• f the generalized coordinates qj:

δWc =

N′′

∑

i=1

f c

i · δr i = n

∑

j=1

Fc

j δqj = F c · δq

(82) Non Conservative Forces The N ′ non conservative forces f nc

i

do a work equal to δWnc =

N′

∑

i=1

f nc

i · δr i = N′

∑

i=1

f nc

i ·

[ n ∑

j=1

∂r i ∂qj δqj ] =

n

∑

j=1

[ N′ ∑

i=1

f nc

i · ∂r i

∂qj ]

• Fnc

j

δqj (83) This last expression highlights the so called generalized non conservative forces Fnc

j

and allows to transform the virtual work from a function of the positions r to a function of the generalized coordinates qj: δWnc =

N′′

∑

i=1

f nc

i · δr i = n

∑

j=1

Fnc

j δqj = F nc · δq

(84) In conclusion, only two types of forces will do work in a system subject to homoge- neous constraint

• the j-th generalized conservative force:

Fc

j (q) = − N′′

∑

i=1

∇Pi · ∂r i ∂qj

SLIDE 40

40

• the j-th generalized non conservative force:

Fnc

j (q) = N′

∑

i=1

f nc

i · ∂r i

∂qj The generalized force, being the result of a scalar product, will be itself a scalar quantity. In case of torques acting on the system, the generalized forces due to them will give

• rigin to the following generalized torques

T c

j = N′′

∑

i=1

−∇Pi · ∂αi ∂qj ; T nc

j

=

N′

∑

i=1

τ nc

i · ∂αi

∂qj (85) In the following, the symbol used to identify both the generalized forces and the generalized torques will be identical, namely F.

0.7.4 Lagrange Equations with Holonomic Constraints

In a multi-body system subject to holonomic constraints, the formulation of the Lagrange equations may take different forms, from the most general one to that describing particular cases; we will see them all briefly in the following Sections We start with the knowledge of the total co-energy of the system, that is the sum of the co-energies of the N composing parts, as a function of the generalized coordinates and velocities C∗(q, ˙ q) =

N

∑

k=1

C∗

k(q, ˙

q) The Lagrange equations are a set of n equations (one for each generalized coordinates qi) defined as d dt (∂C∗ ∂ ˙ qi ) − ∂C∗ ∂qi = Fi i = 1, . . . , n (86) where Fi = Fc

i + Fnc i

is the i-th generalized force, with a positive sign if applied by the environment to the body, with a negative sign if applied by the body to the environment. Indeed we have seen in the previous Sections that the generalized forces are determined by the virtual work principle: the work can be done by the external environment on the body, or from the body on the external environment; in the two cases the sign will be opposite. It is a convention that the positive sign is that of the work done on the body. Considering relation (81), we can write equations (86) in a slightly different way, decomposing Fi into the two terms Fc

i and Fnc i , and bringing portiamo Fc i on the

left hand side d dt (∂C∗ ∂ ˙ qi ) − ∂C∗ ∂qi − Fc

i = Fnc i

SLIDE 41

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 41 Since (81) is valid, we can write d dt (∂C∗ ∂ ˙ qi ) − ∂C∗ ∂qi + [ N′ ∑

k=1

∇Pk · ∂r k ∂qi ] = Fnc

i

we notice that the term inside the square bracket is equal to the partial derivative ∂P ∂qi ; therefore we can write d dt (∂C∗ ∂ ˙ qi ) − (∂C∗ ∂qi − ∂P ∂qi ) = Fnc

i

Again we notice that ∂P ∂ ˙ qi = 0 so we can write d dt (∂C∗ ∂ ˙ qi − ∂P ∂ ˙ qi ) − (∂C∗ ∂qi − ∂P ∂qi ) = Fnc

i

i = 1, . . . , n (87) The Lagrange state function (or simply the Lagrange function) is defined as the difference between the kinetic co-energy C∗ of the system and its potential energy P: L(q, ˙ q) = C∗(q, ˙ q) − P(q) (88) Considering (89), we can therefore write n differential equations d dt (∂L(q, ˙ q) ∂ ˙ qi ) − ∂L(q, ˙ q) ∂qi = Fnc

i (q)

i = 1, . . . , n (89) each one relative to the i-th generalized coordinate. The term ∂L ∂ ˙ qi is called generalized momentum and is usually indicated by the symbol µi; the vector of generalized momenta is therefore indicated by µ(q(t)). Dissipative and Friction Forces If we prefer to keep the friction or other dissipative forces f fric

i

separate from the

• ther non conservative forces, considering that the friction represents an energy

dissipated by the body in the form of heat, and therefore appears with the minus sign, the Lagrange equation becomes: d dt (∂L ∂ ˙ qi ) − ∂L ∂qi = Fi − f fric

i

(90)

SLIDE 42

42 In general the friction forces are due to complex physical phenomena involving the interaction between solid/solid or solid/fluid surfaces; tribology is the science that studies in details the friction forces. To keep the discussion simple, we can schematically describe the friction force f fric

i

as a nonlinear function of the relative velocity v between the two contact surfaces

• f the involved bodies.

This nonlinear function is approximately due to three distinct phenomena that takes place at different velocity regimes: when the velocity is zero or very low, the friction is highly nonlinear and is called stiction or Stribeck friction or Stribeck effect, from the name of the German engineer that studied it at the beginning of the XX century. When the velocity increases, the friction is essentially proportional to the velocity; it is the so-called viscous regime. Moreover there is also a constant component that is independent of the sliding velocity, called Coulomb friction from the name of the French scientist that studied friction phenomena at the end of XVIII

• century. At nonzero velocities friction introduces a damping effect on the system

motion. In conclusion, we can model the total friction force as in Figure 17, and write f fric

total = f fric stiction + f fric coulomb + f fric viscous

(91) Figure 17: The three components of the non-conservative friction forces While stiction and Coulomb friction must be explicitly introduced as non-conservative forces, it is a common assumption, as suggested by the English physicist Rayleigh,

SLIDE 43

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 43 to express the viscous dissipative phenomenon as the derivative of a dissipation function, also called Rayleigh function, given by: Di( ˙ q) = 1 2 ˙ q T(βiI ) ˙ q = 1 2βi ∥ ˙ q∥2 (92) where the coefficient βi is the viscous friction coefficient, and ∥ ˙ q∥ is the norm of the relative velocity between the moving body and the surface responsible of the viscous friction effect. Care should be taken not to call this quadratic expression dissipation “energy” since it is only a conventional way to introduce it in the Lagrange, reformulating (90) as follows d dt (∂L ∂ ˙ qi ) − ∂L ∂qi + ∂Di ∂ ˙ qi = Fi (93) Now the term Fi includes only the non conservative forces, with the exclusion of the forces modelled by the dissipation function, that now are given by ∂Di ∂ ˙ qi = βi ∥ ˙ q∥ In particular the forces Fi include the motion forces/torques applied by the actuators and the interaction forces exchanged with the environment. The n equations (93) can be written in matrix form, with a slight abuse of notation, as d dt (∂L ∂ ˙ q ) − ∂L ∂q + ∂D ∂ ˙ q = F (94) To conclude this paragraph, we point out that in the scientific literature the holo- nomic systems in which the forces are solely due to generalized potential functions P(q, ˙ q), are called Lagrangian or Lagrangean systems; when the kinetic co- energy and the potential energy explicitly depend on time, i.e., C∗ = C∗(q, ˙ q, t) and P = P(q, ˙ q, t), the systems are called Hamiltonian systems.

0.7.5 Lagrange Equations with Non Holonomic Constraints

The Lagrangian approach has the advantage to produce the dynamical equations as function only of the generalized coordinates. In this way only the minimum number

• f differential equations, namely n equations, is generated, and this is due to the

fact that the constraint forces do not generate work in correspondence to the virtual displacements compatible with the geometrical constraints. The constraint forces never appear in the equations: the holonomic systems are described by a set of independent generalized coordinates, with no constraints. Conversely, systems where non holonomic constraints are present, cannot be com- pletely described by the independent generalized coordinates alone. The equations

SLIDE 44

44 must be augmented including also the constraint characterization; as a result we

• bserve the appearance of the constraint forces.

If the problem requires it, we can make appear the constraint forces also in holonomic systems; to do this, we consider that:

• 1. the geometric constraints are representable as hyper-surfaces in the χj config-

uration space;

• 2. the geometric constraints are enforced by suitable forces/torques called, as

already known, constraint forces;

• 3. at every instant the constraint forces are orthogonal to the hyper-surfaces

described by the constraints. The holonomic constraints are “translated” into generalized constraint forces, that are added to the non-constraint forces already included into the equations. Un- fortunately these constraint forces cannot be a-priori computed, since they depend

• n the actual motion and, in order to obtain them, we must in advance solve the

differential equations. However we can use a simple method to obtain them, that consists in introducing new artificial coordinates, that represent the constraint violation: they are zero when the constraints are satisfied, otherwise they are positive or negative, according to a suitable sign convention related to how the constraints are violated. These coordinates are called superfluous coordinates and the generalized forces associated with them are exactly the constraint forces. If the set of the generalized coordinates together with the superfluous coordinates is independent, then the resulting dynamical equations will contain the constraint forces; these will appear only in the equations related to the superfluous coordinates. The constraint forces can be computed solving the motion equations and setting to a constant value the superfluous coordinates. Conversely, when non-holonomic constraints are involved, the solution require to introduce the so called Lagrange multipliers. Lagrange Multipliers Let us assume a system, defined by n generalized coordinates qi, constrained by m non-holonomic constraints, where the k-th constraint is written in pfaffian form as a1kdq1 + a2kdq2 + · · · + ankdqn + a0k, k = 1, . . . , m The constraint equation expressed in virtual displacements, recalling that δt = 0, becomes a1kδq1 + a2kδq2 + · · · + ankδqn = 0

SLIDE 45

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 45 written in vector form as ak · δq = 0 (95) where ak = [ a1k a2k · · · ank ]T. This relation states that the virtual displace- ments are orthogonal to the vector ak, that is therefore aligned with the constraint forces, also orthogonal to the constraints, as in point 3) above. Therefore the con- straint forces f c

k and the vector ak are proportional, as in

f c

k = λk(t)ak

and we call the proportionality coefficient the k-th Lagrange multiplier. Since δW =

N

∑

i=1

Fiδqi =

N

∑

i=1

(Fnc

i

+ Fc) δqi =

N

∑

i=1

Fnc

i δqi + N

∑

i=1 m

∑

k=1

λiδqi the total generalized force relative to the coordinate qi will be now obtained as the sum of the external generalized and constraint forces: Fnc

i

+

m

∑

k=1

λiaik The i-th Lagrange equation will then become d dt (∂L ∂ ˙ qi ) − ∂L ∂qi = Fnc

i

+

m

∑

k=1

λiaik (96) Equation (95) together with (96), gives origin to n+m equations in n+m unknowns, among which there are the m Lagrange multipliers λk(t). When these equations are solved and the λk are obtained, one use f c

k = λkak

to compute the constraint forces. The matrix form of the Lagrange equation becomes d dt (∂L ∂ ˙ q ) − ∂L ∂q = F nc + Aλ (97)

• r else, separating friction effects from other non-constraint forces,

d dt (∂L ∂ ˙ q ) − ∂L ∂q + ∂D ∂ ˙ q = F + Aλ (98) where A = [ a1 · · · ak · · · am ] ∈ Rn×m; λ =    λ1 . . . λm    ∈ Rm×1

SLIDE 46

46 This equation, together with the non-holonomic constraint equations AT ˙ q + a0 = 0 (99) form a system of n+m equations in n+m unknowns (q, λ), that must be solved si- multaneously to get the n generalized coordinates q and the m Lagrange multipliers λ. After that, it will be possible to compute the m constraint forces as F c = Aλ. Notice thta this method can still be used to find, if necessary or advisable, the constraint forces associate to the mh holonomic constraints, expressed as ϕk(q, t) = 0, k = 1, . . . , mh building the matrix A and the vector a0 in (99) as follows: A =       ∂ϕ1 ∂q1 . . . ∂ϕmh ∂q1 . . . . . . . . . ∂ϕ1 ∂qn . . . ∂ϕmh ∂qn       and a0 =      ∂ϕ1 ∂t . . . ∂ϕmh ∂t      Example 0.7.2 A small sphere with mass m is constrained to move without friction along an ideal wire that rotates with at a constant velocity ˙ θ around the z axis of a given reference frame; the sphere does not rotate. The wire has a parabolic shape, given by the relation z = kx2 (see Figure 18). Figure 18: Example of a mobile constraint.

SLIDE 47

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 47 Using the cartesian coordinates12, we define two reference frames, one fixed R0 and

• ne mobile with the wire R1, related by the rotation matrix

R0

1 = Rot(k, θ) =

  cθ −sθ sθ cθ 1   The constraint are given by the relation        constraint 1: z(t) = kx2(t) in frame R1 constraint 2: ω(t) =   ˙ θ = cost   in both frames Both are holonomic and this implies we can write the virtual displacements as { δz − 2kxδx = 0 δθ = ωδt = 0 The two constraints imply two Lagrange multipliers, λ1 and λ2; the generalized forces associated to the constraints are:

• 1. a force f c that keeps the sphere on the wire,
• 2. a torque τ c, aligned along the axis k 0, that keeps the constant rotation of the

wire. There are two degrees of freedom and two generalized coordinates: we can choose q = (x, θ)

• r

q ′ = (z, θ) Since we want to compute f c and τ c in R0, and we know that they are aligned as the columns of ak, as in (95), it is more convenient to write the constraint in matrix form as functions of the cartesian coordinates ξ = (x, y, z, α1, α2, θ) and not of the generalized coordinates q; with αi we have denoted two generic orientation angles that are not relevant to the problem, since the only orientation we are interested in is θ, i.e., that relative to the rotation around k 0. In forma matriciale abbiamo perci`

• :

ATδξ = [ aT

1

aT

2

] δξ = 0 ⇒ [ −2kx 1 1 ]           δx δy δz δα1 δα2 δθ           = 0 (100)

12In this case the use of the cylindrical coordinates would have been more appropriate, but, in

• rder to make the development straightforward, we have use the cartesian ones.

SLIDE 48

48 We write the Lagrange equations (97), functions of q, after having computed the different quantities necessary to build L:

• mass coordinates in R1:

p1(t) =   x(t) z(t)   =   x(t) kx2(t)  

• mass coordinates in R0:

p0(t) = R0

1p1(t) =

   x(t)cθ x(t)sθ z(t)    =    x(t)cθ x(t)sθ kx2(t)   

• mass velocity in R0, expressed in the generalized coordinates q:

˙ p0(t) =    ˙ x(t)cθ − x(t)sθ ˙ θ ˙ x(t)sθ + x(t)cθ ˙ θ 2kx(t) ˙ x(t)    =    cθ −x(t)sθ sθ x(t)cθ 2kx(t)    [ ˙ x(t) ˙ θ ] = J(q(t)) ˙ q(t)

• kinetic co-energy of the mass

C∗(q, ˙ q) = 1 2m

• ˙

p0 2 = 1 2m ˙ q T(t)J TJ ˙ q(t) = 1 2m ˙ q T(t) [1 + 4k2x2(t) x2(t) ] ˙ q(t) = 1 2m [ ˙ x2(t)(1 + 4k2x2(t)) + x2(t) ˙ θ2(t) ] (101)

• Potential energy of the mass

P(q) = −mgz(t) = −mgkx2(t)

• Lagrange function

L = C∗ − P = 1 2m ( 4k2x2(t) ˙ x2(t) + x2(t) ˙ θ2(t) ) + mgkx2(t)

• generalized forces: no generalized forces are applied, but only constraint forces,

that will be computed once the Lagrange equations (97) are solved

SLIDE 49

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 49

• constraint forces: there are two constraints and two generalized constraint

forces, see (100) f c

1 = λ1a1 = λ1

[ −2kx 1 ]T f c

2 = λ2a2 = λ2

[ 1 ]T Equazione 1) relative to q1 = x(t): d dt ( ∂L ∂ ˙ q1 ) − ∂L ∂q1 = F1 + ∑6

k=1 λ1a1k

d dt (∂L ∂ ˙ x ) − ∂L ∂x = −2kxλ1 d dt [1 2m(2 ˙ x + 8k2x2 ˙ x) ] − [1 2m(8k2x ˙ x2 + 2x ˙ θ2) ] − 2mgkx = −2kxλ1 hence m¨ x(1 + 4k2x2) + 4mk2x ˙ x2 − x ˙ θ2 = 2kx(mg − λ1) (102) Equazione 2) relative to q2 = θ(t): d dt ( ∂L ∂ ˙ q2 ) − ∂L ∂q2 = F2 + ∑6

k=1 λ2a2k

d dt (1 22mx2 ˙ θ ) = λ2 m ( x2¨ θ + 2x ˙ x ˙ θ ) = λ2 (103) since ˙ θ = cost, it follows that ¨ θ = 0, and hence 2mx ˙ x ˙ θ = λ2 Solving these two nonlinear differential equations we obtain also the values for λi; as a matter of fact λ2 can be directly obtained from the last equation.

• 0.7.6

Characterization of the Resulting Equations

Having defined the set of generalized coordinates q(t) ∈ Rn and of the generalized velocities ˙ q(t) ∈ Rn, assuming all holonomic constraints, we have seen that the Lagrange approach generates n differential equations d dt (∂L ∂ ˙ qi ) − ∂L ∂qi = Fi i = 1, . . . , n

SLIDE 50

50 The term d dt (∂L ∂ ˙ qi ) will generate in the equations terms that are proportional to the second time derivatives ¨ q(t), and similarly there will be terms proportional to the first time derivatives ˙ q(t) and to the generalized coordinates q(t). In conclusion, collecting the n equations in one vector equation d dt (∂L ∂ ˙ q ) − ∂L ∂q + ∂D ∂ ˙ q = F (104) this can be linear or non linear: if the equations are linear or if we consider small perturbations around some equilibrium point, we will have a linear formulation that in the general case can be expressed as a second order differential vector equation A1¨ q(t) + A2 ˙ q(t) + A3q(t) = F (105) where F is the generalized force vector. The three matrices Ai ∈ Rn×n can have different physical interpretations in relation to the treated problem; some problems will produce matrices whose elements are function of the generalized coordinates or velocities, in other cases the matrices will be constant. Gravitational terms, when present, will appear on the left hand side as a column vector g(q(t)) added to A3q(t). In alternative, we can include these terms in the generalized forces F. Recalling from (??) that a generic A ∈ Rn×n may be factored as A = As + Aa where As is a symmetric matrix As = AT

s and Aa is a skew-symmetric matrix

Aa = −AT

a .

Therefore equation (105) can be rewritten as M ¨ q(t) + (D + G) ˙ q(t) + (K + H )q(t) = F (106) where M = M T mass or inertia matrix, positive definite symmetric D = DT viscous damping matrix, symmetric G = −GT gyroscopic matrix, skew-symmetric K = K T stiffness (elasticity) matrix, symmetric H = −H T circulatory matrix (due to constrained damping), skew-symmetric (107) Some of these matrices can be function of the generalized coordinates or velocities, as in the industrial robot kinematic chains, where we have configuration-dependent mass matrices M = M (q(t)). Not all the above terms are always present in the dynamical equations of multibody systems; the most common cases are:

SLIDE 51

0.7. LAGRANGE EQUATIONS IN MECHANICAL SYSTEMS 51

• 1. Systems where the structure of (106) is

M ¨ q(t) + Kq(t) = F (108) whre M e K are symmetrical and positive definite, are called conservative

• systems. They are called “conservative” despite the fact that not all the

forces F are due to conservative phenomena; indeed in the left hand side of the equation there are no forces associated to dissipating energy elements, and consequently the total energy is conserved inside the system.

• 2. Systems where

M ¨ q(t) + G ˙ q(t) + Kq(t) = F (109) They are still conservative, in the sense explained above (the internal energy is conserved), but they are properly called gyroscopic conservative sys- tems or undamped gyroscopic systems. This type of dynamical behavior is associated to systems where there is an important rotation motion, as in gyroscopes or satellites.

• 3. Systems where

M ¨ q(t) + D ˙ q(t) + Kq(t) = F (110) where M , D and K are, other than symmetric are also positive definite, are called damped non gyroscopic systems and sometimes damped con- servative systems, also if the internal energy is not conserved, since it is dissipated through the viscous friction. Systems with symmetric positive def- inite matrices are sometimes called passive systems; we will see in Section ?? a more detailed definition of this class of systems.

• 4. Systems where

M ¨ q(t) + (K + H )q(t) = F (111) are called circulatory systems. The circulatory forces modeled by the cir- culatory matrix H are not easily described; they appear in the dynamics of rotating mechanical or electrical machines (rotodynamics). The circulatory matrix contains non conservative terms linked with the internal damping of rotating elements and with the damping of the fluid surrounding the rotor [?], and can cause instability.

• 5. Systems with the general form

M ¨ q(t) + D ˙ q(t) + (K + H )q(t) = F are often the result of modeling rotating shafts and structures, where, beside the usual dissipative phenomena, an additional constraint damping func- tion exists.

SLIDE 52

52 In kinematics chains composed by n rigid links, as for example the industrial robots, the nonlinear Lagrange equation takes the following form: M (q(t))¨ q(t) + C(q(t), ˙ q(t)) ˙ q(t) + g(q(t)) = F (112) where M (q(t)) is the mass or inertia matrix of the robot, depending on the instantaneous configuration q(t), C(q(t), ˙ q(t)) ˙ q(t) represents the Coriolis torques and g(q(t)) represents the gravitational torques.

0.8 Numerical Algorithms for Direct and Inverse Dynamics

The Newton-Euler equations, presented in Section 0.5, are vector equations that must include the constraint forces/torques between the rigid bodies of the multi- body system; when the number of these bodies increases, the equations become more and more hard to solve by hand. Conversely, the Lagrange equations are simpler to write, since they require only the knowledge of the positions and velocities of the bodies, and furthermore, they emphasize some energy aspects that are useful for the design of a special class of controllers, the so-called passive controllers, i.e., those that preserve the passivity property of the controlled system (see Section ??). However, from a numerical point of view, the Newton-Euler equations, when written in a recursive form, are the basis for numerical algorithms with a computational complexity that is lower than that of the algorithms based on the Lagrange approach. To show an example of this recursive formulation, we refer to a robotic structure, that will be discussed without presenting the details that allow to write its final equation. The model consists of n equations; for each i-th robot arm we write the following differential equation, whose general form was already presented in (112))

n

∑

j=1

Mij(q)¨ qj+

n

∑

j=1 n

∑

k=1

cijk(q) ˙ qj ˙ qk + gi(q) = τi i = 1, . . . , n (113) where

• qi(t), ˙

qi(t) and ¨ qi(t) are the generalized coordinates, velocities and accelera- tions;

• τi are the motor torques applied to the joints;
• the terms Mij, cijk and gi model the inertial, centrifugal and Coriolis torques,

as well the gravitational ones. These terms depend on the configuration given by q.

SLIDE 53

0.8. NUMERICAL ALGORITHMS FOR DIRECT AND INVERSE DYNAMICS53 Given the input torques τi(t) it is possible to integrate the differential equations (113) and find a solution for the generalized coordinates qi(t); this problem is the so-called direct dynamics problem. The inverse dynamics problem is that of finding the torques τi(t) from the generalized coordinates qi(t) and their derivatives ˙ qi(t) and ¨ qi(t). The direct dynamics is useful to simulate the behavior of the system given the input torques, while the inverse dynamics allows to compute the control torques to achieve the desired qi(t), ˙ qi(t) and ¨ qi(t). The number of operations required to compute these quantities increase rapidly with

• n. It is possible to show that the number of products is proportional to n4, that is,

the complexity is O(n4). For a particular robot where n = 6, the Lagrange approach requires at every step the computation of 66.271 multiplications. There exists a number of algorithms that reduce the computational complexity of the problem; the most cited is due to Luh, Walker and Paul, that is detailed in [?]. It consists in a two-steps recursive procedure, where, using the kinematics equations, the positions, linear and angular velocities and accelerations of the center-of-mass

• f every joint are computed, starting from the base to the most external one.

At this point, taking into consideration also the generalized forces applied by the environment to the various parts of the structure and the gravitational effects, a backward recursion in implemented, that allows to compute the joint torques. In particular, indicating the following quantities

• ωℓ

k the angular velocity, αℓ k the angular acceleration and aℓ k the linear acceler-

ation of the center-of-mass of the k-th arm in the reference frame Rℓ attached to the ℓ-th arm;

• Rk

k−1 is the rotation matrix from Rk−1 to Rk;

• d ℓ

k−1,k is the translation from the origin of Rk−1 to the origin of Rk, represented

in Rℓ;

• d ℓ

k−1,kc is the translation from the origin of Rk−1 to the center-of-mass of the

k-t arm, represented in Rℓ;

• f ℓ

k−1,k are the constraint forces applied from the (k − 1)-th arm to the k-th

arm, represented in Rℓ;

• τ ℓ

k−1,k are the constraint torques applied from the (k − 1)-th arm to the k-th

arm, represented in Rℓ;

• F ℓ

k,c is the total equivalent force applied to the k-th arm, having its line of

action across the center-of-mass of the k-th arm, represented in Rℓ;

SLIDE 54

54

• nℓ

kc is the total torque on the k-th arm, that is equal to the momentum with

respect to the center-of-mass of the k-th arm, represented in Rℓ;

• Γ ℓ

k/c is the inertia matrix of the k-th arm with respect to the center-of-mass

• f the k-th arm, represented in Rℓ;
• mk is the total mass of the k-th arm;
• k =

[ 1 ]T; Recalling that d k

k,kc = d k k−1,kc − d k k−1,k

and that d k

k,kc, d k k−1,kc, d k k−1,k are constant in Rk, we have this recursive algorithm,

valid for a robot with all rotating joint: Initialization We set ω0

0 = 0, α0 0 = 0 e a0 0 = −g 0, where g 0 is the gravitational acceleration.

Forward recursion (k = 1, . . . , 6) ωk

k = Rk k−1(ωk−1 k−1 + ˙

θkk) (114) αk

k = Rk k−1(αk−1 k−1 + S(ωk−1 k−1) ˙

θkk + ¨ θkk) (115) ak

k = Rk k−1ak−1 k−1 + S(αk k)d k k−1,k + S(ωk k)S(ωk k)d k k−1,k

(116) Motion equations relative to the center-of-mass (k = 1, . . . , 6) F k

kc = mk[ak k + S(αk k)d k k,kcS(ωk k)S(ωk k)d k k,kc]

(117) nk

kc = Γ k k/cαk k + S(ωk k)Γ k k/cωk k

(118) Backward recursion (k = 5, . . . , 0) f k

k−1,k = Rk+1 k

f k+1

k,k+1 + F k kc

(119) τ k

k−1,k = Rk+1 k

τ k+1

k,k+1 + S(d k k−1,k)Rk+1 k

f k+1

k,k+1

(120) + S(d k

k−1,kc)F k kc + nk kc

(121) Motor torque computation τk = k TRk−1

k

τ k

k−1,k

This algorithm requires 117n − 24 products and 103n − 21 sums; for a robot with n = 6 arms the total is 678 products and 597 sums, largely smaller than those required with the Lagrangian approach.

SLIDE 55

Bibliography

[1] S.H. Crandall, D.C.Karnopp, jr. E.F Kurtz, and D.C. Pridmore-Brown. Dynam- ics of mechanical and Electromechanical Systems. McGraw-Hill, 1968. [2] J. Meisel. Principles of Electromechanical-Energy Conversion. Robert E. Krieger Publishing Company, 1984. 55