Multibody Dynamics B
wb1413
 (Spring 2009)

Delft University of Technology 

Description: After the first courses in dynamics a student can often deal well with the dynamics of one rigid body. In this course we will cover a systematic approach to the generation and solution of equations of motion for mechanical systems consisting of multiple interconnected rigid bodies, the so-called multibody systems. This course differs from ``advanced dynamics'', which mostly covers theoretical results about classes of idealized systems (e.g. Hamiltonian systems), in that the goal here is to find the motions of relatively realistic models of systems (including, for example, motors, dissipation and contact constraints).

The main approach is the reduction of the constraint Newton-Euler equations using the principle of virtual work and the principle of D'Alembert. However the relation of this approach to Lagrange's equations, Kane's equations, and to the differential algebraic equations approach will also be covered. The course will start with 2D examples and then move on to 3D systems. Rotations will be represented with matrices, various Euler angles, and Euler parameters. Various computer solution techniques for the ODEs or DAEs will also be covered. 

Goal: By the end of the course students will be competent at finding the motions of linked rigid body systems in two and three dimensions including systems with various kinematic constraints (sliding, hinges and rolling, closed kinematic chains). Collisional interactions will be considered in a unified manner for all the different ways of formulating the equations of motion.

Homework: There will be weekly homework assignments and a final project. The homework is due a week after hand out and will be graded. Hand in your homework at the start of class at the front. The homework is strictly individual but in doing the homework I encourage you to work together. To get credit, on every homework assignment please do the following things:

The graded work can be picked up at the end of class or at Guido's office.

Grading: Total course grading is 70 % homework and 30 % final project. Homework is the average of the weekly homework assignments. Final project is the grading of the written report on the final project (strictly individual!).

News
- July 6: Final Grades can be found at: wb1413Spring2009.htm Please check this and contact me if things are not ok.
- May 28: Attendance to the Mechanics Colloquium on Thursday June 4, 12:45-13:30, room D, will give you one bonus point on the final project.
- Apr 16: Answer to puzzler6.m.
- Mar 19: Attendance to the Mechanics Colloquium on Wednesday April 8, 12:45-13:30 room B will give you one bonus point on the final project.
-Mar 12: Question for Puzzler5.txt
-Mar 12: Answer to last weeks puzzler: puzzler4.txt
-Mar 5: Answer to last weeks puzzler: puzzler3.txt , and did you get the Bonus?
-Feb 18: Answer to last weeks puzzler: Stanley from the Stanford Racing Team, see DARPA 2005 Grand Challenge.
-Feb 2: NO LECTURE on Feb 5, I'm sick wit the flu.

Hand-Outs
- The course contents.
- Preliminary course notes, first chapters MultibodyDynamicsB.pdf
- Reprint from Lanczos1949.
- My lecture notes on numerical integration NumericalIntegrationCh7.pdf
- My lecture notes on the coordinate projection method:wb1413_diktaatCh8.pdf.
- My notes on 3D rotation and a conference paper on this topic DETC2006Paper99307.pdf.

Homework
- Assignment 1, due Thu Feb 19, 2009, 15:45 h.
- Assignment 2, due Thu Feb 26, 2009, 15:45 h.
- Assignment 3, due Thu Mar 05, 2009, 15:45 h.
- Assignment 4, due Thu Mar 12, 2009, 15:45 h.
- Assignment 5, due Thu Mar 19, 2009, 15:45 h.
- Assignment 6, due Thu Apr 09, 2009, 15:45 h.
- Assignment 7, due Thu Apr 16, 2009, 15:45 h.
- Assignment 8, due Thu Apr 23, 2009, 15:45 h.
- Assignment 9, due Thu May 14, 2009, 15:45 h.
- Assignment 10, due Thu May 28, 2009, 15:45 h.

Final Project
Pick one of the final projects and hand-in your report on paper Mon June 15, 2009 the latest to me, Arend Schwab, or my paper mailbox.
- Final Project 1, Simulating the motion of an ejection seat reverse bungee.
- Final Project 2, Dynamic analysis of the motion of a bicycle, and some notes on the bicycle project.

