Sticky Collisions (perfectly inelastic)

Question

There is an awesome simulation https://mathematica.stackexchange.com/a/124926/87086 which shows N particles bouncing in a box (with elastic collisions), made by @Feyre https://mathematica.stackexchange.com/users/7312/feyre. Since then, they have moved on from commenting on StackExchange, and I was wondering how one could modify such an example's WhenEvent(s) such that when any number of particles collide, instead of ricocheting off of one another, they "stick" together and continue on as clump, with momentum:

m*v1 + m*v2 + ...m*v_n = n*m*v3

(where n is the number of particles in the clump)

Once stuck, I would like them to remain in contact with one another (maintaining their clumpiness) and be treated as a new body which is made up of the n particles that have collided (not a single particle with mass n*m and velocity v3).

Is there any way to go about doing this with WhenEvents at all? I have been at it for a while, and am at a loss as of now. Any insight would be greatly appreciated.

I have included the code below, but upon following the first link listed in this post you will find the simulation that @Feyre produced.

data = Table[{RandomReal[{-0.85, 0.85}, {3}], 
RandomReal[{-1, 1}, {3}], 
RGBColor[RandomReal[{-1, 1}], RandomReal[{-1, 1}], 
 RandomReal[{-1, 1}]]}, {i, 10}];
n = Length[data];
positions = Transpose[Table[data[[i, 1]], {i, n}]];
velocities = Transpose[Table[data[[i, 2]], {i, n}]];
(Setting up stuff for NDSolve[])
{xt, yt, 
   zt} = {ToExpression[Table["x" <> ToString[i], {i, n}]], 
   ToExpression[Table["y" <> ToString[i], {i, n}]], 
   ToExpression[Table["z" <> ToString[i], {i, n}]]};
{xm, ym, zm} = {Through[xt[t]], Through[yt[t]], Through[zt[t]]};
{xz, yz, zz} = {Through[xt[0]], Through[yt[0]], Through[zt[0]]};
yt = ToExpression[Table["y" <> ToString[i], {i, n}]];
zt = ToExpression[Table["z" <> ToString[i], {i, n}]];
rm = Flatten[
    Table[If[i != j, 
      Sqrt[(xm[[j]] - xm[[i]])^2 + (ym[[j]] - ym[[i]])^2 + (zm[[j]] - 
           zm[[i]])^2]], {i, n}, {j, n}]] /. Null -> Sequence[];
(No friction or gravity:)
xf = Thread[D[D[xm, t], t] == ConstantArray[0, n]];
yf = Thread[D[D[ym, t], t] == ConstantArray[0, n]];
zf = Thread[D[D[zm, t], t] == ConstantArray[0, n]];
(Friction and gravity:)
xf = Thread[D[D[xm, t], t] == -0.01 D[xm, t]];
yf = Thread[D[D[ym, t], t] == -0.01 D[ym, t]];
zf = Thread[D[D[zm, t], t] == -0.01 D[zm, t] - ConstantArray[0.5, n]];
(The final equations for the differential equation:)
pos = {Thread[xz == positions[[1]]], Thread[yz == positions[[2]]], 
   Thread[zz == positions[[3]]]};
vel = {Thread[D[xm, t] == velocities[[1]]] /. t -> 0, 
   Thread[D[ym, t] == velocities[[2]]] /. t -> 0, 
   Thread[D[zm, t] == velocities[[3]]] /. t -> 0};
col = Table[data[[i, 3]], {i, n}];
we = {};
(Events)
Table[AppendTo[we, 
   WhenEvent[
    Sqrt[(xm[[i]] - xm[[j]])^2 + (ym[[i]] - ym[[j]])^2 + (zm[[i]] - 
           zm[[j]])^2] <= 0.2 // Evaluate, {
      D[xm[[i]], t] -> 0.5 (D[xm[[i]], t] + D[xm[[j]], t]), 
      D[xm[[j]], t] -> 0.5 (D[xm[[i]], t] + D[xm[[j]], t]), 
      D[xm[[i]], t] -> 0.5 (D[ym[[i]], t] + D[ym[[j]], t]), 
      D[xm[[j]], t] -> 0.5 (D[ym[[i]], t] + D[ym[[j]], t]), 
      D[xm[[i]], t] -> 0.5 (D[zm[[i]], t] + D[zm[[j]], t]), 
      D[xm[[j]], t] -> 0.5 (D[zm[[i]], t] + D[ym[[j]], t])} // 
     Evaluate]], {i, n - 1}, {j, i + 1, n}];
