Transcript ME451 Kinematics and Dynamics of Machine Systems

ME451
Kinematics and Dynamics
of Machine Systems
Dynamics of Planar Systems
Tuesday, April 16, 2009
Constraint Reaction Forces – 6.6
© Dan Negrut, 2009
ME451, UW-Madison
Quote of the day:
Nobody goes there no more; it’s too crowded.
Yogi Berra
Before we get started…

On Tuesday…


No class, I’m out of town
Exam coming up on April 23

Review session on Wd evening, 7:15 pm




Look at an example
Answer your questions
We’ll meet in 3126ME
I’m open to having the Thursday exam in the evening


Thursday 7:30 – 9:00 PM would be a possibility
If we shift the time of the exam we will use the regular Thursday
lecture to cover one more example and review the material
2
Reaction Forces: The Basic Idea

Recall the partitioning of the total force acting on our mechanical system

Applying a variational approach (principle of virtual work) we ended up
with this equation of motion

After jumping through hoops, we ended up with this:

It’s easy to see that
3
The Important Observation

IMPORTANT OBSERVATION:

Actually, you don’t care for the “generalized” QC flavor of the reaction force,
but rather you want the actual force represented in the Cartesian global
reference frame



You’d like to have Fx, Fy, and a torque T that is due to the constraint
You report these quantities as they would act at a point P
The strategy:

Look for a force (the classical, non-generalized flavor), that when acting on
the body would lead to a generalized force equal to QC
4
The Nuts and Bolts



There is a joint acting
between Pi and Pj and we
are after finding the reaction
forces/torques Fi and Ti, as
well as Fj and Tj
Figure is similar to Figure
6.6.1 out of the textbook
Textbook covers topic well (pp. 234), I’m only modifying one thing:


The book expresses the reaction force/torque Fi and Ti in a body-fixed
reference frame
attached at point Pi
I didn’t see a good reason to do it that way

Instead, start by deriving in global reference frame OXY and then
multiply by AT to take it to the body-fixed reference frame
5
The Main Result
(Expression of reaction force/torque in a joint)

Suppose that two bodies i and j
are connected by a joint, and that
the equation that describes that
joint, which depends on the
position and orientation of the two
bodies, is

Suppose that the Lagrange multiplier associated with this joint is

Then, the presence of this joint in the mechanism will lead at point P on
body i to the presence of the following reaction force and torque:
6
Comments
(Expression of reaction force/torque in a joint)

Note that there is a Lagrange multiplier associated with each constraint equation

Number of Lagrange multipliers in mechanism is equal to number of constraints

Each Lagrange multiplier produces (leads to) a reaction force/torque combo

Therefore, to each constraint equation corresponds a reaction force/torque
combo that throughout the time evolution of the mechanism “enforces” the
satisfaction of the constraint that it is associated with


Since each constraint equation acts between two bodies i and j, there will also
be a Fj/Tj combo associated with each constraint, acting on body j


Example: the revolute joint brings along a set of two kinematic constraints and
therefore there will be two Lagrange multipliers associated with this joint
According to Newton’s third law, they oppose Fi and Ti, respectively
Note that you apply the same approach when you are dealing with driving
constraints (instead of kinematic constraints)

You will get the force and/or torque required to impose that driving constraint
7
Reaction Forces
~ Remember This ~


As soon as you have a joint
(constraint), you have a Lagrange
multiplier 

As soon as you have a Lagrange
multiplier you have a reaction
force/torque:
The expression of  for all the usual
joints is known, so a boiler plate
approach renders the value of the
reaction force in all these joints
Just in case you want another form for the torque T above, note that
8
Example 6.6.1: Reaction force in
Revolute Joint of a Simple Pendulum
Pendulum driven by motion:
1) Find the reaction force in the
revolute joint that connects
pendulum to ground at point O
2) Express the reaction force in
the
reference frame
9
End Constraint Reaction Forces
10
ME451
Kinematics and Dynamics
of Machine Systems
Dynamics of Planar Systems
April 30, 2009
Elements of the Numerical Solution for Ordinary Differential Equations
Chapter 7
Before we get started…

Exam coming up in one week (Tu, Dec. 4, 9:30 AM)





Exam draws on Dynamics Analysis topic (material covered since first
exam)
Take home component emailed to you over the weekend
Rough guidelines provided for what’s expected from you in terms of
the preliminary report on the Final Project
Review session at 7:15 PM on Monday, Dec. 3, in this room
Last Time

Discussed about initial conditions (ICs)



Second order differential equations needs ICs for positions and
velocities
These ICs should be consistent with the joint constraints associated
with your mechanism
Learned how to compute the constraint reaction force that appears in
the joints that connect bodies in a mechanism
12
Today’s Lecture…

Looking at algorithms to find an approximation of the solution
of the equations that govern the time evolution of a dynamic
system




Finding an exact solution within pen/paper framework impossible even
for the swinging motion of a pendulum in gravitational field
We need to resort to numerical methods (algorithms) to produce an
approximation of the solution
We’ll continue this discussion on Th when we focus on an ME451
specific method
Chris Hubert from Oshkosh Truck will talk about use of
kinematics/dynamics analysis in their work


Chris is with Modeling and Simulation group at Oshkosh Truck
Analysis goes beyond kinematics and dynamics (CFD, FEA, etc.)
13
Numerical Method
(also called Numerical Algorithm, or simply Algorithm)

Represents a recipe, a succession of steps that one takes to find an
approximation the solution of a problem that otherwise does not admit an
analytical solution



Analytical solution: sometimes called “closed form” or “exact” solution
The approximate solution obtained with the numerical method is also called
“numerical solution”
Examples:

Evaluate the integral

Solve the equation

Solve the differential equation that governs time evolution of simple pendulum

Many, many others (actually very seldom can you find the exact solution of a
problem…)
14
Where/How are Numerical Methods Used?