Office Hours
Arend L. Schwab:     Monday,    14-16 h., room 4B-1-24, phone: 015 278 2701, a.l.schwab@tudelft.nl.
Guido Delhaes:         Wednesday, 9-12 h., room 5A-1-55, phone: 015 278 6509, g.m.j.delhaes@tudelft.nl

Log

Week 1

Thursday, Feb 5, 2009 15:45-17:30 u, room B.

NO LECTURE, I'm sick wit the flu.

Week 2

Thursday, Feb 12, 2009 15:45-17:30 u, room B.

Short introduction to Multibody Dynamics. Showed a video with an example of real dynamical systems.
- A stable monocycle. (look under media->movie working models 1956-1965)
Talked about Newton and Euler, and the rigid body concept. Posed the equations of motion for a rigid body in 2D, the so-called Newton-Euler equations. Derived the complete Differential and Algebraic Equations for the double pendulum motion. Handed out assignment 1, this homework is due in one week, Feb 19, 2009, 15:45 h..

Week 3

Thursday, Feb 19, 2009, 15:45-17:30 u, room B.

Solved the accelerations together with the joint forces of the double pendulum for the upright configuration at zero speed. Discussed in dept the results qualitatively. Did two Euler numerical integration steps to find the motion in time. Showed this motion  and discussed the matlab function to calculate the accelerations and the constraint forces daehwa1.m and the script to integrate the eqns of motion hwa1.m . Talked about the discrepancy between the calculated result (things fly apart!) and real life. Started with the systematic approach for deriving the equations of motion for a system of rigid bodies with constraints. Applied the principle of virtual power and included the D'Alembert forces to come to the equations of motion instead of static equilibrium. Talked about virtual velocities which fulfill the constraints. Introduced the Langrange multipliers lambda to include the constraints on the virtual velocities. Derived the equations of motion f-M*xdd=D'*lambda. Added the constraints on the accelerations as derived from the constraints on the coordinates diff twice with respect to time i.e. D*xdd+D2*xd*xd=0. Came up with a full set, n+m, of linear equations to solve for the n unknown accelerations of cm's of the bodies together with the m langrange multiliers lambda. Most of you looked confused, more on this next week. Handed out assignment 2, this is due in one week, Feb 26.

Week 4

Thursday, Feb 26, 2009, 15:45-17:30 u, room B.

Recapitulation of last week's lecture. Showed that the eqn's of motions for the constrained system derived by means of virtual power and Lagrange multipliers indeed express external force equilibrium by setting xdd to zero and looking at f=D'*lambda. Interpreted lambda as the constraint forces by drawing for the double pendulum in the upright position the four equilibrium force vectors. Handed out a  reprint from Lanczos1949 together with assignment 3 which is on this weeks lecture and the reprint. The homework is due in one week, Mar 05.

Week 5

Thursday, Mar 05, 2009, 15:45-17:30 u, room B.

Added Passive and Active elements to system by adding the virtual power term of these elements to the virtual power expression. Derived the equations of motion with the inclusion of these elements. As an example I have added a spring to body one of the double pendulum and considered the effect on the equations of motion. Impacts. Started with the example of two point masses impacting. Discussed the limit case where the speeds change step-wise and introduced the impulse as the limit case of a finite force-time integral. Discussed Newton's impact law from the source, The Principia (1686) and talked about the coefficient of restitution e. Introduced the contact conditions and formulated Newton's impact law in terms of the relative contact velocities and solved for the velocities after impact together with the contact impulse. Generalized the impacts to impacts in multibody systems. Introduced the contact conditions and formulated Newton's impact law in terms of the relative contact velocities. Introduced the applied impulses and the constraint and contact impulses and finally derived the impact equations for a constrained multibody system. Handed out homework assignment 4 which is due in one week, March 12.

Week 6

Thursday, March 12, 2009, 15:45-17:30 u, room B.

