A Quick Tutorial On Multibody Dynamics

A Quick Tutorial on Multibody DynamicsC. Karen LiuSumit JainSchool of Interactive ComputingGeorgia Institute of Technologyi

Contents1 Introduction22 Lagrangian Dynamics33 Review: Newton-Euler equations64 Rigid Body Dynamics: Lagrange’s equations85 Articulated Rigid Body Dynamics5.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .5.2 Cartesian and generalized velocities . . . . . . . . . . . . . . . . . . . . . . .5.3 Equations of motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .131314156 Conversion between Cartesian and Generalized Coordinates6.1 Velocity conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .6.2 Force conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .1717197 Recursive Inverse Dynamics7.1 Dynamics in the local frame . . . . . . . . . . . . . . . . . . . . . . . . . . .7.2 Pass 1: Compute velocity and acceleration . . . . . . . . . . . . . . . . . . .7.3 Pass 2: Compute force and torque . . . . . . . . . . . . . . . . . . . . . . . .212122241

1IntroductionIf you have not read the excellent SIGGRAPH course notes on physics-based animationby Witkin and Baraff, you can stop reading further right now. Go look for those notes athttp://www.cs.cmu.edu/ baraff/sigcourse/ and come back when you fully understandeverything in those notes.If you are still reading this document, you probably fit the following profile. You are a computer scientist with no mechanical engineering background and minimal training in physicsin high school but you are seriously interested in physics-based character animation. Youhave read Witkin and Baraff’s SIGGRAPH course notes a few times but don’t know whereto go from simulating rigid bodies to human figures. You have played with some commercialphysics engines like ODE (Open Dynamic Engine), PhysX, Havok, or Bullet, but you wishto simulate human behaviors more interesting than ragdoll effects.Physics-based character animation consists of two parts: simulation and control. This document focuses on the simulation part. It’s quite likely that you do not need to understandhow the underlying simulation works if your control algorithm is simple enough. However,complex human behaviors often require sophisticated controllers that exploit the dynamicsof a multibody system. A good understanding of multibody dynamics is paramount fordesigning effective controllers.There are many ways to learn multibody dynamics. Reading a textbook on this topic ortaking a course from the mechanical engineering department will both do the job. However,if you only want to learn the minimal set of multibody dynamics necessary to jump startyour research in physics-based character animation, this document might be what you arelooking for. In particular, this document attempts to answer the following questions. I know how to derive the equations of motion for one rigid body and I have seen peopleuse the following equations for articulated rigid bodies, but I don’t know how they arederived.M (q)q̈ C(q, q̇) Q I have seen the Euler-Lagrange equation in the following form before, but I don’t knowhow it is related to the equations of motion above. d Ti Ti Q 0dt q̇ q I use generalized coordinates to compute the control forces, how do I convert themto Cartesian forces such that I can use simulators like ODE, PhysX, or Bullet whichrepresent rigid bodies in the maximal coordinates? I heard that inverse dynamics can be computed efficiently using a recursive formulation.How does that work?2

2Lagrangian DynamicsArticulated human motions can be described by a set of dynamic equations of motion ofmultibody systems. Since the direct application of Newton’s second law becomes difficultwhen a complex articulated rigid body system is considered, we use Lagrange’s equationsderived from D’Alembert’s principle to describe the dynamics of motion. To simplify themath, let’s temporarily imagine that the entire human skeleton consists of a collection ofparticles {r1 , r2 , . . . , rnp }. Each particle, ri , is defined by Cartesian coordinates that describethe translation with respective to the world coordinates. We can represent ri by a set ofgeneralized coordinates that indicate the joint configuration of the human skeleton:ri ri (q1 , q2 , . . . , qnj , t)(1)where t is the time and qj is a joint degree of freedom (DOF) in the skeleton. Each qj is afunction of time but we assume that ri is not an explicit function of time.The virtual displacement δri refers to an infinitesimal change in the system coordinates suchthat the constraint remains satisfied. In the context of human skeleton, the system coordinates are the generalized coordinates qj and the constraint manifold lies in the Cartesianspace. The virtual displacement δri is a tangent vector to the constraint manifold at a fixedtime, written asX riδri δqj(2) qjjWe can now write the virtual work done by a force fi acting on particle ri asfi · δri fi ·X rij qjδqj XQij δqj Qi · δq(3)j T rifi is defined as the component of the generalized force associated withwhere Qij qjcoordinate qj . In vector form, Qi is the generalized force corresponding to the Cartesianforce fi with the relation Qi JiT fi , where Ji is the Jacobian matrix with the j th column ridefined as q.jFrom D’Alembert’s principle, we know that the sum of the differences between the forcesacting on a system and the inertial force of the system along any virtual displacementconsistent with the constraints of the system, is zero. Therefore, the virtual work at ri canbe written asX riδWi fi · δri µi r̈i · δri µi r̈i ·δqj(4) qjjwhere µi is the infinitesimal mass associated with ri . The component of inertial force associated with qj can be written as d rid ri ri µi ṙi ·(5) µi ṙi ·µi r̈i · qjdt qjdt qj3

Now let us consider the velocity of ri in terms of joint velocity q̇jṙi X rij qj(6)q̇jfrom which we derive the following two identities: ṙi q̇j ṙi qj ri qjX 2 rid ri q̇k qj qkdt qjk Using these two identities, we rewrite Equation (5) as 1 T1 Td ri µi ṙ ṙiµi ṙ ṙi µi r̈i · qjdt q̇j 2 i qj 2 i(7)(8)(9)We can denote the kinetic energy of ri as1Ti µṙTi ṙi ,2(10)and rewrite Equation (9) as ridµi r̈i · qjdt Ti q̇j Ti qj(11)Combining the definition of generalized force (Equation (3)), D’Alembert’s principle (Equation (4)), and the generalized inertial force (Equation (11)), we arrive at the following equation: d Ti Ti δqj Qij δqj(12)dt q̇j qjIf the set of generalized coordinates qj is linearly independent, Equation (12) leads to Lagrangian equation: Tid Ti Qij 0(13)dt q̇j qjEquations of Motion in Vector Form. Equation (13) is the equation of motion for onegeneralized coordinate in a multibody system. We can combine nj scalar equations into thefamiliar vector formM (q)q̈ C(q, q̇) Q(14)where M (q) is the mass matrix, C(q, q̇) is the Coriolis and centrifugal term of the equationof motion, and Q is the vector of generalized forces for all the degrees of freedom (DOFs) inthe system. M only depends on q and C depends quadratically on q̇.4

How do we derive M and C from Equation (13)? Let us go back to the velocity of oneparticle ri :X riṙi q̇j Ji (q)q̇(15) qjjwhere Ji denotes the Jacobian of ri . By summing up all the particles in the system, thekinetic energy of the system can then be expressed asT XiTi X1i2µṙTi ṙi X11µJiT Ji )q̇ q̇T M (q)q̇µ(Ji q̇)T (Ji q̇) q̇T (222iX1iwhere we define the mass matrix, M (q) mass matrix in Equation (14).Pi(16)µJiT Ji , and will shortly show it is indeed theFrom Equation (16), we can derive the derivative terms to construct the Lagrange’s equation(Equation (13)): T1d T M q̈ Ṁ q̇ q̇Tdt q̇ q2 M q Tq̇ M q̈ C(q, q̇)(17)Comparing Equation (17) to Equation (14), we confirm that the mass matrix is identical inboth equations. C is the Coriolis and centrifugal term in Equation (14) and is defined as TC Ṁ q̇ 21 Mq̇q̇. q, which implies that theNote. In the second term of C, we introduce tensor notation M q M Mthj element of the tensor q is the matrix qj . Note that, in general, the quantity with Mq̇ is not equal qP)k Mq̇ or k (Mq k , qj qjthnotationvectorIn contrast, the jto Ṁ . This is because, the j th column of the matrixwhere the notation (A)j denotes the jP (M )column of the matrix Ṁ is k qk j q k .th Mq̇ qis thecolumn of the matrix A.Once we know how to compute the mass matrix, Coriolis and centrifugal terms, and generalized forces, we can compute the acceleration in generalized coordinates, q̈, for forwarddynamics. Conversely, if we are given q̈ from a motion sequence, we can use these equationsof motion to derive generalized forces for inverse dynamics.The above formulation is convenient for a system consisting of finite number of mass points.However, for a dynamic system that consists of rigid bodies, there are infinitely many pointscontained in each rigid body making the above formulation intractable. In the following twosections, we view a rigid body as a continuum and derive compact equations of motions inboth Cartesian coordinates and generalized coordinates.5

3Review: Newton-Euler equationsThis section reviews Newton-Euler equations for rigid body dynamics. The derivation of massmatrix M (q) and Coriolis and centrifugal term C(q, q̇) for a rigid body will be presented inthe next section. If you are familiar with Newton-Euler equations, you can skip this sectionand continue to the next. However, many math notations used in Witkin and Baraff’s coursenotes are also reviewed in this section, such as linear momentum, angular momentum, skewsymmetric matrix and its properties.math notations and definitions used in Witkin and Baraff’s course notes, such as linearmomentum, angular momentum, and skew-symmetric matrices and their properties, youcan safely skip this section.To derive Newton-Euler equations, we begin with the momenta of the rigid body whosemass, position of the center of mass (COM), orientation, linear velocity of the COM, andangular velocity are m, x, R, v, and ω respectively (these definitions are the same as arefound in Witkin and Baraff’s course notes). The linear momentum P is computed as:XXXP Pi µṙi µ(v ω r′i )iii mv(18)Pwhere r′i ri x. Because i µr′i 0 (property of the COM), the second term vanishes.The angular momentum L about the COM is computed as:XXL Li r′i Pii Xiµr′i (v ω r′i )i 0 Xµ[r′i ][ω]r′i iXi! µ[r′i ][r′i ] ω(19)The notation [a]b denotes the cross product a b with [a] being the skew-symmetric matrixcorresponding to the vector a: 0 a3 a20 a1 (20)[a] a3 a2 a10Therefore the following identities hold: [a]b [b]a and [a]T [a].Now Precall the inertia tensor about the COM defined in Witkin and Baraff’s course notes:′′ ′TIc Pi µ((r′Ti ri )I3 ri ri ), where I3 is the 3 3 identity matrix. We can easily show that′′Ic i µ[ri ][ri ] by verifying the identity [a][a] (aT a)I3 aaT . As a result, we writethe angular momentum of a rigid body as:L Ic ω6(21)

where the inertia tensor can be written as Ic RI0 RT . R is the rotation matrix corresponding to the orientation of the body and I0 is the constant inertia tensor defined at zerorotation. From Witkin and Baraff’s course notes, we also learned that the angular velocityin the skew-symmetric form is related to the rotation matrix R as [ω] ṘRT .With these definitions, we can derive the equations of motion for a rigid body. The equationscorresponding to the linear force can be evaluated as:f ṗ mv̇(22)The equations corresponding to the torque can be evaluated as:τ L̇ (Ic ω) Ic ω̇ (RI0 RT )ω Ic ω̇ ṘI0 RT ω RI0 ṘT ωIc ω̇ ṘRT Ic ω Ic (ṘRT )T ωIc ω̇ [ω]Ic ω Ic [ω]ω (Using the identity [ω]T [ω])Ic ω̇ ω Ic ω(23)Combining Equation (22) and Equation (23), we arrive at the Newton-Euler equations: f0v̇mI3 0(24) τω Ic ωω̇0 Ic7