Posted: February 12, 2013
This post was inspired by a blog post made by Matt Sottile.
His blog post discussed implementing a dynamic system model of an interesting experiment from the article Pendulum waves: A demonstration of wave motion using pendula. He chose to implement his model in Haskell. I only came to know about this because of a conversation between Matthis Thorade and Matt Sottile about the suitability of using Modelica for building such a model.
The experiment was so fascinating to watch, I decided that I just had to recreate it in Modelica. This example turns out to be not only a great way to demonstrate how well suited Modelica is for such problems, but also how the array capabilities in particular make building such a model quite straightforward. After I published the models to GitHub and wrote this post, Dietmar Winkler sent me a pull request with improved parameters that led to a nicer visualization. Those changes have been incorporated in the video at the end of this post.
Rather than reading a detailed account of the physics behind this experiment (the links above will give you plenty of that), just watch this video and you'll get the idea:
The key thing here is that the lengths of the strings that the balls are suspended from differ in precisely the way necessary to give them interesting relative natural frequencies. Specifically, the length of each pendulum needs to follow this relationship:
\[ l_i = g \left( \frac{T}{2 \pi (X+n-i)} \right)^2 \]
where \(n\) is the number of pendula, \(g\) is the Earth's gravitational constant, \(T\) is 54 seconds and \(X\) is 60 (again, see Matt's blog post for all the details).
Note, I've neglected damping so the system should (mathematically speaking) return to its initial conditions after 54 seconds.
OK, so hopefully you have an understanding of the physical experiment now. The question is then, "How would we model this in Modelica?" The answer comes in two parts.
We start by modeling a single pendulum. As usual, the source code for this is post is hosted on GitHub. There, we find that the source code for an individual pendulum looks like this:
within HarmonicMotion;
model Pendulum "A single individual pendulum"
parameter Modelica.SIunits.Position x "Axial position of this pendulum";
parameter Modelica.SIunits.Angle phi "Initial angle";
parameter Modelica.SIunits.Length L "String length";
parameter Modelica.SIunits.Diameter d=0.01 "Reference diameter for visualization";
parameter Modelica.SIunits.Mass m "Mass of mass point";
Modelica.Mechanics.MultiBody.Parts.Fixed ground(r={0,0,x}, animation=false);
Modelica.Mechanics.MultiBody.Parts.PointMass ball(m=m, sphereDiameter=5*d);
Modelica.Mechanics.MultiBody.Parts.BodyCylinder string(
density=0, r={0,L,0}, diameter=d);
Modelica.Mechanics.MultiBody.Joints.Revolute revolute(
phi(fixed=true, start=phi),
cylinderDiameter=d/2,
animation=false);
equation
connect(string.frame_a, ball.frame_a);
connect(revolute.frame_b, ground.frame_b);
connect(revolute.frame_a, string.frame_b);
end Pendulum;
The four key parameters here are x, phi, L and d. The x
parameter represents the position of this pendulum along the axis of
rotation and we will use this parameter to position instances of this
pendulum one after another. The phi parameter represents the
initial angle of the pendulum and is fed into the Revolute joint as
an initial condition. The L parameter is used to scale the
BodyCylinder model which represents a mass-less string from which the
ball hangs. Finally, d is a representative length scale used to
determine the physical sizes of the various 3D objects in the scene.
All of these parameters are used to specify the characteristics of
the components that make up the pendulum (via modifications).
This is where things get interesting because in order to create a collection of pendula all swinging with different natural frequencies, we'll need to use some tricks in Modelica for dealing with arrays.
The actual model of our system is simpler than the pendulum:
within HarmonicMotion;
model System "A system of pendula"
import Modelica.Constants.g_n;
import Modelica.Constants.pi;
parameter Integer n=10 "Number of pendula";
parameter Modelica.SIunits.Time T = 54;
parameter Modelica.SIunits.Time X = 60;
parameter Modelica.SIunits.Angle phi0 = 0.5;
parameter Modelica.SIunits.Position x[n] = linspace(0,n-1,n);
parameter Modelica.SIunits.Length lengths[n] = { g_n*(T/(2*pi*(X+(n-i))))^2 for i in 1:n};
Pendulum pendulum[n](x=x, each m=1, each phi=phi0, L=lengths);
inner Modelica.Mechanics.MultiBody.World world;
end System;
The majority of the model is related to specifying the parameters of the model. There are really only two component declarations. Let's go through these one at a time.
The first thing we need to do is specify the number of pendula we want. This is pretty straightforward:
parameter Integer n=10 "Number of pendula";
We also need to specify a couple of the characteristics of the dynamic system as well as the initial angle we will release the pendula from. These are also pretty straightforward:
parameter Modelica.SIunits.Time T = 54;
parameter Modelica.SIunits.Time X = 60;
parameter Modelica.SIunits.Angle phi0 = 0.5;
Now things get interesting. We now need to specify the axial position of the pendula along the axis of rotation and the lengths of each pendulum. Since these values are different for each pendulum, we need to create an array of these values, one for each pendulum:
parameter Modelica.SIunits.Position x[n] = linspace(0,n-1,n);
parameter Modelica.SIunits.Length lengths[n] = { g_n*(T/(2*pi*(X+(n-i))))^2 for i in 1:n};
We can see in creating the lengths array that we have used array
comprehensions in Modelica to define an expression for each element in
the array parameterized by the pendulum number, i.
Once we've done all this, the rest of the model is only two more lines:
Pendulum pendulum[n](x=x, each m=1, each phi=phi0, L=lengths);
inner Modelica.Mechanics.MultiBody.World world;
The declaration of the inner World element is required by the
multi-body dynamics library in Modelica to specify some parameters of
the 3D system. The creation of the pendula is simply the declaration
of the array pendulum[n]. When specifying the parameters of the
pendulum, the modifications come in two flavors. The first, qualified
by the each keyword provides an expression which is then used to
modify each element in the array. In other words, all elements of
the array will have the same value for this parameter. The other
flavor is when each element in the array should have a different
value. In this case, we leave off the each qualifier and pass an
array of data. The elements of such an array (e.g., lengths) will
then be applied element-wise to the underlying component being
modified, pendulums in this case.
That's it! That is all the Modelica code we need.
Before demonstrating the results, let's take a second to comprehend
the size of the underlying problem. I've used the
Modelica.Mechanics.Multibody library. As a result, this is a
complete 3D model. This means it computes the 3D position and 3D
forces and torques for every physical element in the model. I could
trivially do all kinds of interesting things like incline the angle of
the axis of rotation or change the pivots of pendula into universal
joints instead of revolute joints.
I compiled this model with Dymola, which provides some interesting statistics about the model:
Original Model
Number of components: 299
Variables: 2855
Constants: 20 (20 scalars)
Parameters: 1051 (2003 scalars)
Unknowns: 1784 (5149 scalars)
Differentiated variables: 230 scalars
Equations: 1515
Nontrivial : 986
Translated Model
Constants: 4552 scalars
Free parameters: 233 scalars
Parameter depending: 692 scalars
Continuous time states: 20 scalars
Time-varying variables: 410 scalars
Alias variables: 1285 scalars
Assumed default initial conditions: 10
Number of mixed real/discrete systems of equations: 0
Sizes of linear systems of equations: {7, 7, 7, 7, 7, 7, 7, 7, 7, 7}
Sizes after manipulation of the linear systems: {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
Sizes of nonlinear systems of equations: { }
Sizes after manipulation of the nonlinear systems: { }
Number of numerical Jacobians: 0
So to be clear, this system has 5149 scalar unknowns. Now, if you've read Matt Sottile's blog post, you will immediately wonder why I'm approaching the problem this way. Clearly, you can simulate this system with *far fewer equations and unknowns.
Well, the answer is simple. It was actually easier to do it this way. I already had all the 3D primitives and it was pretty easy to just drag, drop and connect everything together with the existing components. I certainly could have done it using the simpler 1D rotational pendulum model, but then I'd have to deal with adding the animation after the fact.
But the bigger question is, what did this cost me, if anything? Yes, the system is much, much larger. But building it was quick. Compiling the model was quick. The actual compilation time (this includes flattening, equation sorting, index reduction, code generation and code compilation) took about 1 second in my estimation. Finally, the simulations were quick. Setting extremely tight relative tolerances, \(10^{-12}\), so that the end condition after 54 seconds would closely match the initial conditions, the simulation of the entire system took only 0.445 seconds. The results include 6000 result points and all the physical signals and animation data. This means I can plot things like the reaction forces in the joints or the tension in the string in addition to the position of the pendulum.
My point here is that symbolic manipulation of the kind that most Modelica tools do means really fast simulation code. So it usually doesn't make sense for engineers to spending much time or effort up front trying to simplify their system because the tools will do that for them automatically anyway. In fact, in the end Dymola reduced the problem down to only 20 states, exactly how many you'd get if you formulated it as a 1D rotational problem in the first place.
If you have a Modelica tool at your disposal, you can simulate the model yourself since the source code is available on GitHub. But for those who just want to see an animation of the simulation results, here they are:
Note how the system returns to almost exactly the initial conditions? There will invariably be some numerical artifact that prevents 100% of the energy from being conserved (at least with the integration schemes I'm using). But given all the motion, it is remarkably close at the end. You can loop this video and never really notice the differences.
When I started this, I just thought it would be a cool model to build since it is fascinating to watch the motion. But in the end, I recognized that this model is a great way to showcase both the benefits of symbolic manipulation and some of the nice array features in Modelica.
Share your thoughts
comments powered by Disqus