Modelling an inelastic, rough, constrained collision

Question

Understood situations:

a) Inelastic, rough collision of free spheres

In such a collision, two coefficients are used. The coefficient of restitution in the normal direction (the ratio $c_N$ of relative normal velocities before and after, between $0$ and $1$) determines elasticity, with $1$ being perfectly elastic and $0$ perfectly inelastic. The coefficient of restitution in the tangential direction ($c_T$, the same for tangential velocities) determines smoothness, with $1$ being perfectly smooth and $-1$ perfectly rough.

Both may be somewhat dependent on the impact parameters (and not intrinsic to the spheres), but they (together with momentum and angular momentum conservation laws) define the collision result.

A further discussion can be found in this PDF from page 14 onwards.

Important source for understanding perfect roughness and the relation of roughness to energy conservation.

b) Elastic, smooth collision with constraints

One of the spheres is constrained to move along a circle and simplified to a point mass. The other sphere is free and impacts it at an oblique angle. We have three degrees of freedom and three conservation laws (two for angular momentum and one for energy). From these, we can determine the post-collision state of the system.

The problem:

Generalizing situation b) to inelasticity and roughness cannot be done through coefficients of restitution. A perfectly smooth collision need not always have $c_T = 1$ and a perfectly elastic collision need not have $c_N = -1$. A simple counterexample is a perfectly elastic and perfectly smooth collision of a very light orbiting point mass and a stationary free sphere. The point mass will bounce back with little change to the speed of the free sphere, effectively a $c_T$ of $-1$, even though the collision is smooth.

How can such a collision be characterised instead?

Insights:

There are two effects at hand: The normal bounce and the tangent bounce. If we disregard one, the other behaves in line with coefficients of restitution (i.e. a COR of $-1$ is a perfect bounce and a COR of $1$ is nothing changing).

In the general situation, we have four unknowns: The two speeds of the free marble (the most convenient coordinate system here being the speed in the normal direction and the speed in the tangent direction), the angular velocity of the point, and the rotational angular velocity of the marble. Angular momentum around the center of orbit (of the constrained point) is conserved and all forces act through the contact point, so the angular momentum of just the free marble around the contact point is also conserved. This gives us two equations.

There are three known solutions that conserve energy: A "collision" where all parameters stay the same (total pierce), a perfectly smooth collision, and a perfectly rough collision where the normal velocity of the free marble stays the same (just the tangent component of the collision).

By analogy with nonconstrained collisions, there should be a perfectly rough, perfectly elastic collision that conserves energy, but I cannot find it. Simply finding the differences of pre- and post-collision speeds and rotation speeds for both $c_T = -1$ and $c_N = -1$ and adding them together with the original values (adding the impulses) leads to the total energy changing (in either direction depending on the setup).

$c_T$ and $c_N$ are also the subject of a Q&A I posted with a simpler setup and more details on the results.

Answer 1

Sigh, I wrote this answer when there was more details in the question... I agree with JAlex that the simplest way to account for conversation of momentum is to use a single impulse $J$ to represent the collision.

$$\hat{v_t} = \frac{J_t}{m_2} + v_t$$ $$\hat{v_n} = \frac{J_n}{m_2} + v_n$$ $$\hat{\omega_S} = -\frac{J_t}{L} + \omega_S$$ $$\hat{\omega} = \frac{J \bullet r_1}{m_1(r_1 \bullet r_1)} + \omega$$

As a slight change in notation, $L$ is the angular inertia of the sphere.

To ensure that the collision doesn't result in penetration:

$$\hat{v_n} \leq -\hat{\omega} \, \sin(\alpha)$$

This results in a deformed half plane (topologically) constraint on the impulse.

To ensure that energy is conserved:

$$m_2(\hat{v_t}^2+\hat{v_n}^2)+L\hat{\omega_S}^2 + m_1 r_1^2 \hat{\omega}^2 \leq m_2(v_t^2+v_n^2)+L\,\omega_S^2 + m_1 r_1^2 \omega^2 $$

This is a deformed disk constraint.

The Intersection of these two constraints define the area of valid collisions. All the stuff about slip velocity and friction is really just to get a better estimate of what the tangential impulse will be, but unusual internal geometry / structure or external constraints can violate those rules.

In particular, because of the constraint on the point, you can have the relative slip-velocity double, or double in the opposite direction, even with no tangential impulse (frictionless).

If you want to model very high friction but elastic (aka super ball) collisions, then you need to define the deformation model you're going to use in order to get a single defined answer, rather than a valid range.