Powerful and inexpensive computers have revolutionized the use of
numerical methods and their impact
 Simulation of a car crash in minute detail
 Formation of galaxies
 Folding of a protein
 Finding the electron distribution around nuclei in a nanostructure

Numerical methods enable the concept of “simulation-based
engineering”
 You use computer simulation to understand how your mechanism
(mechanical system) behaves, how it can be modified and controlled
15
Numerical Methods in ME451

In regards to ME451, one would use numerical method to solve the dynamics problem
(the resulting set of differential equations that capture Newton’s second law)

The particular class of numerical methods used to solve differential equations is
typically called “numerical integrators”, or “integration formulas”

A numerical integrator generates a numerical solution at discrete time points
(also called grid points, station points, nodes)

This is in fact just like in Kinematics, where the solution is computed on a time grid

Different numerical integrators generate different solutions, but the solutions are
typically very close together, and [hopefully] closed to the actual solution of our
problem

Remember: In 99% of the cases, the use of numerical integrators is the only
alternative for solving complicated systems described by non-linear differential
equations
16
Numerical Integration
~Basic Concepts~


So what’s the problem?




Initial Value Problem:
(IVP)
You are looking for a function y(t) that depends on time (changes in time), whose
time derivative is equal to a function f(t,y) that is given to you (see IVP above)
In other words, I give you the derivative of a function, can you tell me what the
function is?
Remember that both y0 and the function f are given to you. You want to find y(t).
In ME451, the best you can hope for is to find an approximation of the unknown
function y(t) at a sequence of discrete points (as many of them as you wish)

The numerical algorithm produces an approximation of the value of the unknown
function y(t) at the each grid point. That is, the numerical algorithm produces y(t1),
y(t2), y(t3), etc.
17
Relation to ME451




When carrying out Dynamics Analysis, what you can
compute is the acceleration of each part in the model
Acceleration represents the second time derivative of your
coordinates
Somewhat oversimplifying the problem to make the point
across, in ME451 you get the second time derivate
This represents a second order differential equation since
it has two time derivatives taken on the position q
18
How do you go about solving an Initial Value Problem?
A look at:
Euler’s Method
Predictor-Corrector Method
Runge-Kutta Method
19
Numerical Integration:
Euler’s Method

Find solution of this
Initial Value Problem:
 y(t )  f ( y, t )

 y(t0 )  y0
The idea: at each grid point tk, turn the
differential problem into an algebraic
problem by approximating the value of
the time derivative:
This step is called “discretization”. It transforms
the problem from a continuous ODE problem
into a discrete algebraic problem
Euler’s Method (t is the step size):
yk 1  yk  t f (tk , y k )
20
y  10 y  0
Example:
y (0)  1
- Integrate 5 steps using Euler’s Method
- Compare to exact solution
Exact solution:
y(t) = e-10t
f(t,y) = -10y (note no explicit dependency on time t for f)
k=0
k=1
k=2
k=3
k=4
k=5
y0 = 1.0
y1 = y0+f(t0,y0)t
y2 = y1+f(t1,y1)t
y3 = y2+f(t2,y2)t
y4 = y3+f(t3,y3)t
y5 = y4+f(t4,y4)t
= 1.0
+ (-10*1.0 )*0.01 = 0.9
= 0.9
+ (-10*.9
)*0.01 = 0.81
= 0.81 + (-10*.81 )*0.01 = 0.729
= 0.729 + (-10*.729 )*0.01 = 0.6561
= 0.6561 + (-10*.6561)*0.01 = 0.5905
Solution:
y(0) =1.0000
y(0.01)=0.9048
y(0.02)=0.8187
y(0.03)=0.7408
y(0.04)=0.6703
y(0.05)=0.6065
21
Predictor-Corrector Methods

Instead of computing f(y,t) at one point, you can
average over multiple points (two in this case)
t
yk 1  yk   f (tk , y k )  f (tk 1 , y k 1 ) 
2

To implement, must predict yk+1 using forward Euler Method
yk( p1)  yk  t f (tk , y k )
yk 1  yk 
t
 f (tk , y k )  f (tk 1 , yk( p1) ) 
2
The predictor
The corrector
22
Example
- Integrate in Matlab for 1 second using Euler’s Method
- Compare to exact solution
y  10 y  0
y (0)  1
Matlab
y0=1;
% Euler’s Method
dt=.01;
t=0:.01:1.;
yh=y0;
for i=2:length(t)
f=-10*yh(i-1);
yh(i)=yh(i-1)+f*dt;
end
% Exact solution
y=y0*exp(-10*t);
% Plot and compare
plot(t,y,'b-',t,yh,'ro');
1
0.8
0.6
X

0.4
0.2
0
0
0.2
0.4
0.6
0.8
1
Time (s)
23
Euler Method:
~ Effect of Step Size ~


 y  0.1y  sin t

 y (0)  0
Solve using step sizes t=0.1, 1 and 5 sec
dt=5 sec
Matlab
% Exact solution
t=0:.01:50;
y=.99*exp(-t/10)+.995*sin(t-1.47);
% Numerical solution
dt=.1;
th=0:dt:50;
yh=0;
for i=2:length(th)
f=-yh(i-1)/10+sin(th(i));
yh(i)=yh(i-1)+f*dt;
end
plot(t,y,'b-',th,yh,'ro-');
dt=1 sec
dt=0.1 sec
6
4
2
0
-2
-4
-6
0
10
20
30
Conclusion: If you use large step-sizes t, the ACCURACY of the
solution is very poor (you can’t be too aggressive with size of t)
40
50
24
MATLAB Support for solving IVP
25
Ordinary Differential Equations
(Initial Value Problem)


 y  f (t , y )
An ODE + initial value: 
 y (t0 )  y0
