What are the general solutions to a hard sphere collision? [duplicate]

Question

Surely someone has found the solutions to the hard sphere collisions (in $n$ dimensions) of two bodies of mass $m_1$ and $m_2$, respectively--that is the resultant velocities (or momenta) of the two bodies after a hard sphere collision given $\vec{p}_1$ and $\vec{p}_2$. I have a solution (both for elastic and inelastic collisions), but I wonder whether or not it's correct.

FYI, here are my solutions for elastic, hard-sphere collisions:

$$ \vec{p}_{1f} = \vec{p}_1 + 2\frac{\left(m_1\vec{p}_2 - m_2\vec{p}_1\right)\circ \hat{r}}{m_1 + m_2}\hat{r}, \\ \vec{p}_{2f} = \vec{p}_2 - 2\frac{\left(m_1\vec{p}_2 - m_2\vec{p}_1\right)\circ \hat{r}}{m_1 + m_2}\hat{r}. $$

wher $\hat{r}$ is in the direction of the "hard sphere" (i.e. the unit vector connecting the two centers of the two spheres). You can see that my formulation trivially conserved momentum--the question is whether or not it conserves energy as well.

If it helps, here is my full explanation. For these purposes, you should skip to the section on n-dimensional cases (page 6). I wrote this as a reference to myself for doing physics simulations since I could not find any reference on hard sphere collisions (other than references setting up the problem which give no solutions).

Here is an abbreviated work on finding the final momenta from the elastic case:

Just as in the $1$D case, we add and subtract the same momentum from each particle--except now we are adding and subtracting a vector: $\vec{p} = p\hat{r}$, which means that $|\vec{p}| = p$ (such that $\vec{p}_{1f} = \vec{p}_1 + p\hat{r}$ and $\vec{p}_{2f} = \vec{p}_2 - p\hat{r}$). This gives the final energy as:

\begin{align} E_f = \frac{|\vec{p}_1 + p\hat{r}|^2}{2m_1} + \frac{|\vec{p_2} - p\hat{r}|^2} {2m_2} \end{align}

We now write the final energy in terms of the magnitude of $p$ and the dot product between $\hat{r}$ and $\vec{p}_1$ and likewise for $\vec{p}_2$. Note that if you choose the direction for $\vec{p}$ to act a priori (like in the hard-sphere case), the dot products $\vec{p}_1\cdot \hat{r}$ and $\vec{p}_2\cdot\hat{r}$ will be known:

\begin{align*} E_f = \frac{p_1^2 + p^2 + 2p\hat{r}\circ\vec{p}_1}{2m_1} + \frac{p_2^2 + p^2 - 2p\hat{r}\circ\vec{p}_2}{2m_2} \end{align*}

As usual the elastic case greatly simplifies:

\begin{align*} \require{cancel} \cancel{\frac{p_1^2}{2m_1}} + \cancel{\frac{p_2^2}{2m_1}}= \frac{\cancel{p_1^2} + p^2 + 2p\hat{r}\circ\vec{p}_1}{2m_1} + \frac{\cancel{p_2^2} + p^2 - 2p\hat{r}\circ\vec{p}_2}{2m_2} \\ \frac{m_1 + m_2}{2m_1m_2}p^2 - 2\frac{m_1 \vec{p}_2\circ\hat{r} - m_2\vec{p}_1\circ\hat{r}}{2m_1m_2}p = 0 \end{align*}

This leads to a trivial solution of $p = 0$ (no collision occurs) and:

\begin{align} \frac{m_1 + m_2}{2m_1m_2}p - 2\frac{m_1 \vec{p}_2\circ\hat{r} - m_2\vec{p}_1\circ\hat{r}}{2m_1m_2} = 0 \\ p = 2\frac{m_1 \vec{p}_2\circ\hat{r} - m_2\vec{p}_1\circ\hat{r}}{m_1 + m_2} \end{align}

Hence my solutions of:

