Restitution and friction at a contact¶

Restitution and friction at a contact, given the unit normal n out of the surface.

The normal part of the relative velocity, v_n = v . n, reversed and scaled by the restitution e (Newton's law of restitution); the tangential part v_t = v - v_n n reduced by the friction coefficient mu times the normal impulse j_n = -(1 + e) v_n, never past zero (Coulomb's law as an impulse); a body already leaving the surface (v_n >= 0) untouched. v . n and |v_t| are written component by component, the same law. Source: Brach, Mechanical Impact Dynamics, ch. 3; Mirtich, Impulse-based Dynamic Simulation of Rigid Body Systems (PhD thesis, 1996), ch. 3. Checked against a reference in numpy, which has its own tests.

Source

Brach, Mechanical Impact Dynamics, ch. 3; Mirtich, Impulse-based Dynamic Simulation of Rigid Body Systems (PhD thesis, 1996), ch. 3.

Proof

curve, invariant, limit.

contact.trait

namespace lusine.integrator;

vector<float<m/s>, 3> contact(vector<float<m/s>, 3> velocity, vector<float, 3> normal, float restitution,
                              float friction)
{
  let vn = velocity.x * normal.x + velocity.y * normal.y + velocity.z * normal.z;
  if (vn >= 0.0 * vn) {
    return velocity;
  } else {
    let vt = velocity - vn * normal;
    let jn = -(1.0 + restitution) * vn;
    let speed = sqrt(vt.x * vt.x + vt.y * vt.y + vt.z * vt.z);
    if (speed <= friction * jn) {
      return -restitution * vn * normal;
    } else {
      return vt * ((speed - friction * jn) / speed) - restitution * vn * normal;
    }
  }
}

The normal velocity reversed and scaled by e; the slide reduced by mu times the normal impulse, never past stopping; a body leaving the surface untouched.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
from lusine.kernels._si import G0
EPS = float(np.finfo(np.float32).eps)   # Trait's float is binary32

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import trait_run from lusine.kernels._si import G0 EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32

The kernel¶

The law as the library states it, and as this page runs it.

In [2]:

print(trait_run.source('integrator/contact.trait'))

print(trait_run.source('integrator/contact.trait'))

// lusine.integrator.contact: restitution and friction at a contact, given the unit normal n out of the surface. The
// normal part of the relative velocity, v_n = v . n, reversed and scaled by the restitution e (Newton's law of
// restitution); the tangential part v_t = v - v_n n reduced by the friction coefficient mu times the normal impulse
// j_n = -(1 + e) v_n, never past zero (Coulomb's law as an impulse); a body already leaving the surface (v_n >= 0)
// untouched. v . n and |v_t| are written component by component, the same law.
// Source: Brach, Mechanical Impact Dynamics, ch. 3; Mirtich, Impulse-based Dynamic Simulation of Rigid Body Systems
// (PhD thesis, 1996), ch. 3.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.integrator;

vector<float<m/s>, 3> contact(vector<float<m/s>, 3> velocity, vector<float, 3> normal, float restitution,
                              float friction)
{
  let vn = velocity.x * normal.x + velocity.y * normal.y + velocity.z * normal.z;
  if (vn >= 0.0 * vn) {
    return velocity;
  } else {
    let vt = velocity - vn * normal;
    let jn = -(1.0 + restitution) * vn;
    let speed = sqrt(vt.x * vt.x + vt.y * vt.y + vt.z * vt.z);
    if (speed <= friction * jn) {
      return -restitution * vn * normal;
    } else {
      return vt * ((speed - friction * jn) / speed) - restitution * vn * normal;
    }
  }
}

The program¶

Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.

In [3]:

print(trait_run.law_program('integrator/contact.trait', 'contact'))

print(trait_run.law_program('integrator/contact.trait', 'contact'))

namespace probe;

lusine.version version = { { number = 1 } };

struct item {
  host vector<float<m/s>, 3> velocity;
  host vector<float, 3> normal;
  host float restitution;
  host float friction;
  vector<float<m/s>, 3> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.integrator.contact(items.velocity, items.normal, items.restitution, items.friction); }

Curve: each bounce reaches e² of the last¶