Use ode45 for non-stiff IVPs and ode23t for stiff IVPs
(concept of “stiffness” discussed shortly)
[t,y] = ode45(odefun,tspan,y0,options)
function dydt = odefun(t,y)
[initialtime
Initialvlue
finaltime]
• Use odeset to define options parameter
26
IVP Example (MATLAB at work):
function dydt = myfunc(t,y)
dydt=zeros(2,1);
dydt(1)=y(2);
dydt(2)=(1-y(1)^2)*y(2)-y(1);
»
[t,y]=ode45('myfunc',[0 20],[2;0])
3
27
Note:
Help on odeset to set options
for more accuracy and other
useful utilities like drawing
results during solving.
2
1
0
-1
-2
-3
0
2
4
6
8
10
12
14
16
18
20
The concept of stiff differential equations,
and how to solve the corresponding IVP
28
Example: IVP
y  100 y  0 
100 t

y
(
t
)

e

y (0)  1 
- Integrate 5 steps using forward Euler formula: t=0.002, t=0.01, t=0.03
- Compare the errors between numerical and analytical solutions (Algorithm Error)
Algorithm Error
when t=0.002:
0
0.01873075307798
0.03032004603564
0.03681163609403
0.03972896411722
0.04019944117144
Algorithm Error
when t=0.01:
0
0.36787944117144
0.13533528323661
0.04978706836786
0.01831563888873
0.00673794699909
Algorithm Error
when t=0.03:
0
2.04978706836786
-3.99752124782333
8.00012340980409
-15.99999385578765
32.00000030590232
29
y  100 y  0 
100 t

y
(
t
)

e

y (0)  1 
Example:
5000
1
Numerical Solution
0.9
0.8
0.7
0
0.6
0.5
0.4
-5000
0.3
0.2
0.1
0
-10000
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
Forward Euler
Analytical Solution
(t=0.03)
30
Concept of stiff IVP’s

IVP’s for which forward Euler doesn’t work well (see example)
 In general, the entire class of so called explicit formulas doesn’t work


Forward Euler, Runge-Kutta (RK23, RK45), DOPRI5, Adams-Bashforth, etc.
Stiff IVP’s require a different class of integration formulas
 Implicit formulas

Example: backward Euler
yk 1  yk  f (tk 1 , y k 1 )t
31
Explicit vs. Implicit Formulas
(look at Euler family)

Initial Value Problem
 y  f (t , y )

 y (t0 )  y0
• Forward Euler
yk 1  yk  t yk
yk  f ( yk , tk )
yk 1  yk  f (tk , y k )t
• Backward Euler

yk 1  yk  t yk+1
  f ( yk 1 , tk 1 )
yk+1
yk 1  yk  f (tk 1 , y k 1 )t
32
y  100 y  0 
100 t

y
(
t
)

e

y (0)  1 
Example:
Exact Solution
1
5000
Solution
Approximation
Numerical Solution
0.9
0.8
0.7
0
0.6
0.5
0.4
-5000
0.3
0.2
0.1
0
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
-10000
0
0.05
Backward Euler and
Analytical Solution
0.1
0.15
0.2
0.25
0.3
0.35
0.4
Forward Euler
(t=0.03)
33
Other Popular Algorithms for
Stiff IVPs
The family of BDF methods (Backward-Difference Formulas):

BDF of 1st order:

yn1  yn  hyn+1

BDF of 2nd order:
4
1
2

yn+1= yn  yn-1  hyn+1
3
3
3

BDF of 3rd order:
yn+1=

BDF of 4th order:
yn+1=
48
36
16
3
12

yn 
yn-1 
yn-2 
yn-3  hyn+1
25
25
25
25
25

BDF of 5th order:
yn+1=
300
300
200
75
12
60

yn 
yn-1 
yn-2 
yn-3 
yn-4 
hyn+1
137
137
137
137
137
137
18
9
2
6

yn  yn-1  yn-2  hyn+1
11
11
11
11
34
The Two Key Attributes of a
Numerical Integrator

Two attributes are relevant when considering a
numerical integrator for finding an approximation of
the solution of an IVP

The STABILITY of the numerical integrator

The ACCURACY of the numerical integrator
35
Numerical Integration Formula:
The STABILITY Attribute

The stability question:

How big can I choose the integration step-size t and be safe?

Tough question answered in a Numerical Analysis class

Different integration formulas, have different stability regions

You’d like to use an integration formula with large stability region:
 Example:

Backward Euler, BDF methods, Newmark, etc.
Why not always use these methods with large stability region?
 There
is no free lunch: these methods are implicit methods that require the
solution of an algebra problem at each step (we’ll see this on Th)
36
Numerical Integration Formula :
The ACCURACY Attribute

The accuracy question:

How accurate is the formula that I’m using?

If I start decreasing t , how will the accuracy of the numerical solution
improve?



Tough question answered in a Numerical Analysis class
Examples:

Forward and Backward Euler: accuracy O(t )

RK45: accuracy O(t 4)
Why not always use methods with high accuracy order?
 There
is no free lunch: these methods usually have very small stability regions
 Therefore,
you are limited to very small values of t
37
ODE solvers in MATLAB
Solver
Problem
Type
Order of
Accuracy
ode45
Nonstiff
Medium
ode23
Nonstiff
Low
ode113
Nonstiff
Low to high
ode15s
Stiff
Low to
medium
ode23s
Stiff
Low
If using crude error tolerances to solve stiff systems and
the mass matrix is constant
ode23t
Moderately
stiff
Low
For moderately stiff problems is you need a solution
without numerical damping
ode23tb
Stiff
Low
If using crude error tolerances to solve stiff systems
When to use
Most of the time. This should be the first solver tried
For problems with crude error tolerances or for solving
moderately stiff problems.
For problems with stringent error tolerances or for solving
computationally intensive problems
If ods45 is slow because the problem is stiff
38
End Numerical Methods at Large
Beginning Numerical Methods for ME451
39
ME451
Kinematics and Dynamics
of Machine Systems
Dynamics of Planar Systems
November 29, 2007
Elements of the Numerical Solution of the Dynamics Problem
Chapter 7
Before we get started…

