# Rigid body dynamics: modelling a polygon bounce off ground

I'm currently making a physics simulator, but I'm having some trouble making a polygon bounce off the floor.
I know that collisions are normally modeled as described below by ja72, but I did it somewhat differently.

What I'm trying to do:
When the polygon touches (penetrates) the ground, I want to apply a linear impulse $$\overrightarrow{P}$$ at the contact point such that the energy $$E'$$ of the polygon after the bounce is some elasticity constant $$\epsilon\in[0,1]$$ times the energy $$E$$ of the polygon before the bounce. $$E'=\epsilon E$$

My calculations:
Consider a cartesian coordinate system where the floor is $$y=0$$. Then $$\overrightarrow{P}(0,p)$$ where we expect $$p$$ to be positive. Define:

• $$m$$: the mass of the polygon
• $$I$$: the rotational inertia
• $$\overrightarrow{v}$$: the initial velocity
• $$\overrightarrow{v'}$$: the velocity after the collision
• $$\overrightarrow{\omega}$$: the initial rotational velocity (a vector parallel to the $$z$$-axis)
• $$\overrightarrow{\omega'}$$: the rotational velocity after the collision
• $$\overrightarrow{r}$$: a vector going from the polygon's center of mass to the contact point

Furthermore, for a vector $$\overrightarrow{a}$$, we denote its length by $$|\overrightarrow{a}|$$ and we can index it as $$a[0]$$ or $$a[1]$$ to get the $$x$$- and $$y$$-component (scalars, not vectors), respectively.

$$m$$, $$I$$, $$\overrightarrow{v}$$, $$\overrightarrow{\omega}$$ and $$r$$ are given, and we want to calculate the magnitude $$p$$ of the impulse $$\overrightarrow{P}$$. Because $$\overrightarrow{v'}=\overrightarrow{v}+\frac{\overrightarrow{P}}{m}$$ and $$\overrightarrow{\omega'} = \overrightarrow{\omega} + \frac{\overrightarrow{r}\times\overrightarrow{P}}{I}$$, we can calculate: $$\frac{1}{2}m\cdot\overrightarrow{v'}^2 + \frac{1}{2}I\cdot\overrightarrow{\omega'}^2 = E' = \epsilon E = \epsilon\cdot\left(\frac{1}{2}m\cdot\overrightarrow{v}^2 + \frac{1}{2}I\cdot\overrightarrow{\omega}^2\right)$$ $$\Leftrightarrow m\cdot\left(\overrightarrow{v}+\frac{\overrightarrow{P}}{m}\right)^2 + I\cdot\left(\overrightarrow{\omega}+\frac{\overrightarrow{r}\times\overrightarrow{P}}{I}\right)^2 = \epsilon\cdot\left(m\cdot\overrightarrow{v}^2 + I\cdot\overrightarrow{\omega}^2\right)$$ $$\Leftrightarrow m\cdot\overrightarrow{v}^2 + 2\overrightarrow{v}\cdot\overrightarrow{P} + \frac{\overrightarrow{P}^2}{m} + I\cdot\overrightarrow{\omega}^2 + 2\overrightarrow{\omega}\cdot(\overrightarrow{r}\times\overrightarrow{P}) + \frac{(\overrightarrow{r}\times\overrightarrow{P})^2}{I} = \epsilon\cdot\left(m\cdot\overrightarrow{v}^2 + I\cdot\overrightarrow{\omega}^2\right)$$

And now, we can use the fact that $$\overrightarrow{P}(0,p)$$ to calculate the dot and cross products: $$0 = (1-\epsilon)\cdot\left(m\cdot\overrightarrow{v}^2 + I\cdot\overrightarrow{\omega}^2\right) + 2p\cdot v[1] + \frac{p^2}{m} + 2p\cdot r[0]\cdot|\overrightarrow{\omega}| + \frac{(p\cdot r[0])^2}{I}$$ $$\Leftrightarrow \left[\frac{1}{m}+\frac{r[0]^2}{I}\right]\cdot p^2 + \left[2v[1]+2r[0]\cdot|\overrightarrow{\omega}|\right]\cdot p + (1-\epsilon)\cdot\left[m\cdot\overrightarrow{v}^2 + I\cdot\overrightarrow{\omega}^2\right]$$

This is a quadratic equation in $$p$$. If we set:

• $$A = \frac{1}{m}+\frac{r[0]^2}{I}$$
• $$B = 2v[1]+2r[0]\cdot|\overrightarrow{\omega}|$$
• $$C = (1-\epsilon)\cdot\left[m\cdot\overrightarrow{v}^2 + I\cdot\overrightarrow{\omega}^2\right]$$
• We can calculate the determinant $$D = B^2-4AC$$.

Now, it is clear that $$p=\frac{-B\pm\sqrt{D}}{2A}$$. Actually, $$p=\frac{-B+\sqrt{D}}{2A}$$, as $$B$$ must be negative to enforce a collision, and the sign of $$v[1]$$ must be reversed.

When I run this in a simulation, it looks very natural. However, it always glitches when $$D<0$$. I've tried to resolve this problem in countlessly many ways $$(*)$$, but it never worked out well.

Can anyone check whether my method and calculations are correct? And if not, how can I model this collision instead?

Here is a snippet of my code (in Python):

class POLYGON:
def bounce(self):
#Checks whether the polygon touches the ground and makes it rebound if necessary

low_nodes = []    #makes a list of the nodes with y<0
for node in self.nodes:
if node[1] <= 0:
low_nodes.append(node)

if len(low_nodes) >= 1:
speed = self.speed
self.translate([0, -2*lowest_node[1])    #lifts the lowest node above the ground

for node in low_nodes:
r = node - self.pos    #calculates r

A = 1/self.mass + r[0]**2/self.rot_inertia
B = 2*self.speed[1] + 2*self.rot_speed*r[0]
C = (1-restituence_constant) * (self.mass*np.linalg.norm(self.speed)**2 + self.rot_inertia*self.rot_speed**2)
D = B**2 - 4*A*C

if D<0:
impulse = [0, -B/(2*A)]
else:
impulse = [0, (-B-np.sqrt(D)) /(2*A) /len(low_nodes)]
self.applyLinearImpulse(impulse, node)

$$(*)$$
For example:
using the absolute value of $$D$$
regarding $$p$$ as a complex number and applying the impulse $$\overrightarrow{P}(\Im(p), \Re(p))$$
setting $$p=\frac{-B}{2A}$$ if $$D<0$$
...

• Usually, the bounce speed is reduced by the coefficient of restitution and not the energy, – ja72 Jul 30 '19 at 21:25
• Yes, I know. However, I found my objects bouncing unnaturally when I applied the impulse P(0, -(1+e) * self.velocity[1] * self.mass), as no rotational characteristics are concerned. I thought this would be a better way of modelling the collision, as the kinetic energy keeps decreasing. – Jonas De Schouwer Jul 30 '19 at 21:48
• So instead of the ratio of speeds, you specify the ratio of energy. But how is the contact direction enter into the calculation? If you have a contact that was largely tangential, it seems to me, your method would slow it down like it enters a bath of oil instead of bouncing off at a specific direction. – ja72 Jul 31 '19 at 13:07
• When r[0]^2/I is close to 0 (almost always), then A~=1/m. Let's say A=1/m. Then p = (-B+sqrt(D)) / (2A) >= -v_tan[1] / (1/m) >= -m*v_tan[1], so the object bounces off. The problem occurs when D<0 or when r[0]^2/I ís big (happens very rarely). I have not (yet) found a solution to this issue. – Jonas De Schouwer Jul 31 '19 at 13:32
• @ggcg we use the fact that P is perpendicular to the ground to evaluate the cross products (with respect to the magnitude p). – Jonas De Schouwer Aug 13 '19 at 17:04

The physics of contact can be found at many online resources, so I am just going to summarize below.

The impulse $$\vec{P}$$ needs to decomposed into the magnitude $$p$$ and the known contact normal direction $$\vec{n}$$, such that $$\vec{P} = p\, \vec{n}$$. The impulse magnitude is found from the law of collision that states that the relative velocity at the contact point and along the contact normal after the collision is a fraction of the relative impact velocity

$$v_{\rm relative}' = -\epsilon \, v_{\rm impact}$$

mathematically the above is expressed as follows when an object impacts the ground (single body impact).

$$\vec{n} \cdot (\vec{v}'+\vec{\omega}' \times \vec{r}) = -\epsilon \;\; \vec{n} \cdot ( \vec{v} + \vec{\omega} \times \vec{r})$$

where $$\vec{v}$$ and $$\vec{v}'$$ are expressed at the center of mass before and after the impact, and the rest according to the author's post.

State the equations of motion & kinematics to get equation 8-18 from the linked paper that gives the impact magnitude

$$p = \frac{ -(1+\epsilon)\; \vec{n}\cdot \vec{v} }{ \frac{1}{m} + \vec{n} \cdot \left( I^{-1} (\vec{r} \times \vec{n}) \right) \times \vec{r} }$$

The final velocity will be adjusted then accordingly

\begin{aligned} \vec{P} & = p \,\vec{n} \\ \vec{v}' & = \vec{v} + \tfrac{1}{m} \vec{P} \\ \vec{\omega}' & = \vec{\omega} + I^{-1} \left( \vec{r} \times \vec{P} \right) \end{aligned}

I have also linked an answer to a similar question in 2D.

• Thanks for the link. I had already found this formula, but without explanation. So I thought that v_rel was the relative velocity between (the centers of mass of) A and B, instead of the relative linear velocity of the contact points. – Jonas De Schouwer Jul 31 '19 at 10:42
• I'm not going to accept this answer (yet), though, as I would still like to know whether my method is sound and how to improve it. – Jonas De Schouwer Jul 31 '19 at 12:30
• @JonasDeSchouwer - thank you for pointing out the inconsistency in nomenclature in my post. I have corrected it to have $\vec{v}$ specified at the center of mass (per your post). – ja72 Jul 31 '19 at 12:47
• @JonasDeSchouwer - I have seen the same effect where an object falling flat barely bounces, but when slightly rotated it bounces vigorously. I think the culprit is a succession of rapid contacts on alternating corners, each sapping some energy. Each contact counter-rotating the object, causing the next contact on the opposite corner. I don't have a good fix, but this can form the basis of a great new question to be posted in Physics. Note that robust modeling of contacts is still a subject of ongoing research. There are still novel answers to be found there. – ja72 Jul 31 '19 at 14:25
• I've asked this question in Game Development (instead of Physics). You might want to check it out. – Jonas De Schouwer Jul 31 '19 at 15:14