If you want to assume that the slip/sliding velocity will reduce from a non-zero value to zero during the collision then it doesn't make sense to try to conserve energy, as there must have been rubbing to reduce that slip velocity which must result in frictional losses.

Deformation Modeling

On possible deformation model is fully elastic:

$$ F_t = -k_t \, x_t $$ $$ F_n = -k_n \, x_n $$

Where $x$ is the contact displacement, and $k$ represents the stiffness of the material.

Geometric constraints: $$\frac{d \, x_n}{d\,t} = v_n + \omega\, r_1 \, \sin(\alpha) $$ $$\frac{d \, x_t}{d\,t} = v_t - \omega\, r_1 \, \cos(\alpha) - r_2 \, \omega_S $$

Equations of motion:

$$\frac{d \, v_t}{d\,t} = \frac{F_n}{m_2} $$ $$\frac{d \, v_n}{d\,t} = \frac{F_t}{m_2} $$ $$\frac{d \, \omega_S}{d\,t} = \frac{F_t}{L} $$ $$\frac{d \, \omega}{d \, t} = \frac{-F_n \sin(\alpha)-F_t \cos(\alpha)}{m_1 \, r_1}$$

Then if we initialize $x$ to zero we can integrate until $x_n$ is once again zero, and at that point we'll have our new velocities. Note that even though there's no damping in these equations it still doesn't guarantee that there isn't energy lost. If $x_t$ doesn't reach zero at the same time that $x_n$ does then there will be energy "lost" that's stored in the tangential stiffness when the collision ends.

This same integration could be done with damping terms added to the force equations to model less elastic collisions.

Answer 2

A few notes first. I disagree with handling friction with a coefficient of restitution. It makes more sense to me with first calculating first the tangential impulse $J_e$ that would be needed for the parts not to slip past each other (in addition to the normal impulse $J_n$ due to bounce) and then capping the magnitude to a value such that $|J_e| \leq \mu | J_n |$.

These causes three possible cases

• Frictionaless makes $J_e = 0$ always
• Frictional caps tangential impulse using the friction coefficient $\mu$ such that $|J_e| \leq \mu |J_n|$
• Rough acts like friction is infinite and does not modify the tangential impulse $J_e$, leaving it what it needs to be for parts to stick to each other tangentially (no slip).

As far as handling impulses to constrained bodies, it make the situation a little harder because you have a kinematic constraint you need to enforce as well as the law of contacts. You quickly find yourself assembling large complex vector expressions and trying to solve them using linear algebra.

I feel like the proper handling of impulses on constrained bodies is a subject to a separate question, with the subset of 1DOF, 2DOF or 3DOF planar bodies. In general, this is two step process

• First calculate the impulse(s) on the contact from the reduced mass of each body $m^{-1}_{i}$ and the impact speed $v_{\rm imp}$ with resulting impulse $$ J_n = (1+\epsilon) \tfrac{1}{ m_1^{-1} + m_2^{-1}} v_{\rm imp}$$ Something similar happens in the tangential direction to calculate $J_e$ but with different reduced mass, since the direction where the impulse is applied is different from the normal case.

• Then calculate the change in joint speeds from inverse kinematics. Part of the applied impulse goes into joint rotation, and the remaining goes into a reaction impulse at the joint.

  So in a free body diagram kind of sense each body is subject to a contact normal impulse $J_n$, a contact frictional impulse $J_e$ and one or more constraint reaction impulses $G_j$, resulting in a change in the joint degree of freedom speeds $\Delta \dot{q}_{3-j}$.

  _{The subscript $j$ iterates through the number of constraints on the joint, and $3-j$ iterates through the number of degrees of freedom.}

The exact details details with such problems are quite complex, and in the case of articulated bodies subject to contacts, part of ongoing research in graduate level robotics or mechanics research (engineering/physics/computer science all have parts in this).

Answer 3

For posterity I am adding another answer, showing specifically how to handle a single constrained rigid body (pinned) under the influence of an impulse along a specified direction.

Consider a pined body with the center of mass at $\boldsymbol{c}$ and an impulse $J$ along the direction $\boldsymbol{n}$ at position $\boldsymbol{r}_{A}$. The pin is at the origin G, has direction $\boldsymbol{z}$ and speed $\dot{\theta}$. The center of mass is at point C, the contact at point A and the pin at the origin. The mass is m and the mass moment of inertia at the center of mass is $\mathcal{I}_{C}$.

Kinematics