Exam coming up in one week (Tu, Dec. 4, 9:30 AM)





Exam covers Dynamics (material covered since first exam)
Take home component emailed to you: due at the time of the exam
Preliminary report on Final Project due
Review session at 6:00 PM on Monday, Dec. 3, in this room
Last Time

Discussed how numerical algorithms (methods) play an essential part
in finding the time evolution of a mechanism (or any dynamic system
for that matter)



Solution of IVP found by either explicit or implicit numerical integration
formulas
Key Idea: Differential problem transformed into Algebraic problem
through discretization (use of discretization formulas)
Two important attributes associated with a numerical integration
algorithm:


Stability
Accuracy
41
Summary of the Lagrange form of the
Constrained Equations of Motion

Equations of Motion:

Position Constraint Equations:
The Most
Important Slide
of ME451

Velocity Constraint Equations:

Acceleration Constraint Equations:
42
What’s special about ME451 problems?

Looking at the previous slide there are several things that make the
ME451 dynamics problem challenging:

The problem is not in standard form

Furthermore, our problem is not a first order Ordinary Differential Equation
(ODE) problem


Rather, it’s a second order ODE problem, due to form of the equations of
motion (contain the second time derivative of the positions)
After all, it’s not even an ODE problem

The unknown function q(t); that is, the position of the mechanism, is the
solution of a second order ODE problem (first equation previous slide) but
it must also satisfy a set of kinematic constraints at position, velocity, and
acceleration levels, which are formulated as a bunch of algebraic equations



To conclude, you have to solve a set of differential-algebraic equations (DAEs)
DAEs are much tougher to solve than IVPs
This lecture is about using a numerical method to solve DAEs of multibody dynamics
43
Example: Find the time evolution of the
pendulum

Simple Pendulum:



Mass 20 kg
Length L=2 m
Force acting at tip of pendulum


Although not shown in the picture,
assume RSDA element in revolute joint




F = 30 sin(2 t) [N]
k = 45 [Nm/rad] & 0=3/2
c = 10 [N/s]
ICs: hanging down, starting from rest
Stages of the procedure (three):



Stage 1: Derive constrained equations of motion
Stage 2: Indicate initial conditions (ICs)
Stage 3: Apply numerical integration algorithm to
discretize DAE problem and turn into algebraic problem
44
Example: Simple Pendulum
Stage 1: Deriving EOM (see also posted solution)
Stage 2: Indicate Initial Conditions (ICs)
45
Detour: Algorithm for Resolving
Dynamics of Mechanisms

If you have the EOM and ICs, how do you go
about solving the problem?

This is a research topic in itself

We’ll present almost the simplest algorithm possible

It is based on Newmark’s integration formulas

That is, we are going to use Newmark’s formulas to discretize
our differential problem
46
Solution Strategy: Important Slide

This slide explains the strategy through which the numerical solution; i.e., an
approximation of the actual solution of the dynamics problem, is produced

Step 1: two integration formulas (Newmark in our case) are used to express
the positions and velocities as functions of accelerations


These are also called “discretization formulas”
Step 2: everywhere in the constrained equations of motion, the positions and
velocities are replaced using the discretization formulas and expressed in
terms of the acceleration

This is the most important step, since through this “discretization” the differential
problem is transformed into an algebraic problem


The algebraic problem, which effectively amounts to solving a nonlinear system, is
approached using Newton’s method (so again, we need to produce a Jacobian)
Step 3: solve a nonlinear system to find the acceleration and the Lagrange
multipliers
47
Overview Newmark Integration Formulas

Newmark method (N.M. Newmark – 1957)

One step method designed to integrate directly second order equations
of motion:

The “discretization” formulas that relate position to acceleration and
velocity to acceleration are:

The goal is to find the positions, velocities, accelerations and Lagrange
multipliers on a grid of time points
48
Newmark (Cntd.)




Newmark Method
 1st Order
 Very good stability properties
 Choose =0.3025, and =0.6 (these are two
parameters that control the behavior of the method)
If we write the equation of motion at each time tn+1 one gets
Now is the time to replace
and
with the
discretization formulas (see previous slide)
You end up with an algebraic problem in which the unknown
is exclusively the acceleration
49
The Problem at Hand

Our rigid multibody dynamics problem is slightly more complicated than the Linear
Finite Element problem used to introduce Newmark’s discretization formulas

More complicated since we have some constraints that the solution must satisfy

We also have to deal with the Lagrange multipliers that come along with these constraints (from a
physical perspective remember that they help enforce satisfaction of the constraints)
Linear Finite Element
dynamics problem
Nonlinear multibody dynamics problem.
Newmark algorithm still works as before,
problem is slightly messier…
50
Quantities we are interested in

Generalized accelerations:

Generalized velocities:

Generalized positions:

Reaction Forces:
All these quantities are functions of time (they change in time)
51
Stage 3: The Discretization
of the Constrained Equations of Motion

The Discretized Equations Solved at each Time-Step tn+1:

Above you see
and
functions of the accelerations:
52
, but remember that they are
Again, these are Newmark’s formulas
that express the generalized position
and velocity as functions of
generalized accelerations
The Discretization
of the Constrained Equations of Motion (Cntd.)

Remarks on the discretization and its outcome:

Your unknowns are the accelerations and the Lagrange multipliers

The number of unknowns is equal to the number of equations

The equations that you have to solve now are algebraic and nonlinear

Differential problem has been transformed into an algebraic one


The new problem: find acceleration and Lagrange multipliers that satisfy
You have to use Newton’s method


This calls for the Jacobian of the nonlinear system of equations (chain rule will be
used to simplify calculations)
This looks exactly like what we had to do when we dealt with the Kinematics analysis
of a mechanism (there we solved (q,t)=0 to get the positions q)
53
Solving the Nonlinear System =0
~ Newton Method ~

Based on Newmark Integration formulas, the Jacobian is calculated as:

Corrections Computed as :
54
Note: subscripts “n+1”
dropped to keep
presentation simple
Simple Pendulum Example:
Looking at Stage 3
55
End Numerical Methods
Begin Formulating EOM in NonCartesian Generalized Coordinates
56
ME451
Kinematics and Dynamics
of Machine Systems
Dynamics of Planar Systems
December 6, 2007
EOM in non-Cartesian Reference Frames
~ not in textbook~
Before we get started…

HW due next Th


HW emailed to you, has two main components
 Formulating eqs of motion following example presented in class
 Find solution of EOM using MATLAB’s ode45.
Next Tu:

Director of Engineering Learning Center to come to class to discuss
with you about ME451 and my teaching


I will have to leave the room when this happens
Last Time

Exam

Hope to get it graded over the weekend and return next Tu
58
A Look Ahead

Cover of Kinematics and Dynamics of 2D mechanisms
has been completed

Moving on to discuss some topics related to material
covered

Formulating the equations of motion using non-Cartesian
generalized coordinates (today)

Equilibrium and Inverse Dynamics problems (Tuesday)

Rotation matrix A in 3D kinematics (Thursday, last lecture of
semester)
59
Today lecture’s Question:

You have a mechanism and are
interested in finding its time
evolution

What if you want to express the
equations of motion (EOM) in a set
of generalized coordinates that is
not Cartesian?

Example: Simple Pendulum
 Why not use the angle  to express
the time evolution of the pendulum?
60
Non-Cartesian GCs: Further comments

Benefits of working with non-Cartesian GCs:

You might be able to formulate the problem as an ODE problem as
opposed to a DAE problem as it is almost always the case when you use
Cartesian coordinates


For ODE problems you don’t have the complications deriving from the use of
Newmark’s formulas, and using any classical numerical integration formula
will do
The dimension of the problem in general is reduced

Typically, you can reduce it all the way to having a number of differential
equations equal to the number of degrees of freedom associated with your
mechanism

Example:
 For the simple pendulum you’d have one second order ODE
 For the double pendulum you’d have two second order ODEs
61
Non-Cartesian GCs: Further comments
(Cntd.)

Drawbacks of working with non-Cartesian GCs:

You might lose the ability to compute the reaction forces since
you eliminate their presence from the problem
 This goes back to trading the DAE for an ODE problem (there are
no Lagrange multipliers associated with an ODE problem…)

The expression that the ODE problem assumes is tremendously
messy
 Ex: Slider-Crank
62
EOM for Slider Crank:
Where k1, k2, and k3 are defined as:
Nomenclature & Notation
OLD
NEW
64
Example: Simple Pendulum

Scenario 1:


Note that there are no constraints on w
Scenario 2:
65
Taxonomy

Based on how we select the generalized coordinates w,
we can be in one of TWO cases:

1) There are no constraints that the new GCs must satisfy
 These GCs are called Lagrangian GCs
 They are equal in number to the number of degrees of freedom
 See Scenario 1 on previous slide, the angle  was not subject to
any constraints

2) There are some constraints that the new GCs must satisfy
 See Scenario 2 on previous slide
 w must satisfy in this case the nonlinear algebraic equations
66
The case of Lagrangian Gen. Coordinates



Here w can change in an arbitrary fashion, no
constraints that need to be satisfied
How is a virtual displacement in q related to a virtual
displacement in w?
Go back to variational form of the equations of motion:
67
The case of Lagrangian Gen. Coordinates
(Cntd.)

Note that the acceleration in q and w are related:

Simple manipulation lead to

From here, since w is arbitrary (no constraints acting on w),
68
Example: Simple Pendulum




Simple Pendulum:
Mass 20 kg
Length L=2 m
Force acting at tip of pendulum




F = 30 sin(2 t) [N]
ICs: hanging down, starting from rest
FIND EQUATION OF MOTION IN TERMS OF 
NOTE:

This is a Lagrangian generalized coordinates scenario
69
The case of non-Lagrangian and
non-Cartesian Gen. Coordinates
(this is pretty thick…)
70
The case of non-Lagrangian and nonCartesian Gen. Coordinates



Here w cannot change in an arbitrary fashion, it must satisfy
constraints
How is a virtual displacement in q related to a virtual
displacement in w?
It can be proved that if w represents a consistent virtual
displacement, that is,
… then by taking q = ww you got yourself a consistent virtual
displacement q, that is,
71
The case of non-Lagrangian and nonCartesian Gen. Coordinates

Proof of previous result is skipped
but it relies on the following
observation:

If w is a consistent configuration, that is
… then q= (w) is also a consistent
configuration, that is
72
The case of non-Lagrangian and nonCartesian Gen. Coordinates

Note: Here w must satisfy at all times

Variational form of the equations of motion:

The above holds only when w is consistent, that is

Applying Lagrange’s theorem, one ends up with the following:
73
What GCs do people use?

ADAMS uses Cartesian GCs, as does the second most widely used
commercial package out of US

Most widely used commercial packages out of Asia, RecurDyn, uses nonCartesian GCs

For open tree topologies they use Lagrangian GCs

Most widely used commercial packages out of Europe, SimPack, uses nonCartesian GCs

My take on this:

I prefer Cartesian GCs



EOMs are larger, but sparser
EOMs are easy to formulate, I let the computer deal with their solution
Where are non-Cartesian GCs widely used?


In robotics
When you have open tree topologies

In these cases you can actually easily select a set of Lagrangian GCs and end up with an ODE
problem for your EOM.
ME451
Kinematics and Dynamics
of Machine Systems
Dynamics of Planar Systems
December 11, 2007
Inverse Dynamics 6.4
Equilibrium Analysis 6.5
Before we get started…

HW due Th

Today: Director of Engineering Learning Center to come to class to
discuss with you about ME451 (last 30 minutes of the class)

You can have two HWs dropped (instead of one) when you fill out survey put
together by the Engineering Learning Center




Provides a snapshot of my teaching in this class. It takes 10 minutes to fill out
Helps with organizing & teaching ME451 better next time
I emailed details about survey yesterday. Very grateful if you fill out.
Last Time

Discussed about formulating the equations of motion using a set of
generalized coordinates (GCs) different than Cartesian GCs


The idea is that you use precisely the same ingredients that you use in the
Cartesian case
Additionally, you need only one thing – a “bridging” function that tells you how
the Cartesian GCs are obtained in terms of the new GCs that you prefer to use
Cartesian GCs
New GCs
76
Exam 2: Grade Distribution
100%
90%
80%
70%
60%
1 2 3

Exam 2 Grades:


Average: 88 %
Solution posted online
4 5 6 7 8
9 10 11 12 13 14 15 16 17 18 19 20 21 22
Inverse Dynamics: The idea

First of all, what does dynamics analysis mean?



In *inverse* dynamics, the situation is quite the opposite:


You apply some forces/torques on a mechanical system and look at how the
configuration of the mechanism changes in time
How it moves also depends on the ICs associated with that mechanical system
You specify a motion of the mechanical system and you are interested in finding
out the set of forces/torques that were actually applied to the mechanical system
to lead to this motion
When is *inverse* dynamics useful?

It’s useful in controls. For instance in controlling the motion of a robot: you know
how you want this robot to move, but you need to figure out what joint torques you
should apply to make it move the way it should
78
Inverse Dynamics: The math

When can one talk about Inverse Dynamics?

Given a mechanical system, a prerequisite for Inverse Dynamics is that the
number of degrees of freedom associated with the system is zero

You have as many generalized coordinates as constraints (THIS IS KEY)

This effectively makes the problem a Kinematics problem

The two stages of the Inverse Dynamics analysis

First solve for accelerations (recall the acceleration equation):

Next you solve for the reaction forces:
79
Inverse Dynamics: Closure

Are we done once we computed the reaction forces?

Yes, because among the forces you computed, you get all
the forces/torques that are necessary to impose the driving
constraints D that you imposed on the system
Here constraint D acts between body i
and some other body. Reaction forces
are computed as “felt” by body i

This gives you the forces/torques that you need to apply to
get the prescribed motion
80
Example: Inverse Dynamics

Compute torque that electrical motor
applies to open handicapped door like

Free angle of the spring:

Mass m = 30
Mass Moment of Inertia J’ = 2.5
Spring/damping coefficients:
K=8
C=1
All units are SI.



81
82
End Inverse Dynamics
Beginning Equilibrium Analysis
83
Equilibrium Analysis: The Idea

A mechanical system is said to be in equilibrium if the
following conditions hold:

Equivalently, the system is at rest, with zero acceleration

So what does it take to be in this state of equilibrium?



You need to be in a certain configuration q
The applied forces should assume certain values
(OK, tell me more about what “certain” means…)
84
Equilibrium Analysis: The Math

Equations of Motion:

Position Constraint Equations:

Velocity Constraint Equations:

Acceleration Constraint Equations:
85
Equilibrium Analysis: The Math
(Cntd.)

To conclude, one needs a
configuration q and the applied
forces QA should be such that

How can you go about finding such a
configuration?

Approach 1 (dumb, but powerful)


Approach 2 (OK, but you need a good
starting point)


Add damping in a system and watch it
move till it stops
Simply solve the nonlinear system to
find q and QA
Approach 3 (not that common)


Cast it as an optimization problem
Works for conservative systems only
86
Example: Pendulum Equilibrium

Free angle of the spring:

Spring constant: k=25

Mass m = 10

Length L=1

All units are SI.
87
88
ME451
Kinematics and Dynamics
of Machine Systems
Dynamics of Planar Systems
December 13, 2007
Rotation in 3D
Before we get started…

Last class of the semester, no more HW ;-)

Final Project presentation coming up on Dec. 22



One presentation per group
Each team member gets to talk and then answers questions at
end of presentation
Final Project Score




Project and PPT presentation packaging & looks - 10%
Project content – 65 %
Your presentation of your contribution to the project – 25%
 Each of you will have to stand up and make a presentation
 The 25% is meant to reflect your contribution to the project
Email was sent to you regarding scheduling of the presentations
90
The Problem at Hand…

In Kinematics, the starting point is to position and orient a body in space

In 2D Kinematics,




The position of a body is given by x and y coordinates of the origin of a
REFERENCE FRAME attached to the rigid
The orientation of that body is then described by the orientation angle 
You end up with three generalized coordinates the describe the
position/orientation
In 3D Kinematics,


For position you need x, y, and z, no big difference other than that
However, for orientation you need three angles



Orientation in 3D is quite a different animal and it leads to some major differences
You end up with six generalized coordinates the describe the position/orientation
THIS LECTURE IS A VERY QUICK OVERVIEW OF HOW PEOPLE GO
ABOUT DESCRIBING (SPECIFYING) THE ORIENTATION IN 3D
KINEMATICS
91
Specifying Orientation in 3D

There are several ways to describe how a rigid
body is oriented in space

Euler angles (probably the most common)

Rodriguez Parameters

Quaternions
 Euler Parameters, as a sub-case

Bryant angles
92
Purpose of the Orientation Matrix A

The matrix A is what you need to take
the representation of a vector in a
local RF and then find its
representation in the global RF

Recall that when we deal with s`P, we
have a representation of a vector in the
local RF

If you want to get the location of P in
the global RF,

THIS RELATION STAYS THE SAME
FOR 3D KINEMATICS
How is A obtained?

Suppose that you have a global RF, OXYZ, and some
other RF, O`x`y`z`.

Given a vector s` represented in O`x`y`z`, you’d like to know
what it’s representation of the vector is in the global OXYZ
reference frame
This is what the matrix A does…

It turns out that




The first column of A is the representation in the global reference
frame of the unit vector along the O`x` axis
The second column of A is the representation in the global
reference frame of the unit vector along the O`y` axis
The third column of A is the representation in the global reference
frame of the unit vector along the O`z` axis
94
How is A obtained?
(Cntd.)

Example: Derive the orientation matrix for the 2D case

Steps:



Consider unit vector in the x` direction and represent in global RF
Consider unit vector in the y` direction and represent in global RF
Assemble matrix A by using the two columns determined above
95
How is A obtained?
(Cntd.)

It all hinges upon representing a vector, which is a
quantity/entity independent of a reference frame, in two
different reference frame:
This is the important step: the unit
vectors i’ and j’, when expressed in the
new reference frame, become ex and ey
This is our buddy A…

NOTE: if you express any other vector
, you end up with the same matrix
A. In other words, it is an invariant associated with switching from one
reference frame to a different one (and it’s not specific to
or
above):
96
Z
Euler Angles

y
Euler Angles:


A set of three angles used to
describe the orientation of a
reference frame in 3D space
Draws on the following
observation (intuitive but
mathematically deep/rich, not
going to get into details):


You can align the global
reference frame to any
arbitrary reference frame
through a sequence of THREE
rotation operations
z
Y
x
X
The sequence of three angle rotations that we’ll be working with is
about axes Z then X then Z again (sometimes called 3-1-3 sequence)
 The three angles are denoted by , , and , respectively, and called the
Euler angles
http://prt.fernuni-hagen.de/lehre/KURSE/PRT001/course_main_sh3d/node10.html
97
Relationship between Euler Angles
and Orientation Matrix A

Why do we keep talking about both the orientation
matrix A and the Euler angles , ,  in parallel?


Because these angles will help us get the orientation matrix
A of the local reference frame with respect to the global
reference frame
In other words, if I choose any vector b` represented in the
local reference frame, it’s representation (“image”) in the
global reference frame will be
98
Expressing A using Euler Angles
(Part 1 of 3)




Recall that I have a 3-1-3 rotation sequence.
The last sequence, is a rotation of angle  about z`` to get
O`x`y`z`.
Therefore, the rotation matrix that relates O`x`y`z` to O``x``y``z`` is
In other words, if I have a vector represented as a` in O`x`y`z`, it
will be represented in O``x``y``z`` as
99
Expressing A using Euler Angles
(Part 2 of 3)


Just dealt with the last “3” in the 3-1-3 rotation sequence (angle )
Focus next on the “1” rotation in the 3-1-3 rotation sequence



This the rotation of angle 
The rotation of angle  about the axis O```x```, takes
O```x```y```z``` into O``x``y``z``
In other words, if I have a vector represented in O``x``y``z`` as
a``, the same vector will be represented in O```x```y```z``` as
100
Expressing A using Euler Angles
(Part 3 of 3)


Just dealt with the “1” in the 3-1-3 rotation sequence (angle )
Focus next on the first “3” rotation in the 3-1-3 rotation sequence



This is the rotation of angle 
The rotation of angle  about the axis OX, takes OXYZ into
O```x```y```z```
In other words, if I have a vector represented in O```x```y```z``` as
a```, the same vector will be represented in OXYZ as
101
Expressing A using Euler Angles
~ Putting it Together ~


Z
A Z-X-Z (or 3-1-3) rotation
sequence takes OXYZ to
O`x`y`z`
Recall that

Given a`, you get a`` as

Given a`` you get a``` as
102
y
z
Y
x

Given a``` you get a as
X

Therefore,
Expressing A using Euler Angles
~ Putting it Together ~

Using the expression of the matrices A1, A2, A3, one gets for the
expression of the orientation matrix A

Carry out multiplications to get

Note: it’s easy to see that AT A=A AT=I
103
Expressing A using Euler Angles
~ Critical Considerations ~

Some questions that should be asked
1.
2.
Should the order (sequence)  -  -  be observed?
In general, how many parameters do you actually need to capture
an arbitrary orientation of a rigid body?
3.
For any arbitrarily oriented reference frame, can you always find a
set of 3-1-3 Euler angles to relate the global and local reference
frame?
4.
Suppose that you can always find one set of  -  -  angles. Are
they unique? Do you run into any type of singularity issue when
you have an infinite number of 3-1-3 rotations that can relate the
global and local reference frame?
104
Expressing A using Euler Angles
~ Critical Considerations ~

1.
2.
3.
4.
Answers to previous four questions:
The order of the rotation should be observed since it’s important. If you change
the order this will land you in a completely different configuration
Euler’s theorem says that you need three parameters to represent an arbitrary
orientation of any rigid body
Yes, you can always find a set of angles  -  -  that take the global reference
frame into another reference frame and the transformation matrix is A given on the
previous slide (the really ugly 3X3 matrix)

In fact the situation can become quite bad, when you end up having not one
but an infinite number of rotations to lead to a good matrix A (see below)
The only time when you run into a problem is when the angle  turns out to be
zero (the “1” rotation in the 3-1-3 sequence is zero). In this case it’s not clear how
to split the rotations about the Z axes (the first and the last rotations “3”).
105
Expressing A using Euler Angles
~ Critical Considerations ~

Example:

Take =20, =0, and =60. Compute A
Take =60, =0, and =20. Compute A
Take =40, =0, and =40. Compute A

In all cases you will end up with the same A