$$ \vec{p}_{1f} = \vec{p}_1 + 2\frac{\left(m_1\vec{p}_2 - m_2\vec{p}_1\right)\circ \hat{r}}{m_1 + m_2}\hat{r}, \\ \vec{p}_{2f} = \vec{p}_2 - 2\frac{\left(m_1\vec{p}_2 - m_2\vec{p}_1\right)\circ \hat{r}}{m_1 + m_2}\hat{r}. $$

Answer 1

Finding the missing equations

Coordinate transforms just complicate the issue. The heart of the matter is that in n dimensions you have n degrees of freedom for the velocities of the COM of each sphere, and you only have n momentum conservation equations plus one energy conservation equation. That means you need an additional n-1 equations to solve the problem.

The additional equations come from the fact that if you assume a frictionless impact*, then the impulse/force imparted to each object must be normal to the surface. This gives a constraint on the direction of the impulse. Along with the magnitude, this direction would give us n degrees of freedom. So without it, it gives us n-1: exactly the number we need.

Solving for final velocities

Since we have spheres, the direction normal to the impact surface will be the vector that points from one COM to the other. Let's call this normal vector $\hat{j}$.

Now the velocities of balls A and B can be given as $V_{a1}$ and $V_{b1}$ before the collision and $V_{a2}$ and $V_{b2}$ after the collision.

Let's call the impulse that is imparted on ball A $J$. Then we have:

$$J \bullet\hat{j} =|J|$$ $$ V_{a1} + \frac{J}{M_a}=V_{a2}$$ $$ V_{b1} - \frac{J}{M_b}=V_{b2}$$

Then we can write the conservation of energy (multiplied by 2) as: $$M_a V_{a1} \bullet V_{a1} + M_b V_{b1} \bullet V_{b1} = M_a V_{a2} \bullet V_{a2} + M_b V_{b2} \bullet V_{b2}$$

Plugging in for our velocities after to collision gives:

$$M_a V_{a1} \bullet V_{a1} + M_b V_{b1} \bullet V_{b1} = M_a \left( V_{a1} + \frac{J}{M_a}\right) \bullet \left( V_{a1} + \frac{J}{M_a}\right) + M_b \left( V_{b1} - \frac{J}{M_b}\right) \bullet \left( V_{b1} - \frac{J}{M_b}\right)$$

$$M_a V_{a1} \bullet V_{a1} + M_b V_{b1} \bullet V_{b1} = M_a V_{a1}\bullet V_{a1} + 2 J\bullet V_{a1} + \frac{J\bullet J}{M_a} + M_b V_{b1}\bullet V_{b1} - 2 J\bullet V_{b1} + \frac{J\bullet J}{M_b}$$

$$0 = 2 J\bullet V_{a1} + \frac{J\bullet J}{M_a} - 2 J\bullet V_{b1} + \frac{J\bullet J}{M_b}$$

dividing by the magnitude of $J$ yields:

$$0 = 2 \hat{j} \bullet (V_{a1}-V_{b1}) + |J|\left(\frac1{M_a} + \frac1{M_b}\right)$$

$$|J| = 2 \frac{\hat{j} \bullet (V_{a1}-V_{b1})}{\frac1{M_a} + \frac1{M_b}}$$

Final Velocities

$$V_{a2}=V_{a1}+2 \frac{\hat{j} \bullet (V_{a1}-V_{b1})}{M_a\left(\frac1{M_a} + \frac1{M_b}\right)}$$ $$V_{b2}=V_{b1}+2 \frac{\hat{j} \bullet (V_{b1}-V_{a1})}{M_b\left(\frac1{M_a} + \frac1{M_b}\right)}$$

Finding Direction of Impulse

Let's call the vector pointing from the center of ball B to the center of ball A $R$. We know that before the impact $R=R_0+R't$ where $R_0$ is the difference in the balls' initial positions $P_a-P_b$ and $R'$ is the difference in the balls' initial velocities $V_{a1}-V_{b1}$.

