Pair collisions

We can find the collision rule between two particles by considering the conservation of momentum, and conservation of energy. We take all particles masses, $m=1$.

Momentum conservation gives $ \textbf{v}_1+ \textbf{v}_2 = \textbf{v'}_1+ \textbf{v'}_2 $, for the collision between two particles $1$ and $2$ of equal mass. Energy conservation implies that the kinetic energy $K= \frac{\textbf{v}_1^2}{2} + \frac{\textbf{v}_2^2}{2}$, is also the same before and after collision.

This implies $$\textbf{v}_1' =\textbf{v}_1 -\hat{ \textbf{r}} (\hat {\textbf{ r}} \cdot ( \textbf{v}_1 -\textbf{ v}_2) )$$, $$\textbf{v}_2' =\textbf{v}_2 +\hat{ \textbf{r}} (\hat {\textbf{ r}} \cdot ( \textbf{v}_1 -\textbf{ v}_2) )$$ where $\hat{\textbf{r}}$ is a unit vector in the direction of the vector joining the centres of the particles.

When coding the collisions it is a good crosscheck on your code to verify these two conservation laws. Python has two function for this. isclose() to check if two floating numbers are close. assert stops the program if the comparison fails. We are working with floating point numbers. You should never test for exact equality. Round-off errors are always present in any calculation.

import math
a=1+ 1.e-12
b=1
print("a b\t",a,b)

# this compares a,b to within 1.e-9. they are very nearly the same
# nothing is printed
assert math.isclose(a,b)
# this makes a and b different by a larger amount, the assert will fail and print and error
a=1+1.e-8
print("a b\t",a,b)
print("this next line will give an error\n")
assert math.isclose(a,b)

a b	 1.000000000001 1
a b	 1.00000001 1
this next line will give an error

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[1], line 13
     11 print("a b\t",a,b)
     12 print("this next line will give an error\n")
---> 13 assert math.isclose(a,b)

AssertionError: 

Useful vector expressions in numpy

Calculating the norm of all vectors at same time

import numpy as np

position = np.array([[0.25, 0.25], [0.75, 0.25], [0.25, 0.75], [0.75, 0.75]])
i, j = np.triu_indices(position.shape[0], k=1)
rij = position[i]-position[j]  # set of all 6 separation vectors
print('rij=\n',rij, '\n')
rij_sq = (rij**2).sum(1)
print("squared separations", rij_sq)
C=np.array([20, -30])

print('rij-C=\n',rij-C, '\n')

print ('rij*C', rij*C, '\n')

rij=
 [[-0.5  0. ]
 [ 0.  -0.5]
 [-0.5 -0.5]
 [ 0.5 -0.5]
 [ 0.  -0.5]
 [-0.5  0. ]] 

squared separations [0.25 0.25 0.5  0.5  0.25 0.25]
rij-C=
 [[-20.5  30. ]
 [-20.   29.5]
 [-20.5  29.5]
 [-19.5  29.5]
 [-20.   29.5]
 [-20.5  30. ]] 

rij*C [[-10.  -0.]
 [  0.  15.]
 [-10.  15.]
 [ 10.  15.]
 [  0.  15.]
 [-10.  -0.]] 

Collision times

Particles are moving at constant speed so, for instance the $x$-coordinate of a particle is $$x=x_0 + v_x t$$ as a function of time. Show that if the particles have diameter $\sigma$ the contact time is the solution to a quadratic equation.

Collisions between particles with different masses (ignore for first reading)

Changes in velocity are along the radial vector $\hat {\bf r}$:

$$\Delta {\bf v}_1 = a \hat {\bf r}$$$$\Delta {\bf v}_2 = b \hat {\bf r}$$

Momentum conservation implies: $$m_1 a + m_2 b =0 $$ so we can write $$a = \alpha m_2 $$ and $$b=-\alpha m_1 $$

Consider now the conservation of kinetic energy

$$m_1 v_1^2 +m_2 v_2^2 $$

is constant or $$ m_1 ({\bf v}_1 + \alpha m_2\hat{\bf r})^2 + m_2 ({\bf v}_2 - \alpha m_1\hat{\bf r})^2 $$ is constant, implying $$ 2 \alpha m_1 m_2 \hat{\bf r}.\cdot ({\bf v_1-v_2}) + (m_1 m_2^2 + m_2 m_1^2)\alpha^2=0 $$ or

$$ 2 \hat{\bf r}.\cdot ({\bf v_1-v_2}) + (m_2 + m_1)\alpha = 0 $$

giving: $$\Delta {\bf v}_1 = -\hat {\bf r} \frac{2 m_2}{m_1+m_2} ( \hat {\bf r} \cdot ({\bf v_1-v_2}) ) $$ $$\Delta {\bf v}_2 = +\hat {\bf r} \frac{2 m_1}{m_1+m_2} ( \hat {\bf r} \cdot ({\bf v_1-v_2}) ) $$