This is called the singularity associated with the Euler
angles, and there is no way around it
That’s why people have been attempting to use
different generalized coordinates to determine in a
more robust fashion the orientation of a body with
respect to a global reference frame
106
Angular Velocity of a Body
107
Preliminary Discussion

The concept of skew-symmetric matrix (for 3X3 matrices)



Matrix D is skew-symmetric if it satisfy the condition
Note that diagonal entries of a skew-symmetric matrix are zero
Moreover, the skew-symmetric matrix will look like

Note that all it takes to assemble the matrix D is a set of three values:
d1, d 2, and d3

In general you’d store only the vector d = [ d1 d2 d3]T and quickly
generate D on an as needed basis


For this reason, d is called the “generator” of D
If you have the generator vector d, the matrix D is denoted by “tilde d”:
108
The KEY Slide

Note that since A is orthonormal,

Take a time derivative to eventually get

This means that the matrix product at the left is a skew-symmetric matrix
 Therefore there is a generator vector  that can be used to represent this
Identity matrix
skew-symmetric matrix:


By definition, the generator vector  is called the angular velocity of
the body on which the local reference frame is attached
NOTE:  is an attribute of the body, and not of the reference frame
attached to the body (it’s easy to see this…)
109
Angular Velocity for 3-1-3
Euler Sequence

When using the 3-1-3 Euler angles, based on the expression
of A, one gets:

110
What you
additionally know is that:
Angular Velocity for 3-1-3 Euler Sequence


If you carry out the horrific matrix multiplications you end up
with the following expression for the generator vector :
Equivalently, you can write this in matrix form like
…where the array  is defined as
111
Angular Velocity for 3-1-3 Euler Sequence

The previous slide says that you can compute the time derivatives
of the Euler angles given  by solving the linear system
Note that det(B)=sin, which yet again indicates that when the second
rotation (the “1” in the 3-1-3 sequence) is a integer multiple of , the
matrix B is singular and you are in trouble yet again (singular
configuration)


112
To avoid this situation people introduced a set of four “Euler parameters”
which always work, except that you have some redundancy in
information (three parameters are necessary, yet you have four Euler
parameters ) leads to one additional constraint equation)
 This is ME751 material…
Title
113
Title
114
Title
115
Title
116
Title
117
Title
118
Example 2.4.1

Express the position of P as a function of the two angles 1 and 2
11
Virtual Displacement: Examples
12
Example 3.2.3: Slider Crank
(HW, due date: TBD)

For slider crank in Figure 3.1.3, a motion is prescribed on 1





Find the equations of motion (use Cartesian coordinates as in Figure)
Write MATLAB or C code to find the time evolution of the mechanism
Plot the position, velocity, and acceleration of points O2 and P1
Is the motion of O2 in the x direction sinusoidal? Explain why (don’t use
the plots, but look at the equations)
Is the motion of P1 in the y direction sinusoidal? Explain why (don’t use
the plots, but look at the equations)
121
ME451
Kinematics and Dynamics
of Machine Systems
Absolute Constraints 3.2
Relative Constraints 3.3
September 25, 2007
© Dan Negrut, 2009
ME451, UW-Madison
Before we get started…

Next lecture:
 We’ll continue coverage of Section 3.3
 HW: MATLAB problem (to be emailed to you this morning), 3.1.1, 3.1.2,
3.1.4, 3.3.2, 3.3.5

Last Time

We discussed some very key concepts:


In Kinematics, it all boils down to being able to identify the set of constraints
associated with your mechanism
If you have these constraints (expressed as (q,t)=0) the time evolution of
the kinematic system is easy to determine:




Perform position analysis (challenging)
Perform velocity analysis (simple)
Perform acceleration analysis (OK)
If you keep these concepts in mind the first half of the semester will
make sense, otherwise you won’t be seeing the forest for the leaves
123
Important Slide of Last Lecture

The most critical point in this course is being able to produce the set of constraints that
are present in the physical system

The three stages of Kinematics Analysis: position analysis, velocity analysis, and
acceleration analysis they each follow *very* similar recipes for finding for each body of
the mechanism its position, velocity, and acceleration, respectively

ALL STAGES RELY ON THE CONCEPT OF JACOBIAN MATRIX:

q – the partial derivative of the constraints wrt the generalized coordinates

ALL STAGES REQUIRE THE SOLUTION OF A SYSTEM OF EQUATIONS

WHAT IS DIFFERENT BETWEEN THE THREE STAGES IS THE EXPRESSION OF
THE RIGHT-SIDE OF THE LINEAR EQUATION, “b”
124
Example

Consider the slider-crank in Fig. P3.3.2. Come up with the set of
kinematic constraint equations to kinematically model this mechanism

Length of crank and connecting rod is L1 and L2
125
ME451
Kinematics and Dynamics
of Machine Systems
Position, Velocity, and Acc. Analysis 3.6
October 16, 2007
Before we get started…

HW due next Tu:

Redo Problem 3.5.5 in ADAMS


Are you getting the same results as in the solution posted at Learn@UW?
Last Time
 ADAMS demo

Prior to that (last Tu):


Finished driving constraints
Started to talk about the Kinematic Analysis solution procedure:

Step A: Identify *all* physical joints and drivers present in the system
Step B: Identify the corresponding constraint equations (q,t)
Step C: Solve for the Position as a function of time (q is needed)
Step D: Solve for the Velocities as a function of time ( is needed)

Step E: Solve for the Accelerations as a function of time ( is needed)



12
Grade Distribution: ME451 Exam 1
Exam 1
Grades:
100
95
90



Max: 100
Min: 72.5
Average: 91.45
85
80
75
70
65
60
1

2
3
4
5
6
7
8
9
10
11 12
13 14 15
16 17
18 19 20
Further remarks:




Each bar in plot above represents a student (22 students in class)
With few exceptions, almost everybody received a full 50% for the
take home component
Difficulties taking partial derivatives (the bonus problem)
Overall, the results were good
21 22