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:
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$
...
