Classical Mechanics: Problem Set 2

This illustrates what happens when you flip a book about its unstable axis (basically end over end). Had you flipped the book about any of the other two axes, you would have found that the dynamics are stable. What we find in this problem set is that the axis with the 'medium' moment of inertia is unstable. In fact, the conservation of energy equations indicate that the angular velocity vector must be constrained to lie on the surface of an ellipsoid.

1. Axis of symmetry

If a body has a plane of symmetry passing through the centroid, show that the normal to the plane of symmetry is a principal axis of the inertia matrix at the centroid.

Spoiler

2. Another route to the angular velocity vector

Suppose that $(\mb{e}_1, \mb{e}_2, \mb{e}_3)$ is a time-dependent orthonormal basis. Explain why $\dot{\mb{e}}_1 \cdot \mb{e}_1 = 0$ and give two similar identities. Deduce from these that for some $\alpha, \beta, \gamma, \lambda, \mu, \nu$, \[ \dot{\mb{e}}_1 = \beta \mb{e}_3 - \gamma \mb{e}_2, \quad \dot{\mb{e}}_2 = \lambda \mb{e}_1 - \alpha \mb{e}_3, \quad \dot{\mb{e}}_3 = \mu \mb{e}_2 - \nu \mb{e}_1. \]

Now show that $\dot{\mb{e}_1} \cdot \mb{e}_2 + \dot{\mb{e}}_2 \cdot \mb{e}_1 = 0$, and deduce from this and similar identities that $\lambda = \mu$, $\mu = \alpha$, $\beta = \nu$. Deduce that there exists $\omega$ such that for all $i$, $\dot{\mb{e}}_i = \omega \times \mb{e}_i$.

Spoiler

3. Angular velocity in a rotating frame

Using the same definitions of $\mb{e}_i$, $\hat{\mb{e}}_i$, $H_{ij}$, and $\omega$ as in lectures, show that the components of $\omega$ in the $\hat{\mb{e}}_i$ frame are given by \[ \frac{1}{2} \epsilon_{ijk} (H^T \dot{H})_{jk}. \]

Spoiler

4. Special solutions of Euler's equations

Consider the system of Euler's equations: \begin{align*} A\dot{\omega}_1 &= (B - C) \omega_2 \omega_3 \\ B\dot{\omega}_2 &= (C - A) \omega_3 \omega_1 \\ C\dot{\omega}_3 &= (A - B) \omega_1 \omega_2. \end{align*}

Show that there are always solutions of the form $(\omega_1, \omega_2, \omega_3) = (\Omega, 0, 0)$ or $(0, \Omega, 0)$ or $(0, 0, \Omega)$, for a constant $\Omega$. If $A < B < C$, show that two of these motions are stable, and one is unstable.

Spoiler

5. Stability of free rotations

Again from Euler's equations, consider the case when $A < B$, $C = A + B$, and the body is set in motion with $\omega_2 = 0$ and $\sqrt{A + B} \omega_3 = \sqrt{B-A}\omega_1$.

Using the constants $T$ and $L^2$ of the motion, show that \[ 2BT = L^2, \qquad 2CT - L^2 = AB(\omega_1^2 + \omega_2^2. \]

Hence show that \[ A^2 (\dot{\omega}_1)^2 = (B - C)^2 \omega_2^2 \omega_3^2 = A^2 \left(\frac{B-A}{A+B}\right)^2\omega_1^2 \left(\frac{2T}{B} - \omega_1^2\right), \]

and find a solution of the form $\omega_1 = c_1 \text{sech}(c_2 t)$. What happens to $\omega_2$ as $t\to \infty$.

Spoiler

Matlab code for a flipped book

The following Matlab code solves the Euler equations to derive the rotations of a book. The trick is that we need to use the values of $\omega_1$, $\omega_2$, $\omega_3$ to obtain the rotation matrix. To be more specific, we solve Euler's equations simultaneously with \begin{align*} \dot{e}_1 &= \omega \times \mb{e}_1 \\ \dot{e}_2 &= \omega \times \mb{e}_2 \\ \dot{e}_3 &= \omega \times \mb{e}_3 \end{align*}

which indicate the directions of the principle axes. Once we have obtained these values, we normalize them, and then we define the transformation matrix by \[ A = [\mb{e}_1 \ \mb{e}_2 \ \mb{e}_3]. \]

Let the original location of a point in the book be $\mb{r}(0) = x_1 \mb{i} + x_2 \mb{j} + x_3\mb{k}$. This is then mapped to \[ \mb{r}(t) = A\mb{r}(0). \]

1
% A simplified code for the rotation of a flipped book
% Written April 2013
% PHT, Oxford 2013
function rigidrotsimp2
 
    close all
 
    % Principle moments of inertia
    A = 1; B = 2; C = 3;
 
    % Flip about x with small perturbation
    %winit = [0.99; sqrt(1-0.99^2); 0];
    % Flip about y with small perturbation
    winit = [sqrt(1-0.99^2); 0.99; 0];
 
    % Kinetic energy
    T = 1/2*(A*winit(1)^2 + B*winit(2)^2 + C*winit(3)^2);
 
    % This is the initial orientation of the book
    e10 = [1 0 0]; e20 = [0 1 0]; e30 = [0 0 1];
 
    % We for the three components of angular velocity and three vectors for
    % the principle axes (i.e. a system of 12 variables)
    solinit = [winit; e10(:); e20(:); e30(:)];
    [t Y] = ode45(@(t,W)eulerde(t, W, A, B, C), [0 50], solinit);
 
    % Assign variables to the solution vector
    w1 = Y(:,1); w2 = Y(:,2); w3 = Y(:,3);
    e1 = Y(:,4:6); e2 = Y(:,7:9); e3 = Y(:,10:12);
 
    % Define the faces of the book
    fm = [1 2 6 5; 2 3 7 6; 3 4 8 7; 4 1 5 8; 1 2 3 4; 5 6 7 8];
    % Define the corners of the book
    vm = [-A -B -C; A -B -C; A B -C; -A B -C;
          -A -B C; A -B C; A B C; -A B C];
 
    figure(1); clf(1)
    for j = 1:length(t)
        clf(1)
        title(['Time = ', num2str(t(j))]);
 
        %Normalize the new principle axes
        e1(j,:) = e1(j,:)/norm(e1(j,:));
        e2(j,:) = e2(j,:)/norm(e2(j,:));
        e3(j,:) = e3(j,:)/norm(e3(j,:));
 
        % Find a rotation matrix that tells us how (i, j, k) are mapped to
        % the current configuration 
        M = rotmat(e1(j,:), e2(j,:), e3(j,:));
 
        plot3([0 A*e1(j,1)], [0 A*e1(j,2)], [0 A*e1(j,3)], 'b');
        hold on
        plot3([0 B*e2(j,1)], [0 B*e2(j,2)], [0 B*e2(j,3)], 'r');
        plot3([0 C*e3(j,1)], [0 C*e3(j,2)], [0 C*e3(j,3)], 'g');
 
        % The new vertices are given by M times the original vertex
        % locations 
        newvm = (M*vm.').';
        patch('Vertices', newvm, 'Faces', fm, ...
                            'FaceVertexCData', hsv(6), ...
                            'FaceColor', 'flat', 'FaceAlpha', 0.7);
        hold off
        axis([-4 4 -4 4 -4 4]);
        axis square;
        box on;
        title(['Time = ', num2str(t(j))]);
        drawnow
        pause(0.1)
    end
end
 
function F = eulerde(t, Y, A, B, C)
 
    w1 = Y(1);
    w2 = Y(2);
    w3 = Y(3);
    e1 = Y(4:6);
    e2 = Y(7:9);
    e3 = Y(10:12);
 
    % Euler equations for w1, w2, w3
    w1p = (B-C)/A*w2*w3;
    w2p = (C-A)/B*w3*w1;
    w3p = (A-B)/C*w1*w2;
 
    % Evolution of unit vectors, using dr/dt = w x r
    e1p = w3*e2 - w2*e3;
    e2p = w1*e3 - w3*e1;
    e3p = w2*e1 - w1*e2;
 
    F = [w1p; w2p; w3p; e1p(:); e2p(:); e3p(:)];
end
 
function M = rotmat(e1, e2, e3)
    % Transformation matrix is 
    M = [e1(:) e2(:) e3(:)];
 
end