0
$\begingroup$

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}]```
$\endgroup$
6
  • $\begingroup$ Should the particles alos stuck at contact with the boundary? $\endgroup$ Jul 22 at 6:23
  • $\begingroup$ Radius of the bodies is .15? $\endgroup$ Jul 22 at 8:24
  • $\begingroup$ They should bounce off all boundaries, but can be stopped at the bottom. Radius of sphere looks like 1 $\endgroup$
    – Zach
    Jul 22 at 14:52
  • 1
    $\begingroup$ 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. $\endgroup$ Jul 22 at 15:43
  • 1
    $\begingroup$ For me the main problem of the simulation is the contur of the new clumped masses. Have to think about ... $\endgroup$ Jul 23 at 8:13

1 Answer 1

2
$\begingroup$

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[# . #] &[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]] -> V1], Thread[vi[[j]] -> 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.

enter image description here

Hope it helps!

$\endgroup$
7
  • 2
    $\begingroup$ When we set SeedRandom[1]; we can find two balls (green and blue)merge into one solid. $\endgroup$
    – cvgmt
    Jul 24 at 14:59
  • $\begingroup$ @cvgmt Thanks, initial position data has to be improved $\endgroup$ Jul 24 at 16:01
  • $\begingroup$ @UlrichNeumann Das ist die Antwort. Vielen Dank!!! $\endgroup$
    – Zach
    Jul 25 at 18:58
  • $\begingroup$ @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. $\endgroup$
    – Zach
    Jul 25 at 22:04
  • $\begingroup$ @zak20 Thanks. Perhaps you can define contact forces (NDSolve, DiscreteVariables ->...) which appear after the collision. Have to think about it $\endgroup$ Jul 26 at 6:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.