At the moment of impact: $$\hat{j}=\frac{R}{|R|}$$ and $$|R|=r_a+r_b$$ Where $r_a$ and $r_b$ are the radii of the balls.

We can transform this second equation to find t:

$$R\bullet R=(r_a+r_b)^2$$ $$R'\bullet R't^2+2R_0\bullet R't+R_0\bullet R_0-(r_a+r_b)^2=0$$ $$t=-\frac{R_0\bullet R'+\sqrt{(R_0\bullet R')^2-R'\bullet R'(R_0\bullet R_0-(r_a+r_b)^2)}}{R'\bullet R'}$$

If this time is negative or imaginary, then the collision will not occur. Otherwise $R$ at the moment of impact can be found by plugging in t:

$$R=R_0-R'\frac{R_0\bullet R'+\sqrt{(R_0\bullet R')^2-R'\bullet R'(R_0\bullet R_0-(r_a+r_b)^2)}}{R'\bullet R'}$$

Wrap Up

So given initial positions $P_a$ and $P_b$, initial velocities $V_{a1}$ and $V_{b1}$, and ball radii $r_a$ and $r_b$ the final velocities $V_{a2}$ and $V_{b2}$ can be found be following the following procedure:

Calculate the relative positions and velocities: $$R_0=P_a-P_b$$ $$R'=V_{a1}-V_{b1}$$

Calculate impact time: $$t=-\frac{R_0\bullet R'+\sqrt{(R_0\bullet R')^2-R'\bullet R'(R_0\bullet R_0-(r_a+r_b)^2)}}{R'\bullet R'}$$ If t is negative or imaginary, then no collision occurs, and the initial velocities are the final velocities. Otherwise:

Calculate the impulse direction: $$\hat{j}=\frac{R_0+R't}{|R_0+R't|}$$ Calculate the impulse: $$J=2 \hat{j}\frac{\hat{j} \bullet R'}{\frac1{M_a} + \frac1{M_b}}$$ Calculate the final velocities: $$V_{a2} = V_{a1} + \frac{J}{M_a}$$ $$V_{b2} = V_{b1} - \frac{J}{M_b}$$

*The collision must be frictionless if we're not considering moments of inertia as any friction would apply torque to the ball which would transfer energy and rotational inertia to the ball. The problem becomes slightly more complicated in that case as one has to assume how the friction will cause the balls to interact, which is not trivial as many of the simple assumptions conflict with energy conservation.

Answer 2

Here is a general figure of an hard spheres collision drawn in the center of mass of the mass $m_2$ before the collision. The black dot is attached to this frame.

To solve the problem, you need to observe

• Conservation of energy: $m_1v_1^2=m_1(v'_1)^2+m_2(v'_2)^2$.
• Conservation of momentum: $m_1\vec v_1=m_1\vec v'_1+m_2\vec v'_2$
• Conservation of torque with respect to the center of mass of $m_2$: $m_1v_1b=m_1v'_1b'$.

I have drawn a 2D problem, so you have four unknowns (two for $\vec v'_1$ and two for $\vec v'_2$). And there are four equations (one for energy and torque, two for momentum). So the problem is solved.

In three dimensions, you have $6$ unkowns, energy conservation provides one equation, momentum conservation provides $3$ equations and torque conservation provides two.

Answer 3

This thread is a little old, but I thought I would chip in in case anyone comes looking for this answer. In short, Tom-Tom is exactly correct in his analysis of the problem. However, I am going to break this down in steps for a full 3D case and try to explain some of the points there was contention about (e.g. angular momentum).

First of all, here are the assumptions we are working with:

• Hard sphere collision is a technical term meaning the following: 1) collisions are elastic (which means kinetic energy is conserved); 2) the spheres have a finite radius (as someone else has pointed out, if they don't this is a point-like collision.)
• While one could work with a simulation in which we accounted for rigid body rotation, I will assume that the "spin" is zero. Spin adds a good deal of complexity to the computation.

I'm not going to actually solve each of the terms (it gets pretty messy), but I am going to outline how the process goes. For completeness, I'm not going to assume that we know the bodies will collide. We shall see as part of the analysis that this determination falls out pretty easily.

Here is the procedure:

1) First move to the rest frame of one of the particles. Now we are comparing a point and a 3D line, which greatly simplifies the geometry, rather than two 3D lines.

2) Recall the fundamental 3D geometry proof that a line and a point in 3D space can always be contained in a single plane. We can, in principle, do a coordinate transformation (rotation in $\theta$ and $\phi$ in spherical coordinates around the particle at the origin).

