Message Boards Message Boards

12 Replies
11 Total Likes
View groups...
Share this post:

Simulating mechanics of a cylinder rolling on a parabola

Posted 11 years ago

I'm trying to simulate a cylinder (or a sphere) of radius R rolling without slipping on a parabola of equation y = ax^2 to study the relation between the frequency of the oscillations and the parameters (R and a). Can I do that with SystemModeler ?

POSTED BY: Jef Barthe
12 Replies
Posted 11 years ago
Pretty much any system that you can express as a system of ordinary differential equations you can simulate in SystemModeler.

For a simple system like this one you can alternatively use Mathematica.

Hope this helps,
POSTED BY: Peter Fleck
Easier and fun if you try the cycloid tautochrone at first. Radius of curvature is bisected by x-axis, it is not too complicated. 

The property of same time period  for any boundary condition ( found by Christian  Huygens) can be verified. hth.
Look at the RollingWheel example in Modelica.Mechanics.MultiBody.Examples.Elementary package. Its probably a bit overkill for you since it is a full-fledged 3D mechanical model, but it show the principle of how to add constraints such as a wheel rolling on a surface. In this example the constraint is simply z=0, but that can be easily adapted.
POSTED BY: Peter Aronsson
@Aronsson, is there a way to change the profile of the ground? 
POSTED BY: Shenghui Yang
Posted 11 years ago
Hello Jef,

( I have edited this post. My beginners attempt was pretty clumsy -- I've since learned how to do better.)

This can be done in Mathematica using Lagrangian mechanics. The code below (something I once did for fun) solves a similar problem using that method. It considers a roller coaster car on a parabolic track with a simple damping. The Lagrangian is T-V, with T = kinetic energy and V = potential energy. It uses the VariationalMethods package to obtain the equations of motion left hand side from the Lagrangian. The method of Lagrange multipliers (one is needed) is used to enforce the constraint that the car follows the track. That produces forces of constraint which appear on the right hand side of the equations of motion, together with the damping force.

For your problem, there is an additional term required in the kinetic energy Iw^2, where I is the moment of inertia (related to your R) and w is the radian velocity. Of course, the radian velocity is related to the linear velocity of the car by the no-slip condition, so that wR = v, with v the magnitude of the velocity of the car.

It's worth noticing that once you have the equations of motion in symbolic form, you see a damped harmonic oscillator, with the oscillation frequency in the x and y directions equal and in phase. So you know a great deal about the system without actually performing a numerical simulation.

You did not indicate where you are in your learning. These methods are at the level of an advanced undergrad in physics. And if your not quite there yet may not be appropriate for now. But you clearly have an interest, so I thought you might enjoy the response.

Kind regards,


Remove symbols and load the VariationalMethods package
SetOptions[Plot, PlotRange -> All];

The constraints to follow the parabola will appear as factors to a Lagrange multiplier
constraint0 = y[t] - x[t]^2;
constraint = y[t] == x[t]^2;
constraintd1 = D[constraint, {t, 1}];
constraintd2 = D[constraint, {t, 2}];
aX = D[constraint0, x[t]];
aY = D[constraint0, y[t]];

Use the variational derivative to obtain the Euler-Lagrange equations. The Lagrange multiplier and the damping factor B appear on the RHS of the equations of motion.
L = m (x'[t]^2 + y'[t]^2)/2 - m g y[t];
eqX = VariationalD[L, x[t],
    t] == -\[Lambda][t] aX + \[CapitalBeta] x'[t];
eqY = VariationalD[L, y[t],
    t] == -\[Lambda][t] aY + \[CapitalBeta] y'[t];

Use NDSolve to solve for the trajectory and the multipliers
 tmax = 50;
 equations = {constraintd2, eqX, eqY};
 values = {g -> 1, m -> 1, \[CapitalBeta] -> 0.2};
 {x, y, \[Lambda]} = NDSolveValue[
     equations /. values,
     x[0] == -1, x'[0] == 0, y[0] == 1, y'[0] == 0
     }, {x, y, \[Lambda]}, {t, 0, tmax}];
 maxTime = x[[1, 1, 2]];

Plot x, y, and \ individualy
xy\[Lambda]Plot =
Plot[{x[t], y[t], \[Lambda][t]}, {t, 0, maxTime},
  PlotLegends -> {"X", "Y",
    "\[Lambda]"}]; Export["xylPlot.png", xy\[Lambda]Plot];

The multipler \ with its factor is the force of constraint acting on the car
fPlot = Plot[\[Lambda][t] {aX, aY} // Evaluate, {t, 0, maxTime},
  PlotLegends -> {"Fx", "Fy"}]; Export["fxyPlot.png", fPlot];

Make an arrow function to plot a vector representation of the force of constraint. "a" is just a scaling factor
arrow[t_, a_] =
  Arrow[{{x[t], y[t]}, {x[t], y[t]} +
     a {aX \[Lambda][t], aY \[Lambda][t]}}];

