Dynamics - IITK - Indian Institute of Technology Kanpur

Document technical information

Format pdf
Size 823.7 kB
First found May 22, 2018

Document content analysis

Category Also themed
Language
English
Type
not defined
Concepts
no text concepts found

Persons

Organizations

Places

Transcript

Intermediate Dynamics
in about 100 pages
Anindya Chatterjee
Department of Mechanical Engineering
Indian Institute of Technology
Kanpur 208016
This version: November 4, 2014
2
Preface
These notes started with a course at the Indian Institute of Science (Bangalore) on the dynamics and control
of mechanical systems, of which I taught the “dynamics” part several times. I have since moved to IIT Kanpur
via IIT Kharagpur.
There are many excellent but fat books on dynamics. My notes, at least, are brief.
I hope you will find my treatment of rotation matrices easy and useful (“rotate the body, not the axes”).
You might enjoy my discussion of nonholonomic constraints in Lagrangian mechanics (“where do the λ’s come
from?”). If you like something else, please send me email.
Many students initially struggle with simple dynamics problems. It is like riding a bicycle. It looks easy
until you try, and later you wonder what the fuss was about; but a journey lies in between.
Several excellent people tried to teach me mechanics. Progress was slow. Professor A. C. Gomes of St.
Xavier’s College, Kolkata, taught me well in my impressionable teens (1983-85). Professor Andy Ruina of
Cornell University, my last formal teacher, helped a lot (1993-96). I took a long time to understand what I
describe in these notes. But then dynamics is a great subject developed by great people. I am just a messenger,
trying to bring the story to you.
Several students, friends, and colleagues have helped with various parts of these notes (typing, figures,
general discussion). I am embarrassed to say that I did not write down their names at the time, and so may
miss someone now. People I remember who helped are Sai Jagan Mohan, Venkatesh, Pradeep Mahadevan,
Pradeep Gudla, Pankaj Wahi and Sovan Das. If you helped and your name is not here, please email me and
I will add it with an apology.
I know there are errors in these notes. Please email me if you find some.
Kanpur
[email protected]
November 4, 2014
3
4
Contents
1 Preliminaries
1.1 Vectors . . . . . . . . . .
1.2 Dynamics of a point mass
1.3 A system of point masses
1.4 The laws of dynamics . .
1.5 Exercises . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
7
7
8
12
13
16
2 Relative Motion
2.1 A note on reference frames . . . . . . .
2.2 Derivative of a vector in a moving frame
2.3 The rolling cone . . . . . . . . . . . . .
2.4 Addition of angular velocities . . . . . .
2.5 Acceleration in a moving frame . . . . .
2.6 Exercises . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
17
17
17
19
22
23
27
. . . .
. . . .
. . . .
. . . .
body
. . . .
. . . .
. . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
31
31
32
33
34
36
37
37
38
4 Rotations and Euler Angles
4.1 Rotations do not commute . . . . . . . . . . . . . . . . . .
4.2 Some basic mathematical facts . . . . . . . . . . . . . . .
4.3 Euler’s theorem and the rotation matrix . . . . . . . . . .
4.4 Obtaining the rotation matrix . . . . . . . . . . . . . . . .
4.5 Successive rotations . . . . . . . . . . . . . . . . . . . . .
4.5.1 Infinitesimal Rotations . . . . . . . . . . . . . . . .
4.5.2 Angular velocity . . . . . . . . . . . . . . . . . . .
4.5.3 The integral of angular velocity is not meaningful .
4.5.4 Rotations upto second order in infinitesimals . . .
4.5.5 Angular velocity is not an exact derivative . . . . .
4.6 Euler angles . . . . . . . . . . . . . . . . . . . . . . . . . .
4.7 Angular velocity, again . . . . . . . . . . . . . . . . . . . .
4.8 Alternative form for the rotation matrix . . . . . . . . . .
4.9 The moment of inertia matrix Icm . . . . . . . . . . . . .
4.10 Angular accelerations . . . . . . . . . . . . . . . . . . . .
4.11 Rotated coordinate systems . . . . . . . . . . . . . . . . .
4.12 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
41
41
41
42
43
44
44
44
44
45
45
46
46
46
47
48
49
49
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3 Momentum Balance for Rigid Bodies
3.1 Linear momentum balance . . . . . . .
3.2 Angular momentum . . . . . . . . . .
3.3 Matrices and vectors . . . . . . . . . .
3.4 Angular momentum of a rigid body . .
3.5 Angular momentum balance for a rigid
3.6 Planar or 2D problems . . . . . . . . .
3.7 The kinetic energy of a rigid body . .
3.8 Exercises . . . . . . . . . . . . . . . .
5
CONTENTS
6
5 Lagrange’s Equations
5.1 Virtual work for a system of particles . . . . . . .
5.1.1 Applied versus constraint forces . . . . . .
5.1.2 Net zero virtual work of constraint forces
5.2 Generalized coordinates and degrees of freedom .
5.3 Constraints . . . . . . . . . . . . . . . . . . . . .
5.4 Holonomic systems . . . . . . . . . . . . . . . . .
5.5 Nonholonomic systems . . . . . . . . . . . . . . .
5.6 The calculus of variations . . . . . . . . . . . . .
5.7 Hamilton’s principle . . . . . . . . . . . . . . . .
5.8 Exercises . . . . . . . . . . . . . . . . . . . . . .
6 The
6.1
6.2
6.3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
51
51
51
52
53
53
55
60
62
64
65
Rolling Coin
67
General comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
Lagrange’s equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
Newton-Euler equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
7 Linear Vibrations
7.1 A single degree of freedom system . . . . . . . .
7.2 Normal modes . . . . . . . . . . . . . . . . . . .
7.3 The generalized eigenvalue problem for vibrations
7.4 Forced vibrations with damping . . . . . . . . . .
7.5 Approximations via Lagrange’s equations . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
71
71
72
73
75
75
8 Some Problems Involving Single Bodies
8.1 Cylinder rolling down an incline . . . . .
8.2 Hemispherical shell on a table . . . . . .
8.3 Torque-free rigid body . . . . . . . . . .
8.3.1 Stability of pure spin . . . . . . .
8.3.2 Axisymmetric bodies . . . . . . .
8.4 Wobbling disk using Euler angles . . . .
8.5 Rigid body on an axle . . . . . . . . . .
8.6 Exercises . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
79
79
80
81
82
83
84
86
87
9 Friction
9.1 Coulomb friction . . . . . . . . . . .
9.2 A spring block system . . . . . . . .
9.3 On distinguishing between μs and μk
9.4 Vibration damping . . . . . . . . . .
9.5 Resonance . . . . . . . . . . . . . . .
9.6 Sliding on moving surfaces . . . . . .
9.7 Rolling friction . . . . . . . . . . . .
9.8 Planar sliding . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
89
. 89
. 90
. 94
. 94
. 95
. 98
. 103
. 103
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Chapter 1
Preliminaries
The material we will cover is at a beginning postgraduate level. However, there are many preliminary topics
that are covered well in some undergraduate courses, and not well in others. Some of these essential topics
are presented briefly in this chapter.
1.1
Vectors
Vectors and tensors can be thought of as physical quantities that exist independently of whether or not we
use them in any way. But they are important in the teaching of dynamics for two reasons. They help to clear
away clutter (on the blackboard and in our minds) due to the compactness with which they can express the
physical ideas we work with. And they have a transparent correspondence with the matrix calculations that
lie at the heart of much scientific work with computers. We will set ourselves the minimal goal of using vectors
and matrices (in place of tensors), which suffices for the material covered here.
Notation: We denote vectors with lowercase letters, either boldface or underlined. For example,
a = a = ax î + ay ĵ + az k̂,
where ax , ay , and az are components along the unit vectors î, ĵ, and k̂ respectively. The lowercase letters
with hats represent unit vectors.
k
j
O
i
Figure 1.1: Unit vectors define a coordinate system.
The dot product of two vectors a and b is given by
a · b = ax bx + ay by + az bz .
7
CHAPTER 1. PRELIMINARIES
8
The cross product is given by
î
a × b = ax
bx
ĵ
ay
by
k̂
az
bz
.
This is equal to (ay bz − az by )î − (ax bz − az bx )ĵ + (ax by − ay bx )k̂.
1.2
Dynamics of a point mass
The reader may have some intuitive idea about what is meant by force and by mass. We will not define these
quantities. The reader may also have some idea about what a frame of reference is, but we will provide a
working definition here. A frame of reference may be thought of as a video camera that records motions of
objects. If the camera itself moves around and rotates, it may give an altered view of the motions of the
objects under study. When we say an object moves in a certain way, we imply that this is the motion as
recorded by some specific video camera. Two video cameras that do not move relative to each other record
identical measurements of the motion of any object (i.e., their observations agree). In textbooks, the “video
camera” is sometimes called an observer. Some authors say a frame of reference is a rigid body. Some say it
is a set of three orthonormal unit vectors.
The idea of a point mass is used in mechanics in a circular but self-consistent way. We take a point mass to
mean any body whose overall translation alone is of interest. Rotations are either unimportant or negligible
in the present context. (The analysis of rotational motions, which is necessary to decide whether rotations are
negligible or not, is developed using the idea of a system of point masses: circular, but self-consistent.)
If F is the net vector force on a particle of mass m then its acceleration a in an inertial frame of reference
satisfies
F = ma.
(1.1)
This is true for any force and any point mass.
This leaves us with the question of what an inertial frame is. Unless we know that we have an inertial
frame to begin with, we cannot be sure that Eq. 1.1 applies. Given a frame, we can experimentally check if
Eq. 1.1 is satisfied to acceptable accuracy. If it is, then the frame is inertial for the purpose at hand. For
many calculations, such as involving the responses of cars to braking forces, or the trajectories of tennis balls
across distances of several metres, the surface of the earth is acceptable as an inertial frame. However, for
monsoon winds which change direction as they cross the equator, the earth is a rotating frame (non-inertial).
One incontrovertible point is that any frame which moves with constant velocity and no rotation with respect
to an inertial frame is also itself an inertial frame. Strictly speaking, therefore, Eq. 1.1 actually provides a
definition of an inertial frame, and not a description of the behaviour of masses acted upon by forces.
These philosophical issues apart, many problems involving single point masses can be solved in a straightforward fashion, on assuming that we have an inertial frame. The answers obtained are usually very accurate
and useful. Two examples follow.
Notation: The position vector of point A with respect to point B, or from B to A, is written as rA/B . By
itself, rA represents rA/O , where O is a point fixed in some frame of reference that should be clear from the
context. Usually, O will be the origin of some coordinate system fixed to an inertial reference frame, such as
stationary ground.
Problem: A point mass m is suspended from three strings A, B and C as shown in figure 1.2. The system
accelerates upwards with an acceleration a. There is gravity. Find the tensions in the strings. It is given that
1.2. DYNAMICS OF A POINT MASS
9
(in some units of length)
rP/O
=
2 î + 2 ĵ + 3 k̂,
rQ/O
=
3 k̂,
rR/O
=
4 ĵ + 3 k̂,
rS/O
= î + 2 ĵ + 0.5 k̂.
Z
Q
string
P
string A
B
g
R
string C
O
X
s
Y
Figure 1.2: A point mass suspended using strings accelerates upwards.
FB
FA
FC
mass = m
mg
Figure 1.3: Free body diagram of the mass.
Solution: We begin with a free body diagram, which is a sketch of the system (here, a point mass) showing
all the external forces (here, from the strings and from gravity) and moments (here, there are none). See figure
1.3.
CHAPTER 1. PRELIMINARIES
10
Various needed position vectors may be found as
rP/S = rP/O − rS/O ,
etc. By Eq. 1.1,
The unit vectors along FA , FB
F = ma = +mg − [FA + FB + FC ].
rZ/S
rP/S rQ/S
,
and
respectively. We obtain
and FC are
|rP/S | |rQ/S |
|rZ/S |
ma k̂ = −mg k̂ + [0.3714FA − 0.2981FB − 0.2981FC ] î + [−0.5963FB − 0.5963FC ] ĵ
(1.2)
(1.3)
+[0.9285FA + 0.7454FB + 0.7454FC ] k̂.
From the above vector equation, we can extract three scalar equations by considering individual components
along î, ĵ and k̂. These may be solved for FA , FB and FC . If m = 1 kg, a = 1 m/s2 and g = 10 m/s2 , for
example, we find FA = 5.9229 N, FB = FC = 3.6897 N.
Problem: See figure 1.4. A bead of mass m slides on a helical coil of diameter D and pitch h. Gravity acts
downwards. Find the bead’s (i) downward acceleration and (ii) steady state speed, respectively, assuming (i)
the coil to be frictionless and (ii) that Coulomb friction acts between the coil and the bead.
^
k
h
^j
^i
m
D
Figure 1.4: Bead sliding down a coil.
Solution:
Case (i) No friction. We seek the downward acceleration of the bead.
In this case, there is no energy dissipation. We can, for this simple problem, state and use energy conservation1 as follows: kinetic energy + potential energy = total energy = constant.
Let the speed of the bead along the wire be v, and its height at any instant be z. Then
1
mv 2 + mgz = constant.
2
Differentiating,
mv v̇ + mg ż = 0.
(1.4)
From figure 1.5,
ż = −v sin α, so z̈ = −v̇ sin α.
1 We include an example of “energy methods” here, and will formalize them in the form of a law by the end of this chapter.
But we will not use them much in the multi-degree-of-freedom problems that follow.
1.2. DYNAMICS OF A POINT MASS
11
motion of bead
h
α
πD
Figure 1.5: The triangle shows one unwrapped turn of the coil, relating motion along the wire with change in
height of the bead.
Thus, from Eq. 1.4,
v̇ = g sin α.
Note that the bead also has a centripetal acceleration, but we do not seek that here. Here we want z̈, which
from the above is
gh2
z̈ = −g sin2 α = − 2
.
h + π 2 D2
Case (ii) With Coulomb friction. We seek the bead’s steady speed (assuming there is one).
At steady sliding speed, the bead also has steady downward speed. In an inertial reference frame that
moves vertically down with this same speed, the bead executes uniform circular motion. In that frame (and
therefore in any inertial frame), the particle has only radially inward acceleration.
^k
m
Fw
v cos α
α
^et
^er (out of page)
et
(a)
(b)
mg
(c)
Figure 1.6: (a) Coordinate system. (b) Horizontal component of velocity. (c) Free body diagram.
We now define some unit vectors shown in figure 1.6(a). Unit vector êt is along the wire in the direction
of motion of the bead; k̂ is vertical; and êr is perpendicular to k̂ and êt (i.e., out of the page at the instant
shown).
The horizontal component of the velocity has magnitude v cos α, as shown in figure 1.6(b). The (centripetal)
acceleration of the bead at steady speed is therefore
a=
where R = D/2.
−v 2 cos2 α
êr ,
R
(1.5)
CHAPTER 1. PRELIMINARIES
12
The free body diagram of the bead is shown in figure 1.6(c). There is a vector force Fw from the wire, in
addition to the weight acting downwards. By Eq. 1.1,
mv 2 cos2 α
êr .
R
The force from the wire itself may be split into a component along the wire and one perpendicular to it, as in
Fw − mg k̂ = −
Fw = FN − Ft êt ,
where we have used boldface for FN to indicate its vector nature, and no boldface in Ft to indicate that it is
a scalar magnitude; note also that the force along the wire opposes the motion. We then have
FN − Ft êt − mg k̂ = −
mv 2 cos2 α
êr .
R
(1.6)
Taking dot products on both sides with êt , and noting that FN · êt = 0, êr · êt = 0, and k̂ · êt = − sin α, we get
Ft = mg sin α.
From Eq. 1.6 we then obtain
|FN | = mg 1 − sin α
2
2
2
+ (mg sin α cos α) +
mv 2 cos2 α
R
2 12
.
We have not so far used Coulomb friction; the results so far apply for any friction law that causes a steady
sliding speed.
Now, Coulomb friction implies
(1.7)
Ft = μ|Fn |.
Substituting the expressions obtained above,
2 12
2
mv 2 cos2 α
2
2
mg sin α = μ mg 1 − sin α + (mg sin α cos α) +
.
R
The steady speed v is then solved for as (substituting for sin α as well)
√
1 1
g 2
√
h + 4π 2 R2 4 h2 − μ2 4π 2 R2 4 .
v=
2π Rμ
1.3
(1.8)
(1.9)
A system of point masses
We now consider a system of several (say, N ) point masses. The following definitions and their consequences
are standard. The center of mass of a system of particles, rcm , is given by
N
mi ri
,
(1.10)
rcm = i=1
mtot
where ri is the position vector of the ith point mass and the total mass
mtot =
N
mi .
i=1
The velocity and acceleration of the center of mass follow by differentiation:
N
i=1 mi vi
,
vcm =
mtot
N
acm
=
i=1
mi a i
mtot
,
where vi and ai are the velocity and acceleration, respectively, of the ith point mass.
1.4. THE LAWS OF DYNAMICS
1.4
13
The laws of dynamics
A free body diagram can be drawn for any system.
Free body diagrams are so useful in mechanics that I like to state this as a law (after my teacher Andy Ruina).
A “system” here includes a collection of one or more parts, of one or more objects, which may each be
solid or fluid or a mixture thereof, stationary or in motion, in equilibrium or accelerating. All external forces
and moments acting on the parts of the system must appear in the diagram; this includes forces from the rest
of the objects at the imaginary cut sections where some parts have been removed or “freed” for inclusion in
the free body diagram. Internal forces, that act between parts of the system, must not be shown in the free
body diagram.
A few examples of systems and their free body diagrams are given in figure 1.7 (a) through (c). In each
(a)
R2
Wrod
(b)
R
M1 1
Wpart 1
M2
Wpart 2
M
(c)
V
R
Figure 1.7: (a) A bomb exploding under water. (b) A rod resting on two blocks. (c) A cantilever beam with
roller supports and a distributed load.
case, the system is shown on the left and a free body diagram is shown on the right.
In figure 1.7(a), a bomb explodes under water. We can draw a free body diagram of a portion of water
including the bomb. The bomb is breaking up inside the water, and solids, liquids and gases are all included
in the system: this is allowed. As a modeling choice, we may choose to ignore the effect of gravity while
analyzing this explosion; accordingly, weights are not shown in the free body diagram. The forces between
the pieces of the bomb casing and the water, for example, are internal forces and not shown in the free body
diagram. The forces acting from the exterior fluids on to the surface of the portion of water included in the
free body diagram are shown. We may choose to model these forces as purely normal, as suggested in the
diagram.
In figure 1.7(b), we consider a rod resting on two blocks. We decide to draw a free body diagram for a
CHAPTER 1. PRELIMINARIES
14
system consisting of portions of the two blocks as well as the rod, as shown using a dashed line. The weight
of the rod appears in the free body diagram, acting through the center of mass of the rod. The weights of
the portions of the blocks considered in the system appear as well, acting through the respective centers of
mass of these portions. The net forces and moments on the portions of the blocks included, from the portions
removed, are also included. The points at which these forces act must be chosen and fixed at this time. For
any such choice of points, there will be a set of forces and moments that represent the net effect of contact
with the cut-off portions, and changing these points will generally change the moments as well.
In figure 1.7(c), we consider a cantilever beam with roller supports and a distributed load. A part of the
beam is selected, as shown, and a free body diagram drawn for it. The distributed load acting on the beam
to the left of the cut section does not appear in the free body diagram; and that acting to the right does
appear. The reaction force from the first roller does not appear, and that from the second roller does. The
roller reaction is shown as vertical, in recognition of the physical nature of the roller contact.
The information that the roller contact force should be vertical cannot be deduced from a free body
diagram without some physical assumptions or knowledge about the contact. For example, if we draw a free
body diagram of the roller alone (not presented here), and put vector forces and moments at the top and
bottom contact points, then the roller is physically not distinguished from, say, a welded or built in support.
The key assumption here needs to be that there are no contact moments. Incorporating that, we find from
static equilibrium conditions that the contact forces must themselves be vertical. This information can then
be fed back to the free body diagram of the portion of the beam under consideration.
The above examples of free body diagrams illustrate basic principles. But practice is needed even for such
apparently simple tasks as drawing free body diagrams (like swimming, or riding a bicycle). Inexperienced
readers are encouraged to consult a good textbook on undergraduate level engineering mechanics.
Once a correct free body diagram has been drawn, application of the laws of momentum balance (which
follow below) is routine. In this sense, the solution of dynamics problems using this approach stands or falls
based on whether a useful system has been identified for purposes of drawing a free body diagram, and whether
that free body diagram has been drawn correctly.
Linear momentum balance
The sum of all forces appearing in a free body diagram equals the total mass of the objects in the diagram
times the acceleration of the center of mass.
(1.11)
Fext = mtot acm .
Examples of using this law for the special case of single point masses have been given above. We will
address more complex problems later.
Angular momentum balance
The net moment (due to the action of all the forces and moments in the free body diagram) about any point
C in space is
M/C =
ri/C × mi ai .
(1.12)
We have not considered problems involving angular momentum balance yet, but will.
Energy balance
For any system with total energy E and net power input Ẇin ,
dE
− Ẇin = 0.
dt
The above can also be integrated over some time interval to get
(1.13)
ΔE − ΔWin = 0.
A simple example of using this law, for the special case of Ẇin = 0, has been considered above. This
law, though powerful, works best for systems with one degree of freedom and will not be used much in what
follows.
1.4. THE LAWS OF DYNAMICS
15
Problems with energy balance
We mention that though energy conservation seems like a clear enough principle, its application can have
pitfalls for the unwary. In an academic setting such pitfalls are more entertaining than dangerous, and so I
give below an example that I like.
Consider a uniform chain of length L and total mass m lying coiled on the ground. One end of the chain
is held and moved upwards with a steady velocity v. The net gain in energy of the chain at the instant when
its other end finally leaves the ground equals the sum of its kinetic and potential energies, i.e.,
Net energy gain of chain =
1
1
mv 2 + mgL.
2
2
(1.14)
Now let us analyze the mechanics of lifting the chain. Assume the ground to be frictionless. A free body
diagram of the partially lifted chain is shown in figure 1.8. At the instant shown, a length h of the chain has
F
mg
h
R
Figure 1.8: Partially lifted chain.
been lifted off the ground. The forces shown in the free body diagram are the pulling force F , the net weight
mg, and the ground reaction R. We are interested only in the vertical motions, and so the horizontal positions
of the forces (their lines of action) are not important.
It is clear that the center of mass of the chain is at a height
hh
2
L2 = h .
m
2L
m
ycm =
The acceleration of the center of mass is then
acm =
v2
ḣ2 + hḧ
= ,
L
L
because ḣ = v = constant. From Eq. 1.1, we are sure that
F + R − mg = m
v2
.
L
(1.15)
But what is R? Some people would assume, as we do here, that
R = mg
i.e., the weight of the chain still on the ground.
L−h
,
L
(1.16)
CHAPTER 1. PRELIMINARIES
16
This assumption seems reasonable, though it is conceivable that (for example) internal longitudinal stress
waves might be set up in the lifted portion of the chain, which might then change R. We accept this assumption
here, as many might.
Accepting Eq. 1.16, we find from Eq. 1.15
L−h
v2
+ mg − mg
.
L
L
The net work done by the lifting force is now easily calculated as
L
1
F dh = mv 2 + mgL,
2
0
F =m
which is more than the energy gained by the chain (see Eq. 1.14).
There is no real trick in the foregoing calculation. Under the conditions assumed in this problem energy
simply cannot be conserved, in the sense that some dissipation or transformation into other forms must occur.
Equation 1.13 merely says that the difference in energy must have gone into other forms, which might be heat,
sound, or possibly even residual longitudinal vibrations in the chain.
What if our assumption about R was wrong? Suppose R was a little greater than we assumed, because the
ground was pushing up on a chain a little more than we thought? Then the force F would be smaller, and less
energy would be wasted. Would the ground then be doing work on the chain to offset dissipation in the chain,
or would there simply be less dissipation in the chain to start with? Equation 1.13 cannot answer this question.
Better understanding of the microscopic material behavior of the chain in the infinitesimal region where it
instantaneously experiences infinitely large accelerations may provide an answer. Other “explanations” are
possible. But that is not the point. The point is that Equation 1.13, though correct, can be tricky to use on
unfamiliar problems.
1.5
Exercises
1. Let a = î + ĵ + k̂. It is given that b · k̂ = 1, that b · a = 2, and that (a × b) · î = 3. Find b.
Answer: 3î − 2ĵ + k̂.
2. Find a unit vector perpendicular to both î + 2.2ĵ and 1.2î − 0.9ĵ + 0.4k̂.
Answer: ± (0.2398î − 0.1090ĵ − 0.9647k̂).
3. A truck of mass 12,000 kg, traveling at 10 m/s, is brought to a halt using a constant braking force F over
a time duration of 2 seconds. What is F ?
Answer: 60,000 N.
4. A system consists of 4 point masses, each of mass 1 kg. They are moving at 1 m/s each, in the directions
North, North-East, East, and East respectively. What is the velocity of the center of mass?
Answer: Letting î be along East and ĵ along North, v cm = (2.707î + 1.707ĵ)/4.
5. A chain of 5 blocks, each of mass m, moves frictionlessly along a straight horizontal line. Successive masses
are connected to each other using 4 springs of constants k12 , k23 , k34 and k45 . At some instant, the extensions
in these 4 springs are x12 , x23 , x34 and x45 . A single external force F acts on the second mass in the direction
of its motion; no external forces act on the other masses. What is the acceleration of the center of mass of
this system?
Answer: F/4m.
6. A car moves along a straight line. The center of mass of the car is 0.7 m above the ground, and 1.6 m
ahead of the rear wheel. The distance between the front and rear wheels is 2.2 m. The front wheel brakes have
failed, but the rear brakes are working. The coefficient of friction between the road and tires is 0.8. What is
the maximum achievable deceleration? Ignore the inertia of the wheels, and assume the front wheels rotate
freely.
Chapter 2
Relative Motion
2.1
A note on reference frames
Reference frames have been discussed earlier. A set of three mutually perpendicular unit vectors will be called
a coordinate system. Coordinate systems are used in a reference frame in order to quantify positions and
velocities of objects. Any coordinate system can be chosen to describe the motion of an object as seen in
any reference frame. Thus, the concept of reference frame is fundamentally different from that of coordinate
system. However, if we choose and rigidly fix a unique coordinate system to each reference frame, then each
such coordinate system is equivalent to its corresponding reference frame.
2.2
Derivative of a vector in a moving frame
Let a frame xyz rotate with angular velocity ω with respect to a fixed (stationary) frame XY Z. Let a vector
A be fixed in xyz. Then A is not fixed in XY Z, and its rate of change as seen from (or by) that frame is
dA
=ω × A
(2.1)
dt XY Z
To see why this is true, we must first define angular velocity.
Here we will use an intuitive definition, which will later be given in a form better suited to actual calculation.
In an infinitesimal time period Δt, it is intuitively clear that reference frame xyz (like any other rigid body)
undergoes an infinitesimal rotation. That rotation involves some infinitesimal angle (say Δθ), and occurs
about some unit vector n̂. We define
Δθ
n̂.
(2.2)
ω = lim
Δt→0 Δt
See figure 2.1(a). The initial position of A is marked by points OP, and the infinitesimally rotated position
by points OQ. If the rotation Δθ was continued further about the same axis, the tip of A would describe a
circle whose plane is indicated by grey shading in the figure. The projection of A onto this plane lies along a
unit vector which we denote by λ̂. A third unit vector ê is then defined as n̂ × λ̂.
Let the angle between A and n̂ be φ. See figure 2.1(b). L is the length of arc between points P and Q of
figure 2.1(a), and is given by
L = |A| sin φ Δθ,
and so the infinitesimal change in A is
ΔA = |A| sin φ Δθ ê.
(2.3)
It is clear now that the time derivative of A as seen in frame XY Z is in the direction of ê, and given by
dA
Δθ
= |A| sin φ lim
ê
(2.4)
Δt→0 Δt
dt XY Z
17
CHAPTER 2. RELATIVE MOTION
18
Z
e
P
Δθ
Δθ
Q
As
e^
L
in φ
^λ
n^
A
_
O
Y
X
(a)
(b)
Figure 2.1: (a) Vector A is fixed in frame xyz which rotates by an angle Δθ in time Δt about the unit vector
n̂. (b) Projection of vector A on to the plane normal to n̂.
Since ê = n̂ × λ̂, we have
dA
dt
= lim
XY Z
Δt→0
Δθ
n̂ × |A| sin φλ̂ = ω × A,
Δt
because the component of A along n̂ does not contribute to the cross product. This proves Eq. 2.1.
We now consider the case when A is not fixed in xyz, even as xyz rotates relative to XY Z. In this case
the total infinitesimal change of A during Δt is the sum of (i) the change of A as occurring in xyz, and (ii)
the change in the infinitesimally changed A due to the rotation of xyz.
An example may help here. Imagine a clock that faces east. The second hand points to 3 at some instant
of time. Over the next 5 seconds, the clock rotates about a vertical axis so that it ends up facing south east.
During these same 5 seconds, the second hand moves from 3 to 4. What is the net change in the “vector”
represented by the second hand? It is clearly the same as we would get if we let the clock be stationary
for 5 seconds (during which interval the second hand keeps moving), and then very rapidly rotate it to face
south east (during which interval the second hand has no time to move). For this latter case, however, we see
that the net change is simply the vector sum of two changes: (i) the change in the second hand from 3 to 4
considering the clock to be fixed in the initial configuration (facing east), and (ii) the change in the second
hand obtained by holding it fixed, pointing to 4, as we rotate the clock to face south east. This way of splitting
the net change into two vector parts is correct even if the two changes are non-infinitesimal.
So, in an infinitesimal time interval Δt, we think of the change in A as occurring in two parts. The first
part occurs smoothly over time Δt, during which xyz does not rotate. This part of the net change is then
dA
Δt.
ΔA1 =
dt xyz
The vector A has now changed to A + ΔA1 . The angle φ has changed to, say, φ + Δφ. The new tangential
unit vector, in place of ê, is now the slightly different unit vector ê + Δê. The second part of the change in A
is now calculated. Time stands still, and the vector A + ΔA1 is instantaneously rotated (along with xyz) by
Δθ about n̂. The second part of the change is (by Eq. 2.3)
:
ΔA2 = |A + ΔA1 | sin(φ + Δφ) Δθ (ê + Δê) = |A| sin φ Δθ ê + second order
infinitesimals.
0
2.3. THE ROLLING CONE
19
Now, having computed the change in A over a time interval of Δt, we can write
dA
dA
ΔA1 + ΔA2
=
= lim
+ ω × A.
dt XY Z Δt→0
Δt
dt xyz
(2.5)
Equation 2.5 will be used often in what follows, and the reader is encouraged to spend some time thinking it
through.
2.3
The rolling cone
Let us begin with a problem and see what is needed to solve it.
Problem: See figure 2.2. A right circular cone of base circle radius R and slant height L rolls without slip
on the ground. The line of contact with the ground rotates about the vertical at an angular rate Ω. Find the
angular velocity and angular acceleration of the cone at the instant when the line of contact coincides with
the Y axis.
Z
Ω
R
n
Y
O
L
X
Figure 2.2: A cone of base radius R and slant height L rolls without slipping on the ground. The line of
contact rotates with an angular velocity Ω about the Z axis.
Solution: To solve this problem, we will consider more than one rotating reference frame.
The angular velocity of a body B relative to (or with respect to, or as seen by) some frame F1 is the vector
sum of the angular velocity of the body relative to any other frame F2 and the angular velocity of F2 itself
relative to F1 . We write this statement as
ω B/F1 = ω B/F2 + ω F2 /F1 ,
(2.6)
use it for now, and will prove it later.
Choose a frame xyz rotating about the Z axis as shown in figure 2.3(a), such that the y axis is always
along the instantaneous line of contact. At the instant of interest, the y and Y axes coincide. From figure
2.3(b), we see that
sin β =
R
.
L
CHAPTER 2. RELATIVE MOTION
20
Z
Ω
z
R
n
O
y
O’
L
Y
P
n
h
x
R
β
X
L
(a)
(b)
Figure 2.3: (a) Frame XY Z is stationary and frame xyz rotates about the Z axis with angular rate Ω. The
origins of axes attached to the two frames, O and O , are coincident but shown slightly displaced for clarity.
(b) Projection of the cone on the Y Z plane.
What is the angular velocity of the cone relative to xyz? Imagine sitting on a turntable at the origin,
always looking out along the line of contact: this is how the video camera of frame xyz sees the cone. So,
relative to the xyz frame the cone simply rotates about n̂, the unit vector along the axis of symmetry of the
cone:
ω cone/xyz
= λ n̂,
where λ is an as yet undetermined magnitude. From Eq. 2.6,
ω cone/xyz
= λ n̂ + Ω k̂
We can now proceed along one of a few different directions.
Method 1 (intuitive):
Many people enjoy intuitive approaches. Unfortunately, intuition is mysterious, and cannot be reliably
strengthened by systematic methods. A typical complex engineering problem will be solved with a disciplined
and systematic attack, not with flashes of intuition interspersed with indefinite periods of waiting. This may
be one reason why engineering is called not a subject but a discipline. So, although we appreciate intuition
when we find it, we must deemphasize it in teaching. Here, we allow an exception.
When the cone revolves once around the Z axis, the distance covered on the ground by the tip of the line
of contact is 2πL. The cone rolls without slip. The observer of frame xyz sees the cone spinning as the ground
rotates. So, during the same time interval, the cone spins through m revolutions such that
2πRm = 2πL,
i.e., m = L/R. Therefore,
L
Ω.
R
In the above the minus sign comes from visualization, and represents a trap for the careless.
λ=−
2.3. THE ROLLING CONE
21
From Eq. 2.6,
ω cone/XY Z
= −
From figure 2.3(b),
√
n̂ = cos β ĵ + sin β k̂ =
L
Ωn̂ + Ωk̂
R
L2 − R2
R
ĵ + k̂.
L
L
This gives
√
L2 − R2
Ω ĵ.
ω cone/XY Z = −
R
The negative sign says the angular velocity vector points towards O.
(2.7)
Method 2 (also intuitive):
Since the cone rolls without slip, all the points on the line of contact have zero velocity. Then by our
previous definition of angular velocity, the line of contact is the instantaneous axis of rotation of the cone; and
ω cone/XY Z = ω ĵ.
ω cone/XY Z · k̂ = λn̂ + Ωk̂ · k̂ = 0,
This means
which gives
L
Ω
R
as above, with the minus sign coming without need for visualization.
λ=−
Method 3:
Consider the material point P (shown in figure 2.3(a)) at the tip of the line of contact. Here by “material
point” P we mean the point on the surface of the cone, which leaves the ground as the cone rolls.
The position vector of P is rP/O . At the instant of interest, P has zero velocity in frame XY Z (no slip).
Note that rP/O is rigidly fixed to the cone. Hence, using Eq. 2.5,
drP/O
dt
*
drP/O × rP/O +
,
dt
cone
XY Z
= 0 = ω cone/XY Z
0
which gives
ω cone/XY Z × ĵ = 0,
which in turn means ω cone/XY Z is along ĵ, giving
ω cone/XY Z · k̂ = 0
as before.
We have now calculated the angular velocity of the cone in three different ways. Of these, the third way
is least intuitive, most systematic, and most recommended.
To calculate the angular acceleration of the cone let us start with
ω cone/XY Z = ω cone/xyz + ω xyz/XY Z ,
where both terms on the right hand side are now known.
(2.8)
CHAPTER 2. RELATIVE MOTION
22
We write the angular acceleration of the cone as
d
ω
αcone/XY Z =
dt cone/XY Z XY Z
:0
d
d
ω cone/xyz
ω
=
+
Z
dt
dt xyz/XY
XY Z
XY Z
where the second term on the right hand side is zero because the line of contact rotates with constant Ω.
Moreover,
d ω cone/xyz
dt
XY Z
=
:0
d + ω
ω
xyz/XY Z × ω cone/xyz
dtcone/xyz xyz
=
Ωk̂ × λn̂
=
(Ωλ) k̂ × n̂,
where the first term on the right hand side of the first equation is zero because Ω, and therefore λ (as defined
above), is constant. Substituting for λ and n̂,
αcone/XY Z =
2.4
√
Ω2 L2 − R2
î.
R
Addition of angular velocities
See figure 2.4. A rigid body B rotates with an angular velocity ω B/F2 relative to frame F2 . Frame F2 itself
has an angular velocity ω F2 /F1 relative to frame F1 . Equation 2.6 claims
ω B/F1
= ω B/F2 + ω F2 /F1 .
We now prove it.
Z
F2
ωB/ F
z
2
F1
P
O
y
x
Y
G
X
Figure 2.4: A rigid body B with an angular velocity specified relative to frame F2 . Frame F2 itself has an
angular velocity relative to frame F1 . P is an arbitrary point in B.
2.5. ACCELERATION IN A MOVING FRAME
23
Let P be an arbitrary material point in body B. Let the point O in the body be fixed at the origin of xyz.
Then
drP/O
dt
*0
drP/O
=
+ ω B/F1 × rP/O
B
dt
F1
where the first term on the right hand side is zero since P is fixed in B. The same quantity can be calculated
using the frame F2 as well. Then
drP/O
drP/O
=
+ ω F2 /F1 × rP/O ,
dt
dt
F1
F2
*0
drP/O
drP/O
=
+ ω B/F2 × rP/O ,
where in turn
B
dt
F2
dt
where again the first term on the right hand side is zero because P is fixed in B.
Equating the two identical quantities,
ω B/F 1 × rP/O = ω B/F 2 + ω F 2/F 1 × rP/O ,
or
ω B/F 2 + ω F 2/F 1 − ω B/F 1 × rP/O = 0.
Since P is arbitrary, we must have
ω B/F 1 = ω B/F 2 + ω F 2/F 1 .
(2.9)
In the above we have used the following fact: if a vector A satisfies the condition A × B = 0 for all vectors
B, then A = 0. Why this is true is easy to see. If A has, say, a nonzero î component, then we can choose
B = ĵ to obtain a nonzero k̂ component in A × B. Since A × B is zero for all choices of B, however, it must
be true that A’s î component is zero. Similar arguments apply for A’s ĵ and k̂ components.
2.5
Acceleration in a moving frame
Consider a frame xyz rotating with angular velocity ω and angular acceleration α relative to a fixed frame
XY Z. Consider a moving point P , as shown in figure 2.5.
In what follows, imagine that the position of the origin of xyz, i.e., R, is easily observed from XY Z, but
that the position of P is more conveniently described using the vector ρ in xyz. The velocity of P is then
understood to be
dρ
dρ
dr
dR
dR
=
+
=
+
+ ω × ρ.
vP =
dt XY Z
dt XY Z
dt XY Z
dt XY Z
dt xyz
The acceleration of P is obtained by differentiation as follows:
2 d dr
d R
d d dρ
aP =
ω
=
+
+
×
ρ
dt dt XY Z XY Z
dt2 XY Z
dt dt xyz
dt
XY Z
XY Z
2 2
dρ
dρ
d ρ
d R
dω
=
+
+ω×
+
×ρ+ ω×
+ω×ρ
dt2 XY Z
dt2
dt xyz
dt XY Z
dt xyz
xyz
2 dρ
d2 ρ
d R
=
+
+ 2ω ×
+α×ρ+ω× ω×ρ .
dt2 XY Z
dt2
dt xyz
xyz
CHAPTER 2. RELATIVE MOTION
24
α
Z
ω
z
P
ρ
r
y
x
R
Y
X
Figure 2.5: Frame xyz rotates with angular velocity ω and angular acceleration α relative to a fixed frame
XY Z. Point P moves in an arbitrary way.
We can write the above as
r̈XY Z = R̈XY Z + ρ̈xyz + 2 ω × ρ̇xyz + ω̇ × ρ + ω × ω × ρ ,
or still more briefly,
r̈ = R̈ + ρ̈ + 2ω × ρ̇ + ω̇ × ρ + ω × ω × ρ.
(2.10)
In the above, there are several implicit notational conventions. Derivatives of r and R are taken in XY Z,
derivatives of ρ are taken in xyz, and ω × ω × ρ is taken to mean ω × (ω × ρ), not (ω × ω) × ρ. Equation 2.10
is often called the five term acceleration formula, and will be useful in what follows.
In order to understand the five term acceleration formula it may be helpful to consider some situations in
which only one or two of the terms are nonzero.
α
B
A
ω
B
A
Figure 2.6: A gramophone record rotates with an angular velocity ω about the axis BB.
1. In the case of a body undergoing pure translation, we can attach xyz to the body. Then, for any point
P on the body, only the R̈ term of Eq. 2.10 is nonzero. If xyz is attached to XY Z instead, then ρ̈ is
nonzero instead of R̈.
2. See figure 2.6. A turntable rotates with constant angular velocity ω about a vertical axis (ignore α and
the soccer goalpost for the moment). We fix the frame xyz to the rotating table with its origin at the
2.5. ACCELERATION IN A MOVING FRAME
25
center of the table. The origin of the fixed frame XY Z is chosen at the same point. If P is any point
on the turntable (except the center), the acceleration of P is given by ω × ω × ρ, where ρ is the position
vector of P . This term, often called the centripetal term, is the only nonzero contributor to Eq. 2.10 in
this case. If the spin rate ω varies, then ω̇ × ρ is nonzero as well. And if the turntable starts from rest
with nonzero ω̇, then ω̇ × ρ is the only nonzero term.
3. Consider the same turntable again, with ω varying, but with a radial groove in the turntable. A point
P moves along this groove with a constant speed. At the instant when it reaches the center of the table
the only nonzero term in Eq. 2.10 is 2ω × ρ̇. If the speed of P is not constant in the groove, then ρ̈ is
nonzero as well.
4. Consider the same turntable yet again. This time the spin rate is constant, but the base of the turntable
is itself being accelerated from rest with an angular acceleration α about a horizontal axis AA, as shown
in the figure. Also, the assembly of two vertical rods and a horizontal one (soccer goalpost) is attached
to the turntable. A point P moves along the top of the goalpost. At the instant when point P crosses
the axis of the turntable, 2ω × ρ̇ and ω̇ × ρ are the only nonzero contributors to Eq. 2.10 (note that here
the ω̇ comes from α).
Problem: A point mass falls under the action of gravity. Air resistance is negligible; and so is the variation
in g with height. Find the acceleration as seen on Earth, taking the Earth’s rotation into account.
Solution: Fix the origin of a non-rotating reference frame XY Z (assumed inertial) to the center of the Earth.
Another reference frame xyz with origin at the same point rotates with the earth. The rotation is at a fixed
rate ω and about the Z axis. So R = 0, and r = ρ (we use the same symbols as in Eq. 2.10).
By Eq. 1.1,
ρ
mr̈ = −mg .
|ρ|
By Eq. 2.10,
r̈ = R̈ + ρ̈ + 2ω × ρ̇ + ω̇ × ρ + ω × ω × ρ.
Since R̈ and ω̇ are zero, we get
ρ̈ = −g
ρ
− 2ω × ρ̇ − ω × ω × ρ,
|ρ| correction terms
which includes two correction terms. Of these, the first is called the Coriolis term (which, e.g., is important
for monsoon winds); and the second affects the apparent value of the acceleration due to gravity. For order of
magnitude estimates, we can use
ω=
2π
rad/s and ρ = radius of the Earth = 6.4 × 106 m.
86400
Since the rotation rate of the Earth is fairly small, the ω × ω × ρ term is small as well (compared to g ≈ 10),
with ω 2 ρ = 0.034 m/s2 . However, for 2ω × ρ̇ to be comparable to this correction, we need something like
2ω|ρ̇| = 0.034 m/s2 , which gives |ρ̇| ≈ 233.7 m/s, which is in excess of 840 km/h.
Problem: See figure 2.7. A turntable rotates at a constant angular rate Ω. A wheel is located on the turntable
as shown in figure 2.7. The wheel spins at a rate of β̇ about an axis fixed relative to the turntable. Find the
acceleration of a point P on the rim of the wheel.
Solution: See figure 2.8. Fix a non-rotating reference frame XY Z with its origin at the center of the turntable.
Consider a frame xyz that rotates with the turntable and also has its origin at the center of the turntable.
Let Q be the center of the wheel. We can write
rOP
= rOQ + rQP .
In what follows, rOP will play the role of ρ in the five term acceleration formula (Eq. 2.10).
CHAPTER 2. RELATIVE MOTION
26
Z
Ω
Y
Z
P
β
Y
O
β
X
X
(a)
(b)
Figure 2.7: (a) A spinning wheel on a rotating turntable. (b) Top view.
Z
z
P
O’
y
Y
O
β
x
X
Figure 2.8: Frame xyz rotates with the turntable, with its origin coinciding with the center of the turntable
and the origin of the fixed frame XY Z.
The velocity of P is
vP
drOP
dt
xyz
drOP
dt
XY Z
drOP
=
+ Ω k̂ × rOP , where
dt
xyz
*0 drOQ drQP
=
+
.
dt
dt
xyz
xyz
=
The canceled term above is zero because Q is fixed in xyz. Since the wheel is rotating at an angular rate
of β̇ relative to xyz, and since vector rQP is fixed to the wheel, application of Eq. 2.5 (using the wheel as the
rotating frame and xyz as the fixed frame) gives
drQP
=
−β̇ î × rQP .
dt
xyz
2.6. EXERCISES
27
The right hand side above plays the role of ρ̇ in the five term acceleration formula (Eq. 2.10).
Taking stock of our progress towards using Eq. 2.10, we write
aP
0
0
*
= R̈
× ρ̇ + ω̇ ×
ρ+ω×ω×ρ
+ ρ̈ + 2ω
= ρ̈ + 2Ωk̂ ×
−β̇ î × rQP + Ωk̂ × Ωk̂ × rOP .
So we still need ρ̈. To find it, we do a separate calculation on the side.
Z
Ω
z’
β
P
y’
Y
O
x’
X
Figure 2.9: A new coordinate system x y z is fixed to the center of the wheel.
Side calculation: We put this separate calculation in a box to avoid confusion. Here, too, we will
use the five term acceleration formula. For ease of understanding, we will use it in terms of the same
symbols as in Eq. 2.10. So, for example, ρ̈ inside this box will not mean the same thing as outside
it. Consider a new frame x y z with origin at the center of the wheel, and rotating with the wheel,
as shown in figure 2.9. Now treating xyz as the fixed frame, the acceleration of the point P with
respect to xyz is given by
0 0
0
0
*
:
2ω×
R̈
ρ̇ + ω̇ ×
ρ + ω × ω × ρ,
+ ρ̈ +
= −β̇ î × −β̇ î × rQP .
Here ρ is the position vector of P as seen in x y z frame. Since the origin of frame x y z is fixed in
xyz, the first term is zero. Since P is fixed in x y z the second and third terms are zero. Since the
wheel is rotating at a constant rate, the fourth term is zero as well, leaving only the last centripetal
term.
Finally, the acceleration of point P is
aP = −β̇ î × −β̇ î × rQP + 2Ωk̂ × −β̇ î × rQP + Ωk̂ × Ωk̂ × rOP .
2.6
Exercises
1. See figure 2.10. A turntable T is placed in a cradle C. The cradle rotates about a point P with angular
velocity Ωk̂. A disk D on the turntable rotates relative to the cradle at a constant rate ω about a vector
CHAPTER 2. RELATIVE MOTION
28
P
Ω
^
j
^i
ω
D
T
C
Figure 2.10: A turntable on a cradle.
directed towards P as shown. At the instant shown, find the angular velocity and angular acceleration of the
disk relative to stationary ground.
2. See figure 2.10 again. Let the center of the disk be O. Let rO/P at the instant of interest be −Lĵ. Consider
a point Q on the disk, such that at the instant of interest rQ/O = Rî. Find the velocity and acceleration of Q
relative to stationary ground.
3. A gramophone record rotates at a constant rate of 78 RPM. The tip of a pin touching the record at an
average radial distance of 10 cm from the centre of the record executes sinusoidal oscillations in this radial
direction with an amplitude of 0.1 mm and a frequency of 1000 Hz, when viewed from a stationary frame.
At the instant of interest, the pin tip is passing through its average position, moving outwards. What is
the velocity and acceleration of the pin tip at this instant, considered relative to the rotating gramophone
record frame? Take the x-axis to lie along the radial direction, y-axis to lie in the plane of the record and
perpendicular to x, and z-axis to be normal to the plane of the record.
4. See figure 2.11. A system of 4 linked rigid bodies is shown schematically. Link P Q rotates about P , relative
Q
R
P
T
S
Figure 2.11: A linkage.
to ground, about a unit vector n̂P fixed relative to ground, at a time varying rate ωP . Link QR rotates about
Q, relative to link P Q, about a unit vector n̂Q fixed relative to link P Q, at a time varying rate ωQ . Link
RS rotates about R, relative to link QR, about a unit vector n̂R fixed relative to link QR, at a time varying
rate ωR . Link ST rotates about S, relative to link RS, about a unit vector n̂S fixed relative to link RS, at a
time varying rate ωS . Write a computer program that, given all other quantities of interest at some instant of
2.6. EXERCISES
29
time, will find the velocities and accelerations of points Q, R, S and T ; as well as the angular velocities and
accelerations of the links.
Partial solution: Start at the bottom.
v P = 0,
aP = 0,
ω P Q = ωP n̂P ,
αP Q = ω̇P n̂P .
Now move up link by link.
Work out the angular velocities and angular accelerations first. Consider link QR. Attach a moving frame
to link QR.
(k)
ω QR = ω QR/P Q + ω P Q = ωQ n̂Q + ω P Q ,
where the superscript “(k)” says a quantity has been determined and is now known. Also,
dω QR/P Q
(k)
(k)
αQR =
+ ω P Q × ω QR/P Q + αP Q = ω̇Q n̂Q + ω P Q × ωQ n̂Q + αP Q .
dt
PQ
Now consider link RS. Attach a moving frame to link QR.
(k)
ω RS = ω RS/QR + ω QR = ωR n̂R + ω QR .
Also
(k)
(k)
αRS = ω̇R n̂R + ω QR × ωR n̂R + αQR ,
and so on up the chain to find all the angular velocities and accelerations.
Now consider the velocities and accelerations of the pivot points. Consider point Q. Attach a moving
frame to P Q.
(k)
(k)
v Q = v P + ω P Q × rQ/P .
Also, towards using the five term acceleration formula, we note that
(k)
r̈ = aP ,
ρ̈ = 0,
ρ̇ = 0,
giving
(k)
(k)
(k)
(k)
aQ = aP + ω P Q × ω P Q × rQ/P + αP Q × rQ/P .
Now consider point R. Attach a moving frame to link QR.
(k)
(k)
v R = v Q + ω QR × rR/Q .
Also,
(k)
(k)
(k)
(k)
aR = aQ + ω QR × ω QR × rR/Q + αQR × rR/Q ,
and so on up the chain.
30
CHAPTER 2. RELATIVE MOTION
Chapter 3
Momentum Balance for Rigid Bodies
3.1
Linear momentum balance
The linear momentum of a system of particles is
L=
mi v i = mtot v cm ,
i
where the velocities are understood to be measured in an inertial frame of reference.
Differentiating.
L̇ = mtot acm or, equivalently, mtot acm .
If F is the net external force acting on this system of particles (for which we refer to a free body diagram),
then (see the laws of dynamics, section 1.4)
F = L̇.
If there are n rigid bodies, we can consider their n distinct masses and centers of mass, using
F =
n
mtot,k acm,k .
k=1
Problem: A spool of radius R and mass m has a string wound around it. At some instant the string is
vertical and pulled upwards by a force T . What is the acceleration of the center of mass of the spool?
R
T
m
Figure 3.1: Spool pulled upwards.
Solution: The spool may spin as it moves, but this is irrelevant. A free body diagram of the spool shows only
a vertical force F acting upwards, and the weight mg acting downwards. So the acceleration is in the vertical
31
CHAPTER 3. MOMENTUM BALANCE FOR RIGID BODIES
32
direction; and is
T − mg
.
m
a =
3.2
Angular momentum
The angular momentum of a system of particles about some point C in space is given by
H /C =
ri/C × mi v i ,
(3.1)
i
where the velocities are understood to be measured in an inertial frame of reference.
Referring to the location and velocity of the center of mass, we write
v i = v i/cm + v cm
and
ri/C = ri/cm + rcm/C ,
whence
H /C
=
ri/cm × mi v i/cm +
i
+
(1)
ri/cm × mi v cm +
i
rcm/C × mi v i/cm
i
(2)
rcm/C × mi v cm
i
(3)
(4)
is seen to be the sum of four terms which we have labeled for individual consideration below.
Term no. 3 in Eq. 3.2 is
ri/cm × mi v cm =
mi ri/cm × v cm since mi is a scalar.
i
i
This term is actually zero. To see why, we start with
ri = ri/cm + rcm ,
multiply each equation by mi , and then sum over i to find
mi r i =
mi ri/cm +
mi rcm .
i
i
( =
mtot
rcm
((
(
0 =
i
i
(, or
mi ri/cm + (
mtot
rcm
((
mi ri/cm .
i
Term no. 2 in Eq. 3.2 is
i
rcm/C × mi v i/cm = rcm/C ×
mi v i/cm
,
i
which is zero as well. To see this, consider
d
mi ri/C =
mi ri/cm +
mi rcm/C
.
mtot rcm/C =
dt
i
i
i
inertial frame
(3.2)
3.3. MATRICES AND VECTORS
33
From the above,
mtot v cm =
mi v i/cm +
i
(=
mtot
v cm
((
(
and so
i
mi v cm ,
i
(,
mi v i/cm + (
mtot
v cm
((
mi v i/cm = 0 .
i
So Eq. 3.2 becomes
H /C =
ri/cm × mi v i/cm +
i
which in turn is
H /C =
rcm/C × mi v cm ,
i
ri/cm × mi v i/cm
+ rcm/C × mtot v cm .
(3.3)
i
In physical terms, the angular momentum of a system of particles about a point C can be broken into two
parts. One consists of the contribution from particles moving relative to the center of mass, and is independent
of the location of C . The other considers the system to be an effective point mass mtot concentrated at the
system’s center of mass, and is independent of the relative motions between particles in the system. Interpreting
these terms in the context of a single rigidy body, the first contribution comes from rigid body rotations while
the second comes from center of mass translations.
So far, we have made no assumptions specific to rigid bodies. The simplification for a single rotating rigid
body is that (using, e.g., Eq. 2.5 with A fixed in a rotating frame attached to the rigid body)
v i/cm = ω × ri/cm .
We will exploit this below, after introducing some matrix notation.
3.3
Matrices and vectors
We first choose and fix a right handed orthonormal co-ordinate system.
Now any vector a is written as
a = ax î + ay ĵ + az k̂.
This can be unambiguously written as
⎧
⎫
⎨ ax ⎬
ay
.
a≡
⎩
⎭
az
Comment: Strictly speaking, matrices of components are not equal to the vectors they represent because
changing coordinate systems changes the components without changing the vectors. However, because we
choose and fix a specific right handed orthonormal coordinate system, vectors are equivalent to their matrices
of components. While retaining the right to change coordinate systems (and thus not allowing vectors and
their matrices of components to be equivalent) has theoretical advantages for more advanced treatments of
mechanics, we find that accepting and using the vector-matrix equivalence helps at this level, at least from a
problem solving and computer programming point of view.
Notation: “a” or “a” represents a vector and “a” represents its matrix of components. We will use these
interchangeably.
The dot of product of two vectors can now be written as
a · b = aT b = bT a,
where the T -superscript denotes matrix transpose.
CHAPTER 3. MOMENTUM BALANCE FOR RIGID BODIES
34
Consider
a×b
(ay bz − az b) î + (−ax bz + az bx ) ĵ + (ax by + ay bx ) k̂
⎫
⎤⎧
⎡
0
−az ay
⎨ bx ⎬
0
−ax ⎦
by
≡ ⎣ az
.
⎩
⎭
−ay ax
0
bz
=
Define the skew symmetric matrix
⎡
0
S(a) = ⎣ az
−ay
−az
0
ax
⎤
ay
−ax ⎦ .
0
(3.4)
Then
a × b ≡ S(a)b = −S(b)a.
3.4
Angular momentum of a rigid body
Recall Eq. 3.3, in which the term within round brackets is
ri/cm × mi v i/cm , where for a rigid body v i/cm = ω × ri/cm ,
i
which in the above matrix notation becomes
T
mi S(ri/cm )S (ri/cm ) ω =
i
T
mi S(ri/cm )S (ri/cm )
ω.
(3.5)
i
The quantity within round brackets in the right hand side above is a 3 × 3 matrix which is independent of
ω, and is a property of a given rigid body in a given configuration. This matrix, called the moment of inertia
matrix about the center of mass, is denoted by Icm .
In the above, if ri/cm and S(ri/cm ) are represented as
⎧
⎫
⎨ x ⎬
y
ri/cm =
⎩
⎭
z
⎡
and
0
S(ri/cm ) = ⎣ z
−y
respectively, then
⎡
S(ri/cm )S T (ri/cm )
−z
0
x
⎤
y
−x ⎦ ,
0
⎤⎡
0 −z y
0 −x ⎦ ⎣
= ⎣ z
−y x
0
⎡ 2
2
y +z
−xy
x2 + z 2
= ⎣ −xy
−xz
−yz
⎤
z −y
0
x ⎦
−x 0
⎤
−xz
−zy ⎦ .
x2 + y 2
0
−z
y
If the rigid body is now seen as constructed from infinitely many point masses distributed in a continuum
(see figure 3.2), then Icm is viewed as a mass-weighted integral of the form
⎤
⎡
Ixx Ixy Ixz
Icm = ⎣ Iyx Iyy Iyz ⎦ ,
Izx Izy Izz
3.4. ANGULAR MOMENTUM OF A RIGID BODY
35
Z
dm
z
y
O
Y
x
X
Figure 3.2: Moment of inertia of an arbitrary rigid body calculated by integration.
with
Ixx
2
=
2
Body volume
Iyy
2
=
2
(x2 + z 2 )ρ dv,
(x + z ) dm =
Body mass
Body volume
Izz
(y 2 + z 2 )ρ dv,
(y + z ) dm =
Body mass
2
(x + y ) dm =
(x2 + y 2 )ρ dv,
Body volume
= −
xy dm = −
xyρ dv,
Body mass
Body volume
= −
xz dm = −
xzρdv,
Body mass
Body volume
= −
yz dm = −
yzρ dv,
=
2
Body mass
Ixy = Iyx
Ixz = Izx
and Iyz = Izy
Body mass
Body volume
where ρ is density (not to be confused with the vector ρ used in the five term acceleration formula). Note
that if the body is rotated, but the coordinate system is rotated identically, then Icm remains the same: we
imagine such rotations, but do not actually conduct them, because we have chosen and fixed a coordinate
system. However, this constancy of Icm is expressed here for further use as follows:
d
Icm
= 0,
(3.6)
dt
Body
where the 0 represents a zero matrix.
We observe that Ixx ≤ Iyy + Izz , Iyy ≤ Ixx + Izz , and Izz ≤ Ixx + Iyy for a general body. As a special case,
consider a flat thin object in the x-y plane (figure 3.3). For this object, taking z = 0 for every mass element,
we find
(3.7)
Iyz = Ixz = 0, and Izz = Ixx + Iyy .
In the above example, if the object is not very thin, but simply symmetrical about the x-y plane, then
Iyz = Ixz = 0, and Izz < Ixx + Iyy .
(3.8)
Comment: Icm , if we do not use the matrix-vector equivalence described in the previous section, should be
written as I cm , a symmetric Cartesian tensor of rank 2. Definitions of tensors and how they behave may, e.g.,
be seen at http://mathworld.wolfram.com/Tensor.html or read in a good book on tensor analysis. Here,
CHAPTER 3. MOMENTUM BALANCE FOR RIGID BODIES
36
Z
Y
X
Figure 3.3: A flat, thin object.
we will write Icm ω ≡ I cm · ω.
Now, Eq. 3.3 becomes
H /C = I cm · ω + rcm/C × mtot v cm .
(3.9)
If we simultaneously consider n rigid bodies, then their net angular momentum about C will be the sum of n
pairs of terms, each of the form in Eq. 3.9.
3.5
Angular momentum balance for a rigid body
Recall Eq. 3.1. Differentiating,
Ḣ /C =
i
:0 ×
m
v
+
ri/C × mi ai ,
v
i
i
i
i
or (using Eq. 1.12)
M /C = Ḣ /C
for any system of particles.
get
What is Ḣ /C for a single rigid body? Differentiating Eq. 3.9 with respect to time in our inertial frame, we
Ḣ /C
:0
v + r
= v cm
×
m
tot
cm
cm/C × mtot acm +
d
[I · ω]
dt cm
.
Inertial frame
In the above,
d
[I · ω]
dt cm
Inertial frame
d
[I · ω]
=
+ ω × I cm · ω
dt cm
Body
dω
= I cm
+ ω × I cm · ω,
dt Body
where we have used Eq. 3.6. Moreover,
dω
dω
:0
=
+
ω
×
ω.
dt Inertial frame
dt Body
3.6. PLANAR OR 2D PROBLEMS
37
Therefore
where α =
3.6
dω
dt
M /C = Ḣ /C = rcm/C × mtot acm + I cm .α + ω × I cm · ω,
(3.10)
.
Inertial frame
Planar or 2D problems
For planar problems, from Eq. 3.8, we know that
⎡
Icm
Ixx
=⎣ ∗
0
∗
Iyy
0
⎤
0
0 ⎦,
Izz
where the “∗” represents some nonzero value. Also,
⎧
⎫
⎨ 0 ⎬
0
ω=
⎩
⎭
ωz
⎧
⎨
⎫
0 ⎬
0
Icm ω =
,
⎩
⎭
Izz ωz
and so
which is parallel to ω. Hence
ω × I cm · ω = 0
in Eq. 3.10.
3.7
The kinetic energy of a rigid body
Kinetic energy does not enter into direct considerations of momentum balance. However, in using Lagrange’s
equations (Chapter 5), we will have use for the kinetic energy of a rigid body.
We begin with a system of particles viewed in an inertial frame. The kinetic energy is
mi mi
vi · vi =
v cm + v i/cm · v cm + v i/cm
KE =
2
2
i
i
=
mi
i
2
v cm · v cm +
: 0 mi
v
mi v
v i/cm
+
·v
.
·
cm
2 i/cm i/cm
i
i
The second term above is zero by Eq. 3.3. We then have
mi
mtot
v cm · v cm +
v
KE =
·v
.
2
2 i/cm i/cm
i
The first term above represents the kinetic energy of overall translation, and is obtained by treating the entire
mass of the system as concentrated at the center of mass; the second term represents the additional kinetic
energy due to relative motion between the individual particles. The above expression is valid for any system
of particles. In particular, for a single rigid body, we have
v i/cm = ω × ri/cm = −ri/cm × ω,
and so (using matrix notation)
KE
=
=
mi T
mtot T
vcm vcm +
S(ri/cm )ω S(ri/cm )ω
2
2
i
mtot T
1 T T
v vcm + ω
mi S (ri/cm )S(ri/cm ) ω.
2 cm
2
i
CHAPTER 3. MOMENTUM BALANCE FOR RIGID BODIES
38
Comparing with Eq. 3.5, we observe that the quantity within brackets is Icm , and so the kinetic energy of a
rigid body is
1
mtot T
v vcm + ω T Icm ω.
(3.11)
KE =
2 cm
2
The first term above is called the translational kinetic energy of the rigid body, and the second term is called
its rotational kinetic energy. By the discussion in section 3.6, the rotational kinetic energy of a planar object
1
restricted to the x-y plane is Izz ωz2 . Many authors denote Izz in such situations using the symbol J.
2
3.8
Exercises
1. Find the moment of inertia matrix about the center of mass of a uniform sphere of mass m and radius R;
a uniform cube of mass m and side a; a uniform right circular cylinder of mass m, radius R and height h; and
a uniform right circular cone of base radius R and height h. The cube has its faces parallel to the coordinate
planes. The cone and cylinder rest with a flat face on the x-y plane.
2. A car is being driven at a constant speed on a circular road so that its rear right tyre traverses a circle
of radius 50 metres. The speed of the car is 72 km/h (you can take this to be the speed of the centre of the
rear right tyre). The radius of the tyre is 30 cm. Model the tire as a planar, circular, rigid object of mass m
and moment of inertia J about the axle; assume that the force on the tyre from the ground has no component
along the path (i.e., no forward component); and find the vertical component of the moment exerted on the
tyre by the axle.
3. See figure 3.4. A uniform thin disk of mass m and radius R spins at a rate ω on frictionless bearings
Ω
ω
Figure 3.4: Wheel on an axle.
mounted on a horizontal shaft as shown in the figure. The shaft itself is rotated about a vertical axis at a rate
Ω. Find the vector moment exerted by the bearing on the disk.
4. See figure 3.5. A uniform circular cylinder of mass m, radius R and height h is mounted on a light, rigid
^
e
cylinder
θ
^i
shaft
Figure 3.5: Skewed cylinder on a shaft.
shaft. The shaft passes through the center of mass of the cylinder. The unit vector along the axis of symmetry
3.8. EXERCISES
39
of the cylinder, ê, makes an angle θ with the shaft as shown. The shaft is mounted on frictionless bearings,
and rotates at an angular rate of ω. Find the kinetic energy of the cylinder. Take î, the unit vector along the
x-axis, to be along the shaft as shown. Let ĵ, the unit vector along the y-axis, be in the vertical direction. At
an instant when ê lies in the x-y plane, find the angular momentum of the cylinder. At the same instant, find
the moment exerted by the shaft on the cylinder. You may like to do your intermediate calculations using the
unit vectors k̂, ê, and (say) λ̂ = k̂ × ê.
40
CHAPTER 3. MOMENTUM BALANCE FOR RIGID BODIES
Chapter 4
Rotations and Euler Angles
4.1
Rotations do not commute
Rotations, although they might have magnitudes and directions, are fundamentally different from translations.
If a body undergoes two successive rotations, the final configuration depends on which rotation occurs first.
For example, if a person standing upright and facing north is first given a quarter rotation about the
vertical axis, and then a quarter rotation about the westerly direction, then he ends up lying on his side,
facing west, with his head pointing north.
Reversing the sequence, if the person standing upright and facing north is first given a quarter rotation
about the westerly direction and then a quarter rotation about the vertical, then he ends up lying on his chest,
facing down, with his head pointing west. An entirely different configuration.
The fact that rotations do not commute (i.e., “rotation 1 followed by rotation 2” = “rotation 2 followed
by rotation 1”) makes them more complicated than translations.
4.2
Some basic mathematical facts
The following will be used below.
1. Standard basis: We call e1 = {1, 0, 0}T , e2 = {0, 1, 0}T , and e3 = {0, 0, 1}T . These represent î, ĵ and k̂
respectively.
2. Linearity: A (real) function f (x) is linear if and only if the following two conditions hold true:
• f (βx) = βf (x) for all real numbers β and all x, and
• f (x + y) = f (x) + f (y) for all x and y.
3. Matrix representation of linear functions: Linear functions on finite dimensional vector spaces (ours
is three dimensional), upon choice of a coordinate system (we have chosen one), are represented by a
matrix multiplication of the form f (x) = Rx for some suitable R that is fixed uniquely for any given
coordinate system but will change if we change the coordinate system (we will not).
4. Positive definite matrices: If R is invertible, then RT R is symmetric and positive definite (SPD): it has
three real strictly positive eigenvalues. If R is singular, then RT R is symmetric and positive semidefinite:
it has three real nonnegative eigenvalues.
5. Identity matrix: I (no confusion with Icm ) will denote the identity matrix. Also, if an SPD matrix A
satisfies the condition rT Ar = rT r for all r, then A = I.
6. Determinants: det(AB) = det(A) det(B); det(AT ) = det(A); and det(A) equals the product of the
eigenvalues of A.
41
CHAPTER 4. ROTATIONS AND EULER ANGLES
42
4.3
Euler’s theorem and the rotation matrix
Euler’s Theorem: The most general displacement of a rigid body with one point fixed is a rotation through
some axis passing through that point.
To prove Euler’s theorem we will use the following result, without proof: The most general displacement
of a rigid body with one axis fixed is a rotation about that axis. (The reader can imagine a door on hinges.
The hinges fix points along a certain line on the door, i.e., the door has one axis fixed.) Let a rigid body be
given an arbitrary displacement that keeps fixed a point O on the body (see figure 4.1). Consider any vector
r whose tail is at O. After displacement, r becomes r (say). We write r = f (r).
Initial configuration
r
r+q
Final configuration
q
r'
βr
βr'
O
Figure 4.1: A rigid body with one point fixed is given an arbitrary displacement.
Claim: f is linear.
To see this, we need two checks.
First, imagine a straight line passing through O. Think of it as a black rod, if you like. On this rod,
imagine two dots: one white, one red. Then the red dot can represent r, while the white dot can represent
βr. After displacement, the rod has a possibly different orientation but is otherwise undeformed, so we see
that f (βr) = βf (r).
Next, imagine a parallelogram with one vertex at O. Let the sides of the parallelogram that meet at O
represent vectors r and q. The diagonal starting at O is the vector r + q. After displacement the parallelogram
has a different orientation but is undeformed, and so f (r + q) = f (r) + f (q). This proves the above claim.
Therefore, r = Rr for some R. Now the body is rigid, so rT r = rT r = rT RT Rr for all r, therefore
R R = I. This means R is an orthogonal matrix, and that RRT = I as well. Since
T
det(RT R) = [det(R)]2 = det(I) = 1
we find det(R) = ±1.
Since the displacement of the body evolves continuously from the reference position where det(R) = 1
(since R is simply I); and since f and hence R is a continuous function of displacement; and since det(R) is
a continuous function of R; we find that det(R) cannot jump discontinuously1 to −1 from its initial value of
+1. Therefore det(R) = 1.
1 For
reflections, the determinant will be −1.
4.4. OBTAINING THE ROTATION MATRIX
43
R is 3 × 3, so it has 3 eigenvalues, of which at least one is real. Let λ be a real eigenvalue of R, and u the
associated eigenvector. Then Ru = λu, and uT RT Ru = uT u = λ2 uT u, therefore λ = ±1.
If R has three real eigenvalues, then the possibilities are (1, 1, 1) and (1, −1, −1), and neither (1, 1, −1) nor
(−1, −1, −1), because det(R) = 1. On the other hand, if R has a single real eigenvalue λ, and two complex
conjugate eigenvalues σ and σ̄, then det(R) = λ |σ|2 = 1 which implies |σ| = 1 and λ = 1. In either case,
λ = 1 is an eigenvalue of R. This means there is a vector u such that Ru = u, i.e., there is an axis in the
body, passing through O, that is fixed during the displacement. The displacement of the body is therefore a
rotation about that axis. This proves Euler’s theorem.
4.4
Obtaining the rotation matrix
How do we obtain R?
λ
n
P
θ
r r
C
C
P
O
θ
e
P
=
λ
n
λ
Figure 4.2: Left: A rigid body with one point O held fixed is rotated through an angle θ about a unit vector
n̂. In the process, a vector r embedded rigidly in the body gets mapped to r . During rotation, the tip of r
traces out a circle with center at C. The unit vector in the direction from C to the tip of r (point P ) is called
λ̂. Right: A view of the same, looking down from the tip of vector n̂ towards its tail, defining a new unit
vector ê.
See figure 4.2. A body is rotated about a unit vector n through an angle θ. It is left as an exercise for the
reader to show, guided by the figure, that a typical vector r gets rotated into a new vector r that is given by
the vector relation
r = n(n · r) + cos θ [r − n(n · r)] − sin θ r × n ,
which is equivalent to the matrix relation
!
r = cos θ I + (1 − cos θ)nnT + sin θ S(n) r ,
so we can now define the matrix corresponding to the rotation explicitly as
R(n, θ) = cos θ I + (1 − cos θ)nnT + sin θ S(n) .
(4.1)
This is sometimes called the axis-angle formula. In computational platforms like Matlab or Maple, which we
have not discussed so far, the formula can be implemented in an m-file or Maple procedure respectively.
CHAPTER 4. ROTATIONS AND EULER ANGLES
44
4.5
Successive rotations
If we have a rotation R1 , r goes to r = R1 r. If we follow this with a rotation R2 , r goes to r = R2 r = R2 R1 r.
By Euler’s theorem, the net displacement is a single rotation, which we now see is simply R2 R1 . Again, since
R2 R1 = R1 R2 in general, rotations do not commute.
4.5.1
Infinitesimal Rotations
By Eq. 4.1, an infinitesimal rotation (ignoring second order quantities) is given by
R(n, θ) = I + θS(n) ,
and therefore two successive infinitesimal rotations are given by
R(n2 , θ2 )R(n1 , θ1 ) = I + θ1 S(n1 ) + θ2 S(n2 ) = R(n1 , θ1 )R(n2 , θ2 ) ,
again ignoring second order quantities. This shows that infinitesimal rotations commute, while finite (nonzero)
rotations in general do not.
4.5.2
Angular velocity
Since infinitesimal rotations commute, they are vectors (unlike finite nonzero rotations, which have magnitude
and direction but do not commute). We therefore anticipate that the rate of change of orientation should be
expressible using a vector, which will be called the angular velocity of the body (recall Eq. 2.2).
T
Let R be a function of time. Since RRT = I, we have ṘRT + RṘT = 0, so ṘRT = − ṘRT , or ṘRT is
a skew symmetric matrix. We define the vector ω by
S(ω) = ṘRT .
(4.2)
Now consider the time-varying vector r , given by r = Rr, with r some fixed vector. Then
r˙ = Ṙr = ṘRT Rr = S(ω)r ,
which is the familiar expression for the velocity of a point in a body with one point fixed at the origin,
v = ω × r ,
from undergraduate level vector mechanics (also Eq. 2.1, with A = r ). Thus, ω ≡ ω as defined in Eq. 4.2 is
the angular velocity of the body.
4.5.3
The integral of angular velocity is not meaningful
If ω = ωx î + ωy ĵ + ωz k̂ is the angular velocity of a body, then it can be interpreted in terms of rotation rates.
In an infinitesimal time interval Δt, the infinitesimal rotations about the x, y and z axes are ωx Δt, ωy Δt and
ωz Δt respectively; these infinitesimal rotations commute, and can be treated as vectors.
The angular velocity can also be integrated over time to give (in matrix notation)
⎧ t2
⎫
⎪
⎪
⎪
⎪
ω
dt
x
⎪
⎪
⎪
⎪
⎪
⎪
t
⎪
⎪
t2
⎨ 1t2
⎬
,
ω dt =
ωy dt
⎪
⎪
t1
t1
⎪
⎪
⎪
⎪
t2
⎪
⎪
⎪
⎪
⎪
⎩
⎭
ωz dt ⎪
t1
but the results have no physical significance as rotations. Although the three components of the integrated
vector look like rotations, it is not clear what sequence to take them in; and as we have observed, when more
than one rotation is involved then the sequence is crucial.
If not directly as rotations, can these integrals at least be viewed as changes in some as yet unknown
functions of some as yet unspecified rotation coordinates? We address this in section 4.5.5.
4.5. SUCCESSIVE ROTATIONS
4.5.4
45
Rotations upto second order in infinitesimals
Up to second order in infinitesimals, by Eq. 4.1, a rotation (ignoring third order quantities) is given by
θ2
θ2
nnT + θS(n) .
R(n, θ) = 1 −
I+
2
2
Now consider two successive infinitesimal rotations about two mutually perpendicular axes n1 and n2 . Again
up to second order in infinitesimals, we find
2
2
2
2
θi2
i=1 θi
ni nTi +
R(n2 , θ2 )R(n1 , θ1 ) = 1 −
θi S(ni ) + θ1 θ2 S(n2 )S(n1 ),
I+
2
2
i=1
i=1
where the only non-commuting contribution is from
θ1 θ2 S(n2 )S(n1 ).
However, this is enough. Since
S(n2 )S(n1 )a ≡ n2 × (n1 × a)
for any vector a, we can let a = n̂1 to see that S(n2 )S(n1 ) = S(n1 )S(n2 ) if n1 and n2 are mutually perpendicular.
Therefore, if second order infinitesimals are considered, then successive rotations about two mutually
perpendicular axes definitely do not commute.
4.5.5
Angular velocity is not an exact derivative
How do we represent rotations? By Euler’s theorem or the axis-angle formula, we know that a unit vector
and an angle can together specify any rotation (i.e., three parameters are needed).
Suppose we choose some set of three parameters, say p = {p1 , p2 , p3 }, that can describe rotations. Suppose
that there exists some vector function f (p) such that
df
= ω,
dt
i.e., suppose that angular velocity is the derivative of some as yet unknown vector function of some as yet
unspecified rotation coordinates. Then we must have, for any Δt > 0 and any ω,
t0 +Δt
t0 +Δt
df =
t0
ω dt.
t0
Let us choose two mutually perpendicular vectors n̂1 and n̂2 , and let
#
n̂1 , t0 ≤ t ≤ t0 + Δt
2 ,
ω=
n̂2 , t0 + Δt
≤
t
≤
t
0 + Δt.
2
i.e., the body rotates about n1 for half the time and about n2 for the other half.
The net change in f is
Δt
Δf = (n̂1 + n̂2 )
.
2
Now if we change the sequence of n̂1 and n̂2 in ω for the above calculation, the change in f remains the same.
However, as demonstrated in section 4.5.4, the net rotations in the two cases are not equivalent. Therefore,
the changes in any locally invertible vector function of any possible choice of rotation coordinates are not the
same either.
This proves that angular velocity is not the derivative of any function of any possible choice of rotation
coordinates.
CHAPTER 4. ROTATIONS AND EULER ANGLES
46
4.6
Euler angles
Enough generalities. We will now parameterize the rotated configuration of the body using 3-1-3 Euler angles,
which are described as follows. First, we rotate the body through an angle φ about e3 (see standard basis,
section 4.2); then through θ about the now-rotated e1 ; and finally, through ψ about the now-twice-rotated e3 .
Defining R1 = R(e3 , φ), R2 = R(R1 e1 , θ) and R3 = R(R2 R1 e3 , ψ), the final rotation matrix Rf is given by
Rf = R3 R2 R1 .
(4.3)
Note that other sequences of Euler angles could be used, such as 1-2-3, 3-2-1, 3-2-3, etc. The only constraint
is that no two successive rotations are allowed to be about the same axis.
Note also that if the second rotation θ = 0, then the first rotation φ and the last rotation ψ are made
about the same axis. Therefore if the final configuration is given, say, then from that information only φ + ψ
can be determined: the individual rotations φ and ψ cannot be uniquely determined. For such configurations,
the choice of coordinates (φ, θ, ψ) as defined above is singular. The singularity in the 3-1-3 Euler sequence
is not peculiar to this particular sequence. All three-parameter descriptions of rotations have such singular
configurations.
Four-parameter descriptions of rotations, along with an added constraint equation, are possible and even
popular; they can also be singularity-free. However, we will only use 3-1-3 Euler angles here.
4.7
Angular velocity, again
Equation 4.2 shows that the angular velocity must be a linear function of the time derivatives of the Euler
angles. The actual expression can be found more easily as follows.
Let us use the reference position of the rigid body to define a reference frame F1 . Let F2 be a frame that is
always rotated about the e3 axis through φ relative to F1 . Let frame F3 be always rotated about the φ-rotated
e1 axis through θ relative to F2 . Finally, let F4 be always rotated about the first-φ-then-θ-rotated e3 axis
through ψ relative to F3 . It is clear that F4 is attached to the moving rigid body.
Now, the angular velocity of the body ω is clearly ωF4 /F1 , and
ωF4 /F1
= ωF4 /F3 + ωF3 /F2 + ωF2 /F1
(4.4)
= R2 R1 e3 ψ̇ + R1 e1 θ̇ + e3 φ̇ ,
(4.5)
which is
ω=
e3
R1 e1
⎧
⎫
! ⎨ φ̇ ⎬
R2 R1 e3
,
θ̇
⎩
⎭
ψ̇
which we in turn abbreviate to
ω = B Θ̇ .
(4.6)
If the second rotation θ is not kπ for some integer k, then the matrix B is invertible2 .
4.8
Alternative form for the rotation matrix
We begin with an observation (a theorem). See figure 4.3. Suppose there is a vertical rod, a rod welded to
it at an arbitrary angle, and a can through which the second rod passes. Imagine that we rotate this system
about the positive z axis by an angle θ; then rotate the can about the rotated-second-rod axis by an angle φ,
and then rotate the system back about the unrotated positive z-axis by an angle −θ. The net effect of this
rotation on the can is obviously the same as if there were no rotations about the vertical axis at all, and the
can was merely rotated by φ about the second rod.
2 The second column of B, R e , is perpendicular to e regardless of what φ is. R R e is simply R e because R e = e
1 1
3
2 1 3
2 3
1 3
3
(the first rotation is about e3 ). Since the rotation R2 is about R1 e1 , to which e3 is perpendicular, R2 e3 is perpendicular to R1 e1
as well. Finally, if the rotation angle θ of R2 is 0 or π, then R2 e3 equals e3 or −e3 respectively and B is singular; otherwise, e3
and R2 e3 are not parallel to each other, and also each perpendicular to R1 e1 , so B is invertible.
4.9. THE MOMENT OF INERTIA MATRIX ICM
47
(1) θ
(2) φ
(3) θ
Figure 4.3: Rotation theorem schematic.
The mathematical statement of this theorem is: for any rotation R̃ and any angle φ̃,
R̃T R(R̃n, φ̃)R̃ = R(n, φ̃) .
Using this result, we can derive an alternative form for the rotation matrix in Eq. 4.3, as follows.
R3 = R(R2 R1 e3 , ψ) ,
and so
(R2 R1 )T R3 R2 R1 = R(e3 , ψ) ,
or
R3 R2 R1 = R2 R1 R(e3 , ψ) .
But
R2 = R(R1 e1 , θ) ,
so
R1T R2 R1 = R(e1 , θ) ,
or
R2 R1 = R1 R(e1 , θ) = R(e3 , φ)R(e1 , θ) ,
and thus
Rf = R3 R2 R1 = R(e3 , φ)R(e1 , θ)R(e3 , ψ) .
(4.7)
In the right hand side of the above equation, the matrices appear in reverse sequence to what one expects
when one thinks of successive rotations! Moreover, the matrices are themselves somewhat simpler, which
eases symbolic manipulation. The above form of the rotation matrix can also be derived by considering not
rotating rigid bodies but rather the changing coefficients of fixed vectors in rotating coordinate systems (a
more common approach than adopted here). However, we prefer not to discuss rotating coordinate systems
except, for completeness, briefly at the end of this section.
4.9
The moment of inertia matrix Icm
What happens to the central moment of inertia matrix of a rigid body when the body is rotated? To see this,
consider the vector cross product sketched in Fig. 4.4. The vectors a and b, selcted arbitrarily, together define
c = a × b. These three vectors are imagined to be welded together to make a single rigid body. If that rigid
body is given an arbitrary rotation, it is clear that we will have “rotated a” × “rotated b” = “rotated c”.
This fact may be written in matrix notation as
c = S(a)b
CHAPTER 4. ROTATIONS AND EULER ANGLES
48
a
c=a xb
b
Figure 4.4: A vector cross product.
and
Rc = S(Ra)Rb,
or
RS(a)b = S(Ra)Rb.
Since b is arbitrary, this means
RS(a) = S(Ra)R
or
S(Ra) = RS(a)RT .
(4.8)
Since a is arbitrary as well, the above must be true for all a. We can use this identity in our calculation of the
moment of inertia for the rotated body. Recall, from Eq. 3.5, that for the unrotated body
Icm,ref =
mi S(ri/cm,ref )S T (ri/cm,ref ),
i
where the “ref” subscript denotes some reference position.
It follows that for the rotated body,
Icm =
mi S(Rri/cm,ref )S T (Rri/cm,ref ).
i
By Eq. 4.8,
Icm
=
i
= R
$
%T mi RS(ri/cm,ref )RT RS(ri/cm,ref )RT
=
mi RS(ri/cm,ref )S T (ri/cm,ref )RT
i
mi S(ri/cm,ref )S T (ri/cm,ref ) RT
i
= RIcm,ref RT .
(4.9)
This is an important result, useful for expressing angular momentum balance or rotational kinetic energy
in a general rotated configuration.
4.10
Angular accelerations
A relationship may be desired between the angular acceleration α and the second derivatives of the Euler
angles. To obtain this, we reexamine the intermediate rotating frames temporarily introduced in section 4.7.
We are interested in the time derivative of the angular velocity of F4 relative to the stationary frame F1 ,
evaluated with respect to F1 . For this, we use Eq. 2.5 to write
dω Fi /Fj
dω Fi /Fj
=
+ ω Fj /F1 × ω Fi /Fj .
dt
dt
F1
Fj
4.11. ROTATED COORDINATE SYSTEMS
49
Applying this result to Eqs. 4.4 and 4.5, we find
α = B Θ̈ + q ,
(4.10)
where q is quadratic in the angular rates and is given by
q = S(R1 e1 θ̇ + e3 φ̇)R2 R1 e3 ψ̇ + S(e3 φ̇)R1 e1 θ̇ .
Equation 4.10 can be used to solve for Θ̈, once α has been obtained from the Newton-Euler equations.
4.11
Rotated coordinate systems
What happens to the coordinates of a vector if the vector stays fixed while the coordinate axes themselves are
rotated like a rigid body? Any rotation in whatever Euler angle sequence is finally equivalent to a rotation
about some n through some θ, so imagine that the coordinate axes are rotated like a rigid body through θ
about n. The effect on a fixed vector r is that its matrix of components in the unrotated basis, r, changes to
a different matrix of components in the rotated basis, r∗ .
Now the change in the relationship between r and the basis is the same whether we (1) keep r fixed and
rotate the basis about n through θ, or (2) keep the basis fixed and rotate r about n through −θ. Thus,
r∗ = R(n, −θ) r = [R(n, θ)]T r .
We will not use the matrix representation of axis rotations in our treatment.
4.12
Exercises
1. Consider a rigid body with a vector p fixed in it. When in the reference configuration, p = î − ĵ (note: p is
not a unit vector). The body is first rotated by 30 degrees about the x axis. Then it is rotated by 45 degrees
about the new position of the vector p. Finally, it is rotated by 20 degrees about the space-fixed z-axis. Find
the net rotation matrix. Also, find the axis of the net rotation, and the angle of rotation (these quantities are
related to the real eigenvector and the complex eigenvalues of the rotation matrix).
√
√
√
2. A rigid body is rotated through 0.6 radians about the unit vector 1/ 2 î − 1/ 3 ĵ + 1/ 6 k̂. During this
rotation, the origin stays fixed. What is the new position of a point on the body which was at 4î − 3ĵ + k̂
before the rotation? If the moment of inertia matrix of this rigid body, about its center of mass, and in the
initial configuration, was a diagonal matrix with the diagonal elements {1, 1, 2}, what is its moment of inertia
matrix in the rotated configuration?
3. Consider again a rigid body with one point fixed at the origin. The orientation of the body is described
using a time-varying unit vector n̂(t) along with
of θ(t)
a√rotation √
about that unit vector. At some instant of
√
√
√
˙
interest, n̂ = 1/ 2 î − 1/ 3 ĵ + 1/ 6 k̂, n̂ = 1/ 3 î + 1/ 2 ĵ s−1 , θ = 0.6 radians, and θ̇ = 0 radians/s.
What is the angular velocity of the rigid body at that instant?
4. Consider the net rotation matrix R of exercise 1 above. Find a choice of (3,1,3) Euler angles (φ, θ, ψ)
corresponding to this net rotation.
Suggestion: You can begin with the (3,3) element of R to find two tentative choices for θ (between 0 and 2π).
Then the third row gives you two corresponding choices for ψ, while the third column gives you two corresponding choices for φ. These two sets of possible solutions must be verified against the remaining elements
of R.
5. Consider again exercise 1 above. Instead of the given rotations of 30, 45 and 20 degrees, let these rotations
be arbitrary; call them φ, θ and ψ respectively. Write a computer program for finding a set of these angles to
match a given rotation matrix R = I.
Suggestion: You know the net rotation and its axis. Imagine the body rotates with a constant ω about this
axis over unit time. Seek a relation of the form ω = B Θ̇. Solve these ODEs with zero initial conditions and
over a unit time interval. Any ideas for singular configurations?
50
CHAPTER 4. ROTATIONS AND EULER ANGLES
Chapter 5
Lagrange’s Equations
We have so far been discussing Newton-Euler mechanics, where one draws free body diagrams showing external
forces and moments, and then directly uses the equations of linear and angular momentum balance.
A completely different approach which has some powerful advantages for certain problems is that of analytical mechanics, which leads to Lagrange’s equations and beyond. We will stop at Lagrange’s equations.
5.1
Virtual work for a system of particles
We start with the equations of motion written for each of the N individual particles in a system. For the ith
particle, using D’Alembert’s principle:
(5.1)
−mi r̈i + F ai + f i = 0 ,
where the second term on the left hand side is the sum of external forces on the particle, and the third term
represents the resultant of all constraint forces on the particle.
5.1.1
Applied versus constraint forces
The above split of forces into applied forces and constraint forces implies that we are able to distinguish
between these two types of forces. We usually are. For example, if two point masses are connected by a taut
string, then the tension in the string provides a constraint preserving force acting on each point mass, while
the weights of these masses are externally applied forces. For another example, consider a block sliding on a
frictionless inclined plane; the normal reaction from the plane is a constraint force while any other forces (such
as weight) are external applied forces. In general, the constraint forces are not directly known; rather their
net effect on the motion is known, and this effect may in principle be used to calculate what these constraint
forces are (although, as we will see, the aim of Lagrange’s equations is to not have to calculate these constraint
forces at all).
In some cases, things may be less clear. For example, a block sliding on a frictional inclined plane experiences both normal and tangential forces from the plane. It is not clear whether the tangential or frictional
force should be treated as a constraint force. It certainly is not specified in advance, and needs to be calculated
indirectly from the motion of the system: by these qualities, it is like a constraint force. Yet it acts in the
tangential direction, while the constraint is on motion in the normal direction: by this quality, perhaps it is not
a constraint force? As will be seen below, this question will become irrelevant in the Lagrangian formulation,
because such frictional forces will be excluded from the list of constraint forces and will have to be included
as applied forces.
Returning to Eq. 5.1, in any virtual displacement of the system (where we use virtual to mean imagined),
F ai − mi r̈i + f i · δri = 0 for i = 1, 2, · · · , N
and so
N i=1
F ai − mi r̈i + f i · δri = 0 .
51
CHAPTER 5. LAGRANGE’S EQUATIONS
52
5.1.2
Net zero virtual work of constraint forces
Now we assume that the net virtual work of the constraint forces is zero in any virtual displacement that
obeys the constraints on the system. This statement requires attention on two points.
real displacement
virtual displacement
specified motion
Figure 5.1: A block slides on an inclined plane. The plane itself has a specified motion.
First, the virtual displacements are imaginary, and considered to be instantaneous. This means that any
time-dependent constraints on the system are held frozen, along with the time, while we imagine the virtual
displacement. For example, if a block slides on an inclined plane that itself oscillates in a prescribed manner,
then during consideration of virtual displacements at some instant t, the inclined plane is considered held fixed
at its instantaneous location. Thus, all virtual displacements of the block that obey the constraint are parallel
to the inclined plane; while, assuming the plane has nonzero real velocity at instant t, all real infinitesimal
displacements of the block are in fact not parallel to the inclined plane (see figure 5.1).
δx
T
δx
x
T
Figure 5.2: Two masses held together by a taut string move along the x axis.
The second point to consider is that only the net work, and that of only the constraint forces, is zero.
For example, if two point masses held together by a taut string move along the x-axis, then in a virtual
displacement of δx to the right, the constraint forces (string tension T ) do virtual work equal to T δx on the
left side mass, virtual work equal to −T δx on the right side mass, and net zero virtual work (see figure 5.2).
Another example worth considering is that of two rods in a plane, connected at a frictionless hinge. The hinge
connection gives rise to equal and opposite forces on these rods; and the net virtual work done in a constraint
obeying virtual displacement is zero.
The above assumption disallows us from considering the frictional components of contact forces as constraint forces, because they do net virtual work in constraint obeying virtual displacements; these frictional
forces must be included in the formulation as external applied forces. Usually, the normal components will
then be treated as applied forces also.
Accepting the assumption on the basis of its clear utility for at least many frictionless systems, we find
that for constraint-satisfying virtual displacements,
N
(F ai − mi r̈i ) · δri = 0 .
(5.2)
i=1
In the above equation, the individual coefficients F ai − mi r̈i are not necessarily zero because the δri are not
independent (they must satisfy the constraint equations).
The above equation can be systematically exploited to obtain a relatively small number of useful equations
of motion for the system. To do that, we need to use generalized coordinates.
5.2. GENERALIZED COORDINATES AND DEGREES OF FREEDOM
5.2
53
Generalized coordinates and degrees of freedom
For a system of particles with possibly many constraints, we assume there are n quantities or numbers qi ,
i = 1, 2, · · · , n, which we call generalized coordinates, such that the equations
rj = rj (q1 , q2 , · · · , qn , t)
(5.3)
are invertible. By this we mean that given the q’s we can find each r and vice versa.
For example, if a bead moves on a circular wire fixed in the x-y plane, then its position vector r has two
nonzero components along x and y directions respectively. For this one-particle system, possible choices of
generalized coordinates are x and y themselves (with n = 2), or an angular position θ along the circular wire
(with n = 1), or an arc length s along the wire (with n = 1), and so on.
The smallest possible number n for which such generalized coordinates can be found is called the number
of degrees of freedom of the system. A few examples follow.
The bead on the circular wire has n = 1 degree of freedom. A bead on a circular wire whose vertical
diameter coincides with the z-axis, and which spins about the z-axis at an externally specified rate, also has
only 1 degree of freedom.
A point mass constrained to move on the surface of a sphere has 2 degrees of freedom. A double pendulum,
i.e., a system where one pendulum is suspended from another which is itself suspended from a fixed support,
also has 2 degrees of freedom.
A block sliding on a table has 3 degrees of freedom (two for translation and one for rotation). A rigid body
with one point attached to ground using a ball and socket joint also has 3 degrees of freedom (say, three Euler
angles).
A block that slides on a cart, with the cart itself free to move along a fixed path, has 4 degrees of freedom
(one for the cart, measuring distance along the path; and three for the block’s motions relative to the cart).
Two point masses in space held together by a taut string have 5 degrees of freedom (say, three components
of the position vector of the first mass; and two angles to describe the orientation of the string).
A single unconstrained rigid body has 6 degrees of freedom (e.g., three for the location of the center of
mass; and three Euler angles). Two point masses in space, connected by a spring, also have 6 degrees of
freedom.
One of the great advantages of using generalized coordinates is that they allow approximate modeling. For
example, a vibrating beam has infinitely many points in it; and each of these points has three translational
degrees of freedom. Thus, the beam has infinitely many degrees of freedom. We can assume as an approximation, however, that originally plane sections normal to the undeformed neutral axis remain plane and normal
to the neutral axis even in the deformed state. We can further assume that the deformed shape of the neutral
axis can be described to sufficient accuracy by a time dependent cubic polynomial with two free coefficients,
i.e.,
w(x, t) ≈ a0 (t) x2 + a1 (t) x3 ,
where x is the position of a point on the initially undeformed neutral axis, and w(x, t) is the time dependent
lateral deflection of that point. The above approximation gives rise to a 2 degree of freedom system, in which
the coefficients a0 and a1 are generalized coordinates.
5.3
Constraints
We may encounter systems where the generalized coordinates must satisfy some constraints during the motion.
For example, for a bead on a circular wire in the x-y plane, if we use x and y as generalized coordinates, then
these coordinates must satisfy the constraint equation
x2 + y 2 = R 2 ,
where R is the radius of the circle. Such constraints will always arise when we use more generalized coordinates
(here, two) than the system has degrees of freedom (here, one). For another example, consider a point mass
constrained to move on the surface of a sphere, where we use 3 generalized coordinates x, y and z (components
of the position vector of the point mass). Now the constraint equation will be of the form
x2 + y 2 + z 2 = R2 .
CHAPTER 5. LAGRANGE’S EQUATIONS
54
If the radius of the sphere is changed by some external agency to equal some specified function of time, then
the above constraint equation will become
x2 + y 2 + z 2 = R2 (t).
From the above, it is clear that a useful and reasonably broad class of constraint equations can be expressed
in the form
(5.4)
g(q1 , q2 , q3 , · · · , qn , t) = g(q, t) = 0,
where we use the abbreviated notation
q = {q1 , q2 , q3 , · · · , qn }T .
Such constraints are called holonomic. (If t or an explicit function thereof appears in the constraint equation,
it is called rheonomic. Otherwise, it is called scleronomic. These two terms will not be important in our
study.)
In principle, at least, every holonomic constraint can be used to solve for one of the generalized coordinates
in terms of the others; this generalized coordinate may then be eliminated from our consideration, giving a
system with fewer generalized coordinates. For systems which have only holonomic constraints, this process
may in principle be continued until we have a minimal set of generalized coordinates and no remaining
constraints. The study of unconstrained, n degree of freedom systems is an important part of analytical
mechanics.
However, many important constraints are not holonomic. A point mass constrained to be outside a sphere
may satisfy a constraint of the form
x2 + y 2 + z 2 ≥ R2 ,
which does not match Eq. 5.4. We will not consider such inequality constraints, except to say here that they
can be handled by dividing the motion into phases where the constraint is either active (with equality) or
inactive (with strict inequality), along with additional rules for transitions between these two phases.
In systems with friction, we may be unable to eliminate the constraint and reduce the number of generalized
coordinates because the constraint force cannot be eliminated in the calculation of net virtual work. For these
systems, we may need to carry the constraint equation along even though it is holonomic.
A particular type of nonholonomic constraint important for mechanical systems places restrictions on the
rates of change of the generalized coordinates. We consider the case where m possibly nonholonomic constraint
equations (m < n) are of the form
n
aij (q, t)q̇j + ait (q, t) = 0 for i = 1, 2, . . . , m ,
(5.5)
aij (q, t)dqj + ait (q, t)dt = 0 for i = 1, 2, . . . , m .
(5.6)
j=1
which is equivalent to
n
j=1
Finally, in a virtual displacement (where time is frozen), we will have
n
aij (q, t)δqj = 0 for i = 1, 2, . . . , m .
(5.7)
j=1
Obviously, differentiating Eq. 5.4 can give a constraint of the form of Eq. 5.5. Such an apparentlynonholonomic constraint is actually a holonomic constraint in disguise, and can be integrated to find the
holonomic constraint, which can then be used to eliminate a generalized coordinate (at least in principle).
However, there are nonholonomic constraints in the form of Eq. 5.5 which are not obtainable by differentiating
holonomic constraint equations, and cannot be used to eliminate generalized coordinates.
To see this, we proceed with an example called the “skate” constraint. See figure 5.3. The figure shows
a single elliptical flat object sliding on a plane. There is a microscopic knife edge embedded in the object at
5.4. HOLONOMIC SYSTEMS
55
e^
n^
θ
P
x
y
Figure 5.3: The skate constraint.
point P . This knife edge is shown aligned with the major axis of the ellipse, which is along unit vector ê. The
skate constraint is
v P · n̂ = sin θ ẋ − cos θ ẏ = 0.
Note that rotation about point P is allowed.
To show that the skate constraint is not merely the differentiated form of some holonomic constraint, we
will argue by contradiction. Accordingly, let us assume that the constraint is indeed a differentiated holonomic
constraint; then, of the three generalized coordinates x, y and θ, it must be possible to eliminate one in terms
of the other two. However, by first rotating about P so that ê is aligned with the x axis; and then changing
x by an arbitrary amount; and then rotating back about P to the original orientation, we see that arbitrary
values of x are possible for any given values of y and θ. Thus, x cannot be eliminated. By a similar argument,
y cannot be eliminated either. And freedom to rotate about P , of course, means that θ cannot be eliminated.
Since no generalized coordinate can be eliminated, our initial assumption must be false. The skate constraint
is not a differentiated holonomic constraint.
An important point in the above argument is that the constraint equation is a statement of geometrical
or kinematic restrictions only. Mass distributions and the laws of dynamics do not enter into considerations
of whether or not a given constraint is holonomic. Only the geometrical restrictions on the motion can be
probed to establish that a given constraint is nonholonomic.
Nonholonomic constraints of the form of Eq. 5.5 restrict the instantaneously accessible set of velocities
without managing to restrict the globally accessible set of geometrical configurations. This provides a physical
interpretation of such nonholonomic constraints.
In mechanical systems, there are two main nonholonomic constraints (of the form of Eq. 5.5) that are well
understood: one is the skate constraint described above, and the other is the constraint of rolling without slip
(in, e.g., a rolling coin). Other constraints on velocities will not be considered here.
5.4
Holonomic systems
We can now seek the equations of motion of systems with only holonomic constraints, such that all extra
generalized coordinates have been eliminated, and we have a minimal set of them with no further constraints
on them.
Recall Eq. 5.2. Now we move from r’s to q’s using some standard calculations. Using Eq. 5.3, we can write
δri =
n
∂ri
δqj ,
∂q
j
j=1
CHAPTER 5. LAGRANGE’S EQUATIONS
56
where we note that time is frozen and so time derivatives of r(q, t) do not appear. This gives
N
(F ai − mi r̈i ) ·
i=1
n
∂ri
δqj = 0.
∂qj
j=1
(5.8)
The generalized force Qj corresponding to qj is now defined. Consider in Eq. 5.8 the quantity
N
i=1
⎡
⎣F ai
·
n
∂f
j=1
⎤
i
∂qj
δqj ⎦ =
N n
F ai
·
i=1 j=1
where
Qj =
∂f i
∂qj
N
i=1
δqj =
N
n
j=1
F ai ·
i=1
F ai
·
∂f i
∂qj
δqj =
n
Qj δqj ,
j=1
∂ri
.
∂qj
(5.9)
Qj is called the generalized force corresponding to the j th degree of freedom. Though the formal definition
of Qj involves a sum over all N particles in the system (where N might be very large), for most practical
problems Qj is in fact easy to calculate. For example, consider figure 5.4, which shows a frictionlessly hinged
rod acted upon by a force at the middle as well as a spring at the free end; gravity is neglected. For simplicity,
F
L
θ
40
L /2
k
Figure 5.4: A hinged bar with a spring support and an external force.
the deflection θ is supposed to be small so that all motions are vertical. The constraint force at the frictionless
hinge does no work because its point of application does not move. The virtual displacement of the midpoint
L
of the rod is Lδθ/2 in the vertical direction, so the virtual work done by the applied force F is F sin(40◦ ) · δθ.
2
At the free end, given some small displacement θ, the deflection is Lθ and the spring force is therefore kLθ
acting downwards. In a further infinitesimal virtual displacement, the motion of the free end is Lδθ upwards,
and so the virtual work done by the spring force is −KL2 θδθ. On all other particles that make up the rod,
the only forces acting are constraint forces, whose virtual work need not be included because their net virtual
work is zero. Thus, the net virtual work done by the external applied forces (in which we include the force F
and the spring force) is
L
F sin(40◦ ) − KL2 θ δθ.
2
Therefore, comparing with Eq. 5.9,
Q=F
L
sin(40◦ ) − KL2 θ.
2
If gravity is included and the net weight of the rod (assumed uniform) is W , then the additional virtual work
done by gravity is
L
−W δθ,
2
5.4. HOLONOMIC SYSTEMS
57
where we take the weight to act through the center of gravity of the rigid rod (eliminating the need to consider
the individual weights of all the particles that make up the rod). Including the gravity, then, we have
Q=F
L
L
sin(40◦ ) − W − KL2 θ.
2
2
This example shows how computing the generalized forces Qj is usually not as difficult in practice as the
formal definition might suggest.
Moving on, consider in Eq. 5.8 the quantity
⎛
⎞
N
n
∂ri
⎝mi r̈i ·
δqj ⎠ .
∂q
j
i=1
j=1
We can rewrite the above as
N
n
∂r
mi r̈i · i
∂qj
i=1
j=1
δqj ,
where the quantity inside brackets is again rewritten as
N
mi
i=1
dṙi
·
dt
∂ri
∂qj
.
Using the identity
da
d
db
·b=
(a · b) − a · ,
dt
dt
dt
we rewrite the above as
N
mi
i=1
d
dt
ṙi ·
∂ri
∂qj
− ṙi ·
d
dt
∂ri
∂qj
.
(5.10)
Due to the special form of ri (see Eq. 5.3), we have
ṙi =
n ∂r
j=1
i
∂qj
q̇j
+
∂ri
,
∂t
whence
∂ ṙi
∂r
= i,
∂ q̇j
∂qj
where the partial derivative with respect to q̇j is taken while holding fixed all the other q̇’s, as well as all the
q’s and t. Substituting the above in Eq. 5.10 and bringing the mi (independent of time) inside the derivatives,
we obtain
N d
∂ ṙi
d ∂ri
− mi ṙi ·
.
(5.11)
mi ṙi ·
dt
∂ q̇j
dt ∂qj
i=1
We now note that
d
dt
because
∂ri
∂qj
n
∂ 2 ri
∂ 2 ri
q̇k +
,
=
∂qk ∂qj
∂t∂qj
k=1
∂ri
is independent of the q̇’s. Rearranging terms,
∂qj
d
dt
∂ri
∂qj
∂
=
∂qj
n
∂ri
∂r
q̇k + i
∂qk
∂t
k=1
,
CHAPTER 5. LAGRANGE’S EQUATIONS
58
where the partial derivative operator can be pulled out while keeping the q̇k inside because the partial derivative
is to be taken with the q̇’s being held constant anyway. Finally, recognizing the quantity within brackets on
the right hand side, we write
∂ ṙ
d ∂ri
= i.
dt ∂qj
∂qj
Substituting the above into Eq. 5.11, we get
N d
∂ ṙi
∂ ṙi
− mi ṙi ·
,
mi ṙi ·
dt
∂ q̇j
∂qj
i=1
which is rearranged one last time to finally write
⎛
⎞
N
N
n
N
∂r
d
∂
ṙ
∂ ṙ
i
i
⎝mi r̈i ·
δqj ⎠ =
mi ṙi ·
mi ṙi · i .
−
∂q
dt
∂
q̇
∂q
j
j
j
i=1
j=1
i=1
i=1
(5.12)
Now we realize that if the kinetic energy of the system is written as
N
1
T =
i=1
then Eq. 5.12 is simply
2
mi ṙi · ṙi ,
⎞
n
∂ri
d ∂T
∂T
⎝mi r̈i ·
δqj ⎠ =
.
−
∂q
dt
∂
q̇
∂q
j
j
j
i=1
j=1
N
⎛
(5.13)
Putting Eqs. 5.13 and 5.9 into Eq. 5.8, we obtain
n j=1
d
Qj −
dt
∂T
∂ q̇j
∂T
+
∂qj
δqj = 0
for all virtual displacements that obey the constraints on the system. But the system, in this form, has no
constraints! All the n generalized coordinates can be varied arbitrarily. In fact, we may first let δq1 = 0 with
all others zero; then imagine that δq2 = 0 with all others zero; and so on, so that we can conclude
d ∂T
∂T
= Qj , for j = 1, 2, · · · , n.
(5.14)
−
dt ∂ q̇j
∂qj
The above are called Lagrange’s equations of motion.
Suppose now, in addition to our prior assumptions, that the generalized forces may be split up into two
parts,
p
Qj = Qnc
j + Qj ,
where the nc-superscript denotes “non-conservative” and can represent any forces at all, while the p-superscript
denotes “potential” and represents some special types of forces that we have the option of treating differently.
These potential forces are assumed to be derivable from some scalar potential function V (q, q̇, t) through the
equation
∂V
.
Qpj = −
∂qj
The equations of motion then become
d
dt
∂T
∂qj
−
∂
(T − V ) = Qnc
j .
∂qj
5.4. HOLONOMIC SYSTEMS
59
In mechanical systems in particular, such functions V arise whenever there is strain energy (including in
springs) and gravitational potential energy. Finally, if (as will be the case for every mechanical system
considered in our study) V = V (q, t) does not depend on the q̇’s, so that
∂V
= 0,
∂ q̇j
then we define the Lagrangian as
L = T − V,
and Lagrange’s equations become
d
dt
∂L
∂ q̇j
−
∂L
= Qnc
j ,
∂qj
j = 1, 2, · · · , n.
for
Example: Consider a simple pendulum of mass m and length L as shown in figure 5.5. This system has
θ
g
L
m
Figure 5.5: A simple pendulum of length L and end mass m.
one degree of freedom. The angle θ measured from the vertical is chosen as the generalized coordinate. The
kinetic energy of the mass is
1
1
T = mv 2 = m(Lθ̇)2 .
2
2
The potential energy of the mass is
V = −mgL cos θ .
The Lagrangian is
L=T −V =
Lagrange’s equation is
d
dt
1
mL2 θ̇2 + mgL cos θ .
2
∂L
∂ θ̇
−
∂L
= 0,
∂θ
or
mL2 θ̈ + mgL sin θ = 0.
Simplifying,
θ̈ +
g
sin θ = 0.
L
CHAPTER 5. LAGRANGE’S EQUATIONS
60
5.5
Nonholonomic systems
Sometimes we might wish to treat actually dependent q’s as independent coordinates, with added constraint
equations. Added constraint equations are in any case unavoidable in the presence of nonholonomic constraints.
We consider here only the special case where m possibly nonholonomic constraint equations (m < n) are of
the form given by Eqs. 5.5 through 5.7. How, and under what further assumptions, can we write equations of
motion for such systems?
Consider the following strategy:
1. Ignore the constraints, and treat the q’s as independent coordinates.
2. Apply some new (extra) generalized forces Qcj to the system.
3. Choose Qcj so that the motion satisfies the constraints.
The above strategy does not provide a unique solution for the constraint forces Qcj . A simple example shows
why. Consider a block constrained to slide frictionlessly on a horizontal surface. To any system of constraint
forces Qcj that leads to one constraint-obeying solution, we can further add on generalized forces equivalent
to, say, an arbitrary horizontal force on the block; this will lead to a different, but also constraint-obeying,
motion of the block.
So, to items (1) through (3) above, we add on a final assumption:
4. In any virtual displacement that obeys Eq. 5.7, the net virtual work of the forces Qcj is zero, i.e,
n
Qcj δqj = 0 whenever Eq. 5.7 holds.
(5.15)
j=1
Item (4) above may initially seem like a strong assumption. However, for holonomic constraints, we have
assumed exactly this already. For nonholonomic constraints like ideal rolling and the constraint of a skate,
the constraint forces act at the contact point, and do no net virtual work in displacements that satisfy Eq.
5.7. Notice that in our strategy (item (1) above) the system is allowed to have virtual displacements that do
not satisfy Eq. 5.7, and for these virtual displacements the net virtual work of the forces Qcj may be nonzero.
We now proceed as follows. We think of δq [n×1] = {δq1 , δq2 , · · · , δqn }T as a vector in n dimensional space
and rewrite Eq. 5.7 in matrix form as
A[m×n] δq = 0[m×1] ,
(5.16)
where the bracketed subscripts indicate the size of each matrix. At any particular instant of time, t is frozen
(constant) while we consider virtual displacements δq. The actual coordinate vector q is also constant during
this process. Thus, though the matrix A changes as time progresses and q evolves, at each instant of time A
is held constant as we consider a variety of imaginable virtual displacements.
The rows of A, or equivalently, the columns of AT , are m vectors in n dimensional space. The span (or set
of all linear combinations) of these m vectors is a subspace of n dimensional space; let us call it SA . The set of
all vectors that are perpendicular to the rows of A (or the columns of AT ) lie in another subspace orthogonal
to SA , which we call SA⊥ . Any vector in our n dimensional space can be split into two perpendicular or
orthogonal components: one in SA and the other in SA⊥ .
A simple example may help to fix the above abstract ideas in the reader’s mind. Imagine that we are
working in our familiar three dimensional space (i.e., n = 3), and AT has 2 columns. The vectors represented
by these 2 columns are (say) î and î + ĵ. The set of all linear combinations of these two vectors covers the
entire x-y plane, and so SA is the x-y plane. All vectors perpendicular to the columns of AT are along k̂, and
so for this example SA⊥ is the z-axis. Finally, any vector in 3 dimensions can be split into two perpendicular
or orthogonal components: one in the x-y plane and the other along the z axis.
Now any vector δq that satisfies Eq. 5.16 must lie completely in SA⊥ . And any candidate δq chosen from
SA⊥ will satisfy Eq. 5.16.
5.5. NONHOLONOMIC SYSTEMS
61
Let the unknown constraint-enforcing generalized forces be arranged in the column matrix (or vector)
Qc[n×1] = {Qc1 , Qc2 , · · · , Qcn }T .
What can we say about the n dimensional vector Qc ? Simply that
(Qc )T δq = 0 whenever δq is in SA⊥ .
In other words, Qc is perpendicular to every vector in SA⊥ . It must therefore be in SA . In other words, Qc is
a linear combination of the columns of AT (equivalently, the rows of A), i.e.,
Qc = λ[1×m] A ,
where λ = {λ1 , λ2 , · · · , λm } represents as yet undetermined coefficients or parameters (note that there are as
many λk ’s as constraint equations, which is fewer than the number of degrees of freedom). Thus,
Qcj =
m
λk akj (q, t) ,
k=1
and we will have to find out what the λk ’s need to be so that the system obeys the nonholonomic constraints.
Finally, the equations of motion are
d
dt
∂L
∂ q̇j
−
∂L
c
= Qnc
j + Qj
∂qj
n
aij q̇j + ait
= Qnc
j +
m
λk akj ,
for j = 1, 2, · · · , n ,
(5.17)
k=1
=
0,
for i = 1, 2, . . . , m ,
(5.18)
j=1
where in the first equation we recall that (1) all potential forces are included implicitly in L = T − V , (2) all
other externally applied forces not arising from the potential V are included in Qnc
j , as given by Eq. 5.9, and
(3) the constraint forces that the system experiences which make it obey the constraints of Eq. 5.5 (be they
differentiated-holonomic or nonholonomic) are given in terms of m unknown λk ’s.
There are now n + m equations and n + m unknown quantities (n second derivatives and m λk ’s) to be
determined. In analytical work, the constraint equations will usually be differentiated with respect to time so
as to yield equations involving the second derivatives.
The λk ’s are usually equal to the constraint forces or scalar multiples of the same.
Example: Consider a skate on a frictionless horizontal plane with a microscopic knife edge at P as shown in
figure 5.6. The mass of the skate is m. The Izz component of the moment of inertia about the center of mass
G is J. The distance from the center of mass G to the knife edge at P is L. The knife edge itself is aligned
with rP/G and unit vector ê.
This is a 3 degree of freedom system. We choose as generalized coordinates the x and y coordinates of the
center of mass, and angle θ that ê makes with the horizontal.
The knife edge constraint is
v P · n̂ = 0.
Now
v P = v G + ω skate × rP/G = ẋî + ẏ ĵ + θ̇k̂ × Lê = (ẋ − L sin θ θ̇) î + (ẏ + L cos θ θ̇) ĵ,
so the nonholonomic constraint is
sin θ ẋ − cos θ ẏ − Lθ̇ = 0,
(5.19)
sin θ dx − cos θ dy − Ldθ = 0.
(5.20)
or equivalently
Since there is a single scalar constraint equation, there will be only one constraint force parameter λ.
CHAPTER 5. LAGRANGE’S EQUATIONS
62
e^
n^
L
P
θ
G
m,J
y
x
Figure 5.6: A skate of mass m and moment of inertia J about the center of mass.
The kinetic energy of the skate is (see section 3.7)
T =
1
1 2
m ẋ + ẏ 2 + J θ̇2 .
2
2
There is no potential energy term, i.e., V = 0.
Lagrange’s equations (incorporating Eq. 5.20) are
J θ̈ = −λL
mẍ = λ sin θ
mÿ
= −λ cos θ ,
where λ remains to be determined. It may be seen from the equations of motion that λ is the magnitude of a
constraint force that acts on the skate at point P , in the direction of n̂.
The above three equations are augmented with Eq. 5.19, which on differentiation will yield a fourth equation
involving ẍ, ÿ and θ̈.
5.6
The calculus of variations
We will now approach Lagrange’s equations from a different direction.
To begin, consider some unknown function y(x) defined on the unit interval. It is given that y(0) = y(1) = 1.
Of all the possible smooth choices of function y(x), which one minimizes
+
1 * 2
dy
2
H=
+y
dx ?
dx
0
Does the problem even have a meaningful solution? Those who have never encountered the calculus of
variations, and do not know the problem can in fact be solved, might like to think about the problem a little.
First note that if y ≡ 0 then H = 0, but the boundary conditions are violated. On the other hand, if y ≡ 1,
then boundary conditions are met and H = 1. This is not the minimum possible value for H, because a quick
1
calculation shows that for y = 1 − sin πx, H ≈ 0.93 < 1.
7
Can we find an approximate solution? Suppose we assume that
y =1+
5
k=1
ck sin kπx,
5.6. THE CALCULUS OF VARIATIONS
63
for some as yet undetermined coefficients ck , then the boundary conditions are satisfied and (from MAPLE)
H
=
4 c5
1
25 2 2
c1
4 c3
1 2
c5 + 4
+
+ c4 2 +
+
π c5 + 2 π 2 c2 2
2
π
5 π
2
3 π
2
1
9
1
1
1
+ c2 2 + π 2 c3 2 + c1 2 + 8 π 2 c4 2 + π 2 c1 2 + c3 2 .
2
2
2
2
2
1+
Taking partial derivatives with respect to the ck ’s and setting them equal to zero, we get 5 equations; solving
them, we obtain
c1 =
−4
−4
−4
, c2 = 0, c3 =
, c4 = 0, c5 =
,
π (1 + π 2 )
3π (1 + 9 π 2 )
5π (1 + 25 π 2 )
which gives H ≈ 0.924294. The series converges, and the form of the k th term is easy to guess. We could use
more terms and get a better approximation.
But can we get the exact solution by a more direct method? The answer lies in the calculus of variations.
Say we wish to find a function y that satisfies y(a) = ya , y(b) = yb , and minimizes (more correctly, extremizes)
some functional
b
H̄ =
f (x, y, y ) dx.
a
This means that H̄ is insensitive, to first order, to small changes in y. Let
Y (x) = y(x) + η(x),
where y(x) is the as yet unknown minimizing function, η(x) is an arbitrary variation, and is a smallnessenforcing parameter. Note that the variation η must satisfy η(a) = η(b) = 0, because at these points the
values of y are specified.
We now have
b
b
b
∂f
∂f
H̄ =
η + η dx + O(2 ),
f (x, Y, Y ) dx =
f (x, y, y ) dx + ∂y
∂y
a
a
a
where the O(2 ) denotes terms of higher order in .
At minima, as indicated above, we expect functions to be insensitive (up to first order) to arbitrary
variations in the inputs. For H̄ to be a minimum, therefore, we require that
b
∂f
∂f
η + η dx = 0
(5.21)
∂y
∂y
a
for arbitrary η. We now note that
a
so Eq. 5.21 becomes
b
0 b
>
∂f
∂f d
b
η dx = η dx,
y η|a −
∂y ∂y a dx
a
b
#
∂f
d
−
∂y
dx
∂f
∂y ,
η dx = 0
(5.22)
for arbitrary η satisfying η(a) = η(b) = 0.
Now we use something called the fundamental lemma of the calculus of variations (FLCV). Consider a
continuous function g(x) defined on the interval (a, b). If we find that
b
g(x)η(x) dx = 0
a
for all η satisfying η(a) = η(b) = 0, then the FLCV states that g(x) ≡ 0 in the interval (a, b). Why? If
g(x) is nonzero (and, say, positive) at any point inside the interval, then by its continuity it is nonzero (and
CHAPTER 5. LAGRANGE’S EQUATIONS
64
positive) in some small subinterval containing that point. Choosing an η that is nonzero and positive inside
that subinterval and zero everywhere else, we find that that integral is strictly positive. That it may be quite
a small number is irrelevant: it is nonzero. In other words, if g(x) = 0 anywhere inside the interval, then there
is some η for which
b
g(x)η(x) dx = 0.
a
This proves the FLCV.
By the FLCV, then, the fact that Eq. 5.22 holds for arbitrary η implies that
d
∂f
∂f
−
= 0.
∂y
dx ∂y (5.23)
It is easy to show that if f depends on more than one unknown and independently variable function, i.e.,
f = f (x, y1 , y1 , y2 , y2 , · · · ), then Eq. 5.23 holds with y replaced separately and successively by y1 , y2 , · · · .
Applying the above to our originally stated problem, where f = y 2 + y 2 , Eq. 5.23 gives
2y − 2y = 0,
subject to y(0) = y(1) = 1. The solution is (from MAPLE)
−1 + e−1 ex
(e − 1) e−x
y=
+
,
e − e−1
e − e−1
whence
H=2
e−1
≈ 0.924234.
e+1
Comment: Popular problems used to demonstrate the calculus of variations in dynamics textbooks include
finding the shortest path between two points on a plane (a straight line), finding the path in a vertical plane
down which a frictionlessly sliding particle will move fastest (called the brachistochrone problem), and finding
the shape of a chain suspended at two ends from two given supports (the shape is called a catenary). Those
classical problems are worth studying. The reader may also find it amusing to guess the series solution to the
present problem, evaluate it at x = 1/2, and compare the sum to the exact solution given by the calculus of
variations.
5.7
Hamilton’s principle
Equation 5.23 has an obvious relevance to Lagrange’s equations for systems without added constraints and
without non-conservative forces Qnc . For such systems, we now see that the action integral
t2
t1
L dt
has minimal (or more correctly, an extremal) value. This is called the principle of least action and also
Hamilton’s principle.
The principle of least action is important in several ways.
First, it is philosophically interesting and aesthetically pleasing because it provides a different starting
point for obtaining equations of motion for dynamic systems. Instead of accepting F = ma for particles
(Newton’s Second Law) and using the principle of virtual work, we could accept the principle of least action
as the basis for mechanics. More can be made of this than engineering problems warrant, however, because
the principle is not valid for systems with non-conservative forces and/or nonholonomic constraints, while the
usual route to equations of motion can handle such systems. Moreover, as engineers, we trust Newton’s Laws
because they match experiments; and we trust Lagrange’s equations because they come from Newton’s Laws;
and we accept the principle of least action because of its equivalence, via Lagrange’s equations, to Newton’s
5.8. EXERCISES
65
Laws. If the principle of least action did not lead to the same equations of motion as Newton’s Laws, we would
drop aesthetics in favour of empirical truth.
The principle of least action is also important because it opens the door to writing down useful governing
equations for systems involving several different forms of energy, such as mechanical systems interacting with
electrical and magnetic elements. Although important, this aspect is not emphasized in our study, which
focuses on mechanical systems.
Finally, and maybe most importantly in the study of mechanical systems, the principle of least action
and the calculus of variations can together be used to obtain the partial differential equations of motion of
systems with infinitely many degrees of freedom. This aspect is not emphasized in these notes: read a book
on structural vibrations.
5.8
Exercises
1. Write Lagrange’s equation of motion for the simple pendulum (length L, mass m, acceleration of gravity
g) with a force F (t) acting on the point mass. The force acts at an angle of π/4, measured counterclockwise
from the horizontal.
2. A point mass m moves frictionlessly on the x axis. Another point mass 2m moves frictionlessly on the y
axis. These two point masses do not collide at the origin. There is no gravity and no friction. The point masses
are connected to each other by a massless spring of constant k and free length zero. Obtain the equations of
motion, both by Lagrange’s method and the Newton-Euler method (FBD’s and momentum balance).
3. Two point masses labelled A and B, of mass m each, move without dissipation on a horizontal plane. They
are connected by a spring of stiffness k and free length L. The system has a constraint: the velocity of point
mass B is always perpendicular to the position vector from A to B. Use x and y coordinates for each mass,
and Lagrange’s method to obtain the equations of motion.
4. Two point masses labelled A and B, of mass m each, are connected by an inextensible, taut string of length
L. Mass A slides on a frictionless horizontal table. The string passes frictionlessly through a small hole in the
table, and mass B hangs below the table. Use polar coordinates r and θ, measured from the hole, for mass
A; and x and y coordinates for mass B. Write the equations of motion. Simulate using Matlab (let r(0) > 0
and θ̇(0) > 0), and check numerically for conservation of three quantities (what are they?).
5. Two point masses labelled A and B, of mass m and 3m respectively, are connected by an inextensible,
taut string of length L. The string passes over a frictionless peg of negligibly small diameter; the masses hang
down, and oscillate in a vertical plane without colliding or otherwise interefering with each other. Write the
equations of motion, and seek initial conditions such that B does not pull A over and around the peg.
66
CHAPTER 5. LAGRANGE’S EQUATIONS
Chapter 6
The Rolling Coin
The “rolling coin” is an idealized mechanical system where a thin, rigid, uniform circular disk with a sharp
edge rolls without slip on a horizontal surface or ground.
Acquiring the ability to write the equations of motion for this system is something of a milestone in a
student’s progress in dynamics: the system has large rotations in 3D as well as nonholonomic constraints, and
therefore requires much of the material covered so far. For this reason, we give this problem a whole, if short,
chapter to itself.
6.1
General comments
We will use (3,1,3) Euler angles (φ, θ, ψ) to represent the rotation of the coin. We will take the reference
configuration as the one where the coin lies flat on the ground (taken to be the x-y plane). This reference
configuration is physically singular in that there are infinitely many contact points, while the slightest change
in configuration gives a single unique point. This reference configuration is also mathematically singular
because of our use of (3,1,3) Euler angles, which are singular whenever the second rotation angle is zero (in
the reference configuration, all rotations are zero). These two singularities (physical and mathematical) will
not concern us as long as the coin rolls.
We denote the mass of the coin by m, its radius by R (not to be confused with the rotation matrix Rf ),
and know or find by integration that
⎤
⎡
mR2
0
0 ⎥
⎢ 4
⎥
⎢
⎥
⎢
mR2
Icm,ref = ⎢ 0
⎥.
0
⎥
⎢
4
2 ⎦
⎣
mR
0
0
2
We begin by modeling this system with 6 degrees of freedom, along with the constraint that it must
touch the horizontal surface (a holonomic constraint) and the constraint of no slip (two scalar nonholonomic
constraint equations). At a later stage, we will point out how we could eliminate one degree of freedom using
the holonomic constraint, leaving a 5 degree of freedom system with two nonholonomic constraint equations.
The 6 generalized coordinates chosen for now are the x, y, and z coordinates of the center of the disk,
along with the three Euler angles (φ, θ, ψ) used in a (3,1,3) sequence to describe the rotation.
We find the rotation matrix Rf of the coin, in terms of the Euler angles, using Eq. 4.3, as
⎡
⎤
cos ψ cos φ − sin ψ cos θ sin φ − cos ψ sin φ cos θ − sin ψ cos φ
sin θ sin φ
⎢
⎥
⎥
Rf = ⎢
⎣ sin ψ cos θ cos φ + cos ψ sin φ cos ψ cos φ cos θ − sin ψ sin φ − sin θ cos φ ⎦ .
sin ψ sin θ
cos ψ sin θ
cos θ
By Eq. 4.9,
Icm = Rf Icm,ref RfT
67
CHAPTER 6. THE ROLLING COIN
68
⎡
is found to be
Icm =
1 + s2θ s2φ
mR2 ⎢
⎢ cφ sφ c2 − 4
θ
4 ⎣
sθ sφ cθ
cφ sφ c2θ − 4
1 + s2θ c2φ
−sθ cφ cθ
sθ sφ cθ
⎤
⎥
−sθ cφ cθ ⎥
⎦,
2
1 + cθ
where cθ , sθ , cφ and sφ represent cos θ, sin θ, cos φ and sin φ respectively.
The angular velocity ω (or, in matrix notation, ω) of the coin, in terms of the Euler angles and their time
derivatives, is given by Eq. 4.6 as
⎧
⎫
cos φ θ̇ + sin θ sin φ ψ̇ ⎪
⎪
⎪
⎪
⎨
⎬
.
ω=
sin φ θ̇ − sin θ cos φ ψ̇
⎪
⎪
⎪
⎪
⎩
⎭
φ̇ + cos θ ψ̇
We can now proceed in either of two ways.
6.2
Lagrange’s equations
We can write the kinetic energy of the system as
T =
1
m 2
ẋ + ẏ 2 + ż 2 + ω T Icm ω.
2
2
The potential energy is
V = mgz,
and the Lagrangian
L = T − V.
Since we are using 6 degrees of freedom, we must incorporate the holonomic constraint of contact with the
surface in its differentiated form, i.e., we must say that the normal component of the contact point velocity is
zero. This, combined with the no slip condition, simply means that all three components of the contact point
velocity are zero.
^k
G
^
n
ground
C
Figure 6.1: Side view of coin.
At any arbitrary configuration of the coin, its plane intersects the ground along some line CC; looking at
the coin along this line, both the coin and the ground reduce to straight lines. This view is shown in figure 6.1.
The center of the coin is called G, and the instantaneous point of contact is called C. A unit vector normal
to the coin is shown, and denoted by n̂. In the reference configuration, n̂ is aligned with k̂ (which is the same
as e3 in Chapter 4).
Unless the coin is flat (i.e., in the singular configuration), the vector
p = n̂ × k̂
is nonzero, and (in figure 6.1) points out of the page. Notice that we could view the coin along line CC in two
ways; the convention we have adopted is that n̂ emerges from the coin towards the right side.
6.3. NEWTON-EULER EQUATIONS
69
Now we see that
rG/C =
R
p × n̂.
|p|
The velocity of G can be written in two ways:
v G = ẋî + ẏ ĵ + ż k̂ = v C + ω × rG/C .
The nonholonomic constraint equation (in vector form) then turns out to be
ẋî + ẏ ĵ + ż k̂ − ω × rG/C = 0.
For calculations, the above vector equations need to be expressed using matrices. Accordingly,
p = S(Rf e3 )e3 ,
R
rG/C = /
S(p)Rf e3 ,
pT p
and (the nonholonomic constraint equation)
⎧
⎫
⎨ ẋ ⎬
R
ẏ
S(ω)S(p)Rf e3 = 0.
−/
⎩
⎭
pT p
ż
Carrying out the above calculations, we find
rG/C
⎧
⎫
⎨ −R sin φ cos θ ⎬
R cos φ cos θ
=
.
⎩
⎭
R sin θ
The third component above gives
z = R sin θ,
which could have been directly used as a holonomic constraint to eliminate z and reduce the number of degrees
of freedom to 5.
The nonholonomic constraint equations turn out to be
ẋ + R cos φ cos θ φ̇ − R sin θ sin φ θ̇ + R cos φ ψ̇
ẏ + R sin φ cos θ φ̇ + R sin θ cos φ θ̇ + R sin φ ψ̇
ż − R cos θ θ̇
= 0,
= 0,
(6.1)
(6.2)
=
(6.3)
0.
Again, the third constraint equation is obviously a differentiated holonomic constraint.
It only remains to take some derivatives and write down the equations of motion. This is left as an exercise
for the reader.
6.3
Newton-Euler equations
We can also obtain the equations of motion for the rolling coin using the Newton-Euler approach, i.e., by using
the equations of linear and angular momentum balance. To this end, we draw a free body diagram of the coin
(figure 6.2 (right)). There are two forces acting on the coin: its weight acts through its center of mass G, and
an unknown vector force F acts at the point of contact C.
We will use the same 6 degrees of freedom as we used before. Now the acceleration of the center of mass
G, written in matrix notation as
⎧
⎫
⎨ ẍ ⎬
ÿ
aG =
,
⎩
⎭
z̈
CHAPTER 6. THE ROLLING COIN
70
coin
G
ground
C
mg
F
Figure 6.2: Rolling coin and its free body diagram.
is sought in addition to the second derivatives of the Euler angles. The angular acceleration of the coin, α, is
related to the second derivatives of the Euler angles through Eq. 4.10. Linear momentum balance for the coin
gives (in matrix notation)
(6.4)
F − mge3 = maG .
Angular momentum balance about the point C gives
rG/C × (−mg k̂) = rG/C × maG + I cm · α + ω × I cm · ω,
which in matrix notation is
−mgS(rG/C )e3 = mS(rG/C )aG + Icm α + S(ω)Icm ω.
(6.5)
Equations 6.4 and 6.5 represent 6 scalar equations with 9 unknowns (three components each of F , aG and
α). Differentiating Eqs. 6.1 through 6.3 above will give 3 more scalar equations, but introduce 3 new unknowns
(the second derivatives of the Euler angles). Finally, Eq. 4.10 will provide 3 scalar equations that relate these
second derivatives to α. We will then have 12 scalar simultaneous equations in 12 unknowns. Solving them
will yield, among other information, the second derivatives of the generalized coordinates (x, y, z, φ, θ and ψ)
in terms of system parameters and the instantaneous values of the coordinates and their first time derivatives.
Note that if we do not want to find the contact force F , then we can drop Eq. 6.4 and solve just 9 equations
for 9 unknowns. Solution of these equations is left as an exercise for the reader.
Chapter 7
Linear Vibrations
The subject of mechanical vibrations has many practical applications. Much can be achieved with an understanding of linear vibrations, by which we mean that the governing differential equations are linear (the
motions need not be linear, e.g., in torsional oscillations of shafts). This chapter provides a short introduction
to some of the ideas in this subject. If you need to work with vibrations, do read an entire book on the topic.
7.1
A single degree of freedom system
See figure 7.1(a). A point mass m is supported by a spring of constant k. Its displacement x is measured from its
equilibrium position. Gravity is neglected. Let the system be disturbed from its equilibrium position. Figure
7.1(b) shows the force acting on the mass at some later instant (the spring force opposes the displacement of
kx
k
x
m
(b)
(a)
Figure 7.1: (a) A spring mass system. (b) Free body diagram.
the mass). By Eq. 1.1 we have
mẍ = −kx,
or
mẍ + kx = 0.
The solution to the above is
x = A sin(ωn t) + B cos(ωn t),
0
where
ωn =
k
m
and A and B are constants to be determined from the initial conditions of the system. A key point is that for
this undamped and unforced system, all solutions are periodic. Moreover, all these periodic solutions have a
special frequency called the natural frequency of the system, which is determined by system parameters alone
(k and m).
71
CHAPTER 7. LINEAR VIBRATIONS
72
Addition of a damping element (a dashpot) and some external forcing f (t) would change the equation of
motion to
mẍ + cẋ + kx = f (t).
In the above, if c is relatively small (in some suitably nondimensionalized sense), and f (t) is the sum of
possibly many different periodic functions with different frequency components, then the solution will usually
be dominated by responses that may be largely understood using three simple ideas:
1. The steady state response will be seen after a period of initial transient motion. During this transient
motion, superimposed on the steady state response, there will be a slowly decaying oscillation at a
frequency ≈ ωn .
2. Away from resonance, i.e., for parts of the forcing that involve frequencies not close to ωn , the damping
plays an insignificant role in the steady state response.
3. If the forcing involves a significant component at a frequency close to ωn , then the damping term is
important but the other frequency components in the forcing play an insignificant role in the steady
state response.
To see this using examples, consider the special case
ẍ + cẋ + x = F0 sin t + F1 cos 2t.
The general solution (from the symbolic algebra software MAPLE) is
0
+
0
ct *
−
F1 (2c sin 2t − 3 cos 2t) F0 cos t
c2
c2
t + B cos
t
+
.
x = e 2 A sin
1−
1−
−
4
4
9 + 4c2
c
The three qualitative ideas listed above may be seen in the solution for x. The first term shows the damped
transient, the second shows the response to nonresonant forcing, and the third to resonant forcing.
As seen above, the natural frequency of a lightly damped single degree of freedom system directly provides
a lot of insight into the system’s vibrations. We now consider systems with many degrees of freedom, and seek
their natural frequencies.
7.2
Normal modes
Let us begin with figure 7.2. It shows a two mass, three spring system (two degrees of freedom), along with
individual free body diagrams for each mass. The displacement of the first mass is x1 , and that of the second
mass is x2 .
We first seek special motions involving normal modes. In a normal mode, all material points on the structure
execute synchronous harmonic motion at a special (or natural) frequency. The amplitudes of vibration of all
these points maintain specific proportions, giving rise to mode shapes associated with the corresponding natural
frequencies. General motions, it will be seen, can be usefully viewed as superpositions of motions in the normal
modes.
The equations of unforced and undamped motion are:
mẍ1
= −kx1 − k(x1 − x2 ),
(7.1)
mẍ2
= −k(x2 − x1 ) − kx2 .
(7.2)
These can be arranged in the form
M ẍ + Kx = 0,
where the 2 × 2 matrices M and K are symmetric and positive definite (in other problems, K can be positive
semidefinite), and x is the 2 × 1 column matrix {x1 , x2 }T . How the above matrix equation will lead to two
normal modes will be discussed later in mathematical terms. In this section, we provide a more intuitive
discussion.
7.3. THE GENERALIZED EIGENVALUE PROBLEM FOR VIBRATIONS
x1
k
kx1
x2
k
m
k
m
m
m
k(x1- x 2 )
73
kx 2
k(x 2- x1)
Figure 7.2: A two degree of freedom system.
First mode
Second mode
Figure 7.3: Normal modes.
Imagine a motion in which the two masses vibrate in phase, and with the same frequency and amplitude,
such that x1 = x2 in Eqs. 7.1 and 7.2. The spring in the middle plays no role. Such a motion, which
clearly satisfies the equations of motion, is possible here because of the symmetry in the system. This motion
corresponds/
to the first normal mode. Note that the ratio of x1 to x2 is fixed (at 1), and the angular frequency
is fixed (at k/m). The amplitude of the motion is not fixed, however: it could be anything. This is general;
for any normal mode, the ratios of the various displacements to each other, as well as the frequency of the
oscillations, are determined while the amplitude can be anything.
For this system, the second normal mode is also easy to see. We set x1 = −x2 in the equations of motion.
Now the motions of the two masses are mirror images of each other, again made possible by the symmetry in
the /
system. For this normal mode, the ratio of x1 to x2 is fixed (at −1), and the angular frequency is fixed
(at 3k/m).
The two normal modes of the system are sketched in figure 7.3. That there are no other normal modes
can be proved using mathematics, which comes in the next section. For systems without such symmetry, the
basic observations continue to hold. In particular, there is one normal mode per degree of freedom.
7.3
The generalized eigenvalue problem for vibrations
Along the lines of Eqs. 7.1 and 7.2, equations for vibrations of multi degree of freedom systems can be written
as
M ẍ + Kx = 0.
(7.3)
Here the n × n matrix M is called the mass matrix, and is usually1 symmetric and positive definite. The
1 We say “usually” because, if A is any n × n invertible matrix, then clearly AM ẍ + AKx = 0 is a valid set of equations of
motion. However, the new matrix AM need not be symmetric and positive definite. Of course, AM is merely M in disguise.
CHAPTER 7. LINEAR VIBRATIONS
74
n × n matrix K is called the stiffness matrix, and is usually symmetric and positive semidefinite.
Seeking normal modes, we put x = u cos ωt (with u = 0) in Eq. 7.3 to get
Ku = ω 2 M u.
(7.4)
Equation 7.4 is called a generalized eigenvalue problem. The eigenvalues obtained will be the squared
natural frequencies of the system.
For large systems (i.e., large n), it is numerically convenient to tackle the generalized eigenvalue problem
directly. But we are not concerned with that here. Noting that M is positive definite (hence invertible), we
write
−1 M K u = ω 2 u,
which is recognized as a standard eigenvalue problem (of the form Au = λu) with n eigenvalues (possibly
repeated). Provided K is symmetric and positive semidefinite, it can be proved that these eigenvalues are real
and nonnegative. Each such eigenvalue represents a natural frequency and has a corresponding normal mode.
Let ω1 , ω2 , · · · , ωn be the natural frequencies and let u1 , u2 , · · · , un be the corresponding eigenvectors.
Consider
Ku1 = ω12 M u1 ,
Ku2 = ω22 M u2 .
Premultiplying the first equation by uT2 and the second by uT1 we get
uT2 Ku1 = ω12 uT2 M u1
and
uT1 Ku2 = ω22 uT1 M u2 .
T
Noting that the scalar quantity uT1 Ku2 = uT1 Ku2 = uT2 Ku1 because K = K T , and also noting a similar
fact involving M , we subtract the second equation from the first to get
2
ω1 − ω22 uT1 M uT2 = 0.
It follows that if ωi = ωj , then uTi M uj = 0. It can be shown that if ωi = ωj (though i = j), then
the eigenvectors can be chosen such that uTi M uj = 0. Moreover, the positive definiteness of M guarantees
uTi M ui > 0, so each ui can be scaled such that uTi M ui = 1 for i = 1, 2, · · · , n. Thus, the eigenvectors are (or,
in the presence of repeated eigenvalues, can be chosen such that they are) orthonormal in a mass-weighted
sense.
How about a stiffness-weighted sense? There are orthogonality properties involving the stiffness matrix as
well. These follow directly from
Kui = ωi2 M ui ,
giving
and
uTi Kui = ωi2 ,
for
uTj Kui = 0,
i = 1, 2, · · · , n,
for
i = j.
We can define the n × n matrix of eigenvectors, also called the mass weighted modal matrix, as
Φ = [u1 , u2 , · · · , un ] ,
where each ui is an n × 1 column matrix. We also define the diagonal matrix
⎡ 2
⎤
ω1 0 · · ·
0
⎢ 0 ω22 · · ·
0 ⎥
⎢
⎥
Λ=⎢ .
.
.. ⎥ .
.
..
..
⎣ ..
. ⎦
0
0 · · · ωn2
Then the solution to the generalized eigenvalue problem of Eq. 7.4, which lies at the heart of linear vibration
theory, is written as
KΦ = M ΦΛ.
(7.5)
The reader should insert the elements of these matrices and verify the above.
There are systematic ways of deriving the equations of motion so that the M obtained will be symmetric and positive definite.
7.4. FORCED VIBRATIONS WITH DAMPING
7.4
75
Forced vibrations with damping
Consider a linear, n degree of freedom, vibrating system which has some damping and is acted upon by m
forces f1 (t), f2 (t), · · · , fm (t) at some specified locations. The equations of motion will be of the form
M ẍ + C ẋ + Kx = BF (t),
(7.6)
where M , K and x are as defined before; C is an n×n matrix that is often symmetric and positive semidefinite
(exceptions are gyroscopic systems, not discussed here: refer to a vibrations textbook); B is an n × m matrix;
and the m × 1 matrix F (t) = {f1 (t), f2 (t), · · · , fm (t)}T .
Using the modal decomposition of Eq. 7.5, we introduce new variables y defined by
x = Φy.
Inserting the above into Eq. 7.6, we get
M Φÿ + CΦẏ + KΦy = BF (t).
Premultiplying by the invertible matrix ΦT , we incur no loss of information in writing
ΦT M Φÿ + ΦT CΦẏ + ΦT KΦy = ΦT BF (t).
(7.7)
Under some special conditions, which include the special case where C = αM + βK (called proportional
damping), ΦT CΦ is diagonal2 . Assuming that ΦT CΦ = Cmodal is diagonal, with the ith diagonal element
being ci , we write from Eq. 7.7,
ÿ + Cmodal ẏ + Λy = G(t),
where G(t) is a column matrix whose elements are known linear combinations of the given forces f1 (t), f2 (t),
· · · , fm (t).
The equations are now decoupled. The ith equation is
ÿi + ci ẏi + ωi2 yi = gi (t).
Thus, the study of vibrations in n degree of freedom systems, through solving an eigenvalue problem, and
under some assumptions about the nature of damping present, can be reduced to the study of n decoupled
single degree of freedom systems.
7.5
Approximations via Lagrange’s equations
A great advantage of Lagrange’s equations (over Newton-Euler equations) is that they allow systematic approximations.
We illustrate this with an example. Consider an Euler-Bernoulli beam (see figure 7.4) of flexural rigidity
EI, length L and mass per unit length m. There is a spring (spring constant K) attached to it at the free
end. This system has infinitely many natural frequencies. What is the first one?
We will make use of an approximation method called the Rayleigh-Ritz method (sometimes called the
method of assumed modes).
Examine the boundary conditions of the beam. At the left end of the beam, both displacement and slope
are zero at all times. Let the displacment of the beam be given by w(x, t). Then w(0, t) ≡ 0 and wx (0, t) ≡ 0.
At the right end, there are no such restrictions.
The potential energy of the system, neglecting gravity, is (see a book on the strength of materials)
L
EI
kw(L, t)2
wx (x, t)2 dx +
,
V =
2
2
0
−1 i
generally, if C = M n−1
K then ΦT CΦ is diagonal, as may be seen by observing from Eq. 7.5 that
i=0 αi M
2
M −1 KΦ = ΦΛ, that M −1 K Φ = M −1 KM −1 KΦ = M −1 KΦΛ = ΦΛ2 , and so on. Since there are now n free parameters, the
strength of damping in each mode is independently adjustable. However, note that many systems of practical importance do not
fall in this category; for them, ΦT CΦ is not diagonal.
2 More
CHAPTER 7. LINEAR VIBRATIONS
76
E, I, L,m
K
x
Figure 7.4: An Euler-Bernoulli cantilever beam with an end support.
where the first term represents the strain energy in the beam and the second represents that in the spring.
The kinetic energy of the system is (neglecting rotary inertia of beam elements)
L
m
wt (x, t)2 dx.
T =
2L
0
The Lagrangian
L = T − V.
So far, there is no approximation.
Let us now make the approximation
w(x, t) ≈ a1 (t)
x2
x3
+
a
(t)
.
2
L2
L3
In the above approximation, the ai (t) play the role of generalized coordinates, and the functions of x (i.e.,
nondimensionalized versions of x2 and x3 ) are called assumed modes or shape functions. The above approximation is somewhat arbitrary, but not completely so. The shape functions individually satisfy the displacement
and slope boundary conditions at the left end; and they also form the beginning of a power series which (theoretically at least) can describe any reasonable function to high accuracy. There are deeper technical issues
in these approximations which we do not go into.
By the above approximation, we get (from MAPLE)
2
k 2
ȧ2
ȧ1 ȧ2
ȧ2
EI +
+ 1 − 3 6 a2 2 + 6 a1 a2 + 2 a1 2 −
a1 + 2a1 a2 + a2 2 .
L=m
14
6
10
L
2
Writing Lagrange’s equations and collecting terms, we find that for this approximated system
⎡ 1 1 ⎤
⎢ 5 6 ⎥
M = m⎣
⎦,
1 1
6 7
⎤
EI
EI
+
k
6
+
k
4
⎥
⎢ L3
L3
⎥.
K=⎢
⎦
⎣ EI
EI
6 3 + k 12 3 + k
L
L
The eigenvalues are (from MAPLE)
1
2
612 EI + 6 kL3 ± 6 9984 (EI) + 64 EI kL3 + k 2 L6
and
⎡
L3 m
,
whence an estimate for the smallest or first natural frequency is
1
⎞1/2
⎛
2
612 EI + 6 kL3 − 6 9984 (EI) + 64 EI kL3 + k 2 L6
⎠ .
ω1 ≈ ⎝
L3 m
7.5. APPROXIMATIONS VIA LAGRANGE’S EQUATIONS
77
Is the approximation any good? A preliminary check is to put k = 0 to obtain
0
EI
3.533
L3 m
and compare with the known lowest natural frequency of a cantilever beam. From a structural dynamics
textbook3 , this frequency is
0
EI
3.516
,
L3 m
i.e., the error is on the order of half a percent.
The error in estimating the second natural frequency is expected to be higher. From a textbook, the second
frequency is
0
EI
22.0
L3 m
while the above two degree of freedom approximation predicts
0
EI
,
34.8
L3 m
with an error of about 58 percent. We do not view this as a serious failure because we set out to find only the
first natural frequency.
But suppose we want to estimate the first two frequencies? Then the above approximation is clearly not
acceptable. We need to use more shape functions. And we should avoid the powers of x and use somewhat
better behaved shape functions (higher powers of x, such as x99 and x100 , say, are not that different from one
another; approximations based on power series may therefore behave poorly in numerical work).
A reliable way to choose shape functions is as follows. We note that w(x, t) has both zero value as well as
zero slope at x = 0. If we let w(x, t) = (x/L) h(x, t), then h(x, t) is zero at x = 0 but its slope can be nonzero.
Now let
N
mπx x
,
am (t) sin
h(x, t) = a0 (t) +
L m=1
L
where N is an order of approximation that we choose. In the above, we note that the first term allows the
displacement at the right end to be anything (i.e., a0 (t)); all other terms constitute a Fourier sine series, which
can represent any reasonable functions on an interval, with zeroes at the endpoints. Thus, we have posed the
following approximation:
N
mπx x2
x
.
(7.8)
w(x, t) = a0 (t) 2 +
am (t) sin
L
L
L
m=1
Using Eq. 7.8 for a sequence of increasing N will show the second natural frequency converging to its true
value. The eigenvalue problem will have to be solved numerically for fixed parameter values. This calculation
is not given here (the reader may benefit from doing it).
3 R.
W. Clough and J. Penzien, 1982. Dynamics of Structures, McGraw-Hill International Edition, pp. 313.
78
CHAPTER 7. LINEAR VIBRATIONS
Chapter 8
Some Problems Involving Single
Bodies
We have already looked at equations of motion for a coin rolling without slip on a horizontal plane. We now
look at some other problems involving single rigid bodies. These will help develop expertise in using the tools
developed so far.
8.1
Cylinder rolling down an incline
A uniform cylinder of mass m and radius R rolls without slip down an inclined plane that makes an angle θ
with the horizontal. Let the distance traveled by the center of mass down the slope be s, and the angle it rolls
through be φ. The no-slip condition here implies
s = s0 + Rφ,
(8.1)
where s0 gives the initial position of the cylinder (when φ = 0). The above is a holonomic constraint: the
system has only one degree of freedom.
Since Eq. 8.1 is true for all s and φ, it can be differentiated with respect to time to give
ṡ = Rφ̇.
The kinetic energy of the cylinder is
T =
1 2 1 2
1
1 mR2 ṡ2
3
mṡ + J φ̇ = mṡ2 +
= mṡ2 .
2
2
2
2 2 R2
4
The potential energy, taking the position s = 0 as the datum, is
V = −mgs sin θ.
The Lagrangian is
L=T −V =
3 2
mṡ + mgs sin θ.
4
The equation of motion is
3
ms̈ = mg sin θ,
2
or
2
g sin θ.
(8.2)
3
Note that a body sliding down without friction would experience an acceleration of g sin θ. The lower
acceleration here is due to the no-slip condition, which requires gravity to do work in both accelerating the
center of mass as well as increasing the rotation rate of the cylinder.
s̈ =
79
80
CHAPTER 8. SOME PROBLEMS INVOLVING SINGLE BODIES
Another point to note here is that the assumption that it is a uniform cylinder determines J in terms of m
and R, and that in the end the acceleration is independent of both m and R. This conclusion could have been
arrived at from dimensional analysis as well (a powerful and frequently useful tool, but one which is given
little place here). To see how, we first write
s̈ = f0 (m, R, g, θ),
where f0 is an as yet undetermined function of the physical parameters that determine the solution. We can
rewrite this as
s̈ = g f1 (m, R, g, θ),
where f1 is some as yet unknown dimensionless function. Looking for nondimensional combinations of the
dimensional parameters, we write
ma Rb g c = dimensionless,
whence, using M , L and T to denote mass, length and time respectively, we have
M a Lb+c T −2c = M 0 L0 T 0 ,
or
a = b = c = 0.
Thus,
s̈ = gf2 (θ)
for some undetermined function f2 . The sinusoidal dependence on θ and the factor of 2/3 cannot be deduced
from such analysis. The 2/3 would change, for example, if it was a uniform sphere of radius R rolling down
the incline; but the dimensional analysis above would not.
8.2
Hemispherical shell on a table
Consider free oscillations of a hemispherical thin shell of radius R on a table (figure 8.1). We will find the
equations governing small planar oscillations in the absence of (i) slip and (ii) friction.
y
2R
φ
x
Table
Figure 8.1: Hemispherical thin shell on a table.
We begin by locating the center of mass of the shell. Ignoring the thickness of the shell, we find by routine
integration that the height of the center of mass when the shell is in the equilibrium position (see figure, left)
is R/2. Choosing the positive x-axis to the right, the positive y-axis vertically upwards, and the positice z-axis
pointing out of the page (see figure, right), we observe that the moment of inertia matrix in the equilibrium
position is diagonal due to symmetry; and that Izz = J = 5mR2 /12.
Comment: Though the object is not planar, it can have planar dynamics in the x-y plane because rotations
about the z-axis affect neither the third row nor the third column of Icm .
8.3. TORQUE-FREE RIGID BODY
81
We first consider motion without slip. As shown in the figure (right), we use the rotation angle φ as a
generalized coordinate. When there is no slip, the position vector of the center of mass is
R
R
rcm = Rφ − sin φ î + R − cos φ ĵ.
2
2
The velocity of the center of mass is obtained by differentiation; the kinetic energy follows; and the potential
energy is gravitational and can be calculated using the y-component of rcm . This gives the Lagrangian. The
equation of motion then is
1
1
5
mR2 φ̈ − mR2 φ̈ cos(φ) + mR2 φ̇2 sin(φ) + mgR sin(φ) = 0,
3
2
2
which on linearization for small φ becomes
φ̈ +
3g
φ = 0.
4R
We now consider motion without friction. Now the x-coordinate of the center of mass is no longer dependent
on φ, and there are two degrees of freedom. We use x and φ as generalized coordinates. Now
R
rcm = x î + R − cos φ ĵ.
2
The equations of motion are
ẍ = 0
and
13
1
1
1
mR2 φ̈ − mR2 φ̈ cos(2 φ) + mR2 φ̇2 sin(2 φ) + mgR sin(φ) = 0,
24
8
8
2
which on linearization for small φ gives
φ̈ +
6g
φ = 0.
5R
The frequency of small oscillations is lower for the no-slip case. While the potential energy gained for a
given rotation φ is the same in both cases (frictionless and no-slip), the no-slip constraint forces more mass
to be in motion for the same rotation (to use a spring mass analogy, we have the same stiffness but greater
mass).
8.3
Torque-free rigid body
A rigid body moves with no net external torques about its center of mass. Taking angular momentum balance
about its center of mass, we find
[Icm ] · α + ω × [Icm ] · ω = 0
in any inertial frame of reference XY Z, where
α=
dω xyz
dt
.
XY Z
It is convenient to view this classical problem in a non-inertial frame xyz attached to the rigid body. As
discussed earlier,
dω xyz
dω xyz
:0
=
+ ω xyz
×
ω
xyz ,
dt
dt
XY Z
xyz
so we have
[Icm ] ·
dω xyz
dt
xyz
= −ω × [Icm ] · ω.
CHAPTER 8. SOME PROBLEMS INVOLVING SINGLE BODIES
82
In the frame attached to the rigid body, [Icm ] is a constant (this fact was, in fact, used in deriving the
above equation earlier). Moreover, there is a coordinate system (also called xyz) in which [Icm ] is diagonal,
i.e.,
⎤
⎡
0
Ixx 0
[Icm ] ≡ [Icm ]ref ≡ ⎣ 0 Iyy 0 ⎦ .
0
0 Izz
These special coordinate axes are called the principal axes of the body.
Writing our equations in that coordinate system, we have
Ixx ω̇x
= − (Izz − Iyy ) ωy ωz ,
(8.3)
Iyy ω̇y
Izz ω̇z
= − (Ixx − Izz ) ωx ωz ,
= − (Iyy − Ixx ) ωx ωy .
(8.4)
(8.5)
These are called Euler’s equations for a torque-free rigid body. A torque about the center of mass, if present,
could be included on the right hand sides of these equations.
Equations 8.3 through 8.4 show that a torque-free uniform sphere (or, say, cube) has no angular acceleration.
They also show that, for any rigid body, pure spin at any rate about any principal axis is possible (i.e., the
equations are satisfied if any two ω-components are zero and the third is an arbitrary constant).
There are two other aspects of this problem that we consider below.
8.3.1
Stability of pure spin
Let us consider the stability of pure spin motions about the x, y and z axes.
Considering small deviations from pure spin about the x axis, let
ωx = ωx0 + ξx ,
ωy = 0 + ξy ,
ωz = 0 + ξz .
Substituting into Eqs. 8.3 through 8.5 and linearizing, we obtain
Ixx ξ˙x
Iyy ξ˙y
Izz ξ˙z
= 0,
= − (Ixx − Izz ) ωx0 ξz ,
(8.6)
(8.7)
= − (Iyy − Ixx ) ωx0 ξy .
(8.8)
The above can be easily solved, and say that ξx = constant, while ξy and ξz oscillate sinusoidally at the same
frequency. Thus, pure spin about the principal axis of smallest moment of inertia is stable within the linearized
approximation. Nonlinear analysis of the above equations can show that the stability result continues to hold.
Similar stability analyses of pure spin about the y and z axes shows the former to be exponentially unstable
and the latter to be stable in much the same way as about the x axis. Thus, pure spin is stable about the axes
of smallest and largest moment of inertia, but unstable about the axis of intermediate moment of inertia.
A final point. Over long periods of time, minute amounts of unmodeled energy dissipation add up for real
bodies. For such long times and for such real bodies, only pure spin about the axis of largest moment of inertia
is stable. I put it in italics because man has in fact designed satellites to be spin stabilized about the axis of
smallest moment of inertia; those satellites are no more.
There is a nice geometric way to see why pure spin about the axis of smallest moment of inertia is unstable
in the long run. The torque-free rigid body conserves angular momentum H in the XY Z frame. So, in a
coordinate system attached to the xyz frame (or any other frame), the magnitude of H stays constant, i.e.,
2
2
2
ωx2 + Iyy
ωy2 + Izz
ωz2 = constant.
H 2 = Ixx
Over short times, the (rotational) kinetic energy also stays constant. We write
2 T = Ixx ωx2 + Iyy ωy2 + Izz ωz2 = constant.
Each of the above conservation conditions constrains the tip of the vector H to lie on an ellipsoid. That is, since
xyz is non-inertial, H is not conserved in it; the moving tip of the vector H, however, must simultaneously
8.3. TORQUE-FREE RIGID BODY
83
z
T, initial
T, final
H
x
Figure 8.2: Energy and momentum ellipsoids for a torque-free rigid body.
satisfy the above two constraints. It must therefore move along the curves of intersection of these two ellipsoids
(these ellipsoids must intersect or at least touch for a real solution for H to exist).
For identification, let us call the first (H 2 conserved) ellipsoid the H-ellipsoid; and the second (2T conserved) the T -ellipsoid. The situation is depicted in a 2D schematic in figure 8.2. Consider the shapes of these
ellipsoids. The principal axes of these ellipsoids are aligned with the coordinate axes.
Assuming that
Ixx < Iyy < Izz ,
each ellipsoid is longest in the x-direction and shortest in the z-direction. However, since the moment of inertia
ratios are exaggerated by squaring in the H-ellipsoid, its shape is more elongated than that of the T -ellipsoid.
Consider an initial condition where a body spins about the x-axis. This corresponds to the ellipsoid labeled
“T , initial” in the figure touching the H-ellipsoid at a point on the x-axis. Now imagine that over a long period
of time, kinetic energy is slowly lost due to unmodeled effects like minute dissipative vibrations in the body.
Angular momentum continues to be conserved because the system has no external moments acting on it, and
so the H-ellipsoid remains fixed. But the T -ellipsoid begins to shrink (while maintaining a fixed aspect ratio).
Contact occurs no longer on the x axis, but there is now intersection between the two ellipsoids on a small
but growing closed curve encircling the x-axis. Eventually, the T -ellipsoid shrinks enough so that contact with
the H-ellipsoid again occurs at a point, but this time on the z-axis. This configuration is shown labeled “T ,
initial” in the figure. Further energy dissipation is not possible under the assumptions of the analysis, and
subsequent spin continues indefinitely about the z-axis. For further discussion of torque-free motion in the
context of these ellipsoids, the reader is encouraged to consult a classical source (e.g., H. Goldstein, Classical
Mechanics, second ed., Addison-Wesley, 1980).
8.3.2
Axisymmetric bodies
Consider the case
Iyy = Izz = Is = Ixx ,
where the s-subscript denotes symmetry, for the case of small deviations from pure spin about the x-axis.
Letting
Izz − Ixx = Iyy − Ixx = k = 0,
Eqs. 8.6 through 8.8 become
Ixx ξ˙x
Is ξ˙y
Is ξ˙z
=
0,
(8.9)
= k ωx0 ξz ,
(8.10)
= −k ωx0 ξy .
(8.11)
As mentioned above, solutions for ξy and ξz are sinusoidal.
It is helpful to think about the above equations in terms of a uniform disk (or, more roughly, a circular
dinner plate). Imagine such a plate held vertically in the Y -Z plane, and then thrown up into the air with a
spin approximately along the X-axis. Due to imperfect initial conditions, the disk will wobble a little as it
spins. How are the wobble and spin rates related?
CHAPTER 8. SOME PROBLEMS INVOLVING SINGLE BODIES
84
For a disk in the configuration described, we take Ixx = mR2 /2 and Is = mR2 /4. Equations 8.10 and 8.11
become
ξ˙y
ξ˙z
= −ωx0 ξz ,
(8.12)
= ωx0 ξy .
(8.13)
The solution is
ξy = A cos(ωx0 t + φ),
ξz = A sin(ωx0 t + φ),
with A and φ arbitrary constants; we let ξx = 0, because it makes no difference below.
Taking unit vectors î, ĵ and k̂ along the x, y and z axes, respectively, and defining ξ = ξy ĵ + ξz k̂, we note
that ξ lies in the plane of the disk and rotates about the x-axis at an angular rate ωx0 , as sketched in figure
8.3. As also shown using grey shading in the figure, the vectors H, ω, ξ as well as the x-axis all lie in the same
x
_
H
_
ω
z
ωx 0
ξ
_
y
Figure 8.3: A spinning, wobbling disk.
plane. This plane may be thought of as a single rigid body, which we call (say) P . We see that
ω disk = ωx0 î + ξ,
and
ω P/disk = ωx0 î,
so
ω P = 2ωx0 î + ξ.
Since
H = Ixx ωx0 î + Is ξ,
we note that ω P is aligned with H (because Ixx = 2Is ).
All these conclusions are drawn with reference to a moving frame xyz. However, we know in advance that
H is fixed in the inertial frame XY Z because the disk is torque-free. We also know ω P , the absolute angular
velocity of P , and know that it is aligned with H.
It is, finally, clear that the x-axis, being a part of P , describes a cone by going around H with an angular
velocity equal to ω P . Comparing with the angular velocity of the disk itself, and recalling that ξ is small, we
conclude that the wobble rate is twice the spin rate.
8.4
Wobbling disk using Euler angles
The foregoing discussion of the spinning and wobbling disk, based on Euler’s equations, was not developed
along the lines adopted so far in these notes. Results obtained were in terms of the angular velocity, whose
8.4. WOBBLING DISK USING EULER ANGLES
85
relation to configuration is not always trivial (see chapter 4). Interpretation of results was dependent on 3D
visualization, which can be unreliable.
It is conceptually simpler to attack the problem directly using Euler angles. For the disk approximately
in the y-z plane with spin approximately along the x-axis, we use (1,3,1) Euler angles denoted (φ, θ, ψ). That
is, the first rotation is about the x-axis and is called φ; the second rotation is about the rotated body-fixed
z-axis and is called θ; and the third rotation is about the twice-rotated body-fixed x-axis and is called ψ.
We can assume the center of mass of the disk is held stationary by a ball and socket joint. We take
⎤
⎡
2Is 0 0
Icm,ref = ⎣ 0 Is 0 ⎦ .
0
0 Is
There is no translational kinetic energy, and no gravitational potential energy. The Lagrangian equals the
rotational kinetic energy. Using Maple, the equations of motion are:
3 + cos(2θ)
φ̈ + 2 cos(θ)ψ̈ − sin(2θ)φ̇ θ̇ − 2 sin(θ)θ̇ ψ̇
2
sin(2θ) 2
φ̇ + 2 sin(θ)φ̇ ψ̇
θ̈ +
2
ψ̈ + cos(θ)φ̈ − sin(θ)φ̇ θ̇
=
0,
(8.14)
=
0,
(8.15)
=
0.
(8.16)
In the above, we note that any configuration is an equilibrium configuration (i.e., the three angles equal to
three arbitrary constants give a valid solution). This is expected.
Next, we consider a steady motion given by three nonzero constants
θ = θ0 ,
φ̇ = φ̇0 ,
ψ̇ = ψ̇0
(note that θ = 0 gives a singular configuration). The above identically satisfy Eqs. 8.14 and 8.16, while Eq.
8.15 gives
2ψ̇0
.
φ̇0 = −
cos(θ0 )
For small θ0 we have
φ̇0 = −2ψ̇0 .
Moreover, also for small θ0 , we find that the absolute angular velocity of the disk is
ω disk = −ψ̇0 î + O(θ0 ),
where î is along the stationary or inertially fixed X-direction, and where O(θ0 ) represents small terms comparable to θ0 . Since the nominal motion of the disk is pure spin, and the wobble is assumed small, we conclude
that
φ̇0 = 2 × spin rate.
Finally, recall from the definition of (1,3,1) Euler angles that θ0 is the angle made by the body-fixed x axis
with the inertially fixed X-axis. This x-axis describes a cone around the X-axis, producing the wobble we
are interested in here. Consider a unit vector n̂ attached to the disk, such that n̂ref = î. Then the matrix of
components of the rotated n̂, given using the net rotation matrix Rnet , is
⎧
⎫
cos(θ)
⎨
⎬
sin(θ) cos(φ)
n = Rnet nref =
.
⎩
⎭
sin(θ) sin(φ)
From the projection of n̂ on the Y -Z plane we see φ̇0 is the wobble rate, which from the above is twice the
spin rate.
Comment: 3D visualization is hard for many, and solutions that do not rely on it are preferable. For the
wobbling disk, the approach using Euler angles seems better to me for this reason. Numerically obtained
solutions for Eqs. 8.14 through 8.16 could also, as a last resort, be used to produce animated graphics on a
computer (see exercises at the end of the chapter).
CHAPTER 8. SOME PROBLEMS INVOLVING SINGLE BODIES
86
8.5
Rigid body on an axle
Consider an arbitrary rigid body mounted on a massless, rigid axle (figure 8.4). The axle does not pass through
^n
Q
G
P
Figure 8.4: A rigid body on an axle.
the center of mass of the body; is supported on frictionless bearings; and points along a unit vector n̂ that is
not vertical. We seek the equation governing small unforced oscillations about the stable equilibrium position.
The system has only one degree of freedom. We take as generalized coordinate the angle θ through which
the axle rotates abut n̂. We assume that θ = 0 in the stable equilibrium position. In that position, let the
moment of inertia about the center of mass be Icm,ref , and the position vector from P to G be rG/P,ref (both
assumed known).
Although Lagrange’s equations could be written easily for this system, the Newton-Euler approach is
adopted here to illustrate some useful points. A free body diagram of the system (not shown; draw as an
exercise) would contain the following:
1. The weight −mg k̂ acting (downwards) through the center of mass G,
2. Bearing reaction forces F P and F Q acting at points P and Q respectively, and
3. Bearing reaction moments M P and M Q acting at points P and Q respectively.
It is important to note that, since the bearing are frictionless,
M P · n̂ = M Q · n̂ = 0.
The rotation matrix, in terms of earlier notation, is R(n, θ). The angular velocity and acceleration of the
body are
ω = θ̇ n̂ and ω̇ = θ̈ n̂.
Also, in matrix notation,
Icm = R(n, θ) Icm,ref R(n, θ)T , and rG/P = R(n, θ) rG/P,ref .
The absolute acceleration of the center of mass G is
aG = ω̇ × rG/P + ω × ω × rG/P .
Let us consider angular momentum balance about the point P. We have the vector equation
M P + M Q + rQ/P × F Q + rG/P × (−mg k̂) = rG/P × (m aG ) + I cm · ω̇ + ω × I cm · ω.
It is now useful take the dot product of both sides with n̂ (noting that n̂ is along rQ/P ), to obtain
n̂ · rG/P × (−mg k̂) = n̂ · rG/P × (m aG ) + I cm · ω̇ + ω × I cm · ω .
(8.17)
8.6. EXERCISES
87
Since we are interested in small oscillations, we can drop nonlinear terms at this point. Accordingly, noting
that ω is O(θ̇), we drop O(|ω|2 ) terms. Moreover, we note that
aG = O(θ̈),
R(n, θ) = I + θS(n) + O(θ2 ),
Icm = Icm,ref + O(θ),
and
rG/P = rG/P,ref + O(θ).
Equation 8.17 then simplifies to
n̂ · rG/P × (−mg k̂) = n̂ · rG/P,ref × (m aG ) + I cm,ref · ω̇ .
(8.18)
Correct upto first order, we now write (in matrix notation)
rG/P = rG/P,ref + θ S(n)rG/P,ref
and
aG = θ̈ S(n) rG/P,ref .
Equation 8.18 becomes (recall that k̂ ≡ e3 )
−n S(rG/P,ref ) e3 + nT S(e3 ) S(n) rG/P,ref θ =
T
1 T
nT Icm,ref n
n S(rG/P,ref ) S(n) rG/P,ref +
g
mg
θ̈.
By assumption, θ = 0 is an equilibrium configuration of the above system. Thus, the first term on the left
hand side above is zero. This gives
T
1 T
nT Icm,ref n
n S(rG/P,ref ) S(n) rG/P,ref +
n S(e3 ) S(n) rG/P,ref θ =
θ̈.
g
mg
Using properties of the cross product, the above may be rewritten as
1
Icm,ref
S(rG/P,ref ) S T (rG/P,ref ) +
n θ̈.
nT S(e3 ) S T (rG/P,ref ) n θ = nT
g
mg
(8.19)
Equation
8.19
√ may be understood better with an example. Let the axle be in the Y -Z plane, with
√
n̂ = 1/ 2 ĵ + 1/ 2 k̂. Let the rigid body be a uniform disk of mass m and radius R, such that in the stable
equilibrium configuration the disk is flat, in the X-Y plane. Let the axle pass through the circumference of
the disk, i.e., rG/P,ref = Rĵ. The moment of inertia matrix in the reference configuration is
Icm,ref
⎡
⎤
1 0 0
mR2 ⎣
0 1 0 ⎦.
=
4
0 0 2
Equation 8.19 then gives
θ̈ = −
8.6
4g
θ.
7R
Exercises
1. Use Maple to obtain Eqs. 8.14 through 8.16.
2. A right crcular cone of uniform density, with base circle radius R and slant height h, rolls without slip on
an inclined plane making an angle α with the horizontal. Find the equation of motion, and the frequency of
88
CHAPTER 8. SOME PROBLEMS INVOLVING SINGLE BODIES
small oscillations about the stable equilibrium position.
3. Consider Eqs. 8.14 through 8.16. Take arbitrary initial values for the Euler angles, e.g., φ(0) = 1.0,
θ(0) = 1.3, and ψ(0) = −0.4. Let n̂ref = î. Find n̂, i.e., n̂ref after the initial (1,3,1) rotation (φ(0), θ(0), ψ(0)).
Let ω(0) = n̂ plus some small perturbation not parallel to n̂. Integrate the equations numerically, produce
suitable animated graphics on a computer (using, e.g., Matlab’s movie command), and verify visually that
the wobble rate equals twice the spin rate.
4. Use the analysis of section 8.5 to devise an experimental method for obtaining the moment of inertia matrix
for an arbitrary object.
5. An arbitrary rigid body is suspended by attaching a string at an arbitrary point on its surface (the point
should not be the center of mass, nor lie on a principal axis of the moment of inertia). The other end of
the string is attached to a rigid support. The string is extensible and damped: model it using a spring and
dashpot. Simulate this system in Matlab using arbitrary initial conditions; make a movie for visualization as
well. Initially, there is significant energy dissipation. However, does the energy go to zero? Check numerically
for conservation of the vertical component of the angular momentum about the point of attachment on the
rigid support.
Chapter 9
Friction
When two solid objects rub against one another, frictional forces1 may need to be accounted for. Accounting
for friction is dynamics is difficult at two levels. First, accurate mathematical description of the friction forces
is itself difficult (this is a modeling issue). Second, the solution of equations of motion for frictional systems
is also troublesome (this is an algorithmic issue). The commonest friction model is that of Coulomb friction.
9.1
Coulomb friction
Most undergraduate level mechanics courses describe friction as being modeled using a pair of empirical
coefficients, μs and μk , called the coefficients of static and kinetic friction respectively. It is understood that
μs > μk > 0. The idea behind taking μs > μk is that it is harder to get something to begin sliding than to
keep it sliding. This friction model, called the Coulomb friction model, is described as follows.
Consider two objects A and B in contact at a point P in space. Assume that a well defined tangent plane
exists at P (see figure 9.1). The contact force F P is taken to be the force exerted by body A on body B. The
A
B
Figure 9.1: Two bodies in contact.
unit vector n̂ is taken to be normal to the tangent plane, and pointing into body B. The normal component
of the contact force is given by
Fn = F P · n̂.
In the absence of adhesion, which is what we assume here, Fn ≥ 0. The tangential component of the contact
force is then
F T = F P − (F P · n̂) n̂.
1 My
understanding of friction has been heavily influenced by Andy Ruina. Any errors here are mine, of course.
89
CHAPTER 9. FRICTION
90
Now let the velocity of the material point on A that is instantaneously at P be v P,A . Let the velocity of
the corresponding point on B be v P,B . Then the relative velocity at the contact point is
v P,rel = v P,B − v P,A .
We assume that
v P,rel · n̂ = 0.
If this dot product were strictly positive, then the bodies would separate and contact would be terminated.
Conversely, if the dot product were strictly negative, then a collision would occur (not discussed in these
notes). Equality is needed for sustained contact.
In sustained contact, the usual Coulomb friction model has two possibilities:
1. Either v P,rel = 0, and F T ≤ μs Fn ,
2. Or v P,rel = 0, and
F T = −μk Fn
v P,rel
.
v P,rel Comment: The rigidness of a body is a well understood concept. A rigid body has a mass, a center of mass,
a moment of inertia matrix, and (if convenient) a surface; it can be located by its center of mass position
and its orientation. The force-deflection characteristic of a spring, or equivalently, the potential energy stored
in a deformed spring, is also well understood and mathematically described using Hooke’s law. Notice that
the spring characteristic depends on material behavior. Two springs of identical shapes and density distributions could, in principle, have different stiffnesses. The force-deflection behavior of a spring, as described by
F = kx, is called a constitutive relation. Now, studying friction, we observe that the Coulomb friction model
is a different type of constitutive relation describing the force-displacement behavior of contacting surfaces.
Frictional effects, being localized at or near surfaces, are sensitive to environmental conditions like humidity
and surface cleanliness. Day to day variations in measured friction coefficients can be significant. Friction
models are approximate, and usually much less accurate than (say) the force-deflection model for a spring.
The distinction between static and kinetic friction, with μs > μk , is useful for describing certain commonly
observed things like the squeaking of door hinges. Here, we begin with a simpler example.
9.2
A spring block system
See figure 9.2. A block of mass m moves on a horizontal, frictional plane. A spring of stiffness k is attached
to the block. The free end of the spring is moved by an external agent at a constant speed v. The coefficients
v
k
m
immovable floor
μ
vrel
Figure 9.2: Sliding block.
of friction are μs and μk < μs . In the figure, the sketch to the right shows μ against sliding speed: here, μ
should be taken to mean
−f
,
μ=
mg
9.2. A SPRING BLOCK SYSTEM
91
where f is the friction force (measured positive to the right), and mg is the weight of the block.
We will first proceed analytically for a bit; and then do less tedious (though approximate) numerics.
The governing equation is (taking the displacement of the block as x, positive to the right):
mẍ = k(vt − x) + f,
(9.1)
where the friction force f (or rather its horizontal component, taken positive when acting to the right) is given
by the Coulomb friction relations
either
or
ẋ
, ẋ = 0,
|ẋ|
|f | ≤ μs mg,
ẋ = 0.
f = −μk mg
(9.2)
Now suppose that there is some period of sticking, i.e., when ẋ ≡ 0. The spring continues to get compressed
until slip begins. However, at that instant, ẋ = 0. Se we need to add on, for incipient sliding,
f = −μk mg
ẍ
,
|ẍ|
ẍ = 0 but ẋ = 0.
Special cases where ẍ = 0 at incipient sliding can be imagined, but are rare in practice and ignored here. This
gives the friction force f as
⎫
ẋ
⎪
either f = −μk mg , ẋ = 0,
⎪
⎪
⎬
|ẋ|
or
|f | ≤ μs mg,
ẋ = 0, ẍ = 0,
(9.3)
⎪
⎪
ẍ
⎪
or
f = −μk mg , ẋ = 0, ẍ = 0. ⎭
|ẍ|
Equations 9.1 and 9.3, due to the nonsmoothness of friction, are numerically more troublesome than
anything encountered so far in these notes. Systems with many frictional contacts, where each contact can
involve sliding, sticking, or transitions from one to the the other, are genuinely difficult to tackle numerically.
We first seek steady solutions for Eqs. 9.1 and 9.3. Try x = vt − xss , where xss is a constant. This means
ẋ > 0, and so from Eq. 9.3 we have
f = −μk mg.
Equation 9.1 becomes
0 = kxss − μk mg,
whence
xss =
μk mg
.
k
Thus, steady sliding is possible.
We then consider small deviations from steady sliding. Accordingly, we try
x = vt −
μk mg
+ η(t),
k
where η(t) is assumed small enough that |η̇| < v, i.e., ẋ remains strictly positive. Then we still have
f = −μk mg.
Equation 9.1 becomes
mη̈ = −kη,
which has oscillatory solutions of constant amplitude. Any viscous damping in the spring will eliminate these
oscillations, and so we may consider the steady sliding solution to be locally stable: small disturbances die
out, and the solution again approaches steady sliding.
The above two solutions are unaffected by μs . Let us now imagine that the solution begins with zero initial
conditions, i.e., x(0) = 0 and ẋ = 0. Then initially, ẍ = 0 as well. From Eq. 9.1 we have f = −kvt, which for
small enough t is smaller in magnitude than μs mg.
CHAPTER 9. FRICTION
92
At time
μs mg
,
kv
the block begins to slide to the right; at that instant, ẋ = 0 but ẍ > 0, and from Eq. 9.3 we have
t1 =
f = −μk mg.
During the subsequent phase of sliding to the right, Eq. 9.1 becomes
mẍ = k (vt − x) − μk mg.
The above equation governs the motion until such time when ẋ becomes zero again. During this phase of
sliding, the solution is given by
x = vt −
μs g 1 2
μs g 3
v
sin ωt −
− 2 μk + (μs − μk ) cos ωt −
,
ω
ωv
ω
ωv
/
where ω = k/m.
We seek the time t2 when ẋ becomes zero again. Skipping some routine calculations, we find
2ωvg(μk − μs )
μs g = 2 2
sin ωt2 −
.
ωv
ω v + g 2 (μk − μs )2
The right hand side above, being of the form 2ab/(a2 + b2 ), is guaranteed to have magnitude ≤ 1, and so there
is a real solution for t2 . Proceeding in this way, a determined analyst may be able to describe the solution
obtained. We will now retreat to a numerical solution.
We choose some arbitrary parameter values:
m = 1, k = 1, g = 1, μk = 0.3, μs = 0.4, and v = 0.05.
(9.4)
For these parameter values, t1 = 8, and t2 = 12.0689. At t = t2 , we find x = 0.40344, and k(vt − x) = 0.2.
Since the latter is smaller in magnitude than μs mg = 0.4, the block stays stuck for a while. Then the above
cycle of sliding followed by a period of sticking is repeated. Such stick-slip motions will survive in the presence
of a small amount of viscous damping in the spring.
A fully numerical solution of the above problem, or of similar problems with, e.g., other nonlinearities which
make analytical progress impossible, will require careful numerics. Numerical solutions of ODEs will usually
be obtained at some discrete points in time; and the length of some time steps will have to be iteratively
adjusted so that every instant when ẋ either becomes zero or changes from zero coincides with one of these
discrete points.
For less careful numerical work, a simpler approach is possible. We can choose a small approximation
parameter 0 < 1 (where the “
” stands for “much smaller than”), and approximate the friction law as
follows (the choice is somewhat arbitrary):
1
3v
v2
v
− 2 .
(9.5)
+ (μs − μk ) exp
μ(v) = μk tanh
2 2
The above approximation is shown, for μs = 0.4 and μk = 0.3, in figure 9.3, where again we understand that
μ=
−f
.
mg
For large ± v/, the above approximation gives ± 0.3 as desired. The maximum magnitude, attained close to
± 1, is very close to 0.4 (and could be made closer with a minor correction, avoided here for simplicity). With
the above smooth approximation, we can write the equation of motion (Eq. 9.1) as
,
#
1
3ẋ
ẋ2
ẋ
−
.
(9.6)
+ (μs − μk ) exp
mẍ = k(vt − x) − mg μk tanh
2 2 2
9.2. A SPRING BLOCK SYSTEM
93
0.4
0.3
0.2
μ
0.1
0
- 0.1
- 0.2
- 0.3
- 0.4
-8
-6
-4
-2
0
2
4
6
8
v/ε
Figure 9.3: Smooth approximation for friction coefficient.
0.18
1.2
0.16
−1
0.12
0.6
velocity
0.8
velocity
displacement
10
0.14
1
0.1
0.08
0.06
0.4
−3
0.04
0.2
0
−2
10
10
0.02
0
10
20
30
time
0
0
10
20
time
30
5
10
time
15
Figure 9.4: Numerical solution with smooth approximation for friction coefficient.
We use the parameter values of Eq. 9.4, initial conditions x(0) = 0 and ẋ(0) = 0, and choose = 0.001 (smaller
gives greater accuracy, but requires finer step sizes in time). The numerical solution obtained is shown in
figure 9.4. The numerical solution was obtained using Matlab’s ode45 with error tolerances set to 10−8 ; this
routine has adaptive step sizing for error control, which is required for such cavalier calculations. However, the
end results are satisfactory. The displacement plot (left) shows periodic near-sticking and gross sliding. The
velocity plots (middle and right) show the instants when gross slip first starts and first ends. These match the
foregoing exact calculations acceptably well. The periods of slowly increasing velocity, in regimes where the
exact solution has exact sticking, are due to the numerical approximation of the friction law (for extremely
slow speeds, the approximation is more like viscous friction than dry friction).
The match will be better for smaller .
CHAPTER 9. FRICTION
94
9.3
On distinguishing between µs and µk
In what follows, we will not distinguish between μs and μk for two reasons.
The first reason is philosophical. The distinction between μs and μk , while convenient for modeling stickslip oscillations as discussed above, actually arises from what is called a distinguished limit. To understand
how, consider that all bodies are in fact deformable. With this in mind, consider for example the single block
of mass m in the foregoing example (figure 9.2) to actually be two masses of m/2 each, connected by a spring
of high stiffness K. Consider also the friction law to be governed by a small parameter as used in the
smooth approximation above (Eq. 9.5). Now if we first imagine K → ∞ and then imagine → 0, we get
the behavior described above. However, if we first consider → 0 and then consider K → ∞, then we get a
different behavior where sliding starts earlier. If we imagine n masses of m/n each, connected in a chain with
springs each of high stiffness K, then sliding starts earlier still. In this way, the limiting case of a rigid block
with Coulomb friction coefficients μs > μk gives results that depend on the relative rates at which K and approach infinity and zero respectively (a distinguished limit). In addition to these theoretical considerations,
it is also known from careful experiments that the coefficient of friction at any steady sliding speed is actually
a function of that speed (as was assumed through an arbitrary functional form in Eq. 9.5). However, on
sudden changes in the sliding speed, there is some transient behavior in the observed friction coefficient which
has magnitude (and importance) comparable to the change in the steady value. In other words, modeling the
friction coefficient as a function of sliding speed alone ignores important physical phenomena of magnitude
comparable to the variation predicted by the model itself.
The second reason for not distinguishing between μs and μk is practical. A lot of frictional effects can be
approximately modeled and understood using a single friction coefficient (μs = μk ). Moreover, whether we
take μs = μk or not, the model is still only an approximate representation of reality (much more so than,
say, using a linear stiffness for a typical spring). For many problems, nothing is really lost on ignoring the
difference.
For these two reasons, we will not distinguish between μs and μk in the rest of this book.
9.4
Vibration damping
Consider the same frictional spring-block system of figure 9.2, except that now μs = μk = μ > 0 (i.e., μ is now
a constant positive parameter); and that the free end of the spring, instead of being given a steady velocity
v, is now held fixed (i.e., v = 0). Let x(0) = A0 > μmg/k, and ẋ(0) = 0. The initial condition is such that
sliding begins immediately, and x begins to decrease. Sliding eventually stops for some x = −A1 < 0. The
k 2
A0 − A21 . This decrease must be
change in kinetic energy is zero, and the decrease in potential energy is
2
equal to the energy dissipated due to friction, which in this case is μmg (A0 + A1 ), giving
A0 − A1 =
2μmg
.
k
Thus, the interval of sliding results in a constant reduction of oscillation amplitude independent of the amplitude itself. The time duration of sliding is also independent of the oscillation amplitude (it is exactly half the
time period of oscillation in the frictionless case; showing this is left as an exercise).
If A1 is large enough for sliding to start again, then a similar motion as above results in a new amplitude
A2 that satisfies
2μmg
.
A1 − A2 =
k
A key point is that the oscillation amplitude decreases linearly with time, as opposed to exponentially (familiar
from viscous damping). Eventually, sliding stops at a value of displacement that is typically nonzero, but small
enough that sliding does not restart.
An approximate (smooth) solution along the lines of figure 9.3, but with μs = μk = μ (say), could be
based on solving the equation (compare with Eq. 9.6)
3ẋ
ẍ = −x − μ tanh
,
(9.7)
9.5. RESONANCE
95
where the 3 inside the “tanh” has been retained for easier comparison with the foregoing discussion, but could
also be absorbed into the parameter .
Solving the above using Matlab’s ode45 with tight error tolerances, for m = 1, g = 1, k = 1, = 0.001,
μ = 0.6, x(0) = 6.55 and ẋ(0) = 0, gives the results shown in figure 9.5. The solution is shown with a
8
6
displacement, x
4
2
0
−2
−4
−6
0
5
10
15
20
25
30
time, t
Figure 9.5: Solution of Eq. 9.7.
heavy line. By the foregoing discussion, on noting that 2μmg/k = 1.2 here, we have A0 = 6.55, A1 =
−(6.55 − 1.2) = −5.35, A2 = 5.35 − 1.2 = 4.15, A4 = −(4.15 − 1.2) = −2.95, A5 = 2.95 − 1.2 = 1.75 and finally
A6 = −(1.75 − 1.2) = −0.55, where the block remains stuck. Dashed horizontal lines mark these values; it is
seen that the approximate solution matches these to plotting accuracy.
9.5
Resonance
Consider the same frictional spring-block system of figure 9.2, except that now μs = μk = μ > 0 (i.e., μ is now
a constant positive parameter); and that the free end of the spring, instead of being given a steady velocity v,
is given an oscillatory motion, i.e., vt is replaced by U sin ωt. Resonance occurs if ω = 1. In viscously damped,
forced, linear vibrating systems, it is well known that at resonance the amplitude of vibration is inversely
proportional to the damping present. What happens with dry friction?
So consider (compare with Eq. 9.1)
mẍ = k(U sin ωt − x) + f,
(9.8)
with f given by Eq. 9.3, modified as follows because μs = μk
either
or
or
ẋ
, ẋ = 0,
|ẋ|
|f | ≤ μmg,
ẋ = 0, ẍ = 0,
ẍ
f = −μmg , ẋ = 0, ẍ = 0.
|ẍ|
f = −μmg
⎫
⎪
⎪
⎪
⎬
⎪
⎪
⎪
⎭
(9.9)
CHAPTER 9. FRICTION
96
As a particular example, let m = 1, g = 1, k = 1, μ = 0.6, and ω = 1. Is the dry friction able to contain the
resonant growth of the oscillation?
To anticipate the answer through our foregoing smoothened numerical approximation, we consider (letting
= 0.001 and dropping the 3 inside the “tanh”)
ẋ
ẍ = −x + U sin t − 0.6 tanh
.
(9.10)
0.001
The numerical solution obtained with the arbitrarily chosen initial conditions x(0) = 10 and ẋ(0) = 0, with
U = 0.4, is shown in figure 9.6. No resonant growth is observed. Another numerical solution for x(0) = 10
10
8
6
displacement, x
4
2
0
−2
−4
−6
−8
−10
0
10
20
30
40
50
60
time, t
Figure 9.6: Solution of Eq. 9.10 with U = 0.4.
and ẋ(0) = 0, but with U = 1, is shown in figure 9.7. It is clear that there is resonant growth in this case.
These issues may be understood with some qualitative analysis as follows. Imagine that at some instant of
time, somehow or other, x has become large. Let this largeness be represented by some quantity A 1, where
the “” stands for “much greater than”. Then we can change variables to y = x/A, so that y is comparable
to 1. The governing equation for y is (compare with Eq. 9.8)
0
f
k
kU
sin
t+ .
(9.11)
mÿ = −ky +
A
m
A
We see that the last two terms above are small for large A, because both kU and f are bounded while we can
think of A being as large as we like.
Now we use an analogy. If a differentiable function g(x, y) changes infinitesimally due to infinitesimal
changes in x and y, then the individual changes in g due to changes solely in x and solely in y can be added
up to find the total change in g. In calculus notation,
Δg ≈
∂g
∂g
Δx +
Δy.
∂x
∂y
Should such reasoning hold for our system’s response, we could use it to simplify our analysis. Let us assume
that it does. (Such reasoning could be formalized using a perturbation expansion, which we avoid here.)
9.5. RESONANCE
97
10
8
6
displacement, x
4
2
0
−2
−4
−6
−8
0
10
20
30
40
50
60
time, t
Figure 9.7: Solution of Eq. 9.10 with U = 1 (heavy line). The dashed line has slope 0.118 and is drawn for
comparison.
By the above, analysis of the resonance reduces to calculating the growth rate in the amplitude due to the
resonant forcing; calculating the decay rate in the amplitude due to the frictional damping; and subtracting
the latter from the former to find the amplitude growth rate. Having identified this strategy, we can drop the
A-scaling and return to the x-equation.
For resonant forcing, we now consider the equation (using, for simplicity, the same parameter values as
above)
ẍ + x = U sin t,
whose solution is
x=−
Ut
cos t + C0 sin t + C1 cos t.
2
Thus, the amplitude of the cos t component in the solution grows by U π for every period of forcing (which is
2π).
In the same period of time, considering the frictional damping independent of the resonant forcing, we
have found above that the amplitude decreases (after two phases of sliding) by an amount 4μmg/k, which for
our parameters is 4μ = 2.4.
The resultant growth rate of the amplitude is therefore predicted to be U π − 2.4 for every 2π units of
time. In the numerical solutions above, we used U = 0.4 and U = 1. For U = 0.4 the growth rate is negative,
resonance is suppressed, and oscillations cannot stay large (as supported by numerics). For U = 1, the growth
rate is predicted by our simple analysis to be (π − 2.4)/2π = 0.118 per unit time, which agrees well with
numerics (the dashed line has slope 0.118 for comparison).
A comparison with viscous damping is useful here. Consider the damped, forced, single degree of freedom
system given by
ẍ + cẋ + x = F sin ωt.
Resonant dynamics with large amplitudes is expected for ω close to 1. If amplitudes start off at some large
CHAPTER 9. FRICTION
98
value A, then we put y = x/A to get
F
sin ωt.
A
In the above, for fixed F and large enough A, we essentially get the damped and unforced system
ÿ + cẏ + y =
ÿ + cẏ + y = 0,
whose solutions decay exponentially in amplitude until A is not so large any more. Thus, unlike damping
through dry friction, viscous damping is able to restrict resonant amplitudes even if F is fairly large.
9.6
Sliding on moving surfaces
The behavior of objects sliding on surfaces that are themselves moving can be interesting and possibly counterintuitive. We study some examples in this section.
Example 1: Vibration induced regularization of dry friction.
Consider a block on a horizontal table. The table oscillates horizontally with a displacement given by
U sin ωt. Let the displacement of the block with respect to the table be x. Then its absolute displacement is
U sin ωt + x. The equation of motion is
m −ω 2 U sin ωt + ẍ = f,
where f is given by Eq. 9.9. As a particular example, let m=1, g = 1, U = 0.1, ω = 20 and μ = 0.4. Using
our previous smooth approximation for the friction, we write
ẋ
ẍ = 40 sin 20t − 0.4 tanh
.
0.001
Let the initial absolute displacement of the block be zero, i.e., x(0) = 0; and the initial absolute velocity of
the block be 1.5 (in suitable units), so that ẋ(0) + ωU = 1.5, or ẋ(0) = −0.5. The numerical solution for the
velocity (absolute, or relative to ground) for these initial conditions is shown in figure 9.8. It is seen that the
velocity has fluctuations due to the oscillations in the table, but on average decreases roughly exponentially
with time as opposed to linearly with time (which is what would happen if the table was not oscillating).
To understand this behavior, we carry out an approximate analysis. Assume that the velocity of the table,
which is on the order of ωU = 2 in our case, is large compared to the velocity of the block (absolute, or
relative to ground). The block’s velocity itself changes slowly on average, and so let us think of it as some
small quantity v > 0 whose variation can temporarily be ignored.
In an interval when ωU cos ωt > v, the block slides to the left relative to the table, and the friction force
on the block acts to the right. Letting ωt lie between −π/2 and π/2, we find that
cos ωt >
or
|ωt| < cos−1
v
,
ωU
v .
ωU
For small v,
v
π
−
.
2
ωU
Thus, for one cycle of oscillation of the table, the phase of oscillation changes by
ωtmax ≈
π−
2v
ωU
π+
2v
ωU
while sliding is to the left, and by
9.6. SLIDING ON MOVING SURFACES
99
1.6
numerical
averaged
1.4
velocity, dx/dt + ω U cos ω t
1.2
1
0.8
0.6
0.4
0.2
0
0
2
4
6
8
10
12
14
16
18
20
time, t
Figure 9.8: Block sliding on a horizontally oscillating table.
while sliding is to the right. Since the friction force changes direction with direction of sliding, but always has
the same magnitude, the net effect of the friction force is that it acts to the left over a net phase change of
4v
ωU
out of every 2π. The average decelerating force then opposes the average velocity of the block, and has a
magnitude equal to
2μv
.
πωU
The averaged dynamics of the block is therefore
mv̇ = −
2μv
,
πωU
which in our case becomes
v̇ = −0.1273 v.
The solution of the above with v(0) = 1.5 is also shown in figure 9.8. The match is acceptable.
The fast, small, in-plane vibration of the table serves to regularize the dynamics of the sliding block, giving
a smooth, viscous-like friction effect on average. This idea is used in industrial applications. In the home, if
we pour out some sugar from a can, the discontinuous nature of dry frictional contact makes it difficult to
control the flow of sugar; however, small but rapid side to side shaking of the can causes the sugar to flow out
more smoothly, like a fluid on average. Though the average flow of the sugar is perpendicular to the oscillatory
motion of the can, the idea is similar.
Example 2: A vibratory conveyor.
CHAPTER 9. FRICTION
100
In the previous example, what if the oscillation of the table is not symmetrical in the two directions? Let
the displacement of the table be
U × (sin ωt − 0.33 sin 2ωt + 0.08 sin 3ωt).
A plot of
g(ωt) = sin ωt − 0.33 sin 2ωt + 0.08 sin 3ωt
versus ωt is given in figure 9.9, and shows the asymmetry.
1.5
1
g(ωt)
0.5
0
−0.5
−1
−1.5
0
2
4
6
phase, ωt
8
10
12
Figure 9.9: Asymmetrical oscillation of table.
The equation of motion now is
$
%
m −ω 2 U (sin ωt − 4 × 0.33 sin 2ωt + 9 × 0.08 sin 3ωt) + ẍ = f.
(9.12)
As above, we take m=1, g = 1, U = 0.1, ω = 20 and μ = 0.4. Noting that the absolute velocity of the
block at t = 0 is
v(0) = ẋ(0) + 2(1 − 0.66 + 0.24) = ẋ(0) + 1.16,
we let v(0) = 0, i.e., ẋ(0) = −1.16. We also let x(0) = 0. The numerical solution is approximated, as before,
using a smooth approximation. The results are shown in figure 9.10. It is seen that the block acquires an
average steady state forward velocity and stays close to it thereafter. After the previous example, this one is
not surprising. The velocity of the table is to the right for longer periods of time, and to the left for shorter
periods of time; assuming (for simplicity) that the table slips relative to the block at all times, there is a
constant friction force acting to the right for a relatively long time, and the same force acting to the left for a
relatively short time. The result is a net acceleration to the right.
The steady state velocity attained depends on some finer details of the velocity history, and we may estimate
it as follows. The block slides to the right relative to the table, and on average. But the time-averaged friction
force is zero, because (recall that v is the absolute velocity)
v̇ = f.
velocity, dx/dt + 2 (cos(20t) 0.66 cos(40t) + 0.24 cos(60t) )
9.6. SLIDING ON MOVING SURFACES
101
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0
2
4
6
8
10
12
14
16
18
20
time, t
Figure 9.10: Velocity of block under symmetrical oscillation of table.
This means that
ωt=2π
sgn(ẋ)dt = 0,
(9.13)
0
where sgn(ẋ) = 1 if ẋ ≥ 0 and sgn(ẋ) = −1 otherwise (we are ignoring the possibility of extended periods of
time when ẋ ≡ 0). Now consider Eq. 9.12. The right hand side is bounded in magnitude by μmg = 0.4 in our
case. The left hand side has a term on the order of ω 2 U = 40 in our case, i.e., a hundred times bigger. This
means the right hand side could be dropped for some approximation purposes (clearly, dropping it altogether
would change the dynamics completely, but we will not come to that). Dropping the right hand side and
integrating the left with respect to t, we find
ωU (cos ωt − 0.66 cos 2ωt + 0.24 cos 3ωt) + ẋ + C = 0,
where C is an integration constant. Under steady state conditions, ẋ is periodic but has a net nonzero average
value va to the right. This means C = −va . Now Eq. 9.13 implies that va must satisfy the equation
ωt=2π
sgn {va − ωU (cos ωt − 0.66 cos 2ωt + 0.24 cos 3ωt)} dt = 0,
0
which for our parameter values is equivalent to
2π
sgn {va − 2(cos τ − 0.66 cos 2τ + 0.24 cos 3τ )} dτ = 0.
0
Recognizing an idea from elementary statistics, we note that va must be equal to the median value of 2(cos τ −
0.66 cos 2τ + 0.24 cos 3τ ). A quick numerical estimate of the median (from Matlab2 ) is
va = 1.069,
2 >>
t=linspace(0,2*pi,20002); t=t(1:20001); y=cos(t)-0.66*cos(2*t)+0.24*cos(3*t); median(y)
CHAPTER 9. FRICTION
102
which matches satisfactorily with figure 9.10. In deterministic mechanics, it is not as common to encounter
the median as it is to encounter the mean of a time-varying quantity.
Example 3: A frictional harmonic oscillator.
L
x
D
μ N2
μ N1
mg
N1
N2
Figure 9.11: A frictional harmonic oscillator.
See figure 9.11 (top). Two identical rollers separated by a horizontal distance D spin rapidly in opposite
directions. A uniform rod of length L is placed on the rollers. The friction coefficient between the rod and
each roller is μ. The center of the rod is located a horizontal distance x to the right of the midpoint between
the rollers. We assume that the velocity of the rod is always small enough that steady unidirectional slip is
sustained at each roller contact.
A free body diagram of the system is shown in figure 9.11 (bottom). There are shown two normal reaction
forces N1 and N2 at the two roller contacts; opposing friction forces μN1 and μN2 as shown; and the downward
acting weight of the rod.
Since the vertical acceleration of the center of mass is zero,
N1 + N2 = mg.
Since the rate of change of angular momentum of the rod about the center of mass is zero,
N1
D
D
+ x = N2
−x .
2
2
Finally, linear momentum balance in the horizontal direction gives
μN1 − μN2 = mẍ.
Solving the above for N1 , N2 and ẍ, we obtain
ẍ = −
which is the equation for a harmonic oscillator.
2gμ
x,
D
9.7. ROLLING FRICTION
9.7
103
Rolling friction
When a sphere rolls without slip down an inclined plane, there is a friction force at the contact. This is
sometimes mistakenly referred to by beginning students as rolling friction, because in the dynamics of this
system there is both rolling and friction. However, this is the usual frictional interaction at point contacts
between two bodies, and the overall rolling of the sphere plays no role in the contact model (e.g., the Coulomb
friction model). In the usual model for such contact, there is no energy dissipation even though there is a
friction force.
Sometimes we are interested in modeling finer effects. For example, if a sphere rolls over a long distance
on a flat surface, then it eventually comes to a stop although there is no macroscopic dissipative sliding at the
contact. Accurate modeling of this dissipation may be based on recognizing that the sphere and flat surface
actually deform slightly, and possibly inelastically; that contact actually occurs over a region of small but
nonzero size; and that different points in this region may actually be either sticking or sliding governed by
sliding friction models such as, again, the Coulomb friction model. Even within rigid body mechanics, we may
wish to develop a model for such dissipation over relatively long distances of rolling.
To this end, consider a uniform sphere rolling on a horizontal surface (figure 9.12, left). A free body
diagram of the sphere is drawn on the right. Within rigid body mechanics, the sphere does not deform. The
contact force acts directly below the center of mass, and consists of a normal reaction R and a frictional
W
M
F
R
Figure 9.12: Rolling friction.
retarding force F . We have also allowed for a contact moment M , although we do not usually incorporate
contact moments in point contacts between rigid bodies.
The normal reaction R must exactly counteract the weight W , because the vertical acceleration of the
center of mass is zero. The center of mass moves to the right and decelerates, so F must be positive in the
direction shown. However, on applying angular momentum balance, we see that in the absence of M , the
friction force F tends to increase the angular velocity while decreasing the linear velocity, leading to slip. Since
the cylinder rolls without macroscopic slip, absence of M implies absence of F . To incorporate dissipation
in rigid bodies rolling without slip, we must allow for contact moments. Assuming that there is no adhesion
between the contacting bodies, the moment M arises because the normal force R actually acts a little bit to
the right of the center of mass; this is in turn possible if the contact region, though small, has a nonzero size.
Rolling friction in rigid body mechanics is usually small because the contact patch is small. This topic is
not discussed further in these notes.
9.8
Planar sliding
Let us consider a 2D problem. Frictional sliding of flat objects on flat tables turns out to be mathematically
quite complicated.
Disk: Let us begin with a disk. On encountering serious analytical difficulties, we will retreat to a rod.
A flat uniform disk of mass m and radius R slides and rotates on a flat table; the coefficient of friction is
μ. See figure 9.13. Let the center of mass of the disk be G, and let P be a point on the disk located with
CHAPTER 9. FRICTION
104
y
R
r
P
θ
x
G
Figure 9.13: A disk sliding on a table.
respect to stationary x and y axes (fixed to the table) using polar coordinates r and θ. Let the instantaneous
velocity of G be v G = vGx î + vGy ĵ and the angular velocity of the disk be ω k̂. How do v G and ω evolve over
time?
The problem is statically indeterminate because the distribution of normal pressure over the contact area
cannot be found within rigid body dynamics. Let us assume, for simplicity, that the pressure distribution is
a function of radial position alone, i.e., p = p(r) is symmetrical about the center of the disk.
The velocity of point P is
v P = v G + ω k̂ × rP/G = (vGx − ωr sin θ) î + (vGy + ωr cos θ) ĵ.
We will assume that v P is nonzero for almost all points. The friction force per unit area at some typical point
P is then
(vGx − ωr sin θ) î + (vGy + ωr cos θ) ĵ
f (r, θ) = −μp(r) 1
.
2
2
(vGx − ωr sin θ) + (vGy + ωr cos θ)
Linear momentum balance in the plane of motion gives
mv̇ G =
R
2π
f (r, θ) r dθ dr.
0
0
Angular momentum balance about G gives
1
mR2 ω̇ =
2
0
R
2π
3
2
(r cos θî + r sin θĵ) × f (r, θ) · k̂ r dθ dr.
0
The above integrals, when evaluated, yield nonlinear differential equations. Note that the pressure is probably
not purely a function of radius anyway.
Motions of such sliding objects have been the subject of research papers (by, e.g., Suresh Goyal and Andy
Ruina). We will not discuss them further.
Rod: Consider the somewhat simpler case of a rod (e.g., a flat ruler) sliding on a table.
See figure 9.14. A uniform rod of length L slides and rotates on a table with friction coefficient μ. The
center of mass of the rod is at G, and the rod instantaneously makes an angle θ with the x-axis. The angular
velocity of the rod, ω = θ̇, is assumed to be in the positive z-direction, i.e., counterclockwise.
This problem, like the disk, is also statically indeterminate because the distribution of normal reaction
forces from the table cannot be found within rigid body dynamics. For simplicity and illustration, let us
assume that the normal reactions are uniform over the length of the rod.
9.8. PLANAR SLIDING
105
y
P
G
r
θ
x
Figure 9.14: A rod sliding on a table.
Consider an arbitrary point P a distance r along the rod from G. Along the lines of the calculations for
the disk, the velocity of point P is
v P = v G + ω k̂ × rP/G = (vGx − ωr sin θ) î + (vGy + ωr cos θ) ĵ.
Assume that v P is nonzero almost everywhere. The friction force per unit length at some typical point P is
then
mg (vGx − ωr sin θ) î + (vGy + ωr cos θ) ĵ
1
.
f (r) = −μ
L
2
2
(vGx − ωr sin θ) + (vGy + ωr cos θ)
Now linear momentum balance in the plane of motion gives
L/2
mv̇ G =
f (r) dr.
−L/2
Angular momentum balance about G gives
L/2 2
3
1
(r cos θî + r sin θĵ) × f (r) · k̂ dr.
mL2 ω̇ =
12
−L/2
The above integrals can be evaluated (e.g., by Maple) and give the following equations of motion:
,
#
2A3 + A2 − ωL
μg
v̇ Gx =
(9.14)
(2vGx + 2A3 sin θ) ln
+ (A1 − A2 ) sin θ ,
2ωL
2A3 + A1 + ωL
,
#
2A3 + A2 − ωL
μg
(9.15)
v̇ Gy =
(2vGy − 2A3 cos θ) ln
− (A1 − A2 ) cos θ ,
2ωL
2A3 + A1 + ωL
3μg
A4 ,
(9.16)
ω̇ =
2ω 2 L3
where
1
2 + 4v 2 ,
A1 = ω 2 L2 − 4vGx ωL sin θ + 4vGy ωL cos θ + 4vGx
Gy
1
2 + 4v 2 ,
A2 = ω 2 L2 + 4vGx ωL sin θ − 4vGy ωL cos θ + 4vGx
Gy
A3 = vGy cos θ − vGx sin θ,
and
A4
%
$
2
2
ln
= − (2 cos 2θ − 6)vGy
− 8A3 cos θvGy + (2 cos 2θ + 2)vGx
2A3 + A2 − ωL
2A3 + A1 + ωL
+2A3 (A1 − A2 ) − ωL (A1 + A2 ) .
The above nonlinear equations are not studied further here; but they indicate the complicated possibilities for
the disk.

Similar documents

×

Report this document