3) Good, now we have a line and a point in a single plane. It should be obvious that another rotation (in polar angle) can ensure we have the moving particle moving the x-axis only. This reproduces a picture like the first figure Tom-Tom posted. Now we have reduced our 3D problem to the 2D problem Tom-Tom gave a method for!

4) Collision detection: how do we know if they hit? We need to check 2 things: 1) is $R_1 + R_2 > b$, where $R_i$ are the radii and $b$ is the impact parameter; 2) we need to check that the moving particle gets within $\sqrt{(R_1 + R_2)^2 - b^2}$ of the x-axis. That's just the Pythagorean theorem for the triangle the spheres make at contact. If both criteria are met, we have a collision.

5) Now we're all setup, and we work from Tom-Tom's solution. I won't work the rest out, but I will make a couple of notes to clear up confusion.

• Angular momentum: Contrary to Jared's comment above, we cannot in general pick a coordinate system where the angular momentum is zero. Angular momentum will only be zero in two special cases: 1) the particles are at rest with respect to each other; 2) the moving particle is headed directly at the stationary particle. Why? $\vec{L} = \vec{r}\times\vec{p}$. If the particles are at rest, $\vec{p} = 0 \to L = 0$ (obviously we don't have a collision in this case either). If the impact parameter is zero, that means $\vec{r}$ and $\vec{p}$ are lined up with each other. Recall, $|\vec{L}| = |\vec{r}||\vec{p}|Sin(\theta)$, where $\theta$ is the angle between the two vectors. In the $b = 0$ case, $\theta = 0 \to |\vec{L}| = 0$. In all other cases, angular momentum will be non-zero, no matter what coordinate system you choose!

• The radii are important! They will enter into the angular momentum equation. To solve the angular momentum equation, first note that $\vec{r_{2}}' = \vec{v_2}' t$, therefore the angle between momentum and the radial vector will be zero giving the angular momentum of the second body as zero (as Tom-Tom's equation has indicated). They are not parallel for the first sphere. Staring at the geometry for a few minutes should convince you, rather, that $\frac{\vec{v_1}'}{t} = \vec{r_1}' + \vec{\alpha}$, where $\alpha$ is the vector pointing from the center of body 1 to the center of body 2. This vector has magnitude $R_1 + R_2$. We can then use the Pythagorean theorem to break $\alpha$ into components (note: the y component is just b). The x component will be $\sqrt{(R_1 + R_2)^2 - b^2}$. We use the vector addition above to eliminate $\vec{r_i}'$ for $\vec{v_i}'$. Note: t is arbitrary and can be set to 1. Finally, carry out the angular momentum cross product. You will see that the $\sqrt{(R_1 + R_2)^2 - b^2}$ term survives, and will therefore enter into your final equations!

From this point, it's just messy algebra. We have 4 unknowns and 4 equations. Solve those equations, and you'll have the velocities. A small amount of extra work can tell you, e.g. exactly when the collision occurs, which is important for simulations and will tell you how far to translate the spheres with the new velocity after collision (to match simulation timestep).

The very last step, after you've done all that, is to remember to transform back to the system coordinates. To do this, just unroll everything backwards, flipping the signs of angles and translations.

Hope that clears things up for someone that's looking at this problem. This is a classic one of those problems that seems like it should be trivially easy, but turns out to be quite involved to solve in the general case.