$$\begin{aligned}\boldsymbol{\omega} & =\boldsymbol{z}\,\dot{\theta}\\ \boldsymbol{v}_{C} & =-\boldsymbol{c}\times\boldsymbol{\omega}=-\left(\boldsymbol{c}\times\boldsymbol{z}\right)\,\dot{\theta}\\ \boldsymbol{v}_{A} & =-\boldsymbol{r}_{A}\times\boldsymbol{\omega}=-\left(\boldsymbol{r}_{A}\times\boldsymbol{z}\right)\,\dot{\theta} \end{aligned} \tag{1}$$

and similarly for changes in speed, as in $\Delta \boldsymbol{\omega} =\boldsymbol{z}\,\Delta \dot{\theta}$, etc.

Dynamics

The impulse vector $\boldsymbol{J}=\boldsymbol{n}\,J$ at A causes a reaction impulse vector $\boldsymbol{G}$ at the pin, in addition to a reaction moment $\boldsymbol{\tau}_{G}$. The balance of momentum is thus $$\begin{aligned}\boldsymbol{G}+\boldsymbol{n}J & =m\,\Delta\boldsymbol{v}_{C}=-m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\,\Delta\dot{\theta}\\ \boldsymbol{\tau}_{G}-\boldsymbol{c}\times\boldsymbol{G}+\left(\boldsymbol{r}_{A}-\boldsymbol{c}\right)\times\boldsymbol{n}J & =\mathcal{I}_{C}\Delta\boldsymbol{\omega}=-\mathcal{I}_{C}\boldsymbol{z}\,\Delta\dot{\theta} \end{aligned} \tag{2}$$

The above is 6 equations, with 7 unknowns (three impulse reaction forces, three impulse reaction torques and the joint speed).

Joint Condition

To solve the above we need one more equation. This comes from the power relationship of the joint. Namely power though a joint is defined as $P=\Delta\boldsymbol{\omega}\cdot\boldsymbol{\tau}_{G}+\Delta\boldsymbol{v}_{G}\cdot\boldsymbol{G}=0$ or for a pin joint since $\Delta\boldsymbol{v}_{G}=0$ and $\Delta\boldsymbol{\omega=z}\Delta\dot{\theta}$ the power expression is $$ \left(\boldsymbol{z}\cdot\boldsymbol{\tau}_{G}\right)\Delta\dot{\theta}=0 \tag{3}$$

The 7th equation is they key to solving this problem. Substitute $\boldsymbol{G}$ from the first equation in (2) into the second equation and solve for $\boldsymbol{\tau}_{G}$. The use this in the power equation.

$$\begin{aligned}\boldsymbol{G} & =-m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\,\Delta\dot{\theta}-\boldsymbol{n}J\\ \boldsymbol{\tau}_{G} & =-\mathcal{I}_{C}\boldsymbol{z}\,\Delta\dot{\theta}-\left(\boldsymbol{r}_{A}-\boldsymbol{c}\right)\times\boldsymbol{n}J+\boldsymbol{c} \times\boldsymbol{G}\\ & =-\left(\mathcal{I}_{C}\boldsymbol{z}+\boldsymbol{c}\times m\left( \boldsymbol{c}\times\boldsymbol{z}\right)\right)\,\Delta\dot{\theta}-\boldsymbol{r}_{A}\times\boldsymbol{n}J\\ 0 & = \boldsymbol{z}\cdot\left(-\left(\mathcal{I}_{C}\boldsymbol{z}+\boldsymbol{c}\times m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\right)\,\Delta\dot{\theta} - \boldsymbol{r}_{A}\times\boldsymbol{n}J\right)\\ & =-\left(\boldsymbol{z}\cdot\mathcal{I}_{C}\boldsymbol{z} +\boldsymbol{z}\cdot\left(\boldsymbol{c}\times m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\right)\right) \,\Delta\dot{\theta}-\boldsymbol{z}\cdot\left( \boldsymbol{r}_{A}\times\boldsymbol{n}\right)J \end{aligned} \tag{5} $$

with the solution

$$\boxed{\Delta\dot{\theta}=-\left(\frac{\boldsymbol{n}\cdot\left(\boldsymbol{z}\times\boldsymbol{r}_{A}\right)}{\boldsymbol{z}\cdot\mathcal{I}_{C}\boldsymbol{z}+\boldsymbol{z}\cdot\left(\boldsymbol{c}\times m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\right)}\right)\;J} \tag{6}$$

