# Simulating mechanics of a cylinder rolling on a parabola

GROUPS:
 Hi,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 ?ThanksJef
5 years ago
12 Replies
 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,Peter
5 years ago
 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.
5 years ago
 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.
5 years ago
 @Aronsson, is there a way to change the profile of the ground?
5 years ago
 David Keith 5 Votes 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,DavidRemove symbols and load the VariationalMethods packageRemove["Global*"];Needs["VariationalMethods"];SetOptions[Plot, PlotRange -> All];SetDirectory[NotebookDirectory[]];The constraints to follow the parabola will appear as factors to a Lagrange multiplierconstraint0 = 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 \ individualyxy\[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 carfPlot = 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 factorarrow[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[   Graphics[    {     {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]
5 years ago
 David Keith 1 Vote Hello All,My first attempt at a post was pretty clumsy -- I've edited this to make it readble.Best regards,David
5 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 :-)Jef
5 years ago
 Erik Mahieu 2 Votes 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: Manipulate[  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"}],     Dynamic@Show[{        Plot[{a (x^2 - 1), -a Sqrt[1 - x^2]}[[i]], {x, -1, 1},          PlotStyle -> Thick],        Graphics[{          Red,           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
5 years ago
 Moderation Team 1 Vote 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!