|Homework3: Constrained Dynamics
Implement dynamics with a simple constraint.
|(25pt) Problem A: A Bead on a Wire
Simulate the motion of a particle and visualize the result.
The particle is defined in a Cartesian coordinate system p=(x,y), and it
is constrained by a unit circle centered at the origin. Use the initial
position (0, 1), the initial velocity (0, 0), and gravitational force (0,-mg).
Constants are give as m=1Kg, and g=9.8m/sec2.
What I did:
I have created a small program.
program (you need to put qt-mt230nc.dll in
the same directory)
It requires Windows (2000?), openGL, GLVU
(to compile), and the QT runtime dll (included).
You can view the source
code. All the simulation code is contained in the main.cpp. The other
files are used for the user interface.
I chose to use the Lagrange’s equations
of motion, which is overkill for this circle case -- but still very easy.
The concept is to express the state of
the system in a natural manner, and then convert it to the "real world".
For a bead on a circle, the position can be stored as a single parameter
which indicates the angle.
Converting to world co-ordinates x,
is easy, just trigonometry. (See math below.) Then, the Jacobin matrix
J is solved. This matrix is built from the partial derivatives shown. Additionally,
the time derivative of J is found (J-dot).
The equation of motion is then shown. Notice
that I have moved the Mass matrix M over to the force f side
of the equation.
Then, some simple math and simplification
provides us with a solution for the acceleration of the angle variable
Finally, this system can be integrated
with Euler's method. The relevant code is shown below:
static Stopwatch timer;
circle.theta = 3.1415/2;
gravity.Set(0, -9.8, 0);
static float time_simulated = 0;
float time_now = timer.GetTime();
float time_passed = time_now - time_simulated;
circle.mass = pControlPanel->FloatSpinBox_mass->doubleValue();
circle.worldForces += circle.mass * gravity;
float force_amount = pControlPanel->FloatSpinBox_force->doubleValue();
circle.worldForces += Vec3f(-force_amount,0,0);
circle.worldForces += Vec3f(0,force_amount,0);
circle.worldForces += Vec3f(force_amount,0,0);
float timestep = pControlPanel->SpinBox_TimeStep->value()
float dampening = 1-pControlPanel->FloatSpinBox_dampening->doubleValue();
while ((time_simulated < time_now) &&
(timestep > 0))
float acceleration =
circle.theta += circle.thetaVelocity
circle.thetaVelocity += acceleration
circle.thetaVelocity *= dampening;
time_simulated += timestep;
The system performs very well, and is very
stable. The constraint can never be broken. Some dampening is required.