Table[AppendTo[we, 
   WhenEvent[Abs[xm[[i]]] >= 0.9 // Evaluate, 
    D[xm[[i]], t] -> -D[xm[[i]], t] // Evaluate]], {i, n}];
Table[AppendTo[we, 
   WhenEvent[Abs[ym[[i]]] >= 0.9 // Evaluate, 
    D[ym[[i]], t] -> -D[ym[[i]], t] // Evaluate]], {i, n}];
Table[AppendTo[we, 
   WhenEvent[Abs[zm[[i]]] >= 0.9 // Evaluate, 
    D[zm[[i]], t] -> -D[zm[[i]], t] // Evaluate]], {i, n}];
s = NDSolve[Flatten[{xf, yf, zf,
     vel, pos, we}], Flatten[{xt, yt, zt}], {t, 0, 50}, 
   MaxSteps -> ∞];
Manipulate[
 Graphics3D[
  Transpose[{col, 
    Thread[Sphere[
      Evaluate[
       Flatten[Thread[Transpose[{xm, ym, zm} /. t -> a] /. s], 1]], 
      0.1]]}], PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}], {a, 0, 50, 
  0.1}]```

Should the particles alos stuck at contact with the boundary? — Ulrich Neumann, Jul 22 '22 at 06:23
They should bounce off all boundaries, but can be stopped at the bottom. Radius of sphere looks like 1 — Zach, Jul 22 '22 at 14:52
But calculation of data in the range -.85.....85 speaks for radius==0.15. Also the box is in the range of -1...+1, though radius=?=1 isn't plausibel. — Ulrich Neumann, Jul 22 '22 at 15:43
For me the main problem of the simulation is the contur of the new clumped masses. Have to think about ... — Ulrich Neumann, Jul 23 '22 at 08:13

Ulrich Neumann · Accepted Answer · 2022-07-24T14:36:12.513

Your expectation "masses clue together after inelastic collision" isn't correct. Only in the direction of impact the velocities are the same. In the perpendicular direction the velocities don't change!

I adapted your code and introduced a restitution coefficient "stosszahl" which is 1 for elastic impact and 0 for inelastic impact. For more details see wikipedia: Coefficient of restitution

To describe the collision I introduced new vectors position piand velocity vi.

The direction of impact i and j is defined as e=(pi[[i]]-pi[[j]])/Norm[(pi[[i]]-pi[[j]])], e.e==1 . Only in this direction velocity changes.

Coefficient of restitution stosszahl== (vni[[j]]-vni[[i]]).e/(vi[[i]]-vi[[j]]).e (vni: velocity after impact)

data = Table[{RandomReal[{-0.85, 0.85}, {3}], 
    RandomReal[{-1, 1}, {3}], 
    RGBColor[RandomReal[{-1, 1}], RandomReal[{-1, 1}], 
     RandomReal[{-1, 1}]]}, {i, 10}];
damping = 0 .01;
gerd = 10;
stosszahl = 1 ; (* Stoßzahl 0[plastisch]...1[elastisch])
(stosszahl=(V2-V1)/(v1-v2) *)
n = Length[data];
positions = Transpose[Table[data[[i, 1]], {i, n}]];
velocities = Transpose[Table[data[[i, 2]], {i, n}]];
(Setting up stuff for NDSolve[])
{xt, yt, zt} = {ToExpression[Table["x" <> ToString[i], {i, n}]], 
   ToExpression[Table["y" <> ToString[i], {i, n}]], 
   ToExpression[Table["z" <> ToString[i], {i, n}]]};
{xm, ym, zm} = {Through[xt[t]], Through[yt[t]], Through[zt[t]]};
{xz, yz, zz} = {Through[xt[0]], Through[yt[0]], Through[zt[0]]};
yt = ToExpression[Table["y" <> ToString[i], {i, n}]];
zt = ToExpression[Table["z" <> ToString[i], {i, n}]];
rm = Flatten[
    Table[If[i != j, 
      Sqrt[(xm[[j]] - xm[[i]])^2 + (ym[[j]] - ym[[i]])^2 + (zm[[j]] - 
           zm[[i]])^2]], {i, n}, {j, n}]] /. Null -> Sequence[];
(No friction or gravity:)
xf = Thread[D[D[xm, t], t] == ConstantArray[0, n]];
yf = Thread[D[D[ym, t], t] == ConstantArray[0, n]];
zf = Thread[D[D[zm, t], t] == ConstantArray[0, n]];
(Friction and gravity:)
xf = Thread[D[D[xm, t], t] == -damping D[xm, t]];
yf = Thread[D[D[ym, t], t] == -damping D[ym, t]];
zf = Thread[
   D[D[zm, t], t] == -damping D[zm, t] - ConstantArray[gerd, n]];
(The final equations for the differential equation:)
pos = {Thread[xz == positions[[1]]], Thread[yz == positions[[2]]], 
   Thread[zz == positions[[3]]]};
vel = {Thread[D[xm, t] == velocities[[1]]] /. t -> 0, 
   Thread[D[ym, t] == velocities[[2]]] /. t -> 0, 
   Thread[D[zm, t] == velocities[[3]]] /. t -> 0};
pi = Transpose[{xm, ym, zm}]; (* position vectors)
vi = D[pi, t];( velocity vectors *)
col = Table[data[[i, 3]], {i, n}];
we = {};
(Events)
Table[AppendTo[we, 
   WhenEvent[ 
    Evaluate[(pi[[i]] - pi[[j]]) . (pi[[i]] - pi[[j]]) <= 0.2^2],
Evaluate[(e = #/Sqrt[# . #] &amp;[pi[[i]] - pi[[j]]] // Evaluate;(* 
  Stoßrichtung *)  
  V1 = (((vi[[i]] + vi[[j]] - stosszahl (vi[[i]] - vi[[j]]))/2) . 
       e ) e + vi[[i]] - (vi[[i]] . e) e ;
  V2 = (((vi[[i]] + vi[[j]] - stosszahl (vi[[j]] - vi[[i]]))/2) . 
       e ) e + vi[[j]] - (vi[[j]] . e) e ;
   Join[Thread[vi[[i]] -&gt; V1], Thread[vi[[j]] -&gt; V2]] )] ]], {i, 

1, n - 1}, {j, i + 1, n}];
Table[AppendTo[we, 
   WhenEvent[Abs[xm[[i]]] >= 0.9 // Evaluate, 
    D[xm[[i]], t] -> -D[xm[[i]], t] // Evaluate]], {i, n}];
Table[AppendTo[we, 
   WhenEvent[Abs[ym[[i]]] >= 0.9 // Evaluate, 
    D[ym[[i]], t] -> -D[ym[[i]], t] // Evaluate]], {i, n}];
Table[AppendTo[we, 
   WhenEvent[Abs[zm[[i]]] >= 0.9 // Evaluate, 
    D[zm[[i]], t] -> -D[zm[[i]], t] // Evaluate]], {i, n}];
s = NDSolve[Flatten[{xf, yf, zf, vel, pos, we}], 
   Flatten[{xt, yt, zt}], {t, 0, 50}, MaxSteps -> [Infinity]];
Manipulate[
 Graphics3D[
  Transpose[{col, 
    Thread[Sphere[
      Evaluate[
       Flatten[Thread[Transpose[{xm, ym, zm} /. t -> a] /. s], 1]], 
      0.1]]}], PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}], {a, 0, 50, 
  0.1}]

Result shows clustered masses.

Hope it helps!

When we set SeedRandom[1]; we can find two balls (green and blue)merge into one solid. — cvgmt, Jul 24 '22 at 14:59
@UlrichNeumann beautiful work. Is there a condition such that when the particles collide they permanently stay attached? No bouncing off of eachother? Thank you, again. — Zach, Jul 25 '22 at 22:04
@zak20 Thanks. Perhaps you can define contact forces (NDSolve, DiscreteVariables ->...) which appear after the collision. Have to think about it — Ulrich Neumann, Jul 26 '22 at 06:17
@UlrichNeumann I tried defining a discrete variable multiplied by a contact force (set to 0 initially), where upon a collision, the force is activated in my set of equations (by changing the discrete variable to 1). NDsolve doesn't really like the system of equations. Any thoughts on how specifically this could be done? I will create a separate question for this if you have an idea and will mark your answer as correct if you find something. — Zach, Jul 27 '22 at 17:22

Sticky Collisions (perfectly inelastic)

1 Answers1