Message Boards Message Boards

2
|
7025 Views
|
9 Replies
|
10 Total Likes
View groups...
Share
Share this post:

Impact of a ball against a wall

Posted 9 years ago

How do I simulate the impact of a ball with a wall a finite number of times (say 5)?

I could simulate the impact of a pendulum with a wall at x=0 (5 times). I used the Event Locator method for the same. How do I do it similarly for a ball hitting the floor at y=0?

q = {{\[Theta][t]}};
dq = D[q, t];
ddq = D[dq, t];

x = R*Sin[\[Theta][t]];    (*X and Y Coordinates*)
y = -R*Cos[\[Theta][t]];

dx = D[x, t];
dy = D[y, t];


KE = (1/2)*m*(dx^2 + dy^2);
PE = m*g*y;
L = KE - PE;   (*Lagrangian is the difference of Kinetic and \
Potential Energies*)

(*The constraint*)
\[Phi] = x;
d\[Phi] = D[\[Phi], q];

Eq = (Thread[
   D[D[L, dq\[Transpose]], t] - D[L, q\[Transpose]] == 
    0]);    (*Solving the Euler Lagrange Equations*)

ELtemp = Solve[Eq[[1]], Flatten[ddq]];

EL = {\[Theta]''[t] == ELtemp[[1, 1, 2]]};

(*Hamiltonian*)
p = D[L, dq\[Transpose]];
H = p.dq - L;

R = 1;
m = 1;
g = 9.81;

ICnt = 5;(*The impact is to be taken 5 times*)
\[Theta]zero = 
 Pi/2; \[Theta]dzero = 0; Tzero = 0;(*Setting the initial conditions*)
\

\[Theta]1 = {}; Time = {}; P1 = {}; Hset = {};

(*The calculations of the Lagrangian just before and after impact \
with the wall*)
FullSimplify[D[L, dq]*dq - L];
(*This is the condition just before impact and is constrained to hit \
the wall along x=0*)
FullSimplify[
  eq1 = (Thread[(D[L, dq] /. t -> Tplus) - (D[L, dq] /. 
         t -> Tminus) == \[Lambda]*d\[Phi]])];
(*This is the condition just after inpact and there is no constraint*)


FullSimplify[
  eq2 = (Thread[((D[L, dq]*dq - L) /. 
         t -> Tplus) - ((D[L, dq]*dq - L) /. t -> Tminus) == 0])];
