Ellipsoid Rolling

Let $\pb^b_1 = (x_1, y_1, z_1)^\top$ and $\pb_2^b = (x_2, y_2, z_2)^\top$ be the contact point coordinates represented in the body frame of ellipsoid 1 and ellipsoid 2. The superscripts $b$ and $w$ denote the vector is represented in body and world frame, respectively.

If the two ellipsoids are in contact with each other, and have only one contact point, we have

where

Take the first order derivative of the above equations and we obtain


By assuming rolling condition at the contact point, we have

where the second step employs the relation that

Thus, by fixing the velocity of the ellipsoid 1, the translational and angular velocity of ellipsoid 2 is determined by

where $\mathcal{N}(Q_2)$ is the null space of $Q_2$ and $\color{red} \lambdab$ is a $3 \times 1$ free vector corresponding to the 3 degrees of freedom motion of ellipsoid 2 around ellipsoid 1 (rolling in tangential plane $\times 2$ and spinning along normal vector $\times 1$).

To let $\omegab_2^w$ align with $z$-axis, we can choose $\lambdab$ by solving

where $\mathcal{N}_{2,\omega} \in \mathbb{R}^{3\times 3}$ is the first three rows of matrix $\mathcal{N}(Q_2)$.


In addition, subtracting $(2)$ from $(1)$ gives

Since we know that

the derivative of the contact position can be written it in a more compact formulation

The equation above can be used to calculate $\pbdot_1^b,\, \pbdot_2^b$ and then update $\pb_1^b,\, \pb_2^b$ by numerical integration.

Matlab Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
clear
dt = 0.0001;

% ellipsoid 1
a1 = 5; b1 = 4; c1 = 3;
C1 = [
1 / (a1 * a1), 0, 0;
0, 1 / (b1 * b1), 0;
0, 0, 1/ (c1 * c1)];
o1 = [0; 0; 0];
R1 = eye(3);
p1 = [5; 0; 0];
omega1 = [0; 0; 0];
v1 = [0; 0; 0];

% ellipsoid 2
a2 = 5; b2 = 4; c2 = 3;
C2 = [
1 / (a2 * a2), 0, 0;
0, 1 / (b2 * b2), 0;
0, 0, 1/ (c2 * c2)];
o2 = [10; 0; 0];
R2 = eye(3);
p2 = [-5; 0; 0];

figure
for i = 1:100000
V1 = [omega1; v1];
Q1 = [-skew(R1 * p1), eye(3)];
Q2 = [-skew(R2 * p2), eye(3)];
N2 = null(Q2);
lambda = N2(1:3, :) \ [0; 0; 1];
V2 = Q2' * inv(Q2 * Q2') * Q1 * V1 + N2 * lambda;
omega2 = V2(1:3,:);
v2 = V2(4:6,:);
o1 = o1 + v1 * dt;
o2 = o2 + v2 * dt;
R1_dot = skew(omega1) * R1;
R2_dot = skew(omega2) * R2;
R1 = R1 + R1_dot * dt;
R2 = R2 + R2_dot * dt;

n1 = C1 * p1;
n2 = C2 * p2;
A = [n1', zeros(1, 3);
zeros(1, 3), n2';
R1, -R2;
R1*C1, R2*C2];
b = [zeros(5, 1);
-R1_dot * n1 - R2_dot * n2];
p_dot = A \ b;
p1 = p1 + p_dot(1:3, :) * dt;
p2 = p2 + p_dot(4:6, :) * dt;

if mod(i, 1000) == 1
clf
draw(o1, R1, o2, R2, p1)
pause(0.1)
end
end

function draw(c1, R1, c2, R2, p)
[x1, y1, z1] = ellipsoid(c1(1),c1(2),c1(3),5,4,3,30);
S1 = surfl(x1, y1, z1);
axang1 = rotm2axang(R1);
rotate(S1, axang1(1:3), rad2deg(axang1(4)));
hold on

[x2, y2, z2] = ellipsoid(c2(1), c2(2), c2(3),5,4,3,30);
S2 = surfl(x2, y2, z2);
axang2 = rotm2axang(R2);
rotate(S2, axang2(1:3), rad2deg(axang2(4)), c2);
ball(p(1), p(2), p(3), 0.5)
axis equal
end

function ball(x, y, z, sz)
[xs,ys,zs] = sphere;
s1 = surf(sz*xs+x,sz*ys+y,sz*zs+z);
set(s1,'facecolor',[.98 .45 .02]);
end

function X = skew(x)
X=[0 -x(3) x(2) ; x(3) 0 -x(1) ; -x(2) x(1) 0 ];
end