A new topic: Lagrange equations of motions.
In the case of many rigid bodies n>>1 and many constraints m>>1 we have few degrees of freedom dof=n- m. To reduce the number of equations n+m and to eliminate the m constraint forces we want to write our eqn's of motion in terms of these degrees of freedom. For the degrees of freedom we will use Generalized Coordinates. I gave some examples. Joseph-Louis Lagrange [Turin 1736-Paris 1813] was the first to derive the eqn's of motion in terms of these Generalized Coordinates and he wrote it all down in his monumental book Mecanique analitique (1788) by which he became the founder of Multibody System Dynamics. There is a very good English translation of this work from 1997 with a wonderful introduction, look under Lagrange. Lagrange's own introduction is also very noteworthy, read the part about On ne trouvera point de Figures dans cet Ouvrage. Introduced the concept of Energy by integration of Power with respect to time. Used Newton to come to the definition of Kinetic Energy T and Work done by the applied forces W. Introduced the potential energy V as work done by conservative forces, dV/dx=-f. Showed that for these forces we have Energy conservation, T+V=constant. Showed how to derive Newton from T+V=constant (1). Introduced the generalized coordinates q_j, j=1..dof and described the coordinates of the cm of the rigid bodies x_i, i=1..n as functions of these Generalized Coordinates. Showed that the velocities of the cm of the bodies xdot_i are a function of both q_j and qdot_j! Introduced the Generalized forces Q_j by means of power and showed how the applied forces transform to the Generalized forces (2). Combined (1) and (2) and derived the equations of motion in terms of the generalized coordinates, the so-called Langrange Equations (of motion!). As an example I derived the eqn's of motion of a simple model of a container crane in terms of generalized coordinates by means of the Lagrange Equations.
 I showed how to add passive or active elements to the Lagrange eqn's of motion. Talked about some different types of extra elements like springs, dampers, motors and dry friction. Showed how to add a constraint to the system and the extensions to the Lagrange eqn's of motion. I showed the impact eqn's in terms of independent generalized coordinates from the constraint Lagrange eqn's of motion. Please do assignment 5 which is due March 19.

Week 7

Thursday, March 19, 2009, 15:45-17:30 u, room B.

I derived the eqn's of motion for the constraint system by using the principle of virtual power, D'Alembert forces, and transformation of the coordinates of the cm's of the bodies x_i in terms of the generalized coordinates q_j as in x_i=F_i(q_j). I introduced the applied generalized forces Q_j which are dual to qd_j by the Power=Q_j qd_j, as an extra source of Power to the system. I extended the virtual power with this term and derived the equations of motions in terms of the independant coordinates q_j. They read:

F_i,j M_ik F_k,l qdd_l = Q_j + F_i,j(f_i-M_ik F_k,lm qd_l qd_m).

As an example I derived the equations of motion for the simple 2D model of a container overhead crane from last time. Finally, you can study the Matlab files, odehwa6.m and hwa6sym.m from the 2005 homework assignments, which demonstrate the elegance of the method. Handed out assignment 6 which is due Thursday April 9.

We continue after the break on Thursday April 9.

Week 8

Thursday, April 9, 2009, 15:45-17:30 u, room C.