ImpactTemp = 
  Solve[eq1[[1]] && eq2[[1]], Join[{\[Theta]'[Tplus]}, {\[Lambda]}]];

For[
  i = 0, i < ICnt, i++, 
  Initial = {\[Theta][Tzero] == \[Theta]zero, \[Theta]'[
      Tzero] == -\[Theta]dzero};
  sol = NDSolve[Join[EL, Initial], {\[Theta]}, {t, Tzero, 10}, 
    Method -> {"EventLocator", "Event" -> R*Sin[\[Theta][t]], 
      "EventAction" :> 
       Throw[{Tend = 
          t, \[Theta]end = \[Theta][t], \[Theta]dotend = \[Theta]'[
           t]}, "StopIntegration"]}];
  \[Theta]sol = \[Theta][t] /. sol[[1]]; Hsol = H /. sol[[1]];
  \[Theta]1 = Append[\[Theta]1, \[Theta]sol]; 
  Time = Append[Time, Tend];
  Hset = Append[Hset, {Hsol, Tzero < t <= Tend}]; 
  P1 = Append[P1, {\[Theta]sol, Tzero < t <= Tend}];
  \[Theta]zero = \[Theta]end; \[Theta]dzero = \[Theta]dotend;
  Tzero = Tend
  ];
    \[Theta]finalsol = Piecewise[P1];
    Hfinalsol = Piecewise[Hset];

(*Hamiltonian and its Plot*)
(*Hamiltonian and its Plot
p=D[L,Transpose[dq]];
H=p.dq-L/.sol[[1]]/.Hsol[[1]];
Plot[H,{t,0,10}]*)

(*Hamiltonian and its Plot*)
Plot[\[Theta]finalsol, {t, 0, Tend}, 
 AxesLabel -> {t, HoldForm[\[Theta]]}]
Plot[Hfinalsol, {t, 0, Tend}, AxesLabel -> {t, HoldForm[H]}]

(*Plot[xsol,{t,0,10}, AxesLabel\[Rule]{t,x}]
Plot[ysol,{t,0,10}, AxesLabel\[Rule]{t,y}]*)
X[\[Tau]_] := R*Sin[\[Theta]finalsol] /. t -> \[Tau];
Y[\[Tau]_] := -R*Cos[\[Theta]finalsol] /. t -> \[Tau];

Animate[Graphics[{PointSize[0.02], PlotRange -> {{-3, 3}, {-3, 3}}, 
   Line[{{-1.5, 0}, {1.5, 0}}], Line[{{0, 1.5}, {0, -1.5}}], 
   Line[{{0, 0}, {X[\[Tau]], Y[\[Tau]]}}], 
   Point[{X[\[Tau]], Y[\[Tau]]}]}], {\[Tau], 0, Tend}]
Attachments:
POSTED BY: Varun Kulkarni
9 Replies
Posted 9 years ago

No problem. I enjoy this. I'm always learning something new in this forum. If the new question is on another topic, you might consider putting it in a new post. No one on this site minds seeing a lot of questions!

POSTED BY: David Keith
Posted 9 years ago

Hi,

While going through the Wolfram Demonstration Project site, I came across an animation of a Trebuchet, http://demonstrations.wolfram.com/OptimizingTheCounterweightTrebuchet/.

I was wondering if it is possible to animate the Trebuchet in such a way that the projectile of the trebuchet actually leaves the sling and makes impact against a wall at a certain distance say x='some constant'.

Thanks

POSTED BY: Varun Kulkarni
Posted 9 years ago

And for "How do I simulate the impact of a ball with a wall a finite number of times (say 5)?"

This code counts bounces, and stops the integration at the top after the 5th bounce:

eqs = {x''[t] == 0, x'[0] == 10, x[0] == 0, y''[t] == -10., 
   y'[0] == 50, y[0] == 0};

events = {
   WhenEvent[y[t] < 0, {y'[t] -> -.9 y'[t], counter = counter + 1}],
   WhenEvent[y'[t] < 0, If[counter >= 5, "StopIntegration"]]
   };

counter = 0; {xx, yy} = 
 NDSolveValue[{eqs, events}, {x, y}, {t, 0, Infinity}];

ParametricPlot[{xx[t], yy[t]}, {t, 0, xx[[1, 1, 2]]}]

enter image description here

POSTED BY: David Keith
Posted 9 years ago

Hi David,

Is there any way to animate a 2-D body like a rectangle or a square that follows the same motion as the ball above for 5 bounces using elastic impacts so that the height of the bounce is equal everytime?

Thanks.

POSTED BY: Varun Kulkarni
Posted 9 years ago

Sure. First, the height of the bounce was controlled by the -.9 factor you see in the velocity of the rebound. (Physically, the energy at the top is mgh, and at the bounce is 1/2 mv^2, so the relationship is not linear.) We can make the height the same by using -1 for the y velocity change. As to what we draw, the solution just gives us functions for the location. This can be used for any graphic, with Animate or Export used to produce animations. The attached notebook demonstrates this.

Of course, this just moves the object, which still bounces like a ball. It would be more difficult to cause the object shape to be involved in detecting the hit on the surface, and far more difficult to include it in the physics. For example, by really simulating the linear and rotational motions of a cube made of elastic material bouncing chaotically.

Best, David

Attachments:
POSTED BY: David Keith
Posted 9 years ago

Hi,

Works well. Thanks

I had a doubt regarding another question that I saw, would it be ok if I ask.

Hope that I am not bothering you with my trivial queries.

Thanks

POSTED BY: Varun Kulkarni
Posted 9 years ago
q = {{x[t]}, {y[t]}};

dq = D[q, t];
ddq = D[dq, t];


dx = D[x, t];
dy = D[y, t];
ddx = D[dx, t];
ddy = D[dy, t];


KE = 0.5 m (x'[t]^2 + y'[t]^2);
PE = m*g*y[t];

L = KE - PE;

(*The constraint*)
\[Phi] = y;
d\[Phi] = D[\[Phi], q];

Eq = (Thread[D[D[L, dq\[Transpose]], t] - D[L, q\[Transpose]] == 0]);

ELtemp = Solve[Eq[[1]] && Eq[[2]], Flatten[ddq]];

EL = {x''[t] == ELtemp[[1, 1, 2]], y''[t] == ELtemp[[1, 2, 2]]};

p = D[L, dq\[Transpose]];
H = p.dq - L;

m = 1;
g = 9.8;

ICnt = 5;
xzero = 0; xdzero = 3; yzero = 10; ydzero = 0; Tzero = 0;
x1 = {}; y1 = {}; Time = {}; P1 = {}; Hset = {};

(*The calculations of the Lagrangian just before and after impact \
with the wall*)
FullSimplify[D[L, dq]*dq - L];
(*This is the condition just before impact and is constrained to hit \
the wall along x=0*)
FullSimplify[
  eq1 = (Thread[(D[L, dq] /. t -> Tplus) - (D[L, dq] /. 
         t -> Tminus) == \[Lambda]*d\[Phi]])];
(*This is the condition just after inpact and there is no constraint*)


FullSimplify[
  eq2 = (Thread[((D[L, dq]*dq - L) /. 
         t -> Tplus) - ((D[L, dq]*dq - L) /. t -> Tminus) == 0])];
ImpactTemp = 
  Solve[eq1[[1]] && eq2[[1]], 
   Join[{x'[Tplus]}, {y'[Tplus]}, {\[Lambda]}]];

For[i = 0, i < ICnt, i++, 
  InitCon = {\[Theta][Tzero] == xzero && 
     yzero, \[Theta]'[Tzero] == -xdzero && -ydzero};
  sol = NDSolve[Join[EL, Initial], {x, y}, {t, Tzero, 30}, 
    Method -> {"EventLocator", "Event" -> y, 
      "EventAction" :> 
       Throw[{Tend = t, yend = y[t], ydotend = y'[t]}, 
        "StopIntegration"]}];
  xsol = x[t] /. sol[[1]]; ysol = y[t] /. sol[[1]]; 
  Hsol = H /. sol[[1]];
  x1 = Append[x1, xsol]; y1 = Append[y1, ysol]; 
  Time = Append[Time, Tend];
  Hset = Append[Hset, {Hsol, Tzero < t <= Tend}]; 
  P1 = Append[P1, {xsol && ysol, Tzero < t <= Tend}];
  yzero = yend; ydzero = ydotend;
  Tzero = Tend
  ];
    \[Theta]finalsol = Piecewise[P1];
    Hfinalsol = Piecewise[Hset];



(*InitCon={x'[0]\[Equal]3, x[0]==0, y[0]==10, y'[0]==0};


sol=NDSolve[Join[EL,InitCon],{x,y},{t,0,30}] ;
xsol=x[t]/.sol[[1]];
ysol=y[t]/.sol[[1]];*)

(*PLOTS*)
Plot[xsol, {t, 0, 2}, AxesLabel -> {t, x}]
Plot[ysol, {t, 0, 2}, AxesLabel -> {t, y}]
ParametricPlot[{xsol, ysol}, {t, 0, 2}, AxesLabel -> {x, y}]

(*Hamiltonian and its Plot*)
Plot[\[Theta]finalsol, {t, 0, Tend}, 
 AxesLabel -> {t, HoldForm[\[Theta]]}]
Plot[Hfinalsol, {t, 0, Tend}, AxesLabel -> {t, HoldForm[H]}]

(*Plot[xsol,{t,0,10}, AxesLabel\[Rule]{t,x}]
Plot[ysol,{t,0,10}, AxesLabel\[Rule]{t,y}]*)
X[\[Tau]_] := R*Sin[\[Theta]finalsol] /. t -> \[Tau];
Y[\[Tau]_] := -R*Cos[\[Theta]finalsol] /. t -> \[Tau];


(*ANIMATION*)
X[\[Tau]_] := xsol /. \[Theta]finalsol /. t -> \[Tau];
Y[\[Tau]_] := ysol /. \[Theta]finalsol /. t -> \[Tau];

Animate[Graphics[{PointSize[0.05], Point[{ X[\[Tau]], Y[\[Tau]]}], 
   Line[{{-10, 0}, {10, 0}}], Line[{{0, -10}, {0, 10}}]}], {\[Tau], 0,
   2}]

I tried this to get a ball to bounce to 5 times and then stop. But it gives me certain errors like equations are not a quantified system of equalities and some replacement rules are not valid.

Any suggestions on how to proceed?

Thanks Varun

Attachments:
POSTED BY: Varun Kulkarni

How do I do it similarly for a ball hitting the floor at y=0?

You just need kinematics for this. May be something like this? (this was for a HW assignment, so code is not very clean. Did it in about 3 hrs). Most of the code deal with handling the state of the running simulation (reset, stop, step, etc...) if you do not care about all of this, you can remove all that and the code will shrink to 10% of its current size.

enter image description here

Manipulate[
 tick;
 h = 1;
 g = 9.8;
 If[(statex == "RUN" || statex == "STEP") && statex2 == "",
  If[state == "down",
   delS = currentV*delT + 1/2 g delT^2;
   currentV += g *delT;
   currentH -= delS
   ,
   delS = currentV*delT - 1/2 g delT^2;
   currentV -= g *delT;
   currentH += delS
   ];
  distantTravelled += delS
  ];

 gr = Graphics[
   {
    Line[{{-.3, -r/2}, {.3, -r/2}}],
    (*{LightGray,Dashed,Line[{{0,0},{0,1}}]},*)
    {Red, Disk[{0, currentH}, r]}
    },
   PlotRange -> {{-.3, .3}, {-4 r, 1 + 4 r}}, Axes -> {False, True}, 
   ImageSize -> {300, 300}];
  gr = Grid[{
    {Grid[{
       {"height", "speed", "cycle #", "\[CapitalDelta]s", "direction"},
       {padIt2[currentH, {3, 2}], padIt2[currentV, {4, 3}], padIt2[n, {1}], 
        padIt2[delS, {4, 3}], state}
       }, Frame -> All]
     },
    {Grid[
      {
       {Column[{"Theoretical", "total time"}], Column[{"Current", "time"}], 
        Column[{"Theoretical", Column[{"total", "distant"}]}], 
        Column[{"current", "distance"}]},
       {padIt2[tTime, {6, 3}], padIt2[currentT, {6, 3}],
         padIt2[tDistance, {3, 2}], padIt2[totalDist, {6, 3}]
        }
       }, Frame -> All
      ]
     },

    {gr, SpanFromLeft}}];

 If[statex2 == "pass",
  statex2 = ""
  ,
  currentT += delT;
  totalDist += delS;
  If[Abs@distantTravelled >= maxH,
   distantTravelled = 0;

   If[state == "down",
    currentV = e*currentV;
    state = "up";
    n = n + 1;
    maxH = e^(2*n)*h;
    currentH = 0
    ,
    state = "down";
    currentV = 0;
    currentH = maxH
    ]
   ]
  ];

 If[statex == "RUN" && currentT < tTime,
  tick = Not[tick]
  ];
 gr,
 {{tick, False}, None},
 Text@Grid[{
    {Grid[{
       {Button[Text@Style["run", 12], {statex = "RUN"; tick = Not[tick]}, 
         ImageSize -> {50, 40}],
        Button[Text@Style["step", 12], {statex = "STEP"; tick = Not[tick]}, 
         ImageSize -> {50, 40}],
        Button[Text@Style["stop", 12], {statex = "STOP"; tick = Not[tick]}, 
         ImageSize -> {50, 40}],
        Button[
         Text@Style["reset", 12], {statex = "STEP"; currentH = 1; 
          currentT = 0; n = 0;
          distantTravelled = 0; currentV = 0; maxH = 1; state = "down"; 
          totalDist = 0; statex2 = "";
          tDistance = If[e == 1, Infinity, (1 + e^2)/(1 - e^2)];
          tTime = If[e == 1, Infinity, Sqrt[2/9.81]*((1 + e)/(1 - e))];
          tick = Not[tick]}, ImageSize -> {50, 40}]}
       }, Frame -> True, FrameStyle -> Gray
      ], SpanFromLeft},
    {Grid[{

       {"Coefficient of restitution",
        Manipulator[Dynamic[e, {e = #;
            statex = "STEP";
            currentH = 1;
            currentT = 0; n = 0;
            distantTravelled = 0;
            currentV = 0;
            maxH = 1;
            state = "down";
            totalDist = 0;
            tDistance = If[e == 1, Infinity, (1 + e^2)/(1 - e^2)];
            tTime = If[e == 1, Infinity, Sqrt[2/9.81]*((1 + e)/(1 - e))];
            tick = Not[tick];
            statex2 = "pass"} &], {0, 1, .01}, ImageSize -> Tiny], 
        Dynamic[padIt2[e, {2, 2}]],
        SpanFromLeft},

       {"Animation speed",
        Manipulator[Dynamic[delT, {delT = #} &], {0.001, 0.03, 0.001}, 
         ImageSize -> Tiny], Dynamic[padIt2[delT, {3, 3}]],
        SpanFromLeft},

       {"ball size",
        Manipulator[
         Dynamic[r, {r = #; tick = Not[tick], statex2 = "pass"} &], {0.01, 
          0.1, 0.001}, ImageSize -> Tiny], Dynamic[padIt2[r, {3, 3}]],
        SpanFromLeft}

       }, Alignment -> Left, Frame -> True, FrameStyle -> Gray]
     }}],
 {{n, 0}, None},
 {{currentH, 1}, None},
 {{currentT, 0}, None},
 {{state, "down"}, None},
 {{statex, "STEP"}, None},
 {{statex2, ""}, None},
 {{distantTravelled, 0}, None},
 {{currentV, 0}, None},
 {{maxH, 1}, None},
 {{gr, 0}, None},
 {{e, .9}, None},
 {{delT, 0.02}, None},
 {{r, 0.04}, None},
 {{totalDist, 0}, None},
 {{tDistance, (1 + (.9)^2)/(1 - (.9)^2)}, None},
 {{tTime, Sqrt[2/9.81]*((1 + .9)/(1 - .9))}, None},
 TrackedSymbols :> {tick},
 Alignment -> Center,
 SynchronousUpdating -> True,
 SynchronousInitialization -> True,
 FrameMargins -> 1,
 ImageMargins -> 1,
 ControlPlacement -> Left,
 Initialization :>
  {
   padIt1[v_, f_List] := 
    AccountingForm[Chop[v], f, NumberSigns -> {"-", "+"}, 
     NumberPadding -> {"0", "0"}, SignPadding -> True];
   padIt2[v_, f_List] := 
    AccountingForm[Chop[v], f, NumberSigns -> {"", ""}, 
     NumberPadding -> {"0", "0"}, SignPadding -> True]
   }
 ]
POSTED BY: Nasser M. Abbasi
Posted 9 years ago

Works perfectly.

Say if I want to use Event Locator[] like the one I have used for the pendulum and solve an impact equation, do i use it inside Manipulate?

POSTED BY: Varun Kulkarni
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract