1
$\begingroup$

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

    $\endgroup$
    • $\begingroup$ Usually, the bounce speed is reduced by the coefficient of restitution and not the energy, $\endgroup$ – ja72 Jul 30 '19 at 21:25
    • $\begingroup$ 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. $\endgroup$ – Jonas De Schouwer Jul 30 '19 at 21:48
    • $\begingroup$ 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. $\endgroup$ – ja72 Jul 31 '19 at 13:07
    • $\begingroup$ 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. $\endgroup$ – Jonas De Schouwer Jul 31 '19 at 13:32
    • 1
      $\begingroup$ @ggcg we use the fact that P is perpendicular to the ground to evaluate the cross products (with respect to the magnitude p). $\endgroup$ – Jonas De Schouwer Aug 13 '19 at 17:04
    1
    $\begingroup$

    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.

    $\endgroup$
    • 1
      $\begingroup$ 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. $\endgroup$ – Jonas De Schouwer Jul 31 '19 at 10:42
    • $\begingroup$ 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. $\endgroup$ – Jonas De Schouwer Jul 31 '19 at 12:30
    • $\begingroup$ @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). $\endgroup$ – ja72 Jul 31 '19 at 12:47
    • 1
      $\begingroup$ @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. $\endgroup$ – ja72 Jul 31 '19 at 14:25
    • 1
      $\begingroup$ I've asked this question in Game Development (instead of Physics). You might want to check it out. $\endgroup$ – Jonas De Schouwer Jul 31 '19 at 15:14

    Your Answer

    By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

    Not the answer you're looking for? Browse other questions tagged or ask your own question.