Plot an animation of the equations. The arrow represents the force exerted on the car to keep it on the track
 animation = Animate[
     {Thick, Plot[x^2, {x, -1, 1}][[1]]},
     {PointSize[0.025], Point[{x[t], x[t]^2}]},
     {Red, Thick, Arrowheads[0.025], arrow[t, .12]}
     }, PlotRange -> {{-1, 1}, {-.1, 1}}
    ], {t, 0, maxTime, .01},
   DefaultDuration -> 25]

POSTED BY: David Keith
Posted 11 years ago
Hello All,
My first attempt at a post was pretty clumsy -- I've edited this to make it readble.
Best regards,
POSTED BY: David Keith
Posted 11 years ago
I worked on the original question to see the influence of the parameter a in the parabola on the period. Inspired by the previous post of David, I included damping and added an elliptic path as a second choice. A cycloidal path would be very interesting but is more difficult (to me) since there seems to be no cartesian equation of the form y=f to fit inside my code? Here is my attempt so far:
  Quiet@DynamicModule[{g = 9.81, m = 1., deqns, xsol, ysol, times,
     periods, periodNumber},
    deqns = {(*parabolic*){a (1 + x[t]^2) Derivative[2][x][
          t] == -\[Mu] Derivative[1][x][t]/m -
         x[t] (g + 2 a Derivative[1][x][t]^2)},
      (*elliptic*){\[Mu] Derivative[1][x][t]/m + Derivative[2][x][t] ==
         x[t] (-g Sqrt[1 - x[t]^2]/a +
           Derivative[1][x][t]^2/(-1 + x[t]^2))}};
   {xsol, times} =
    Reap@NDSolveValue[{deqns[[i]], x[0] == .999,
       Derivative[1][x][0] == 0, WhenEvent[x[t] == 0, Sow[t]]},
      x, {t, tMax}];
   ysol[t_] := {(*parabolic path*)a (xsol[t]^2 - 1),
      (*elliptic path*)-a Sqrt[1 - xsol[t]^2]}[[i]];
   periodNumber =
    Position[times[[1]], Nearest[times[[1]], tt][[1]]][[1, 1]] - 1;
   periods = Differences[times[[1]]];
   Column[{Row[{"swing# ", periodNumber, ": ",
       If[periodNumber != 0, periods[[periodNumber]], 0], " seconds"}],
        Plot[{a (x^2 - 1), -a Sqrt[1 - x^2]}[[i]], {x, -1, 1},
         PlotStyle -> Thick],
          Disk[.98(*scale to slide inside curve*){xsol[tt],
             ysol[tt]}, .02]}]}, AspectRatio -> Automatic,
       ImageSize -> 400, PlotRange -> {{-1, 1}, {-1.25, 0.05}}]}]],
{{i, 1, "path"}, {1 -> "parabolic", 2 -> "elliptic"}, Setter},
{{a, 1}, .5, 1.25, .001, Appearance -> "Labeled"}, {{\[Mu], .15},
  0, .25, .01, Appearance -> "Labeled"},
Control[{{tt, 0}, 0, .5 tMax, .01, Trigger,
   AnimationRate -> 1}], {{tMax, 100}, None}]
BTW: how can I post the interactive output of this MCode? Erik
POSTED BY: Erik Mahieu
David: that's some very nice code and images - .GIF gives excellent feel of the oscillations - thank you!

Eric: About posting interactive output - this is currently under development. For now we can post only images, animated .GIFs including. Take a look at this post:  Showcasing Manipulate[…] via .GIF animations  - maybe this can be a temporarily solution. And of course you can always publish a Demonstration!
POSTED BY: Moderation Team
Posted 11 years ago
Hi David,

Thank you VERY MUCH for your post. It is extremely helpful. The mechanical problem itself is well within my field of expertise but well beyond my knowledge of Mathematica. The problem I was working on was the impact of the radius of the sphere on the frequency of the undamped oscillation, and the behaviour when this radius is getting close to 1/2a (curvature at the origin). I should be able to derive that from your code :-)

POSTED BY: Jef Barthe
Posted 11 years ago
As you can see from my code, the Reply window where you paste the MCode converts x'' into (x^''). Although x' is converted to Derivative[1]! This must be an error in the reply system?

Anyway, there is no problem with my code if I preventively convert x'' into Derivative[2]. I corrected the code in the original post. I copied, this code and pasted into a nb. Works fine.

But I could not convert it into a GIF although I used Vitaly's function ManToGif. keep on
getting:Import::fmterr: "Cannot import data as Quicktime format\!"
POSTED BY: Erik Mahieu
Eric, thank you for explanations and sending corrected code. We actually think it is best to update the code in your original post - so people will not get confused by copying it into a notebook. We went ahead and updated the original code and left just explanatory text in your last post. We hope you are fine with it.

To avoid corruption of code during posting we recommend first to press the MCode button, and then paste Mathematica code in the appeared gray box.

About ManToGif errors. We think you have elaborate structure of controls and dynamic variables so ManToGif cannot process it. We used another external free program LICEcap to capture control functionality. If you like it you can easily use the program yourself.

POSTED BY: Moderation Team
Posted 11 years ago
Thanks a lot, moderation team! Very useful this dynamic screen capture! Erik
POSTED BY: Erik Mahieu
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract