Impact: two bodies meeting¶
Two bodies meeting, their velocities along the line between their centres afterwards.
Momentum is kept, p = m_a v_a + m_b v_b, and Newton's law of restitution sets the speed they part at, e times the speed they met at: v_a' = (p + m_b e (v_b - v_a)) / (m_a + m_b), v_b' = (p + m_a e (v_a - v_b)) / (m_a + m_b). e = 1 keeps the kinetic energy (an elastic hit), e = 0 leaves both moving together. The restitution is the pair's, data the game passes (its own: the library states no restitution); in 3-D the law is the normal component's, the tangential ones untouched. A body against a fixed surface is lusine.integrator.contact. Source: Brach, Mechanical Impact Dynamics, ch. 2; Goldsmith, Impact, ch. 1. Checked against a reference in numpy, which has its own tests.
- Source
- Brach, Mechanical Impact Dynamics, ch. 2; Goldsmith, Impact, ch. 1.
- Proof
- invariant, limit, curve.
impact.trait
namespace lusine.integrator;
(float<m/s>, float<m/s>) impact(float<kg> mass_a, float<m/s> velocity_a, float<kg> mass_b, float<m/s> velocity_b,
float restitution)
{
let momentum = mass_a * velocity_a + mass_b * velocity_b;
let total = mass_a + mass_b;
return ((momentum + mass_b * restitution * (velocity_b - velocity_a)) / total,
(momentum + mass_a * restitution * (velocity_a - velocity_b)) / total);
}
Two bodies' velocities along the line between their centres after they meet, momentum kept and the speed they part at e times the speed they met at (Newton's law of restitution), v_a' = (p + m_b e (v_b − v_a)) / (m_a + m_b) and v_b' = (p + m_a e (v_a − v_b)) / (m_a + m_b). The restitution is the pair's, which the game passes.
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
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
The kernel¶
The law as the library states it.
print(trait_run.source('integrator/impact.trait'))
// lusine.integrator.impact: two bodies meeting, their velocities along the line between their centres afterwards.
// Momentum is kept, p = m_a v_a + m_b v_b, and Newton's law of restitution sets the speed they part at, e times the
// speed they met at: v_a' = (p + m_b e (v_b - v_a)) / (m_a + m_b), v_b' = (p + m_a e (v_a - v_b)) / (m_a + m_b). e = 1
// keeps the kinetic energy (an elastic hit), e = 0 leaves both moving together. The restitution is the pair's, data
// the game passes (its own: the library states no restitution); in 3-D the law is the normal component's,
// the tangential ones untouched. A body against a fixed surface is lusine.integrator.contact.
// Source: Brach, Mechanical Impact Dynamics, ch. 2; Goldsmith, Impact, ch. 1.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.integrator;
(float<m/s>, float<m/s>) impact(float<kg> mass_a, float<m/s> velocity_a, float<kg> mass_b, float<m/s> velocity_b,
float restitution)
{
let momentum = mass_a * velocity_a + mass_b * velocity_b;
let total = mass_a + mass_b;
return ((momentum + mass_b * restitution * (velocity_b - velocity_a)) / total,
(momentum + mass_a * restitution * (velocity_a - velocity_b)) / total);
}
The program¶
Each section runs the kernel over a lattice of items, one an impact: the host sets every parameter, and main() binds the kernel over the items, its two results to two fields.
print(trait_run.law_program('integrator/impact.trait', 'impact'))
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float<kg> mass_a;
host float<m/s> velocity_a;
host float<kg> mass_b;
host float<m/s> velocity_b;
host float restitution;
float<m/s> result_0;
float<m/s> result_1;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { (items.result_0, items.result_1) = lusine.integrator.impact(items.mass_a, items.velocity_a, items.mass_b, items.velocity_b, items.restitution); }
Invariant: momentum kept, and Newton's restitution¶
Over 500 random impacts, m_a v_a' + m_b v_b' is the momentum before and v_b' − v_a' is e (v_a − v_b), each to binary32's rounding of the terms.
def impact(**kw):
return trait_run.law('integrator/impact.trait', 'impact', **kw)
rng = np.random.default_rng(299)
f32 = lambda x: np.float32(x).astype(float)
ma, mb = f32(rng.uniform(0.1, 5.0, 500)), f32(rng.uniform(0.1, 5.0, 500))
va, vb, e = f32(rng.normal(size=500) * 3), f32(rng.normal(size=500) * 3), f32(rng.uniform(0.0, 1.0, 500))
a, b = impact(mass_a=ma, velocity_a=va, mass_b=mb, velocity_b=vb, restitution=e)
scale = np.abs(ma * va) + np.abs(mb * vb) + (ma + mb) * np.abs(va - vb)
assert np.all(np.abs(ma * a + mb * b - (ma * va + mb * vb)) <= 16 * EPS * scale)
assert np.all(np.abs((b - a) - e * (va - vb)) <= 16 * EPS * scale / np.minimum(ma, mb))
print('momentum and restitution hold over', len(a), 'impacts')
momentum and restitution hold over 500 impacts
Limit: elastic, plastic, and a wall¶
Equal masses with e = 1 swap their velocities, exactly (every value here is exact in binary32); e = 0 leaves both at one speed, the same bits; e = 1 keeps the kinetic energy to binary32's rounding; a body a million times heavier, at rest, is a wall: the light one comes back at −e v_a.
a, b = impact(mass_a=1.25, velocity_a=3.0, mass_b=1.25, velocity_b=-0.5, restitution=1.0)
assert a == -0.5 and b == 3.0, (a, b)
a, b = impact(mass_a=ma, velocity_a=va, mass_b=mb, velocity_b=vb, restitution=0.0)
assert np.array_equal(a, b)
a, b = impact(mass_a=ma, velocity_a=va, mass_b=mb, velocity_b=vb, restitution=1.0)
ke = lambda m, v: 0.5 * m * v * v
before = ke(ma, va) + ke(mb, vb)
assert np.all(np.abs(ke(ma, a) + ke(mb, b) - before) <= 64 * EPS * (before + ke(ma + mb, np.abs(va) + np.abs(vb))))
e = np.array([0.0, 0.5, 0.9])
a, b = impact(mass_a=1.0, velocity_a=2.0, mass_b=1e6, velocity_b=0.0, restitution=e)
assert np.all(np.abs(a + e * 2.0) <= 1e-5 * 2.0) and np.all(np.abs(b) <= 1e-5), (a, b)
print('a wall returns', a)
a wall returns [ 1.99999795e-06 -9.99997020e-01 -1.79999626e+00]
Curve: the energy lost, ½ μ (1 − e²) u²¶
A 1 kg body at 2 m/s meets a 3 kg one at −1 m/s, for e from 0 to 1: the fraction of the closing energy ½ μ u² lost (μ the reduced mass, u the closing speed) is 1 − e², to binary32's rounding of the velocities.
e = np.linspace(0.0, 1.0, 41)
a, b = impact(mass_a=1.0, velocity_a=2.0, mass_b=3.0, velocity_b=-1.0, restitution=e)
mu, u = 1.0 * 3.0 / 4.0, 3.0
lost = (ke(1.0, 2.0) + ke(3.0, -1.0) - ke(1.0, a) - ke(3.0, b)) / (0.5 * mu * u * u)
assert np.all(np.abs(lost - (1 - f32(e) ** 2)) <= 32 * EPS), np.abs(lost - (1 - e ** 2)).max()
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.plot(e, lost, 'o', ms=3, label='the kernel')
ax.plot(e, 1 - e ** 2, label='1 - e²')
ax.set_xlabel('restitution e')
ax.set_ylabel('energy lost / ½ μ u²')
ax.legend()
plt.show()