UPDATE

So, this all got me thinking and I went back to re-examine the problem a bit. In short, the problem turns out to be much more subtle than I originally thought. Here are a couple of new comments:

• The reason you need angular momentum is that angular momentum must be conserved for the simulation to be physical. If you break angular momentum, you break rotational invariance of your system. If you're making a simple game engine, that might be OK, but if you're trying to do physics, that's definitely not allowed. Angular momentum conservation is what makes the radii important.

• I need to retract what I said about spin. Spin is necessary (as I will show below). However, Newtonian gravity cannot be the source of spin. To introduce spin requires having net torque. Recall that for any (conservative) potential of the form $V = V(r)$, the force due to a sphere is equivalent to the force due to a point at its center. Briefly, any torque that would be created on one side of the sphere is canceled by the other side, and only force along the radial line survives. Newtonian gravity is clearly a potential of this form (as is, e.g. Coulomb, Yukakawa, and almost any "real" potential we deal with.)

• So, I worked out a solution that takes into account linear momentum, angular momentum and kinetic energy. I'll only put the x-component of the velocity vector here (the others are simply related by plugging this back into the equations:

$v_{xf} = \frac{\left(\beta ^2-1\right) \sqrt{\frac{v_i^2 \left(\beta ^2 (\mu +2)^2-\mu ^2\right)}{\beta ^2-1}}+v_i \left(\beta ^2 (\mu +2)-1\right)}{\mu +1}$

If you're curious, the following code will get this result in Mathematica:

(* Setup vectors: No Q *)
vi = {v1x, 0, 0};
v1 = {v1xf, v1yf, 0};
v2 = {v2xf, v2yf, 0};
\[Alpha] = {-Sqrt[(R1 + R2)^2 - b^2], b, 0}; (* The vector that connects the two centers at impact *)

r1 = \[Alpha] + 
   v1 \[CapitalDelta]t; (* NOTE: \[CapitalDelta]t should drop out! *)

r2 = v2 \[CapitalDelta]t;

(* Physics *)
Eqns0 = Simplify[{
    m1 (vi.vi) == m1 (v1.v1) + m2 (v2.v2)(*+Q m1(vi.vi)*),
    m1 v1x == m1 v1xf + m2 v2xf,
    0 == m1 v1yf + m2 v2yf,
    m1 v1x b == (m1 (r1 \[Cross] v1) + m2 (r2 \[Cross] v2))[[3]]
    } /. {m1 -> m, m2 -> \[Mu] m, b -> \[Beta] (R1 + R2)},    
\[Beta] < 1 && \[Beta] >= 0 && m > 0 && \[Mu] > 0 && R1 > 0 && 
   R2]; Eqns0 // MatrixForm

Eqns1 = Eqns0 /. Solve[Eqns0[[2]], v2xf] // 
    Flatten; Eqns1 // MatrixForm
Eqns2 = Eqns1 /. Solve[Eqns1[[3]], v2yf] // 
  Flatten; Eqns2 // MatrixForm
