Impulse, Restitution and Friction

Goal

Compute collision response with impulse, and check in numbers and pictures what the coefficient of restitution and friction decide. When this lab is done, you can build the response layer of a physics engine yourself.

Why it matters

The contact time of hard objects is less than 1 millisecond, so in a simulation that runs at 1/60 second you cannot imitate the contact process. That is why a physics engine skips the process and computes only the result. It sets up the relationship between velocities before and after the collision as an equation, and finds and applies the velocity change that satisfies it in one go.

There are two conventions to learn here. One is to store inverse mass. If you set an immovable object to inverse mass 0, the same expression runs as it is with no special handling. The other is not touching anything that is already moving apart. If you leave out this check, you pull back an object that is trying to separate and it sticks to the surface.

Steps

1. Put the toolbox in /root/impulse.
2. /root/impulse/impulse.py — the impulse formula.
3. /root/impulse/out/03-conserve.txt — momentum and energy.
4. /root/impulse/out/bounce.csv, 04-bounce.txt, bounce.png — the bouncing height.
5. /root/impulse/out/05-friction.txt — the sliding box.
6. /root/impulse/out/06-tangent.txt — friction along the tangent.
7. /root/impulse/out/balls.png and 07-balls.txt — five balls in a box.

Notes

• Start a server with nohup python3 -m http.server 8080 -d /root/impulse/out & and open http://localhost:8080/ in the Web preview.
• Common mistake 1: applying the impulse even though the relative velocity is already positive (moving apart). The object sticks to the surface.
• Common mistake 2: not clamping the velocity at 0 in friction. The velocity goes down to negative and the box goes backward.
• Common mistake 3: computing the tangent direction from the velocity after applying the normal impulse. The values change.

Put the drawing toolbox in place

Save /root/impulse/gfxlib.py exactly as in the example, and use /root/impulse/check.py to draw a test pattern and make /root/impulse/out/00-check.png. The pattern is a 64x64 black background with a white (255,255,255) diagonal line from (0,0) to (63,63), and over it a red (255,0,0) horizontal line from (0,32) to (63,32).

From this lab on, you do not rebuild the PNG encoder. We hand you the same code you built by hand in the first lab as a tool — because file formats are not what you learn here.

The lab Pod has no volume, so the files you made in the previous lab are not kept. That is why each lab starts by putting the toolbox in place again.

Create Canvas(w, h, bg), draw the two lines with line(x0, y0, x1, y1, rgb), and then save with write_png(path). Draw the horizontal line later, so that the intersection (32,32) becomes red.

In this lab you use it to draw the graph of the bouncing height and the trajectories of the balls.

The impulse formula

In /root/impulse/impulse.py, make resolve(v1, inv_m1, v2, inv_m2, n, e). Velocities are 2D tuples, and n is the unit vector pointing from 1 to 2. If the normal component of the relative velocity is 0 or greater, do nothing and return them as they are, and if it is negative, find j = -(1+e)*v_rel/(inv_m1+inv_m2), apply v1 -= j*inv_m1*n and v2 += j*inv_m2*n, and return the result as (v1, v2).

v_rel = dot(v2 - v1, n). This value being positive means the two objects are already moving apart, so you must not touch them. If you leave out this check, you pull back an object that is trying to separate and make it stick to the surface.

The reason it takes inverse mass rather than mass is that an immovable object can be expressed naturally as inv_m = 0. If you put in infinity, every division needs special handling.

To check, collide two objects of the same mass head on. With e=1 the velocities are exchanged, and with e=0 both end up with the same velocity.

What is conserved and what disappears

Solve the case where object 1 of mass 2 (velocity (3,0)) and object 2 of mass 1 (velocity (-1,0)) collide along normal (1,0) for e=1, 0.5 and 0, and write eight lines to /root/impulse/out/03-conserve.txt: p_before=, ke_before=, p_e1=, ke_e1=, p_e05=, ke_e05=, p_e0= and ke_e0=. Momentum is the sum of the x components, and kinetic energy is the sum of 0.5*m*v² (six decimal places).

The inverse masses are 0.5 and 1.0 respectively.

Momentum must be the same in all three cases. It is structural, because the two objects receive impulses of equal magnitude and opposite direction. If the values come out different, a sign is wrong.

Kinetic energy is conserved only when e=1, and decreases as e gets smaller. In reality the lost energy goes into sound, heat and deformation.

With e=0, the velocities of the two objects become the same. That means they stick and move together.

The bouncing height shrinks by e squared