A ball dropped from h0 = 2 m on a floor of normal +z, for three restitutions side by side: the free flight between bounces is the page's (v = sqrt(2 g h) down, h = v² / 2g up), the bounce the kernel's. The heights on a log scale fall on lines of slope 2 log e.

In [4]:

e = np.array([0.5, 0.7, 0.9])
up = np.tile([0.0, 0.0, 1.0], (3, 1))
h = [np.full(3, 2.0)]
for _ in range(8):
    v = np.zeros((3, 3))
    v[:, 2] = -np.sqrt(2 * G0 * h[-1])
    after = trait_run.law('integrator/contact.trait', 'contact', velocity=v, normal=up, restitution=e, friction=0.0)
    assert np.all(after[:, :2] == 0.0) and np.all(after[:, 2] > 0)
    h.append(after[:, 2] ** 2 / (2 * G0))
h = np.array(h)
e32 = np.float32(e).astype(float)
assert np.allclose(h, 2.0 * e32 ** (2 * np.arange(9))[:, None], rtol=64 * EPS, atol=0)
fig, ax = plt.subplots(figsize=(7, 3.5))
for i in range(3):
    ax.semilogy(h[:, i], 'o', label=f'e = {e[i]}')
ax.set_xlabel('bounce')
ax.set_ylabel('height (m)')
ax.legend()
plt.show()

e = np.array([0.5, 0.7, 0.9]) up = np.tile([0.0, 0.0, 1.0], (3, 1)) h = [np.full(3, 2.0)] for _ in range(8): v = np.zeros((3, 3)) v[:, 2] = -np.sqrt(2 * G0 * h[-1]) after = trait_run.law('integrator/contact.trait', 'contact', velocity=v, normal=up, restitution=e, friction=0.0) assert np.all(after[:, :2] == 0.0) and np.all(after[:, 2] > 0) h.append(after[:, 2] ** 2 / (2 * G0)) h = np.array(h) e32 = np.float32(e).astype(float) assert np.allclose(h, 2.0 * e32 ** (2 * np.arange(9))[:, None], rtol=64 * EPS, atol=0) fig, ax = plt.subplots(figsize=(7, 3.5)) for i in range(3): ax.semilogy(h[:, i], 'o', label=f'e = {e[i]}') ax.set_xlabel('bounce') ax.set_ylabel('height (m)') ax.legend() plt.show()

Invariant: the normal part is e's, friction never reverses a slide, a body leaving is untouched¶

Two hundred contacts on surfaces of every tilt: the normal velocity after is -e times before, the tangential speed after is between zero and before and in the same direction, all to binary32's rounding; a contact separating comes back exactly as it went in.

In [5]:

rng = np.random.default_rng(2)
n = rng.normal(size=(200, 3))
n /= np.linalg.norm(n, axis=1)[:, None]
v = rng.normal(size=(200, 3)) * 3
v -= np.where(np.sum(v * n, axis=1) > 0, 2 * np.sum(v * n, axis=1), 0)[:, None] * n      # each into its surface
e, mu = rng.uniform(0, 1, 200), rng.uniform(0, 1, 200)
after = trait_run.law('integrator/contact.trait', 'contact', velocity=v, normal=n, restitution=e, friction=mu)
v, n, e = (np.float32(a).astype(float) for a in (v, n, e))
vn, an = np.sum(v * n, axis=1), np.sum(after * n, axis=1)
vt, at = v - vn[:, None] * n, after - an[:, None] * n
scale = 16 * EPS * np.linalg.norm(v, axis=1)
assert np.all(np.abs(an + e * vn) <= scale)
assert np.all(np.linalg.norm(at, axis=1) <= np.linalg.norm(vt, axis=1) + scale)
assert np.all(np.sum(at * vt, axis=1) >= -scale * np.linalg.norm(vt, axis=1))
leave = np.array([[1.0, 0.0, 2.0], [0.5, -3.0, 0.0], [-1.0, 2.0, 0.25]])
same = trait_run.law('integrator/contact.trait', 'contact', velocity=leave, normal=np.tile([0.0, 0.0, 1.0], (3, 1)),
                     restitution=0.5, friction=0.5)
assert np.array_equal(same, leave)
print('normal part -e of before; slides never reversed; a body leaving untouched')