and the response to the impulse $$\begin{aligned}\Delta\boldsymbol{\omega} & =-\boldsymbol{n}\cdot\left(\boldsymbol{z}\times\boldsymbol{r}_{A}\right)\left(\frac{\boldsymbol{z}}{\boldsymbol{z}\cdot\mathcal{I}_{C}\boldsymbol{z}+\boldsymbol{z}\cdot\left(\boldsymbol{c}\times m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\right)}\right)\,J\\ \Delta\boldsymbol{v}_{A} & =\boldsymbol{n}\cdot\left(\boldsymbol{z}\times\boldsymbol{r}_{A}\right)\left(\frac{\boldsymbol{r}_{A}\times\boldsymbol{z}}{\boldsymbol{z}\cdot\mathcal{I}_{C}\boldsymbol{z}+\boldsymbol{z}\cdot\left(\boldsymbol{c}\times m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\right)}\right)\;J\\ \boldsymbol{G} & =-\left(\boldsymbol{n}-\boldsymbol{n}\cdot\left(\boldsymbol{z}\times\boldsymbol{r}_{A}\right)\,\left(\frac{m\left(\boldsymbol{c}\times\boldsymbol{z}\right)}{\boldsymbol{z}\cdot\mathcal{I}_{C}\boldsymbol{z}+\boldsymbol{z}\cdot\left(\boldsymbol{c}\times m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\right)}\right)\right)\,J \end{aligned} \tag{7}$$

Contact Condition

If the impulse J is due to a contact with coefficient of restitution \epsilon then the law of contacts is $$\begin{aligned}\boldsymbol{n}\cdot\left(\boldsymbol{v}_{A}+\Delta\boldsymbol{v}_{A}\right) & =-\epsilon\,\left(\boldsymbol{n}\cdot\boldsymbol{v}_{A}\right)\\ \boldsymbol{n}\cdot\Delta\boldsymbol{v}_{A} & =-\left(1+\epsilon\right)\,\left(\boldsymbol{n}\cdot\boldsymbol{v}_{A}\right)\\ -\left(\frac{\left(\boldsymbol{n}\cdot\left(\boldsymbol{z}\times\boldsymbol{r}_{A}\right)\right)^{2}}{\boldsymbol{z}\cdot\mathcal{I}_{C}\boldsymbol{z}+\boldsymbol{z}\cdot\left(\boldsymbol{c}\times m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\right)}\right)\;J & =-\left(1+\epsilon\right)\,\left(\boldsymbol{n}\cdot\boldsymbol{v}_{A}\right) \end{aligned} \tag{8}$$

Note that the impact speed is $v_{{\rm imp}}=\boldsymbol{n}\cdot\boldsymbol{v}_{A}=\boldsymbol{n}\cdot\left(\boldsymbol{z}\times\boldsymbol{r}_{A}\right)\,\dot{\theta}$, and the impulse is solved with $$\begin{aligned}\boldsymbol{n}\cdot\Delta\boldsymbol{v}_{A} & =-\left(1+\epsilon\right)\,\boldsymbol{n}\cdot\left(\boldsymbol{z}\times\boldsymbol{r}_{A}\right)\,\dot{\theta}\\ \boldsymbol{n}\cdot\left(\boldsymbol{z}\times\boldsymbol{r}_{A}\right)\left(\frac{\boldsymbol{n}\cdot\left(\boldsymbol{r}_{A}\times\boldsymbol{z}\right)}{\boldsymbol{z}\cdot\mathcal{I}_{C}\boldsymbol{z}+\boldsymbol{z}\cdot\left(\boldsymbol{c}\times m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\right)}\right)\;J & =-\left(1+\epsilon\right)\,\boldsymbol{n}\cdot\left(\boldsymbol{z}\times\boldsymbol{r}_{A}\right)\,\dot{\theta}\\ \left(\frac{\boldsymbol{n}\cdot\left(\boldsymbol{r}_{A}\times\boldsymbol{z}\right)}{\boldsymbol{z}\cdot\mathcal{I}_{C}\boldsymbol{z}+\boldsymbol{z}\cdot\left(\boldsymbol{c}\times m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\right)}\right)\;J & =-\left(1+\epsilon\right)\,\dot{\theta} \end{aligned} \tag{9}$$

and the required impulse for the bounce of $$\boxed{J=\left(1+\epsilon\right)\left(\frac{\boldsymbol{z}\cdot\mathcal{I}_{C}\boldsymbol{z}+\boldsymbol{z}\cdot\left(\boldsymbol{c}\times m\left(\boldsymbol{c}\times\boldsymbol{z}\right)\right)}{\boldsymbol{n}\cdot\left(\boldsymbol{r}_{A}\times\boldsymbol{z}\right)}\right)\,\dot{\theta}} \tag{10}$$