Get your own Free Home PageRecently, I derived the 3D dynamics equations with some improvement over the 2D. Basically I added a new physics property "stickiness" 0 to 1 which interpolates between the two alternatives for impact in "Dynamics" (Pestel and Thomson).
3D impact involves solving 12 equations for 12 unknowns (the translational and rotational velocities of both bodies.) Below, you will find a mathematical explanation and I have also written a Java program that implements it. I think you'll find the former more insightful, but the latter more useful. I was able to reduce and optimized the formulas a little more when I wrote them in Java.
Assume two rigid free floating objects (called 1 and 2) impact at the origin; the impact normal is parallel with the X axis. If these assumptions are not correct, then you must translate and rotate your coordinate system to make them true.
Momentum is conserved along each space axis. So we get 3 equations, by multiplying mass and velocity.
The moment of momentum of each body with respect to the origin is conserved. Since angular momentum is 3 dimensional and there are two bodies, this gives 6 equations.
The restitution equation (aka Newton's hypothesis) says that the ratio of the relative velocity in the x direction, before and after collision equals -e, where "e" is the elasticity. This gives 1 equation. (one large ugly equation!)
The remaining 2 equations need to come from our assumptions regarding the change in velocities in the tangent plane. If we assume a perfectly "slippery" impact, then the velocity in the y and z directions are unchanged. (Actually this is four equations, but they are partially redundant with the conservation of momentum equations.)
x, y, z = location of center of mass
v, w = velocity and rotation before impact
c, o = velocity and rotation after impact, these are the 12 unknowns.
I = moment of intertia around each axis
m = mass
s = stickiness [0,1]
E = restitution coeff [0,1]
First we have one for conservation of momentum
m1·vx1 + m2·vx2 = m1·cx1 + m2·cx2
Then we have three for conservation of angular momentum of body 1
wx1·Ix1 + m1·(y1·vz1 - z1·vy1) = ox1·Ix1 + m1·(y1·cz1 - z1·cy1)
wy1·Iy1 + m1·(z1·vx1 - x1·vz1) = oy1·Iy1 + m1·(z1·cx1 - x1·cz1)
wz1·Iz1 + m1·(x1·vy1 - x1·vy1) = oz1·Iz1 + m1·(x1·cy1 - x1·cy1)
Then we have three for conservation of angular momentum of body 2
wx2·Ix2 + m2·(y2·vz2 - z2·vy2) = ox2·Ix2 + m2·(y2·cz2 - z2·cy2)
wy2·Iy2 + m2·(z2·vx2 - x2·vz2) = oy2·Iy2 + m2·(z2·cx2 - x2·cz2)
wz2·Iz2 + m2·(x2·vy2 - x2·vy2) = oz2·Iz2 + m2·(x2·cy2 - x2·cy2)
Then one equation for restitution
(cx1 + y1·oz1 - z1·oy1) - (cx2 + y2·oz2 - z2·oy2) = -1·E·((vx1 + y1·wz1 - z1·wy1) - (vx2 + y2·wz2 - z2·wy2))
Finally we get 4 equations from the stickiness assumptions. There are several options here.
OPTION A. Perfectly slippery, no velocity change in the tangent plane
cy1 = vy1
cz1 = vz1
cy2 = vy2
cz2 = vz2
OPTION B. Perfectly sticky, the tangent velocities become equal for this we must also include conservation of momentum.
cy1 = cy2
cz1 = cz2
m1·vy1 + m2·vy2 = m1·cy1 + m2·cy2
m1·vz1 + m2·vz2 = m1·cz1 + m2·cz2
OPTION C. A parameterized combination of options A and B. "s" stands for the "stickiness" or "friction".
cy1 = (1 - s)·vy1 + s·(m1·vy1 + m2·vy2)/(m1 + m2)
cz1 = (1 - s)·vz1 + s·(m1·vz1 + m2·vz2)/(m1 + m2)
cy2 = (1 - s)·vy2 + s·(m1·vy1 + m2·vy2)/(m1 + m2)
cz2 = (1 - s)·vz2 + s·(m1·vz1 + m2·vz2)/(m1 + m2)
The following is the solution based on OPTION C. I solved for cx1 in terms of the knowns. Using cx1 you can find cx2, and so on backwards to get all 12 unknowns.
cy1 = (1 - s)·vy1 + s·(m1·vy1 + m2·vy2)/(m1 + m2)
cz1 = (1 - s)·vz1 + s·(m1·vz1 + m2·vz2)/(m1 + m2)
cy2 = (1 - s)·vy2 + s·(m1·vy1 + m2·vy2)/(m1 + m2)
cz2 = (1 - s)·vz2 + s·(m1·vz1 + m2·vz2)/(m1 + m2)
ox1 = m1/Ix1·s·((m1·vy1/(m1 + m2) + m2·vy2/(m1 + m2) - vy1)·z1 - (m1·vz1/(m1 + m2) + m2·vz2/(m1 + m2) - vz1)·y1) + wx1
ox2 = m2/Ix2·s·((m1·vy1/(m1 + m2) + m2·vy2/(m1 + m2) - vy2)·z2 - (m1·vz1/(m1 + m2) + m2·vz2/(m1 + m2) - vz2)·y2) + wx2
oy1 = m1/Iy1·((-1·cx1 + vx1)·z1 + (m1·vz1/(m1 + m2) + m2·vz2/(m1 + m2) - vz1)·s·x1) + wy1
oz1 = m1/Iz1·((cx1 - vx1)·y1 + (-1·m1·vy1/(m1 + m2) - m2·vy2/(m1 + m2) + vy1)·s·x1) + wz1
oy2 = m2/Iy2·((-1·cx2 + vx2)·z2 + (m1·vz1/(m1 + m2) + m2·vz2/(m1 + m2) - vz2)·s·x2) + wy2
oz2 = m2/Iz2·((cx2 - vx2)·y2 + (-1·m1·vy1/(m1 + m2) - m2·vy2/(m1 + m2) + vy2)·s·x2) + wz2
cx2 = m1/m2·(-1·cx1 + vx1) + vx2
EE = -1·E·(vx1 - vx2 - wy1·z1 + wy2·z2 + wz1·y1 - wz2·y2)
P = (z1^2/Iy1 + z2^2/Iy2 + y1^2/Iz1 + y2^2/Iz2)·m1·m2
M = m1 + m2
General solution for cx1
(M + P)·cx1 = (EE + wy1·z1 - wy2·z2 - wz1·y1 + wz2·y2)·m2 + m1/M·m2^2·s·((-1·vz1 + vz2)·x1·z1/Iy1 + (-1·vz1 + vz2)·x2·z2/Iy2 + (-1·vy1 + vy2)·x1·y1/Iz1 + (-1·vy1 + vy2)·x2·y2/Iz2) + P·vx1 + m1·vx1 + m2·vx2
Special case s=1
s = 1
(M + P)·cx1 = (EE + wy1·z1 - wy2·z2 - wz1·y1 + wz2·y2)·m2 + m1/M·m2^2·((-1·vz1 + vz2)·x1·z1/Iy1 + (-1·vz1 + vz2)·x2·z2/Iy2 + (-1·vy1 + vy2)·x1·y1/Iz1 + (-1·vy1 + vy2)·x2·y2/Iz2) + P·vx1 + m1·vx1 + m2·vx2
Special case s=0
s = 0
(M + P)·cx1 = (EE + wy1·z1 - wy2·z2 - wz1·y1 + wz2·y2)·m2 + P·vx1 + m1·vx1 + m2·vx2
Special case s=0 and E=1
(M + P)·cx1 = P·vx1 + m1·vx1 - m2·vx1 + 2·m2·(vx2 + wy1·z1 - wy2·z2 - wz1·y1 + wz2·y2)
dimensions -> 2 3 4 5 6 ... n
----------------------------------------------------
conserv. momentum 2 3 4 5 6 n
conserv. ang. mom. 2 6 12 20 36 n(n-1)
restitution 1 1 1 1 1 1
slipperiness 1 2 3 4 5 n-1
A "rotation" is a sine-cosine transformation using two axes. So in N dimensional space there are "N choose 2" or n(n-1)/2 possible rotations. This was rather counter intuitive for me (but then, every thing is in 4D) because in 3D, rotation is represented as a 3D vector, so in 4D shouldn't rotation be represented by a 4D vector? NOT SO! In 4D space, rotation is 6 dimensional! Therefore to compute collisions in 4D you must solve for 20 unknowns! Fortunately, as you can see in the table above, we have exactly 20 equations.
Get your own Free Home Page