(Eqns2[[1]] /. Solve[Eqns2[[4]], v1yf] // Flatten)[[1]]
FullSimplify[
  Solve[%, v1xf], \[Beta] < 1 && \[Beta] >= 0 && m > 0 && \[Mu] > 0 &&
    R1 > 0 && R2] // Flatten 

I have put everything in dimensionless units to make the result more transparent. Specifically, $\mu$ is the ratio $m_1/m_2$, $\beta$ is the ratio $b/(R_1+R_2)$, where $b$ is the impact parameter. Note that the dimensionless parameters allow you to quickly verify that all of the terms on the right hand side do in fact have dimensions of velocity.

But wait! This is where things get interesting. Note that the term in the radical can become negative for real situations, which leads to complex velocities. This can't be allowed. We can plot $\beta$ vs $\mu$ to show the "allowed" region:

This is clearly a problem. It is certainly physically possible for $\beta$ to take any value between 0 and 1, but we get imaginary velocities if $\beta$ is too large compared to $\mu$.

In short, the resolution is that (classical) hard sphere scattering cannot be elastic in general (though it can be in special cases). I'll explain in a moment why that is the case, but first, I'll give a solution that "works." Here, I introduce another dimensionless parameter $Q$, which is the ratio of energy loss over initial kinetic energy. The modified solution is:

$v_{xf} = \frac{\left(\beta ^2-1\right) \sqrt{\frac{v_1^2 \left(\beta ^2 (\mu +2)^2+\mu (\mu (Q-1)+Q)\right)}{\beta ^2-1}}-v_i \left(\beta ^2 (\mu +2)-1\right)}{\mu +1}$

We can now plot $Q$ vs. $\beta$ for varying values of $\mu$:

It's important to note that this only provides a bound on $Q$. It does not tell us what value of $Q$ we should use. In other words, the problem is no longer fully constrained. So, what's going on here and how can we resolve it?

Here's the bottom line: hard sphere scattering is not physical! To have a true hard sphere (i.e. impenetrable) we would require the force at the surface of the spheres to be infinite, i.e. the force law would be a Dirac delta function. The impulse experienced, then, would be: $\Delta p = F \Delta t = (\infty)(0)$. This is undefined. It’s not surprising, then that we might get some strange physics out of hard sphere scattering. Moreover, it’s not surprising the problem would be unconstrained.

We can say a little more though. Up till now, I assumed that we could use the old standby: $\vec{F} = - \nabla V(r)$. But, this rule is only valid for conservative forces. Recall that a conservative force is a path independent force. In other words, if you start at point $A$ and take any path you like returning to point $A$ the potential remains the same. This is clearly not the case for hard sphere scattering. Any path that doesn’t involve a collision is going to be dramatically different than any path that does. In fact, this is a general rule for any truly localized force (note: Yukawa potential is a “localized” force, but in reality, it just has exponential fall-off, so it’s not a localized force in this sense).

In the case of a non-conservative force, we can no longer use $\vec{F} = - \nabla V(r)$. Rather, we must use the Helmholtz decomposition, which states that any force vector field can be decomposed as: $\vec{F} = - \nabla V(r) + \nabla \times \vec{A}$. That is, we must introduce a vector potential.

This gets us back to the notion of spin. A quick calculation will verify that the vector potential will generically introduce a net torque in the scattering process, which will produce spin (of course, for very special $\vec{A}$ this will not happen).

The real problem now is that $\vec{A}$ is completely unconstrained by the problem. It is associated with $Q$, and therefore the internal degrees of freedom of the sphere. But, we have assumed from the beginning that the sphere is a “black box” that we don’t specify physics for. This is an unphysical assumption, and it turns out that it matters. The obvious internal degrees of freedom to account for are the spin vector and the temperature (it’s a fairly natural assumption that collisions between real bodies would effect the heat content of the bodies). The argument above suggests that at least some of that energy should go into spin energy, but probably not all of it. Just how much should go into each of these degrees of freedom depends on the specifics of internal physics.

I think that’s about as far as we can go with the problem. If we want to fully specify the problem we need to create something like an “equation of state” (by analogy) for the spheres which determines how the energy gets distributed in bodies. A reasonable approximation might be to just specify a ratio of temperature to spin, though physical intuition suggests that higher impact parameters should transfer more spin and less temperature, so perhaps this should be a function of $\beta$. A suggestion off the top of my head would be to assume spin as power series in $\beta$, where the series constants could be experimentally determined or estimated by more detailed theoretical models of the internal physics of the bodies.

If you want a physical model of scattering, you need to specify the physics of the bodies. One clean solution that makes sense for a lot of cases (e.g. actual macroscopic rigid bodies like pool balls) would be to use Yukawa scattering. The Yukawa potential is conservative, and was specifically designed to approximate the effects of electrodynamic screening that occurs in macroscopic, (mostly) electrically neutral bodies.