I have two spheres in 3-dimensional space. The mass and initial velocity of each sphere is known. If the spheres collide elastically, how can I calculate their final velocities? The spheres do not necessarily collide head-on.
-
3$\begingroup$ Does this answer your question? Calculating new velocities of $n$-dimensional particles after collision $\endgroup$– John AlexiouCommented Dec 8, 2020 at 23:36
-
$\begingroup$ Are you looking to implement a computational simulation, or an analytical model? $\endgroup$– Alexander KalianCommented Dec 9, 2020 at 2:17
3 Answers
It is trivial to extend this solution to 3 or more dimensions. Here, we will stick to 2D.
Solution:
Let us refer to the initial velocity of the first circle as $\overrightarrow{v}_1$, and the final velocity as $\overrightarrow{v}_1'$. $v_{1x}$ denotes the x component of $\overrightarrow{v}_1$. We will call the masses $m_1$ and $m_2$. Consider a vector from the center of circle 2 to the center of circle 1. We call this vector $\overrightarrow{n}$ since it is normal to both surfaces. The forces in the collision will be exerted along this vector, and the change in momentum will be equal and opposite for the two circles. This tells us that for some constant $c$:
$$\overrightarrow{v}_1' = \overrightarrow{v}_1 + \frac{c}{m_1}\overrightarrow{n}$$ $$\overrightarrow{v}_2' = \overrightarrow{v}_2 - \frac{c}{m_2}\overrightarrow{n}$$
For now, I will simply provide c:
$$c=2\frac{n_x\left(v_{2x}-v_{1x}\right)+n_y\left(v_{2y}-v_{1y}\right)}{\left(n_x^2+n_y^2\right)\left(\frac{1}{m_1}+\frac{1}{m_2}\right)}$$
In higher dimensions:
$$c=2\frac{n_x\left(v_{2x}-v_{1x}\right)+n_y\left(v_{2y}-v_{1y}\right)+\cdots}{\left|\overrightarrow{n}\right|\left(\frac{1}{m_1}+\frac{1}{m_2}\right)}$$
Derivation of $c$:
We can use the elasticity of the collision to solve for $c$. Elasticity tells us that the kinetic energy of the system does not change:
$$m_1v_1'^2+m_2v_2'^2=m_1v_1^2+m_2v_2^2$$
Note: we have dropped the $\frac12$ from the formula $ke = \frac12mv^2$.
The above formula uses magnitudes of velocity vectors, which we will need to calculate.
$$v_1=\sqrt{v_{1x}^2+v_{1y}^2}$$ $$v_2=\sqrt{v_{2x}^2+v_{2y}^2}$$ $$v_1'=\sqrt{\left(v_{1x}+\frac{c}{m_1}n_x\right)^2+\left(v_{1y}+\frac{c}{m_1}n_y\right)^2}$$ $$v_2'=\sqrt{\left(v_{2x}-\frac{c}{m_2}n_x\right)^2+\left(v_{2y}-\frac{c}{m_2}n_y\right)^2}$$
Let's now plug these magnitudes into our kinetic energy equation:
$$m_1\left(\left(v_{1x}+\frac{c}{m_1}n_x\right)^2+\left(v_{1y}+\frac{c}{m_1}n_y\right)^2\right)+m_2\left(\left(v_{2x}-\frac{c}{m_2}n_x\right)^2+\left(v_{2y}-\frac{c}{m_2}n_y\right)^2\right)=m_1\left(v_{1x}^2+v_{1y}^2\right)+m_2\left(v_{2x}^2+v_{2y}^2\right)$$
Evaluate the squares:
$$m_1\left(v_{1x}^2+2\frac{c}{m_1}v_{1x}n_x+\left(\frac{c}{m_1}n_x\right)^2+v_{1y}^2+2\frac{c}{m_1}v_{1y}n_y+\left(\frac{c}{m_1}n_y\right)^2\right)+m_2\left(v_{2x}^2-2\frac{c}{m_2}v_{2x}n_x+\left(\frac{c}{m_2}n_x\right)^2+v_{2y}^2-2\frac{c}{m_2}v_{2y}n_y+\left(\frac{c}{m_2}n_y\right)^2\right)=m_1\left(v_{1x}^2+v_{1y}^2\right)+m_2\left(v_{2x}^2+v_{2y}^2\right)$$
We can drop some terms, including the entire right side:
$$m_1\left(2\frac{c}{m_1}v_{1x}n_x+\left(\frac{c}{m_1}n_x\right)^2+2\frac{c}{m_1}v_{1y}n_y+\left(\frac{c}{m_1}n_y\right)^2\right)+m_2\left(-2\frac{c}{m_2}v_{2x}n_x+\left(\frac{c}{m_2}n_x\right)^2-2\frac{c}{m_2}v_{2y}n_y+\left(\frac{c}{m_2}n_y\right)^2\right)=0$$
We can factor $c$ out of this equation, which tells us that one solution has $c=0$. This is the initial state, and not what we want. So we will remove a factor of $c$. We will also distribute $m_1$ and $m_2$ to cancel the denominators and reorder terms:
$$2v_{1x}n_x+2v_{1y}n_y-2v_{2x}n_x-2v_{2y}n_y+\frac{c}{m_1}n_x^2+\frac{c}{m_1}n_y^2+\frac{c}{m_2}n_x^2+\frac{c}{m_2}n_y^2=0$$
$$c\left(n_x^2+n_y^2\right)\left(\frac{1}{m_1}+\frac{1}{m_2}\right)=2\left(v_{2x}n_x+v_{2y}n_y-v_{1x}n_x-v_{1y}n_y\right)$$
$$c=2\frac{n_x\left(v_{2x}-v_{1x}\right)+n_y\left(v_{2y}-v_{1y}\right)}{\left(n_x^2+n_y^2\right)\left(\frac{1}{m_1}+\frac{1}{m_2}\right)}$$
Thank you to John Alexiou for your answer to a similar question here: https://physics.stackexchange.com/a/220776/281943
-
1$\begingroup$ I have implemented this solution on a little simulation here: liambohl.com/projects/bubbles $\endgroup$ Commented Dec 8, 2020 at 22:14
-
1$\begingroup$ This isn't a simulation, it's a calculation. Perhaps you should rewrite the question title to reflect this. $\endgroup$ Commented Dec 9, 2020 at 0:00
-
$\begingroup$ The summary 'c' requires magnitude of n to be squared... as you indicate elsewhere in the answer. I attempted to make the edit for you but Stack Exchange requires 6 or more characters to be changed for an edit to get accepted $\endgroup$– DAGCommented Jul 6, 2023 at 17:13
-
$\begingroup$ Another thing - the 'c' summary (imo) would benefit from showing the numerator as the dot product of n and the difference of velocities. I would guess (hope?) that most people using this would be utilizing a numerical library, so they'll want to use dot function rather than spelling out the definition of it $\endgroup$– DAGCommented Jul 6, 2023 at 17:22
This linked answer gives the exact solution for colliding spheres in any dimension. But if you want a more general solution, where the contact normal does not go through the center of mass (off-center impulse) then follow the following procedure.
Given two bodies, with mass $m_1$ and $m_2$, and mass moment of inertia 3×3 tensor about each center of mass $\mathbf{I}_1$ and $\mathbf{I}_2$, and velocitity vectors $\boldsymbol{v}_1$ and $\boldsymbol{v}_2$ on their center of mass, in addition to their rotational velocity vectors $\boldsymbol{\omega}_1$ and $\boldsymbol{\omega}_2$.
They contact at a common point and the contact normal direction is $\boldsymbol{n}$. The location of each center of mass relative to the contact point is $\boldsymbol{c}_1$ and $\boldsymbol{c}_2$. Given a coefficient of restitution $\epsilon$, find the change is velocity for each body, $\Delta \boldsymbol{v}_1$ and $\Delta \boldsymbol{v}_2$ as well as $\Delta \boldsymbol{\omega}_1$ and $\Delta \boldsymbol{\omega}_2$.
For spheres you can assume that $\boldsymbol{n} = (\boldsymbol{c}_2 - \boldsymbol{c}_1)/\|\boldsymbol{c}_2 - \boldsymbol{c}_1\|$.
Project the velocity vectors along the contact normal to find the impact speed $$v_{\rm imp} = \boldsymbol{n} \cdot ( \boldsymbol{v}_2 - \boldsymbol{v}_1) \tag{1}$$ The convention here is a negative impact speed means the objects are approaching each other.
Find the reduced mass $m_{\rm eff}$ of the system along the contact normal. $$m_{\rm eff} = \frac{1}{ \frac{1}{m_1} +\boldsymbol{n} \cdot \mathbf{I}_1^{-1} ( \boldsymbol{n} \times \boldsymbol{c}_1) + \frac{1}{m_2} +\boldsymbol{n} \cdot \mathbf{I}_2^{-1} ( \boldsymbol{n} \times \boldsymbol{c}_2)} \tag{2} $$
Find the impulse magnitude $J$ needed to reverse their relative speed and reduce it by the coefficient of restitution $\epsilon$
$$ J = -(1+\epsilon)\,m_{\rm eff}\,v_{\rm imp} \tag{3} $$
Find the change in velocity vectors in both translation and rotation as a result of the impulse acting on equal and opposite terms on each body along $\boldsymbol{n}$ and through the contact point.
$$ \begin{aligned} \Delta \boldsymbol{v}_1 &= - \tfrac{J}{m_1} \boldsymbol{n} & \Delta \boldsymbol{\omega}_1 &= -\mathbf{I}_1^{-1} (\boldsymbol{n} \times \boldsymbol{c}_2) J \\ \Delta \boldsymbol{v}_2 &= + \tfrac{J}{m_2} \boldsymbol{n} & \Delta \boldsymbol{\omega}_2 &= +\mathbf{I}_2^{-1} (\boldsymbol{n} \times \boldsymbol{c}_2 ) J \end{aligned} \tag{4} $$
Update the velocities with $$\begin{aligned} \boldsymbol{v}_1 & \leftarrow \boldsymbol{v}_1 + \Delta \boldsymbol{v}_1 & \boldsymbol{\omega}_1 & \leftarrow \boldsymbol{\omega}_1 + \Delta \boldsymbol{\omega}_1 \\ \boldsymbol{v}_2 & \leftarrow \boldsymbol{v}_2 + \Delta \boldsymbol{v}_2 & \boldsymbol{\omega}_2 & \leftarrow \boldsymbol{\omega}_2 + \Delta \boldsymbol{\omega}_2 \end{aligned} \tag{5}$$
References:*
- Simulation of Rigid Bodies (http://www.cs.cmu.edu/~baraff/sigcourse/notesd1.pdf) and Modeling of Contacts (http://www.cs.cmu.edu/~baraff/sigcourse/notesd2.pdf)
- Off-center dynamics of Rigid Bodies (https://physics.stackexchange.com/a/80449/392)
For a computational simulation, a good solution is to run a molecular dynamics simulation with a simple algorithm such as the Verlet algorithm, which is an approximated implementation of the classical equations of motion (using Taylor expansions to a certain order) which work well computationally. The Verlet algorithm is often used in academic research to create deterministic simulations. By using the Verlet algorithm, you would need your code to calculate the new positions and velocities for each sphere as follows:
$$ x_{t+1} = 2x_t - x_{t-1} + \frac{F_t}{m} \Delta t^2 $$ $$ v_t = \frac{x_{t+1}-x_{t-1}}{2 \Delta t} $$
These equations need to be solved separately for each spatial dimension in your simulation (which can be any arbitrary amount). You can find values such as $x_{t-1}$ by rearranging the two equations above and performing various substitutions (the link gives more detail on this).
As you wish to use elastic collisions, $F_{tot}$ would need to be zero when no collisions occur and resolved via conservation of momentum whenever a collision does occur and split into its $F_x,F_y,F_z$ components (each corresponding to $F_t$ in the equations) by solving for the ratios between distances. $F_t$ can actually be used to incorporate forces from any kind of potential in physics, including chemical potentials, so you need to be specific in defining the how and when the forces apply.
If there is a force $F_{tot}$ acting on a sphere in a collision, $F_x,F_y,F_z$ can be calculated from $\Delta x,\Delta y,\Delta z$ between the two sphere's centres. At a collision, the sphere centres are separated by a distance $2r$ where r is the radius of a sphere (assuming they are of the same size - change if needed). By Pythagoras' theorem: $$2r = \sqrt{(\Delta x)^2 + (\Delta y)^2 + (\Delta z)^2} $$ The ratios between these distance components can be used to determine how $F_{tot}$ should be split among its dimensional components. The closer to $2r$ each distance component is, the closer to $F_{tot}$ its corresponding force component will be.
By implementing the Verlet algorithm and allowing it to evolve with discrete timesteps (taking care to select appropriate values that aren't so large that the spheres "jump" into one another - and consequently it may be useful to consider a "collision" at an accepted range of distances around $2r$), a working simulation can be made.
As this simulation would use 3D coordinates, these can easily be stored and then used to graphically represent the simulation's evolution.