Message Boards Message Boards

Fictitious (and reasonable) solutions to simultaneous 3-body collisions

Posted 1 year ago

In recent explorations how to generate randomness from physical systems consisting of massive disks bouncing around according to typical conservation laws, we've encountered difficult corner cases quite unlike the usual low density limit. The purpose of this memo is to answer two questions: What happens to typical high-school collision formula when two or more collisions happen simultaneously? Can a consistent solution be found by modifying the usual equations?

Assuming that we already know what happens in a two-body elastic collision, let's start with the simplest examples of a simultaneous three-body collision, such as follows:

vel = {
   {v1 Cos[\[Theta]1 - \[Delta]1], v1 Sin[\[Theta]1 - \[Delta]1]},
   {v2 Cos[\[Theta]2 - \[Delta]2], v2 Sin[\[Theta]2 - \[Delta]2]},
   {0, 0}
   };

pos = {
   { -Cos[\[Theta]1 ], -Sin[\[Theta]1]},
   { -Cos[\[Theta]2 ], -Sin[\[Theta]2]},
   {0, 0}
   };

Grid[Partition[Table[
   With[{randomCollision = {
       m1 -> 1, m2 -> 1, ms -> 1,
       \[Theta]1 -> 0,
       \[Theta]2 -> -RandomReal[{Pi/3, Pi}],
       \[Delta]1 -> 0,
       \[Delta]2 -> 0,
       v1 -> 1, v2 -> 1}},
    Graphics[MapThread[{White, EdgeForm[Black],
        Disk[#1, 1/2],
        Gray, Arrow[{#1 - #2, #1}]
        } &, {pos /. randomCollision, vel /. randomCollision}],
     PlotRange -> {{-2, 2}, {-1, 2}}
     ]], {12}], 3]]

simple collisions

A frame of reference has been chosen with the center mass fixed. The impinging masses have equal velocities, which are both directed toward the centroid of the fixed mass. A common first idea is: "Well, what if applying collisions in sequence doesn't produce a noticeable difference?" It's a bit hopeful, though, not entirely wrong. We test this idea by defining a commutator for order of operations as follows:

MomentumDifferential[dx_, v1In_, v2In_, m1In_, m2In_
  ] := Times[Divide[2 m1In*m2In, m1In + m2In],
  Dot[v1In - v2In, dx]*dx]

TwoPartCollisionSymmetryBroken[{dx1_, v1In_}, {dx2_, v2In_}
  ] := Module[{dp1s, dp2s, mag, v1Out, v2Out, vsOut},
  dp1s = MomentumDifferential[dx1, v1In, {0, 0}, m1, ms];
  dp2s = MomentumDifferential[dx2, v2In, dp1s/ms, m2, ms];
  v1Out = v1In - dp1s/m1;
  v2Out = v2In - dp2s/m2;
  vsOut = dp1s/ms + dp2s/ms;
  {v1Out, v2Out, vsOut}
  ]

v1Solution = TwoPartCollisionSymmetryBroken[{
    { -Cos[\[Theta]1 ], -Sin[\[Theta]1]},
    {v1 Cos[\[Theta]1 - \[Delta]1], v1 Sin[\[Theta]1 - \[Delta]1]}
    }, {
    { -Cos[\[Theta]2 ], -Sin[\[Theta]2]},
    {v2 Cos[\[Theta]2 - \[Delta]2], v2 Sin[\[Theta]2 - \[Delta]2]}
    }];

v2Solution = TwoPartCollisionSymmetryBroken[{
     { -Cos[\[Theta]2 ], -Sin[\[Theta]2]},
     {v2 Cos[\[Theta]2 - \[Delta]2], v2 Sin[\[Theta]2 - \[Delta]2]}
     }, {
     { -Cos[\[Theta]1 ], -Sin[\[Theta]1]},
     {v1 Cos[\[Theta]1 - \[Delta]1], v1 Sin[\[Theta]1 - \[Delta]1]}
     }][[{2, 1, 3}]];

commutator = ReplaceAll[Refine[Sqrt[Simplify@ReplaceAll[
      # . # & /@ (v1Solution - v2Solution),
      {m1 -> 1, m2 -> 1, ms -> 1, \[Delta]1 -> 0, \[Delta]2 -> 0,
       \[Theta]1 -> \[Theta] + \[Theta]2, v1 -> 1, 
       v2 -> 1}]], \[Theta] \[Element] Reals],
  Abs -> Identity]

Out[] = {Cos[\[Theta]], Cos[\[Theta]], 2 Cos[\[Theta]] Sin[\[Theta]/2]}

Plot[commutator, {\[Theta], 0, 2 Pi}, PlotRange -> All]

commutator

The resulting commutators for each output velocity vector, $\{\cos[\theta], \cos[\theta], 2\cos[\theta] \sin[\theta/2 \}$, have all zeroes at $\theta = \pm \pi/2$ where the input velocity vectors decompose along orthogonal dimensions. Otherwise there is generally a difference when colliding one particle first and the other second (as compared to the vice versa colliding the other first and the one second). To get a visual sense of this, we can plot alternative output vectors in Red and Blue:

With[{simpleCollision = {m1 -> 1, m2 -> 1, 
    ms -> 1, \[Delta]1 -> 0, \[Delta]2 -> 0, v1 -> 1, v2 -> 1}},
 With[{
   vAsym1 = 
    Expand[v1Solution /. 
       simpleCollision /. {\[Theta]1 -> -Pi/5, \[Theta]2 -> Pi/5}],
   vAsym2 = 
    Expand[v2Solution /. 
       simpleCollision /. {\[Theta]1 -> -Pi/5, \[Theta]2 -> Pi/5}]
   },
  Graphics[MapThread[{White, EdgeForm[Black], Disk[#1, 1/2],
      Red, Arrow[{#1, #1 + #2}],
      Blue, Arrow[{#1, #1 + #3}]
      } &, {pos /. 
       simpleCollision /. {\[Theta]1 -> -Pi/5, \[Theta]2 -> Pi/5},
     vAsym1, vAsym2}]]
  ]]

asym collision

Of course, whichever collision happens first completely stops one disk, and then subsequently reduces the velocity of the other. Considering symmetry, this doesn't seem right at all, and we don't already have a way to continuously vary between solutions to find the mean we're looking for. What can we do?

Let's assume that disks exchange momenta in increments $\pm dp_{1s}$ and $\pm dp_{2s}$ with $s$ standing for "initially still". The output momenta can then be written as: $$p_1 ' = p_1 - k\cdot dp_{1s},$$ $$p_2 ' = p_2 - k\cdot dp_{2s}, $$ $$p_s ' = k\cdot( dp_{1s}+dp_{2s}). $$

Where the constant $k$ can be determined from conservation of energy. By design the sum of momenta is conservative, and because the differentials point toward or away from the center of the still disk, their addition also conserves angular momentum. So this is a valid solution with regard to fundamental principles, but also one where an arbitrary scaling choice has been made.

TwoPartCollisionSymmetric[{dx1_, v1In_}, {dx2_, v2In_}
  ] := Module[{
   dp1s = scale*MomentumDifferential[dx1, v1In, {0, 0}, m1, ms],
   dp2s = scale*MomentumDifferential[dx2, v2In, {0, 0}, m2, ms],
   mag, v1Out, v2Out, vsOut},
  v1Out = v1In - dp1s/m1;
  v2Out = v2In - dp2s/m2;
  vsOut = dp1s/ms + dp2s/ms;
  mag = First[DeleteCases[
     scale /. 
       Solve[Simplify@Subtract[m1 v1In . v1In + m2 v2In . v2In,
           Dot[Dot[#, #] & /@ {m1 v1Out, m2 v2Out, ms vsOut},
            {1/m1, 1/m2, 1/ms}]] == 0, scale
        ] /. Abs -> Identity, 0]];
  {v1Out, v2Out, vsOut} /. {
    scale -> mag, Abs -> Identity
    }]

genVecs1 = Simplify@Expand[
    TwoPartCollisionSymmetric[{
      { -Cos[\[Theta]1 ], -Sin[\[Theta]1]},
      {v1 Cos[\[Theta]1 - \[Delta]1], v1 Sin[\[Theta]1 - \[Delta]1]}
      }, {
      { -Cos[\[Theta]2 ], -Sin[\[Theta]2]},
      {v2 Cos[\[Theta]2 - \[Delta]2], v2 Sin[\[Theta]2 - \[Delta]2]}
      }]];

And we can see by updating the collision plot that this technique produces a much better result:

With[{simpleCollision = {m1 -> 1, m2 -> 1, 
    ms -> 1, \[Delta]1 -> 0, \[Delta]2 -> 0, v1 -> 1, v2 -> 1}},
 With[{
   vAsym1 = 
    Expand[v1Solution /. 
       simpleCollision /. {\[Theta]1 -> -Pi/5, \[Theta]2 -> Pi/5}],
   vAsym2 = 
    Expand[v2Solution /. 
       simpleCollision /. {\[Theta]1 -> -Pi/5, \[Theta]2 -> Pi/5}],
   vSym = 
    Expand[genVecs1 /. 
       simpleCollision /. {\[Theta]1 -> -Pi/5, \[Theta]2 -> Pi/5}]
   },
  Graphics[MapThread[{White, EdgeForm[Black], Disk[#1, 1/2],
      Red, Arrow[{#1, #1 + #2}],
      Blue, Arrow[{#1, #1 + #3}],
      Black, Arrow[{#1, #1 + #4}]
      } &, {pos /. 
       simpleCollision /. {\[Theta]1 -> -Pi/5, \[Theta]2 -> Pi/5},
     vAsym1, vAsym2, vSym}]]
  ]]

symmetric average

We also get seemingly acceptable results if we allow disks to collide with skew vectors:

GraphicsGrid@Partition[Table[With[{randomCollision = {
       m1 -> 1, m2 -> 1, ms -> 1,
       \[Theta]1 -> 0,
       \[Theta]2 -> -RandomReal[{Pi/3, Pi}],
       \[Delta]1 -> RandomReal[{-Pi/4, Pi/4}],
       \[Delta]2 -> RandomReal[{-Pi/4, Pi/4}], v1 -> 1, v2 -> 1}},
    With[{
      vAsym1 = Expand[v1Solution /. randomCollision],
      vAsym2 = Expand[v2Solution /. randomCollision],
      vSym = Expand[genVecs1 /. randomCollision]},
     Graphics[MapThread[{White, EdgeForm[Black],
         Disk[#1, 1/2],
         Gray, Arrow[{#1 - #2, #1}],
         Red, Arrow[{#1, #1 + #3}],
         Blue, Arrow[{#1, #1 + #4}],
         Black, Arrow[{#1, #1 + #5}]
         } &, {pos /. randomCollision, vel /. randomCollision,
        vAsym1, vAsym2, vSym}],
      PlotRange -> {{-2, 2}, {-1, 2}}
      ]]], {16}], 4]

more results

Unfortunately, if we try to add a fourth body, it quickly becomes easy to accidentally violate conservation of angular momentum. So more work remains to be done finding an acceptable, exactly solve-able, semi-realistic convention for dealing with corner cases of hard circle collision simulations.

POSTED BY: Brad Klee

enter image description here -- you have earned Featured Contributor Badge enter image description here Your exceptional post has been selected for our editorial column Staff Picks http://wolfr.am/StaffPicks and Your Profile is now distinguished by a Featured Contributor Badge and is displayed on the Featured Contributor Board. Thank you!

POSTED BY: Moderation Team
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