Up until now we solved for the accelerations. What we really want to know is the positions and velocities of the system as a function of time. You could try to solve the equations of motion, being a set of second order differential equation, algebraically but alas. Most systems are complicated and can not be solved in terms of elementary functions. We therefore turn to numerical integration. A natural (but as will turn out naive) way would be the method as demonstrated at the second lecture, the Euler method. I showed a more refined method, Heun's method. Discussed the accuracy of the methods and introduced the local truncation error e, for Euler e= O(h^2) and for Heun e=O(h^3). I introduced the global truncation error E, being the cumulated error over a fixed timespan t=0..T, and consequently E=(T/h)*e, so for Euler: E=O(h) and for Heun E= O(h^2). In practice the global error E can be estimated from two integrations from t=0..T. First pick a step size h then integrate the diff eqn's with the result at t=T: y_h=Y+C*h^p, where Y is the true answer, C is a constant and p the order of the method. Now half the step size and integrate again over the timespan 0..T this result is: y_h/2=Y+C*(h/2)^p. From these two equations we can solve the unknown Y and C and give a global error estimate on y_h/2: Y = y_h/2 +/- 1/(2^p-1)*(y_h/2- y_h).
Next I addressed the problem of stability of the method, or in other words, how do previous made errors propagate? To investigate the stability of the method we will use a test equation of the form ydot = lambda * y. The solution to this diff eqn is the exponential function, y = exp(lambda*t) and with complex lambda = a + b*i this shows an oscillatory growing or decaying behavior in time y=exp(a*t)*[cos(b*t)+i*sin(b*t)]. As an example of the relation of the test equation with a mechanical system I derived the equations of motion of a single mass-spring- damper system, rewritten them in first order form, substituted the exponential function for a solution and derived the eigenvalue problem. I wrote down the characteristic equation and showed the similarity with the original second order diff eqn of motion. Showed the roots lambda_1,2 and discussed the behavior of the system from these roots. With the test eqn I derived, for the Euler method, the amplification factor C(h*lambda), which is the factor by which truncation errors get propagated. The stability condition |C(h*lambda)|<1 for complex lambda is a unit circle in the complex plane centered at (-1,0). Conclusion: the Euler method is unstable for pure oscillatory systems, and the maximal step size for pure damped systems is h<2/|lambda|. For Heun's method I showed the stability region, an ellipse which semi axis 1 and sqrt(3) centered at (-1,0). From a stability point of view Heun's method behaves the same as Euler's method. Next I showed the famous Runge- Kutta 4-stage method (RK4) from 1901. The one step or local truncation error is e=O(h^5) and the stability region encloses part of the imaginary axis, from 0 to 2.8, whereas the real axis is included from -3 to 0. Conclusion: RK4 is stable for pure oscillatory systems and the restriction on the step size for stable integration is h<2.5/|lambda|.
I revisited global error estimates. I included the local round off error e=O(1). The global round off error is the sum or E=(T/h)*e= O(1/h). In refining our step size we can hit this barrier. The total global error estimate is now E = C1*h^p + c2*1/h. In practice we estimate the global error by taking a number of successive integrations from t_start to t_end with step sizes h_n=T/(2^n) and calculate the difference between two successive solutions at t_end, D_n=|y_n-y_n-1|. The estimated error in the method truncation error region is approximate E_n=1/(2^p-1)*D_n whereas in the round off region we have E_n = D_n. We know that this latter error is the increasing thing at very small h so just plot log(1/(2^p-1)*D_n) versus log(h) and watch were the error goes up. This is the minimal step size h_min with maximal accuracy. Check the slopes of the graph, they should be p right from h_min and -1 left from h_min. 
Handed out assignment 7 which is due Thursday April 19. Note the reference list in the assignment. This list is my top 10 books on numerical integration as applied to dynamics of multibody systems. You could also read the lecture notes of this TUDelft course on Numerical Mathematics Wi3097Wb. Finally a first draft of my lecture notes on numerical integration NumericalIntegrationCh7.pdf.

Week 9

Thursday, April 16, 2009, 15:45-17:30 u, room C.

Finding the coordinates of center of mass for all bodies as a function of the degrees of freedom, x_i=F_i(q_j), is not easy for closed-loop systems. Even for a four-bar mechanism this is not trivial, and to be honest I do not know where to start. I showed a copy from a 1941 paper entitled "Part I - Analysis of Single 4-Bar Linkage" by Guy J. Talbourdet which appeared in Machine Design (May 1941). I will come back to this work. Contrary to this ingenious math work our strategy will be: divide and conquer. By making appropriate cuts in the system we can come up with an open loop system at the cost of more generalized coordinates q_j. For this open loop system its easy to write down x_i=F_i(q_j) but note that the bigger set of q_j is now dependant. This dependency is expressed by the joints we had to cut to open the loops. So we add these cut joints as constraint equations. With these constraints we can write down the DAE and solve for the accelerations of the q_j together with the Lagrange multipliers. Next we do a numerical integration step on the q_j. Now in general these new approximate values for q_j and qdot_j do not comply with the constraint equations. We have to change the generalized coordinates q_j such that the constraints are fulfilled but since we have less constraints then generalized coordinates, we have this Cole Porter problem, or in other words "Anything Goes". If we say we want the smallest change in the gen. coordinates q_j we have to come up with a way to measure, for instance the sum of the squares of the changes in the gen. coordinates. Now we want to minimize this distance where the new coordinates q_j must fulfill the constraints. This is a nonlinear constraint least square problem. We solve this by a Gauss-Newton method. First linearize the constraints e(q) around the coordinates q and try to solve for this increment dq in the coordinates. We now have a linear least square problem:
||dq||=min
e(q)+D(q)*dq=0
where D(q) is the Jacobian matrix of the constraints according to D(q)=de(q)/dq. This linear constraint least square problem can be solved the introduction of the Lagrange multipliers mu as many as we have constraints, resulting in the linear set of equations:
|I D^T||dq|=| 0|
|D 0 ||mu| |-e|
Note the resemblance with our DAE of motion, these equations have the same structure. You can find out more about this when you read on Gauss's principle of the least constraint, for instance in Lanczos. Solving for dq gives us:
dq=-D^T*(D*D^T)^-1*e
The matrix:
D^+=D^T*(D*D^T)^-1
is called the Moore-Penrose pseudo-inverse and gives us the least square solution of the underdetermined linear system of equations with full rank matrix D, according to D*dq=-e.
I made a simple example in class with D=[-1 1 0;0 -1 1] and x=[0.9;1;1] resulting in dx=[0.0666;0-0.0333;-0.0333] and consequently x=[0.9666;0.9666;0.9666]. Note the minimal change in x, f.i. dx=[0.1;0;0] would also do but is larger. In Matlab there is a pseudo inverse command pinv. Note that this uses singular value decomposition to solve the more general problem where D can have dependent rows.
Finally we have to solve for the velocities of the q_j's which fulfill the constraints. Again we use the least square solution method but since the velocity constraints are linear in the velocities this linear least square problem can be solved in one step:
qdot=qdot-D^T*(D*D^T)^-1*edot
By the way, getting rid of the constraint errors is, strangely enough, called Constraint Stabilization. Finally a first draft of my lecture notes on the Coordinate Projection Method wb1413_diktaatCh8.pdf.
Handed out assignment 8 which is due Thursday April 23.

Week 10

Thursday, April 23, 2009, 15:45-17:30 u, room C.

Today I presented the Newton-Euler equations of motion for a 3D rigid body. Talked about the changing inertia matrix Jc in the space fixed coordinate system C-xyz. The inertia matrix Jc' is constant in the body fixed coordinate system C-x'y'z'. Talked about the eigenvalues of Jc' and the principal axis of inertia for a rigid body. Gave the Newton-Euler equations of motion for a 3D rigid body expressed in the body fixed coordinate system C-x'y'z'. Discussed stable and unstable torque-free rotation of a 3D rigid body and demonstrated this phenomena by flipping the blackboard wiper. The same problem arises in the rotary motion of a satellite, a tennis racket, and a remote control (Do try this at home). No homework this work, please finish up late work!

May break, next lecture in one week.

Week 11

Thursday, May 7, 2009, 15:45-17:30 u, room C.

In the displacement of a rigid body B in space we can discriminated between translation and rotation. Translation is easy, just x, y, and z of for instance the centre of mass C. Therefore today I focus on the rotation part. I showed how, after a finite rotation, you can express the position of an arbitrary point p of body B in terms of either the space fixed coordinate system xyz or the body fixed coordinate system x'y'z'. The algebraic components p_i' of the vector p in the body fixed coord. sys. are constant since B is a rigid body. I derived the rotation matrix R_ij as being the transformation of p_j' into p_i as in p_i = R_ij p_j' and showed that this matrix is made out of the dot products of the 6 base vector e_i and e_j' as in R_ij = e_i * e_j'. I showed that the inverse of R is just the transpose, R^-1=R^T, or R is an orthogonal matrix as expressed by R R^T = I. Now R has 9 components and the orthogonality conditions are 6 independent ones (the upper triangle of R R^T = I holds the same equations as the lower triangle). We conclude that we can parameterize the 3D rotation matrix by 3 independent parameters. As an example I derived R and R^-1 in 2D were R only depends on 1 parameter, the angle phi about the z-axis.
One way to parameterize the rotation matrix is by means of the so-called Euler angles: phi, theta, and psi, also known as zxz or 313 angles. These angles are also known as the precession angle, the nutation angle, and the spin angle. The recipe for Euler angles is: Start with the body fixed coordinate system aligned to the space fixed coordinate system. First rotate an angle phi about the z-axis, second rotate an angle theta about the rotated x-axis, third rotate and angle psi about he rotated z-axis. You can materialize the Euler angles by three hinges in series. Draw the hinges as "cans" in series and do not be fussy about the location of the cans which in principal are the same but in practice fan out. First draw a hinge from the fixed space with the rotation axis along the z-axis, second draw at the end of the first hinge a second hinge with the rotation axis along the x-axis, and finally draw at the end of the second hinge a third hinge with the rotation axis along the z-axis. The body fixed coordinate system x'y'z' is located at the end of this third hinge. By considering the intermediate coordinate systems between the three hinges you can derive the rotation matrix R as a composition of three single rotation matrices about principal axis as in x = R x', with R = R(phi) R(theta) R(psi). Note that these Euler angles are no vectors since clearly the matrix multiplication is order dependent, or in other words the individual Euler angles is non-commutative. By the way, there are only 12 meaningfully ways of putting three orthogonal hinges in series, so apart from the zxz Euler angles there are these 11 other ways to parameterize R. Next I derived the expressions for the angular velocities omega in terms of the Euler angles and the time derivatives of these angles by means of differentiation of the identities R^T*R=I. I showed the formal way to derive omega from Rdot*R^T but followed the intuitive way, by adding up the transformed angular velocities of the hinges or "cans". This results in omega, w, in terms of the Euler angles, e, and it's time derivatives edot as in w=A(e)*edot. The rotational state of the system is defined by the Euler angles e and the angular velocities w. The state equitation are the time derivatives of these, so we need edot in terms of e and w. It turns out that the inverse of A is singular at the Euler angle theta=0+k*pi. Which, looking at the cans in this configuration, is obvious since the two rotation axis phi and psi are aligned! One way to circumvent this problem is by redefining theta as theta=theta+pi/2, now the phi and the psi axis are not anymore aligned. A major disadvantage of this method is that the problem now is non-smooth in the state variables and we have to restart our numerical integration procedure. Handed out homework assignment 9 which is due  Thursday May 14.
Finally a pdf of my notes on 3D rotation and a conference paper on this topic DETC2006Paper99307.pdf.

Week 12

Thursday, May 14, 2009, 15:45-17:30 u, room C.

From Euler's theorem on Rotation, "Any rotation in space can be represented as a rotation about a fixed axis at a given angle", I derived the rotation matrix R in terms of the Euler parameters lambda_0=cos(u/2) and lambda=sin(u/2)*h, where h is the axis of rotation with unit length and u the angle of rotation. It turns out that R is bilinear in lambda which is computational advantages. Think of the derivatives dR/dlambda=linear in lambda and d2R/dlambda^2=constant! In the previous lectures we showed that R can be represented by 3 independent parameters. Clearly the 4 Euler parameters are dependent, namely the norm must be one: |lambda|=1. They live on this 4 dimensional unit sphere. In order to derive the expressions for the angular velocities in terms of the Euler parameters and its time derivatives we start again from Rdot*R^T=tilde(omega). After expasnion we get the expressions 2*L*lamdadot=omega, with matrix L linear in lambda. From this we like to solve lambdadot in terms of lambda en omega. This can be done by adding the constraint on the lambdas in a differentiated form, as in 2*lambdadot^T*lambda=0. The solution now becomes lambdadot=1/2*L^T*[0; omega], note that this is without any singularities. This is the main advantage of using Euler parameters. We also could have used quaternion algebra to derive the various expressions. Euler parameters are unit quaternions. In his pursued to multiply triplets Hamilton "invented" the quaternion algebra. You can read all about this algebra in one of our recent conference papers entitled "How to draw Euler angles and Utilize Euler parameters", DETC2006Paper99307.pdf (123 KB).
Finally, I handed out homework assignment 10 (due Thursday May 28)

Thu May 21 NO lecture.

Week 13

Thursday, May 28, 2009, 15:45-17:30 u, room C.

In this last lecture I presented two final projects:

You have to pick one final project. The final project is strictly individual. To get credit on the final project please do the following things:

Hand in your final project on paper at my office during office hours or in my paper mailbox. The due date is Monday June 15, 2009. The graded work can be picked up at my office, again during office hours. There will be no oral exam. The final grade for the course will be 70% on the homework and 30% on the final project.

This was the last lecture, Good Luck with the final project!

Week 14

Thursday, June 4, 2009, 15:45-17:30 u, room C.

NO lecture.