English

not defined

no text concepts found

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 ﬁnd 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, ﬁgures, 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 ﬁnd 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 Inﬁnitesimal 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 inﬁnitesimals . . . 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 brieﬂy 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 scientiﬁc work with computers. We will set ourselves the minimal goal of using vectors and matrices (in place of tensors), which suﬃces 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 deﬁne 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 deﬁne these quantities. The reader may also have some idea about what a frame of reference is, but we will provide a working deﬁnition 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 speciﬁc 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 satisﬁes 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 satisﬁed 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 deﬁnition 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 ﬁxed in some frame of reference that should be clear from the context. Usually, O will be the origin of some coordinate system ﬁxed 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 ﬁgure 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 ﬁgure 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 ﬁnd FA = 5.9229 N, FB = FC = 3.6897 N. Problem: See ﬁgure 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 Diﬀerentiating, mv v̇ + mg ż = 0. (1.4) From ﬁgure 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 deﬁne some unit vectors shown in ﬁgure 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 ﬁgure 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 ﬁgure 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 deﬁnitions 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 diﬀerentiation: 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 ﬂuid 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 ﬁgure 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 ﬁgure 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 eﬀect 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 ﬂuids 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 ﬁgure 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 ﬁxed at this time. For any such choice of points, there will be a set of forces and moments that represent the net eﬀect of contact with the cut-oﬀ portions, and changing these points will generally change the moments as well. In ﬁgure 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 ﬁrst 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 ﬁnd 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 identiﬁed 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 ﬁnally 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 ﬁgure 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 oﬀ 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 ﬁnd 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 diﬀerence 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 oﬀset 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 inﬁnitesimal region where it instantaneously experiences inﬁnitely 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 coeﬃcient 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 diﬀerent from that of coordinate system. However, if we choose and rigidly ﬁx 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 ﬁxed (stationary) frame XY Z. Let a vector A be ﬁxed in xyz. Then A is not ﬁxed 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 ﬁrst deﬁne angular velocity. Here we will use an intuitive deﬁnition, which will later be given in a form better suited to actual calculation. In an inﬁnitesimal time period Δt, it is intuitively clear that reference frame xyz (like any other rigid body) undergoes an inﬁnitesimal rotation. That rotation involves some inﬁnitesimal angle (say Δθ), and occurs about some unit vector n̂. We deﬁne Δθ n̂. (2.2) ω = lim Δt→0 Δt See ﬁgure 2.1(a). The initial position of A is marked by points OP, and the inﬁnitesimally 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 ﬁgure. The projection of A onto this plane lies along a unit vector which we denote by λ̂. A third unit vector ê is then deﬁned as n̂ × λ̂. Let the angle between A and n̂ be φ. See ﬁgure 2.1(b). L is the length of arc between points P and Q of ﬁgure 2.1(a), and is given by L = |A| sin φ Δθ, and so the inﬁnitesimal 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 ﬁxed 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 ﬁxed in xyz, even as xyz rotates relative to XY Z. In this case the total inﬁnitesimal change of A during Δt is the sum of (i) the change of A as occurring in xyz, and (ii) the change in the inﬁnitesimally 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 ﬁxed in the initial conﬁguration (facing east), and (ii) the change in the second hand obtained by holding it ﬁxed, 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-inﬁnitesimal. So, in an inﬁnitesimal time interval Δt, we think of the change in A as occurring in two parts. The ﬁrst 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 diﬀerent 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 inﬁnitesimals. 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 ﬁgure 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 ﬁgure 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 ﬁgure 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 diﬀerent 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 ﬂashes of intuition interspersed with indeﬁnite 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 ﬁnd 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 ﬁgure 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 deﬁnition 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 ﬁgure 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 ﬁxed 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 diﬀerent 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 ﬁrst term on the right hand side of the ﬁrst equation is zero because Ω, and therefore λ (as deﬁned above), is constant. Substituting for λ and n̂, αcone/XY Z = 2.4 √ Ω2 L2 − R2 î. R Addition of angular velocities See ﬁgure 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 speciﬁed 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 ﬁxed at the origin of xyz. Then drP/O dt *0 drP/O = + ω B/F1 × rP/O B dt F1 where the ﬁrst term on the right hand side is zero since P is ﬁxed 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 ﬁrst term on the right hand side is zero because P is ﬁxed 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 satisﬁes 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 ﬁxed frame XY Z. Consider a moving point P , as shown in ﬁgure 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 diﬀerentiation 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 ﬁxed 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 brieﬂy, 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 ﬁve term acceleration formula, and will be useful in what follows. In order to understand the ﬁve 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 ﬁgure 2.6. A turntable rotates with constant angular velocity ω about a vertical axis (ignore α and the soccer goalpost for the moment). We ﬁx 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 ﬁxed 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 ﬁgure. 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 ﬁxed 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 ﬁrst is called the Coriolis term (which, e.g., is important for monsoon winds); and the second aﬀects 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 ﬁgure 2.7. A turntable rotates at a constant angular rate Ω. A wheel is located on the turntable as shown in ﬁgure 2.7. The wheel spins at a rate of β̇ about an axis ﬁxed relative to the turntable. Find the acceleration of a point P on the rim of the wheel. Solution: See ﬁgure 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 ﬁve 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 ﬁxed 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 ﬁxed in xyz. Since the wheel is rotating at an angular rate of β̇ relative to xyz, and since vector rQP is ﬁxed to the wheel, application of Eq. 2.5 (using the wheel as the rotating frame and xyz as the ﬁxed frame) gives drQP = −β̇ î × rQP . dt xyz 2.6. EXERCISES 27 The right hand side above plays the role of ρ̇ in the ﬁve 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 ﬁnd 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 ﬁxed 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 ﬁve 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 ﬁgure 2.9. Now treating xyz as the ﬁxed 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 ﬁxed in xyz, the ﬁrst term is zero. Since P is ﬁxed 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 ﬁgure 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, ﬁnd the angular velocity and angular acceleration of the disk relative to stationary ground. 2. See ﬁgure 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 ﬁgure 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 ﬁxed 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 ﬁxed 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 ﬁxed 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 ﬁxed 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 ﬁnd 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 ﬁrst. 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 ﬁnd 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 ﬁve 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. Diﬀerentiating. 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 ﬁnd 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 eﬀective 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 ﬁrst contribution comes from rigid body rotations while the second comes from center of mass translations. So far, we have made no assumptions speciﬁc to rigid bodies. The simpliﬁcation for a single rotating rigid body is that (using, e.g., Eq. 2.5 with A ﬁxed 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 ﬁrst choose and ﬁx 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 ﬁx a speciﬁc 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 ﬁnd 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 = Deﬁne 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 conﬁguration. 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 inﬁnitely many point masses distributed in a continuum (see ﬁgure 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 ﬁve 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 ﬁxed 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 ﬂat thin object in the x-y plane (ﬁgure 3.3). For this object, taking z = 0 for every mass element, we ﬁnd (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. Deﬁnitions 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 ﬂat, 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. Diﬀerentiating, Ḣ /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? Diﬀerentiating 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 ﬁrst 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 ﬁrst 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 ﬂat 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 ﬁnd the vertical component of the moment exerted on the tyre by the axle. 3. See ﬁgure 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 ﬁgure. 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 ﬁgure 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, ﬁnd the angular momentum of the cylinder. At the same instant, ﬁnd 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 diﬀerent from translations. If a body undergoes two successive rotations, the ﬁnal conﬁguration depends on which rotation occurs ﬁrst. For example, if a person standing upright and facing north is ﬁrst 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 ﬁrst 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 diﬀerent conﬁguration. 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 ﬁnite 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 ﬁxed uniquely for any given coordinate system but will change if we change the coordinate system (we will not). 4. Positive deﬁnite matrices: If R is invertible, then RT R is symmetric and positive deﬁnite (SPD): it has three real strictly positive eigenvalues. If R is singular, then RT R is symmetric and positive semideﬁnite: 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 satisﬁes 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 ﬁxed 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 ﬁxed is a rotation about that axis. (The reader can imagine a door on hinges. The hinges ﬁx points along a certain line on the door, i.e., the door has one axis ﬁxed.) Let a rigid body be given an arbitrary displacement that keeps ﬁxed a point O on the body (see ﬁgure 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 ﬁxed 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 diﬀerent 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 diﬀerent 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 ﬁnd 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 ﬁnd 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 ﬁxed 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 ﬁxed 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, deﬁning a new unit vector ê. See ﬁgure 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 ﬁgure, 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 deﬁne 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-ﬁle 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 inﬁnitesimal rotation (ignoring second order quantities) is given by R(n, θ) = I + θS(n) , and therefore two successive inﬁnitesimal 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 inﬁnitesimal rotations commute, while ﬁnite (nonzero) rotations in general do not. 4.5.2 Angular velocity Since inﬁnitesimal rotations commute, they are vectors (unlike ﬁnite 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 deﬁne the vector ω by S(ω) = ṘRT . (4.2) Now consider the time-varying vector r , given by r = Rr, with r some ﬁxed 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 ﬁxed at the origin, v = ω × r , from undergraduate level vector mechanics (also Eq. 2.1, with A = r ). Thus, ω ≡ ω as deﬁned 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 inﬁnitesimal time interval Δt, the inﬁnitesimal rotations about the x, y and z axes are ωx Δt, ωy Δt and ωz Δt respectively; these inﬁnitesimal 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 signiﬁcance 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 unspeciﬁed 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 inﬁnitesimals, 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 inﬁnitesimal rotations about two mutually perpendicular axes n1 and n2 . Again up to second order in inﬁnitesimals, we ﬁnd 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 inﬁnitesimals are considered, then successive rotations about two mutually perpendicular axes deﬁnitely 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 unspeciﬁed 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 conﬁguration 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 ﬁnally, through ψ about the now-twice-rotated e3 . Deﬁning R1 = R(e3 , φ), R2 = R(R1 e1 , θ) and R3 = R(R2 R1 e3 , ψ), the ﬁnal 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 ﬁrst rotation φ and the last rotation ψ are made about the same axis. Therefore if the ﬁnal conﬁguration is given, say, then from that information only φ + ψ can be determined: the individual rotations φ and ψ cannot be uniquely determined. For such conﬁgurations, the choice of coordinates (φ, θ, ψ) as deﬁned 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 conﬁgurations. 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 deﬁne 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 ﬁrst-φ-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 ﬁgure 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 eﬀect 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 ﬁrst 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 coeﬃcients of ﬁxed vectors in rotating coordinate systems (a more common approach than adopted here). However, we prefer not to discuss rotating coordinate systems except, for completeness, brieﬂy 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 deﬁne 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 conﬁguration. 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 ﬁnd α = 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 ﬁxed while the coordinate axes themselves are rotated like a rigid body? Any rotation in whatever Euler angle sequence is ﬁnally 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 eﬀect on a ﬁxed vector r is that its matrix of components in the unrotated basis, r, changes to a diﬀerent 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 ﬁxed and rotate the basis about n through θ, or (2) keep the basis ﬁxed 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 ﬁxed in it. When in the reference conﬁguration, p = î − ĵ (note: p is not a unit vector). The body is ﬁrst 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-ﬁxed z-axis. Find the net rotation matrix. Also, ﬁnd 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 ﬁxed. 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 conﬁguration, was a diagonal matrix with the diagonal elements {1, 1, 2}, what is its moment of inertia matrix in the rotated conﬁguration? 3. Consider again a rigid body with one point ﬁxed 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 ﬁnd 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 veriﬁed 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 ﬁnding 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 conﬁgurations? 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 diﬀerent 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 eﬀect on the motion is known, and this eﬀect 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 speciﬁed 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 speciﬁed 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 ﬁxed 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 inﬁnitesimal displacements of the block are in fact not parallel to the inclined plane (see ﬁgure 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 ﬁgure 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 ﬁnd that for constraint-satisfying virtual displacements, N (F ai − mi r̈i ) · δri = 0 . (5.2) i=1 In the above equation, the individual coeﬃcients 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 ﬁnd each r and vice versa. For example, if a bead moves on a circular wire ﬁxed 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 speciﬁed 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 ﬁxed 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 ﬁxed 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 ﬁrst 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 inﬁnitely many points in it; and each of these points has three translational degrees of freedom. Thus, the beam has inﬁnitely 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 suﬃcient accuracy by a time dependent cubic polynomial with two free coeﬃcients, 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 deﬂection of that point. The above approximation gives rise to a 2 degree of freedom system, in which the coeﬃcients 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 speciﬁed 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, diﬀerentiating 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 ﬁnd 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 diﬀerentiating holonomic constraint equations, and cannot be used to eliminate generalized coordinates. To see this, we proceed with an example called the “skate” constraint. See ﬁgure 5.3. The ﬁgure shows a single elliptical ﬂat 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 diﬀerentiated form of some holonomic constraint, we will argue by contradiction. Accordingly, let us assume that the constraint is indeed a diﬀerentiated 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 ﬁrst 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 diﬀerentiated 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 conﬁgurations. 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 deﬁned. 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 deﬁnition 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 ﬁgure 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 deﬂection θ 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 deﬂection is Lθ and the spring force is therefore kLθ acting downwards. In a further inﬁnitesimal 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 diﬃcult in practice as the formal deﬁnition 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 ﬁxed 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 ﬁnally 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 ﬁrst 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 diﬀerently. 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 deﬁne 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 ﬁgure 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 satisﬁes 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 diﬀerent, but also constraint-obeying, motion of the block. So, to items (1) through (3) above, we add on a ﬁnal 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 ﬁx 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 satisﬁes 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 coeﬃcients 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 ﬁnd 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 ﬁrst 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 diﬀerentiated-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 diﬀerentiated 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 ﬁgure 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 diﬀerentiation will yield a fourth equation involving ẍ, ÿ and θ̈. 5.6 The calculus of variations We will now approach Lagrange’s equations from a diﬀerent direction. To begin, consider some unknown function y(x) deﬁned 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 ﬁnd 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 coeﬃcients ck , then the boundary conditions are satisﬁed 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 ﬁnd a function y that satisﬁes 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 ﬁrst 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 speciﬁed. 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 ﬁrst 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) deﬁned on the interval (a, b). If we ﬁnd 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 ﬁnd 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 ﬁnding the shortest path between two points on a plane (a straight line), ﬁnding the path in a vertical plane down which a frictionlessly sliding particle will move fastest (called the brachistochrone problem), and ﬁnding 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 ﬁnd 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 diﬀerent 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 diﬀerent 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 diﬀerential equations of motion of systems with inﬁnitely 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 stiﬀness 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 conﬁguration as the one where the coin lies ﬂat on the ground (taken to be the x-y plane). This reference conﬁguration is physically singular in that there are inﬁnitely many contact points, while the slightest change in conﬁguration gives a single unique point. This reference conﬁguration 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 conﬁguration, 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 ﬁnd 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 ﬁnd 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 diﬀerentiated 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 conﬁguration 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 ﬁgure 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 conﬁguration, n̂ is aligned with k̂ (which is the same as e3 in Chapter 4). Unless the coin is ﬂat (i.e., in the singular conﬁguration), the vector p = n̂ × k̂ is nonzero, and (in ﬁgure 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 ﬁnd 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 diﬀerentiated 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 (ﬁgure 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 α). Diﬀerentiating 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 ﬁrst time derivatives. Note that if we do not want to ﬁnd 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 diﬀerential 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 ﬁgure 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 diﬀerent periodic functions with diﬀerent 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 insigniﬁcant role in the steady state response. 3. If the forcing involves a signiﬁcant component at a frequency close to ωn , then the damping term is important but the other frequency components in the forcing play an insigniﬁcant 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 ﬁrst 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 ﬁgure 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 ﬁrst mass is x1 , and that of the second mass is x2 . We ﬁrst 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 speciﬁc 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 deﬁnite (in other problems, K can be positive semideﬁnite), 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 satisﬁes the equations of motion, is possible here because of the symmetry in the system. This motion corresponds/ to the ﬁrst normal mode. Note that the ratio of x1 to x2 is ﬁxed (at 1), and the angular frequency is ﬁxed (at k/m). The amplitude of the motion is not ﬁxed, 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 ﬁxed (at −1), and the angular frequency is ﬁxed (at 3k/m). The two normal modes of the system are sketched in ﬁgure 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 deﬁnite. 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 deﬁnite. Of course, AM is merely M in disguise. CHAPTER 7. LINEAR VIBRATIONS 74 n × n matrix K is called the stiﬀness matrix, and is usually symmetric and positive semideﬁnite. 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 deﬁnite (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 semideﬁnite, 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 ﬁrst 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 ﬁrst 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 deﬁniteness 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 stiﬀness-weighted sense? There are orthogonality properties involving the stiﬀness 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 deﬁne 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 deﬁne 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 deﬁnite. 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 speciﬁed locations. The equations of motion will be of the form M ẍ + C ẋ + Kx = BF (t), (7.6) where M , K and x are as deﬁned before; C is an n×n matrix that is often symmetric and positive semideﬁnite (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 deﬁned 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 ﬁgure 7.4) of ﬂexural 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 inﬁnitely many natural frequencies. What is the ﬁrst 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 ﬁrst 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 ﬁnd 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 ﬁrst 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 ﬁnd only the ﬁrst natural frequency. But suppose we want to estimate the ﬁrst 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 diﬀerent 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 ﬁrst 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 ﬁxed parameter values. This calculation is not given here (the reader may beneﬁt 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 diﬀerentiated 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 ﬁrst 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 (ﬁgure 8.1). We will ﬁnd 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 ﬁnd by routine integration that the height of the center of mass when the shell is in the equilibrium position (see ﬁgure, 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 ﬁgure, 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 aﬀect neither the third row nor the third column of Icm . 8.3. TORQUE-FREE RIGID BODY 81 We ﬁrst consider motion without slip. As shown in the ﬁgure (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 diﬀerentiation; 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 stiﬀness 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 ﬁnd [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 satisﬁed 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 ﬁnal 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 identiﬁcation, let us call the ﬁrst (H 2 conserved) ellipsoid the H-ellipsoid; and the second (2T conserved) the T -ellipsoid. The situation is depicted in a 2D schematic in ﬁgure 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 ﬁgure 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 eﬀects 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 ﬁxed. But the T -ellipsoid begins to shrink (while maintaining a ﬁxed 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 conﬁguration is shown labeled “T , initial” in the ﬁgure. Further energy dissipation is not possible under the assumptions of the analysis, and subsequent spin continues indeﬁnitely 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 conﬁguration 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 diﬀerence below. Taking unit vectors î, ĵ and k̂ along the x, y and z axes, respectively, and deﬁning ξ = ξ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 ﬁgure 8.3. As also shown using grey shading in the ﬁgure, 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 ﬁxed 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, ﬁnally, 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 conﬁguration 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 ﬁrst rotation is about the x-axis and is called φ; the second rotation is about the rotated body-ﬁxed z-axis and is called θ; and the third rotation is about the twice-rotated body-ﬁxed 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 conﬁguration is an equilibrium conﬁguration (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 conﬁguration). 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 ﬁnd that the absolute angular velocity of the disk is ω disk = −ψ̇0 î + O(θ0 ), where î is along the stationary or inertially ﬁxed 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 deﬁnition of (1,3,1) Euler angles that θ0 is the angle made by the body-ﬁxed x axis with the inertially ﬁxed 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 (ﬁgure 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 simpliﬁes to n̂ · rG/P × (−mg k̂) = n̂ · rG/P,ref × (m aG ) + I cm,ref · ω̇ . (8.18) Correct upto ﬁrst 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 conﬁguration of the above system. Thus, the ﬁrst 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 conﬁguration the disk is ﬂat, 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 conﬁguration 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 signiﬁcant 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 diﬃcult at two levels. First, accurate mathematical description of the friction forces is itself diﬃcult (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 coeﬃcients, μs and μk , called the coeﬃcients 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 deﬁned tangent plane exists at P (see ﬁgure 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 inﬂuenced 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-deﬂection 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 diﬀerent stiﬀnesses. The force-deﬂection 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 diﬀerent type of constitutive relation describing the force-displacement behavior of contacting surfaces. Frictional eﬀects, being localized at or near surfaces, are sensitive to environmental conditions like humidity and surface cleanliness. Day to day variations in measured friction coeﬃcients can be signiﬁcant. Friction models are approximate, and usually much less accurate than (say) the force-deﬂection 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 ﬁgure 9.2. A block of mass m moves on a horizontal, frictional plane. A spring of stiﬀness k is attached to the block. The free end of the spring is moved by an external agent at a constant speed v. The coeﬃcients v k m immovable floor μ vrel Figure 9.2: Sliding block. of friction are μs and μk < μs . In the ﬁgure, 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 ﬁrst 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 diﬃcult to tackle numerically. We ﬁrst 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 unaﬀected 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 ﬁnd 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 ﬁnd 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 ﬁgure 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 coeﬃcient. 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 coeﬃcient. 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 ﬁner step sizes in time). The numerical solution obtained is shown in ﬁgure 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 ﬁrst starts and ﬁrst 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 ﬁrst 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 (ﬁgure 9.2) to actually be two masses of m/2 each, connected by a spring of high stiﬀness 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 ﬁrst imagine K → ∞ and then imagine → 0, we get the behavior described above. However, if we ﬁrst consider → 0 and then consider K → ∞, then we get a diﬀerent behavior where sliding starts earlier. If we imagine n masses of m/n each, connected in a chain with springs each of high stiﬀness K, then sliding starts earlier still. In this way, the limiting case of a rigid block with Coulomb friction coeﬃcients μs > μk gives results that depend on the relative rates at which K and approach inﬁnity and zero respectively (a distinguished limit). In addition to these theoretical considerations, it is also known from careful experiments that the coeﬃcient 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 coeﬃcient which has magnitude (and importance) comparable to the change in the steady value. In other words, modeling the friction coeﬃcient 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 eﬀects can be approximately modeled and understood using a single friction coeﬃcient (μ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 stiﬀness for a typical spring). For many problems, nothing is really lost on ignoring the diﬀerence. 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 ﬁgure 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 ﬁxed (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 satisﬁes 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 ﬁgure 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 ﬁgure 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 ﬁnally 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 ﬁgure 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, modiﬁed 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 ﬁgure 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 ﬁgure 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 diﬀerentiable function g(x, y) changes inﬁnitesimally due to inﬁnitesimal 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 ﬁnd 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 ﬁnd the amplitude growth rate. Having identiﬁed 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 oﬀ 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 ﬁxed 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 ﬁgure 9.8. It is seen that the velocity has ﬂuctuations 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 ﬁnd 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 eﬀect 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 ﬁgure 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 eﬀect 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 diﬃcult to control the ﬂow of sugar; however, small but rapid side to side shaking of the can causes the sugar to ﬂow out more smoothly, like a ﬂuid on average. Though the average ﬂow 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 ﬁgure 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 ﬁgure 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 ﬁner 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 ﬁnd ω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 ﬁgure 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 ﬁgure 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 coeﬃcient 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 ﬁgure 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 ﬁner eﬀects. For example, if a sphere rolls over a long distance on a ﬂat 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 ﬂat surface actually deform slightly, and possibly inelastically; that contact actually occurs over a region of small but nonzero size; and that diﬀerent 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 (ﬁgure 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 ﬂat objects on ﬂat tables turns out to be mathematically quite complicated. Disk: Let us begin with a disk. On encountering serious analytical diﬃculties, we will retreat to a rod. A ﬂat uniform disk of mass m and radius R slides and rotates on a ﬂat table; the coeﬃcient of friction is μ. See ﬁgure 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 (ﬁxed 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 diﬀerential 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 ﬂat ruler) sliding on a table. See ﬁgure 9.14. A uniform rod of length L slides and rotates on a table with friction coeﬃcient μ. 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.