rng = np.random.default_rng(2) n = rng.normal(size=(200, 3)) n /= np.linalg.norm(n, axis=1)[:, None] v = rng.normal(size=(200, 3)) * 3 v -= np.where(np.sum(v * n, axis=1) > 0, 2 * np.sum(v * n, axis=1), 0)[:, None] * n # each into its surface e, mu = rng.uniform(0, 1, 200), rng.uniform(0, 1, 200) after = trait_run.law('integrator/contact.trait', 'contact', velocity=v, normal=n, restitution=e, friction=mu) v, n, e = (np.float32(a).astype(float) for a in (v, n, e)) vn, an = np.sum(v * n, axis=1), np.sum(after * n, axis=1) vt, at = v - vn[:, None] * n, after - an[:, None] * n scale = 16 * EPS * np.linalg.norm(v, axis=1) assert np.all(np.abs(an + e * vn) <= scale) assert np.all(np.linalg.norm(at, axis=1) <= np.linalg.norm(vt, axis=1) + scale) assert np.all(np.sum(at * vt, axis=1) >= -scale * np.linalg.norm(vt, axis=1)) leave = np.array([[1.0, 0.0, 2.0], [0.5, -3.0, 0.0], [-1.0, 2.0, 0.25]]) same = trait_run.law('integrator/contact.trait', 'contact', velocity=leave, normal=np.tile([0.0, 0.0, 1.0], (3, 1)), restitution=0.5, friction=0.5) assert np.array_equal(same, leave) print('normal part -e of before; slides never reversed; a body leaving untouched')

normal part -e of before; slides never reversed; a body leaving untouched

Limit: elastic, plastic, frictionless, and a slide stopped¶

On a floor of normal +z, every value here is exact in binary32: e = 1 returns the normal speed whole, e = 0 stops it, mu = 0 leaves the slide alone; a slide of 3 m/s loses mu (1 + e) 2 m/s = 0.6 m/s, and one of 0.5 m/s, less than its share, is stopped, never reversed.

In [6]:

up = np.tile([0.0, 0.0, 1.0], (5, 1))
v = np.array([[1.5, 0.0, -2.0], [1.5, 0.0, -2.0], [1.5, 0.0, -2.0], [3.0, 0.0, -2.0], [0.5, 0.0, -2.0]])
after = trait_run.law('integrator/contact.trait', 'contact', velocity=v, normal=up,
                      restitution=[1.0, 0.0, 0.4, 0.5, 0.5], friction=[0.0, 0.0, 0.0, 0.2, 0.5])
assert np.array_equal(after[0], [1.5, 0.0, 2.0]) and np.array_equal(after[1], [1.5, 0.0, 0.0])
assert after[2, 0] == 1.5 and after[2, 1] == 0.0
assert abs(after[3, 0] - (3.0 - 0.2 * 1.5 * 2.0)) <= 4 * EPS * 3.0 and after[3, 2] == 1.0
assert after[4, 0] == 0.0 and after[4, 2] == 1.0
print(after)

up = np.tile([0.0, 0.0, 1.0], (5, 1)) v = np.array([[1.5, 0.0, -2.0], [1.5, 0.0, -2.0], [1.5, 0.0, -2.0], [3.0, 0.0, -2.0], [0.5, 0.0, -2.0]]) after = trait_run.law('integrator/contact.trait', 'contact', velocity=v, normal=up, restitution=[1.0, 0.0, 0.4, 0.5, 0.5], friction=[0.0, 0.0, 0.0, 0.2, 0.5]) assert np.array_equal(after[0], [1.5, 0.0, 2.0]) and np.array_equal(after[1], [1.5, 0.0, 0.0]) assert after[2, 0] == 1.5 and after[2, 1] == 0.0 assert abs(after[3, 0] - (3.0 - 0.2 * 1.5 * 2.0)) <= 4 * EPS * 3.0 and after[3, 2] == 1.0 assert after[4, 0] == 0.0 and after[4, 2] == 1.0 print(after)

[[1.5        0.         2.        ]
 [1.5        0.         0.        ]
 [1.5        0.         0.80000001]
 [2.4000001  0.         1.        ]
 [0.         0.         1.        ]]