Drop a ball from y=5.0, v=0 and run it for 2400 steps with g=-9.8 and dt=1/240. One step is v += g*dt, y += v*dt, and when y < 0 it bounces back with y = -y, v = -0.7*v. Write the results to /root/impulse/out/bounce.csv with the header step,y,v and 2401 lines, and write the peak heights of the first three bounces to /root/impulse/out/04-bounce.txt as peak1=, peak2= and peak3=. Then, on /root/impulse/out/bounce.png (256x256, black background), plot points of color (120,200,255) at px = int(step*255/2400) and py = 255 - int(y/5.0*255).

The peak height of a bounce is the point where y is greater than or equal to both neighboring values (a local maximum). Pick the first three.

Compute the ratios. Both peak2/peak1 and peak3/peak2 should come out near e² = 0.49. This is because when the velocity shrinks by a factor of e, the height it rises to shrinks by a factor of e squared.

In the graph you get the familiar shape of peaks getting lower and lower. This picture shows the meaning of the coefficient of restitution on a single sheet.

Bouncing back with y = -y is a method that reflects the amount of penetration and returns it. It loses less energy than clamping to 0.

Where does a sliding box stop

Apply Coulomb friction (mu=0.3, g=9.8) to a box sliding at speed 5.0 and run it with dt=1/240. One step is v = max(0, v - mu*g*dt), x += v*dt, and when v becomes 0 it has stopped. Write the time taken and the distance moved to /root/impulse/out/05-friction.txt as three lines, steps=, stop_time= and distance=.

Friction does not scale with velocity but gives a constant deceleration. So the velocity decreases linearly, and you can find the time and distance analytically.

t = v0 / (mu * g)        = 5 / 2.94 = 1.7007 s
d = v0² / (2 * mu * g)   = 25 / 5.88 = 4.2517 m

The values you ran directly should match this expression to within one percent. They are not exactly the same, because the discrete steps cut off the last piece.

Without max(0, ...), the velocity goes down to negative and the box goes backward. It amounts to friction pushing the object along, and this is a common mistake in real engines as well.

Friction is an impulse too

In /root/impulse/tangent.py, make bounce_with_friction(v, inv_m, n, e, mu). It handles one object hitting an immovable surface (inverse mass 0). Find and apply the normal impulse jn = -(1+e)*dot(v,n)/inv_m, then for the tangent direction t = normalize(v - dot(v,n)*n) find jt = -dot(v,t)/inv_m, apply the value clamped to [-mu*jn, +mu*jn], and return the final velocity. Return the three values (v_out, jn, jt), and for the case v=(3,-3), inv_m=1, n=(0,1), e=0.5, mu=0.4, write them to /root/impulse/out/06-tangent.txt as vout=vx,vy, jn= and jt=.

The tangent direction is what remains after subtracting the normal component from the velocity. If that remainder is 0 (a head-on hit), the tangent is undefined, so you must skip friction.

Clamping jt is the whole of Coulomb friction. If the value needed to remove the tangential velocity completely does not exceed mu*jn, the object stops on the spot (static friction), and if it does, only the clamped value is applied and it slides (kinetic friction). One conditional gives both frictions.

In the given case, jn=4.5, jt=-3 before clamping, and the upper limit is 0.4*4.5=1.8, so jt=-1.8. The final velocity is (1.2, 1.5).

Compute the normal and tangent components from the velocity before applying jn. If you mix up the order, the values change.

Five balls in a box

Run five balls of radius 0.5 inside a box of side 10 for 1800 steps with g=-9.8, dt=1/240 and a wall restitution coefficient of 0.9. The initial states are, in order, positions (2,8), (5,9), (8,7), (3,5), (7,4) and velocities (3,0), (-2,1), (0,0), (4,2), (-3,-1). Every 10 steps, plot each ball's position on /root/impulse/out/balls.png (256x256, black background, sx=int(x*25.6), sy=int(255-y*25.6)) in the colors (255,100,100), (100,255,100), (100,150,255), (255,220,80), (220,120,255), and write b1=x,y,vx,vy through b5= and energy_start= and energy_end= to /root/impulse/out/07-balls.txt. The energy is sum(0.5*(vx²+vy²) + 9.8*y).

Wall handling is done per axis separately. If x - r < 0, then x = r and vx = -e*vx, and the other side is the same. y is the same.

The order within a step is applying gravity, updating position, and handling walls. If you change the order, the balls go through the walls.

The energy must decrease. The coefficient of restitution is less than 1, so you lose a little every time a ball hits a wall. If it increased, the part of wall handling that puts the position back is missing or a sign is wrong.

The trajectories remain in the picture. You can see the bouncing as parabolic pieces joined together.