Chapter 6

Rigid Body Dynamics

 

 

 

 

 

6.5 Rotational forces $–$ review of moments exerted by forces and torques

 

You can find a detailed discussion of forces and moments, with lots of examples, in Section 2 of these notes.   Moments and torques don’t come up very often in particle dynamics, but play a very important role in rigid body dynamics.   We therefore review the most important concepts related to torques and moments here.

 

You need to remember, and understand, these ideas:

(1)    A moment is a generalized force that causes an object to rotate (see section 2).

(2)    A force can exert a moment on a rigid body.   The moment of a force (about the origin) is defined as

$M=r\times F$

(3)    In general, a force causes a rigid body to accelerate, and will also induce an angular acceleration (so it influences both translational and rotational motion).

(4)    A ‘torque’ or ‘pure moment’ is a special kind of generalized force that causes an object to rotate, but has no effect on its translational motion.  As an example, a motor shaft (eg the bit on a power-driven screwdriver!) will exert a torque on the object connected to it. 

(5)    A torque or pure moment is a vector quantity $–$ it has magnitude and direction.  The direction indicates the axis associated with its rotational force (following the right hand screw convention); the magnitude represents the intensity of the rotational force. The magnitude of a torque has units of Newton Meters. A moment is often denoted by the symbols shown in the figure

 

6.5.1 Rate of work done by a torque or moment: If a torque $Q=Q_{x}i+Q_{y}j+Q_{z}k$ acts on an object that rotates with angular velocity $\omega$, the rate of work done on the object by Q is

$P=Q\cdot \omega =Q_x{\omega }_x+Q_y{\omega }_y+Q_z{\omega }_z$

 

 

6.5.2 Torsional springs

A solid rod is a good example of a torsional spring. You could take hold of the ends of the rod and twist them, causing one end to rotate relative to the other.   To do this, you would apply a moment or a couple to each end of the rod, with direction parallel to the axis of the rod.   The angle of twist increases with the moment.  Various torsion spring designs used in practice are shown in the picture $–$ the image is from

 http://www.mollificio.lombardo.molle.com/springs/torsion_springs.html

 

More generally, a torsional spring resists rotation, by exerting equal and opposite moments on objects connected to its ends.  For a linear spring the moment is proportional to the angle of rotation applied to the spring. 

 

The figure shows a formal free body diagram for two objects connected by a torsional spring.  If object A is held fixed, and object B is rotated through an angle $\theta$ about an axis parallel to a unit vector n, then the spring exerts a moment

$Q=-\kappa \theta n$

on object B where $\kappa$ is the torsional stiffness of the spring.   Torsional stiffness has units of Nm/radian.

 

The potential energy of the moments exerted by the spring can be determined by computing the work done to twist the spring through an angle $\theta$. 

1. The work done by a moment Q due to twisting through a very small angle $d\theta$ about an axis parallel to a vector n is

$dW=Q\cdot d\theta n$

2. The potential energy is the negative of the total work done by M, i.e.

$V=-{\displaystyle \underset{0}{\overset{\theta }{\int }}Q\cdot d\theta n}=-{\displaystyle \underset{0}{\overset{\theta }{\int }}\left(-\kappa \theta n\right)\cdot d\theta n}={\displaystyle \underset{0}{\overset{\theta }{\int }}\kappa \theta d\theta }=\frac{1}{2}\kappa {\theta }^2$

 

A potential energy cannot usually be defined for most concentrated moments, because rotational motion is itself path dependent (the orientation of an object that is given two successive rotations depends on the order in which the rotations are applied).

 

 

6.6 Dynamics of rigid bodies

 

We predict the position and velocity of a particle by integrating F=ma.    For a rigid body, we need to predict both its position and orientation.   We use the following equations to do this. 

 

The figure shows a rigid body subjected to several forces $F^{(i)}$  and torques (pure moments) $Q^{(i)}$ .  During a representative time interval $t_0<t<t_1$ the forces exert a linear impulse $\Im$ and angular impulse A, and do total work on the rigid body $\Delta W$ .

 

The body has total mass M and mass moment of inertia $I_G$ about the center of mass. 

 

Let $r_G,v_G,a_G$ denote the position, velocity and acceleration of the center of mass, and let $\omega ,\alpha$ denote the angular velocity and acceleration.  

 

The linear and angular momentum (about the origin) of the rigid body follow as $p=M{v}_G$ , $h=M{r}_G\times v_G+I_G\omega$ , and its kinetic energy is $T=\frac{1}{2}M{\left|v_G\right|}^2+\frac{1}{2}\omega \cdot I_G\omega$ .

 

The equations of motion are then

 

 

Force-acceleration relation:  ${{\displaystyle \sum _{i}F}}^{(i)}=M{a}_G$

 

Moment $–$ angular velocity/acceleration relation   ${{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=M{r}_G\times a_G+I_G\alpha }+\omega \times \left[I_G\omega \right]$

 

Force-momentum and impulse-momentum relation: ${{\displaystyle \sum _{i}F}}^{(i)}=\frac{dp}{dt}\Im =p_1-p_0$

 

 

Moment $–$ angular momentum relation: ${{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}}=\frac{dh}{dt}A=h_1-h_0$

 

Power $–$ work $–$ kinetic energy relation ${{\displaystyle \sum _{i}F}}^{(i)}\cdot v^{(i)}+{\displaystyle \sum _j{Q}^{(j)}\cdot \omega }=\frac{dT}{dt}\Delta W=T_1-T_0$

 

 

 

For 2D planar motion we can use the simplified formulas

$\begin{array}{l}{{\displaystyle \sum _{i}F}}^{(i)}=M{a}_G\\ {{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=M{r}_G\times a_G+I_{Gzz}{\alpha }_z}k\end{array}$

Derivations:  It is possible to obtain the equations of motion for a rigid body from Newton’s laws for a particle $–$ the basic idea is to assume that a rigid body consists of an infinite number of particles connected by rigid massless links $–$ but this isn’t really a rigorous proof, because we have to assume that the links are two-force members, and there is no way to prove that this is a realistic description of matter.   Another viewpoint is to accept conservation of linear momentum and angular momentum as two separate physical laws (the linear momentum is just Newton’s law, and the angular momentum equation is sometimes referred to as Euler’s law). We can then ‘prove’ that a rigid body can be represented as a bunch of particles connected by two force members.   We’ll show the first approach here.

 

The figure shows a system of particles connected by rigid massless links.  The length of the link between the ith and jth particle will be denoted by $L_{ij}$ We assume that all the links are two-force members.

 

The particles are subjected to a set of external forces $F^{(i)}$ .   We denote the magnitude of the force in the member connecting the ith and jth particle by $R_{ij}$ (by convention a positive $R_{ij}$ represents an attractive force between the particles).  Note that the $R_{ij}=R_{ji}$ because the two  particles exert equal and opposite forces on each other.    The vector valued force exerted on the ith particle by the jth follows as

$R_{ij}=R_{ij}\frac{r_j-r_i}{L_{ij}}$

(to see this note that $(r_j-r_i)/L_{ij}$ is a unit vector from the ith to the jth particles)

 

We can start the derivation with the force-linear momentum relation for a single particle.   For example, for the ith particle (see section 4 of the notes)

$F^i+{\displaystyle \sum _{j\ne i}{R}_{ij}\frac{(r_j-r_i)}{L_{ij}}}=\frac{d}{dt}{m}_i{v}_i$

Sum this over all particles

 

${\displaystyle \sum _i{F}^i}+{\displaystyle \sum _i{\displaystyle \sum _{j\ne i}{R}_{ij}\frac{(r_j-r_i)}{L_{ij}}}}={\displaystyle \sum _i\frac{d}{dt}{m}_i{r}_i}$

But we know that ${\displaystyle \sum _i{m}_i{v}_i}=M{v}_G$, and since $R_{ij}=R_{ji},L_{ij}=L_{ji}$  the second term on the left hand side is zero.   Therefore

${\displaystyle \sum _i{F}^i}=M\frac{d{v}_G}{dt}=M{a}_G=\frac{dp}{dt}$

Since this is independent of the number of particles, it must also apply to a rigid body.  This shows that the force-momentum and force-acceleration for a rigid body can be derived from Newton’s law for a particle.

 

We can derive the angular momentum relation for a rigid body using the same idea.  For one particle we have the angular momentum equation

$r_i\times F^{(i)}+r_i\times {\displaystyle \sum _{j\ne i}{R}_{ij}\frac{(r_j-r_i)}{L_{ij}}}=r_i\times F^{(i)}+r_i\times {\displaystyle \sum _{j\ne i}{R}_{ij}\frac{r_j}{L_{ij}}}=\frac{d{h}_i}{dt}$

where we have noted that $r_i\times r_i=0$ .   We can sum this over all the particles

${\displaystyle \sum _i{r}_i\times F^{(i)}}+{\displaystyle \sum _i{\displaystyle \sum _{j\ne i}{R}_{ij}\frac{r_i\times r_j}{L_{ij}}}}=\frac{d}{dt}{\displaystyle \sum _i{h}_i}$

The second term here is zero, because $r_i\times r_j=-r_j\times r_i$ and $R_{ij}=R_{ji}{L}_{ij}=L_{ji}$ (just write out the sum term by term for some finite number of particles $–$ eg two $–$ if you don’t see this).   The term on the right hand side is clearly just the total angular momentum of the system.  If we replace some subset of the forces with a statically equivalent torque and force, we obtain the moment-angular momentum equation.

 

 

 

 

6.7 Summary of equations of motion for rigid bodies

 

In this section, we collect together all the important formulas from the preceding sections, and summarize the equations that we use to analyze motion of a rigid body.

 

 

We consider motion of a rigid body that has mass density $\rho$ during some time interval $t_0<t<t_1$ , and define the following quantities:

 

 

6.7.1        Forces, torques, impulse, work, power

 The total force acting on the body ${\displaystyle \sum _i^{}{F}^{(i)}}$ 

  The total linear impulse exerted by forces during the time interval $\Im ={\displaystyle \underset{t_0}{\overset{t_1}{\int }}{\displaystyle \sum _i{F}^{(i)}(t)dt}}$ 

 The total moment (including torques) acting on the body  ${{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}}$

 The tot al angular impulse exerted on the body during the time interval $A={\displaystyle \underset{t_0}{\overset{t_1}{\int }}\left({\displaystyle \sum _i{r}_i\times F^{(i)}(t)+{\displaystyle \sum _j{Q}^{(j)}}}\right)}dt$

 The rate of work done by forces and torques acting on the body  $P={\displaystyle \sum _i{F}^{(i)}\cdot v_i}+{\displaystyle \sum _j{Q}^{(j)}\cdot \omega }$ 

 The total work done by forces and torques on the body during the time interval $W={\displaystyle \underset{t_0}{\overset{t_1}{\int }}P(t)dt}$

 

 

6.7.2 Inertial properties

 The total mass is $M={\displaystyle \underset{V}{\int }\rho dV}$ 

 The position of the center of mass is $r_G=\frac{1}{M}{\displaystyle \underset{V}{\int }r\rho dV}$

 

 The mass moment of inertia about the center of mass

$I_G={\displaystyle \underset{V}{\int }\rho \left[\begin{array}{ccc}{d}_y^2+d_z^2& -d_x{d}_y& -d_x{d}_z\\ -d_x{d}_y& d_x^2+d_z^2& -d_y{d}_z\\ -d_x{d}_z& -d_y{d}_z& d_x^2+d_y^2\end{array}\right]dV}$

where $d=r-r_G$

 

 

 

For a 2D body with mass per unit area $\mu$ we use

 

 The total mass is $M={\displaystyle \underset{A}{\int }\mu dA}$ 

 The position of the center of mass is $r_G=\frac{1}{M}{\displaystyle \underset{A}{\int }r\mu dA}$

 The mass moment of inertia about the center of mass is $I_{Gzz}=\frac{1}{M}{\displaystyle \underset{A}{\int }\mu (d_x^2+d_y^2)dA}$

where $d=r-r_G$

 

 

 

6.7.3 Describing motion

 

 The rotation tensor (matrix) maps the vector connecting two points in a solid before it moves to its position after motion

$r_B-r_A=R(p_B-p_A)$

 The spin tensor is related to R by  $W=\frac{dR}{dt}{R}^T\frac{dR}{dt}=WR$ 

 Rotation through an angle $\theta$ about an axis parallel to a unit vector

$n=n_{x}i+n_{y}j+n_{z}k$  is

$R=\left[\begin{array}{ccc}\cos \theta +(1-\cos \theta )n_x^2& (1-\cos \theta )n_x{n}_y-\sin \theta n_z& (1-\cos \theta )n_x{n}_z+\sin \theta n_y\\ (1-\cos \theta )n_x{n}_y+\sin \theta n_z& \cos \theta +(1-\cos \theta )n_y^2& (1-\cos \theta )n_y{n}_z-\sin \theta n_x\\ (1-\cos \theta )n_x{n}_z-\sin \theta n_y& (1-\cos \theta )n_y{n}_z+\sin \theta n_x& \cos \theta +(1-\cos \theta )n_z^2\end{array}\right]$

 The angular velocity vector $\omega ={\omega }_{x}i+{\omega }_{y}j+{\omega }_{z}k$ is related to W by

$W=\left[\begin{array}{ccc}0& -{\omega }_z& {\omega }_y\\ {\omega }_z& 0& -{\omega }_x\\ -{\omega }_y& {\omega }_x& 0\end{array}\right]$

 

 The angular acceleration vector is $\alpha =\frac{d\omega }{dt}$ 

 

 The velocities of two points A and B in a rotating rigid body are related by

$v_B-v_A=\omega \times \left(r_B-r_A\right)$ 

 

 The accelerations of  A and B are related by

$a_B-a_A=\alpha \times (r_B-r_A)+\omega \times (v_B-v_A)=\alpha \times (r_B-r_A)+\omega \times \left[\omega \times (r_B-r_A)\right]$

 

 

6.7.4 Momentum and Energy

 

 The total linear momentum is $p=M{v}_G$

 

 The angular momentum (about the origin) is $h=r_G\times M{v}_G+I_G\omega$

 

 The total kinetic energy is $T=\frac{1}{2}M{v}_G\cdot v_G+\frac{1}{2}\omega \cdot I_G\omega$

 

For 2D planar problems, we know $\omega ={\omega }_{z}k$ .   In this case, we can use

 The total linear momentum is $p=M{v}_G$

 The total angular momentum (about the origin) is $h=r_G\times M{v}_G+I_{Gzz}{\omega }_{z}k$

 The total kinetic energy is $T=\frac{1}{2}M{\left|v_G\right|}^2+\frac{1}{2}{I}_{Gzz}{\omega }_z^2$

 

 

6.7.5 Conservation laws

 

 Linear momentum        ${{\displaystyle \sum _{i}F}}^{(i)}=\frac{dp}{dt}\Im =p_1-p_0$

 

 Angular momentum  ${{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}}=\frac{dh}{dt}A=h_1-h_0$

 

 Work-Power - Kinetic Energy relation  ${{\displaystyle \sum _{i}F}}^{(i)}\cdot v^{(i)}+{\displaystyle \sum _j{Q}^{(j)}\cdot \omega }=\frac{dT}{dt}\Delta W=T_1-T_0$

 

 Energy equation for a conservative system $\frac{d}{dt}(T+V)=0T_0+V_0=T_1+V_1$ 

 

 

 

6.7.6 Linear and angular momentum equations in terms of accelerations

 

The linear and angular momentum conservation equations can also be expressed in terms of accelerations, angular accelerations, and angular velocities.   The results are

 

$\begin{array}{l}{{\displaystyle \sum _{i}F}}^{(i)}=M{a}_G\\ {{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=M{r}_G\times a_G+I_G\alpha }+\omega \times \left[I_G\omega \right]\end{array}$

 

 

 

For 2D planar motion we can use the simplified formulas

$\begin{array}{l}{{\displaystyle \sum _{i}F}}^{(i)}=M{a}_G\\ {{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=M{r}_G\times a_G+I_{Gzz}{\alpha }_z}k\end{array}$

 

 

6.7.7 Special equations for analyzing bodies that rotate about a stationary point

 

We often want to predict the motion of a system that rotates about a fixed pivot $–$ a pendulum is a simple example.   These problems can be solved using the equations in 6.6.5 and 6.6.6, but can also be solved using a useful short-cut.

 

For an object that rotates about a fixed pivot at the origin:

 

 The total angular momentum (about the origin) is $h=I_O\omega$

 The total kinetic energy is $T=\frac{1}{2}\omega \cdot I_O\omega$

 The equation of rotational motion is ${{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=I_O\alpha }+\omega \times \left[I_O\omega \right]$

Here $I_O$ is the mass moment of inertia about O (calculated, eg, using the parallel axis theorem)

 

For 2D rotation about a fixed point at the origin we can simplify these to

 The total angular momentum (about the origin) is $h=I_{Ozz}{\omega }_{z}k$

 The total kinetic energy is $T=\frac{1}{2}{I}_{Ozz}{\omega }_z^2$

 The equation of rotational motion is ${{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=I_{Ozz}{\alpha }_z}k$

 

Proof:  It is straightforward to show these formulas.  Let’s show the two dimensional version of the kinetic energy formulas as an example. For fixed axis rotation, we can use the rigid body formulas to calculate the velocity of the center of mass (O is stationary and at the origin)

$v_G=\omega \times r_G={\omega }_{z}k\times r_G$

The general formula for kinetic energy can therefore be re-written as

$\begin{array}{l}T=\frac{1}{2}M{v}_G\cdot v_G+\frac{1}{2}{I}_{Gzz}{\omega }_z^2=\frac{1}{2}M{\omega }_z^2(k\times r_G)\cdot (k\times r_G)+\frac{1}{2}{I}_{Gzz}{\omega }_z^2\\ =\frac{1}{2}\left(M{\left|r_G\right|}^2+I_{Gzz}\right){\omega }_z^2=\frac{1}{2}{I}_{Ozz}{\omega }_z^2\end{array}$

The other formulas can be proved with the same method $–$ we simply express the velocity or acceleration of the COM in the general formulas in terms of angular velocity and acceleration, and notice that we can re-arrange the result in terms of the mass moment of inertia about O.

 

 

6.8 Examples of solutions to problems involving motion of rigid bodies

 

The best way to learn how to use the equations in section 6.6 is just to work through a series of examples. 

 

 

6.8.1 Solutions to 2D problems

 

Example 1:  A solid of revolution (eg a cylinder or sphere) with mass M and mass moment of inertia about its COM $I_{Gzz}$  is released from rest at the top of a ramp.  It rolls without slip.   Calculate its velocity at the bottom of the ramp.

 

• The system is conservative, so we can solve the problem using energy conservation.  The energy equation tells us that the sum of kinetic and potential energy of the cylinder is constant: $T_0+V_0=T_1+V_1$
• We can take the datum for potential energy to be the position of the COM at the bottom of the ramp. The initial potential energy is therefore $V_0=Mgh$ ; the final potential energy is zero.
• The initial kinetic energy is zero, because the cylinder is stationary.   The final kinetic energy is  $T=\frac{1}{2}M{v}_x^2+\frac{1}{2}{I}_{Gzz}{\omega }_z^2$ .  
• The energy equation gives  $T_0+V_0=T_1+V_1\Rightarrow Mgh=\frac{1}{2}M{v}_x^2+\frac{1}{2}{I}_{Gzz}{\omega }_z^2$ 

 

• Finally, since the cylinder rolls without slip, we know that $v_x=-R{\omega }_z$ .

 

 

Hence

$\begin{array}{l}2Mgh=M{v}_x^2+\frac{I_{Gzz}}{R^2}{v}_x^2\\ \Rightarrow v_x=\sqrt{\frac{2gh}{1+I_{Gzz}/(M{R}^2)}}\end{array}$

 

This formula predicts that an object with a smaller inertia $I_{Gzz}$ will move faster than an object with a large inertia.   A sphere rolls down the ramp more quickly than a cylinder, for example, and a solid cylinder rolls more quickly than a ring.

 

 

Example 2:  For the problem treated in the preceding section, calculate the critical value of friction coefficient necessary to prevent slip at the contact.

 

If we want to learn about forces, we have to use the linear and angular momentum equations.   This problem can be solved with the 2D formulas in terms of accelerations:

$\begin{array}{l}{{\displaystyle \sum _{i}F}}^{(i)}=M{a}_G\\ {{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=M{r}_G\times a_G+I_{Gzz}{\alpha }_z}k\end{array}$

• The figure shows a free body diagram for the cylinder (or sphere)
• We know that the COM is always a constant height above the ramp, so the acceleration must be parallel to i.   The linear momentum equation gives

$(Mg\sin \alpha -T)i+(N-Mg\cos \alpha )j=M{a}_{Gx}i$

• We can use the angular momentum equation $–$ it is convenient to take moments about the contact point C.   (There are no torques in this problem).

$\begin{array}{l}{{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=M{r}_G\times a_G+I_{Gzz}{\alpha }_z}k\\ \Rightarrow -RMg\sin \alpha k=Rj\times a_{Gx}i+I_{Gzz}{\alpha }_{z}k=-MR{a}_{Gx}k+I_{Gzz}{\alpha }_{z}k\end{array}$

• Finally, we can use the rolling wheel formula for accelerations $a_{Gx}=-R{\alpha }_z$ .
• The preceding results give:

$\begin{array}{l}-RMg\sin \alpha =-MR{a}_{Gx}-I_{Gzz}\frac{a_{Gx}}{R}\\ \Rightarrow a_{Gx}=\frac{MgR\sin \alpha }{MR+(I_{Gzz}/R)}=\frac{g\sin \alpha }{1+I_{Gzz}/(M{R}^2)}\end{array}$

• Finally, substituting back into the i components of (1):

$\begin{array}{l}T=Mg\sin \alpha -M{a}_{Gx}\\ =Mg\sin \alpha -\frac{Mg\sin \alpha }{1+I_{Gzz}/(M{R}^2)}=\frac{I_{Gzz}/(M{R}^2)}{1+I_{Gzz}/(M{R}^2)}Mg\sin \alpha \end{array}$

• The j component of (1) gives $N=Mg\cos \alpha$ 
• For no slip $\left|T\right|\le \mu N\Rightarrow$ $\mu \ge \frac{I_{Gzz}/(M{R}^2)}{1+I_{Gzz}/(M{R}^2)}\tan \alpha$

 

The formula shows that objects with large values of $I_{Gzz}/M{R}^2$ are more likely to slip.   If the inertia is very small, slip will never occur.   A ring will slip on a lower slope than a cylinder, which will slip on a lower slope than a sphere.

 

Example 3:  A vertical mast can be idealized as a slender rod with length L and mass M, which is held in an inverted position by a torsional spring with stiffness $\kappa$ at its base.   Find the equation of motion for the angle $\theta$ in the figure, and hence determine the natural frequency of vibration of the mast.

 

This is a conservative system.  Also, the mast rotates about a fixed point.   We can analyze the problem using energy methods, and use the special formulas for rotation about a fixed point.

 

·         The kinetic energy formula for planar motion is

$T=\frac{1}{2}{I}_{Ozz}{\omega }_z^2$

·         For planar motion we know that

${\omega }_z=\frac{d\theta }{dt}$

·         We can use the parallel axis theorem to calculate the mass moment of inertia of a rod about one end:

$I_{Ozz}=I_{Gzz}+M{d}^2=\frac{1}{12}M{L}^2+M{\left(\frac{L}{2}\right)}^2=\frac{1}{3}M{L}^2$

·         Gravity and the torsional spring both contribute to the total potential energy of the system.  The total potential energy is

$V=Mg\frac{L}{2}\cos \theta +\frac{1}{2}\kappa {\theta }^2$

·         Energy conservation means that

$\begin{array}{l}T+V=const\Rightarrow \frac{d}{dt}\left(T+V\right)=0\\ \Rightarrow \frac{d}{dt}\left(\frac{1}{2}{I}_{Ozz}{\omega }_z^2+MgL\cos \theta +\frac{1}{2}\kappa {\theta }^2\right)=0\\ \Rightarrow I_{Ozz}\frac{d{\omega }_z}{dt}{\omega }_z-MgL\sin \theta \frac{d\theta }{dt}+\kappa \theta \frac{d\theta }{dt}=0\end{array}$

·         Recall that ${\omega }_z=\frac{d\theta }{dt}$ so

$\Rightarrow I_{Ozz}\frac{d^2\theta }{d{t}^2}-MgL\sin \theta +\kappa \theta =0$

·         We assume that $\theta$ is small enough that $\sin \theta \approx \theta$ , so

$\Rightarrow \frac{I_{Ozz}}{\kappa -MgL}\frac{d^2\theta }{d{t}^2}+\theta =0$

·         This is a standard ‘Case I’ undamped vibration EOM, so we can just read off the natural frequency

${\omega }_n=\sqrt{\frac{\kappa -MgL}{I_{Ozz}}}=\sqrt{\frac{3(\kappa -MgL)}{M{L}^2}}$

 

Example 4: A thin uniform disk of radius R, mass m and mass moment of inertia $m{R}^2/2$ is placed on the ground with a positive velocity $v_0$ in the horizontal direction, and a counterclockwise rotational velocity (a backspin) ${\omega }_0$.  The contact between the disk and the ground has friction coefficient $\mu$. The disk initially slips on the ground, and for a suitable range of values of ${\omega }_0$ and $v_0$ its direction of motion may reverse.  The goal of this problem is to calculate the conditions where this reversal will occur.

 

 

General discussion of slipping contacts:  Solving problems with sliding at a contact is always tricky, because we have to draw the friction forces in the correct direction.   Before tackling the example, we will summarize the general rules.  We will consider a wheel as an example, but the rules apply to contact between any object and a stationary surface.   The figure shows a wheel that spins with angular velocity $\omega ={\omega }_{z}k$ while the center moves with speed $v_O=v_{ox}i$ .   The direction of the friction force is determined by the direction of motion of the point on the wheel that instantaneously touches the ground, which can be calculated from the formula

$v_C=(v_{Ox}+{\omega }_{z}R)i$

Friction always acts to try to bring point C to rest $–$ if C is moving to the right, friction acts to the left; if C is moving to the left, friction acts to the right.

 

There are three possible cases:

 

Forward slip: $v_{Ox}+{\omega }_{z}R>0$ Point C moves in the positive i direction over the ground

 

 

 

 

 Slip occurs at the contact,

 We have to use the friction law $T=\mu N$ 

 Point C is moving to the right, so friction must act to the left

 

 

 

 

 

 

 

Pure rolling $v_{Ox}+{\omega }_{z}R=0$.  Point C is stationary.

 

 

 

 No slip occurs at the contact.  

 In this case $\left|T\right|<\mu N$ 

 We can draw the friction force in either direction at the contact (if we choose the wrong direction, our calculations will just tell us that T is negative).  It is usually convenient to choose T to act in the positive i direction, but this is not necessary.

 

 

 

 

Reverse slip: $v_{Ox}+{\omega }_{z}R<0$ Point C moves in the negative i direction over the ground

 

 

 

 

 Slip occurs at the contact,

 We have to use the friction law $T=\mu N$ 

 Point C is moving to the left, so friction must act to the right

 

 

 

 

 

Now we return to the example.

 

4.1 Draw a free body diagram showing the forces acting on the disk just after it hits the ground. 

 

We are given that $v_{x0}$ and ${\omega }_{z0}$ are both positive so we have $v_{Ox}+{\omega }_{z}R>0$.  This is forward slip, and we use the corresponding FBD.

 

 

 

 

 

 

 

4.2   Hence, find formulas for the initial acceleration $a$ and angular acceleration $\alpha$ for the disk, in terms of g , R and $\mu$.  Note that the contact point is slipping.

The equations of motion are

${{\displaystyle \sum _{i}F}}^{(i)}=M{a}_G{{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=M{r}_G\times a_G+I_{Gzz}{\alpha }_z}k$

$\begin{array}{l}-\mu Ni+(N-mg)j=m{a}_{x}i\\ -\mu NRk=I_{Gzz}{\alpha }_{z}k\end{array}$

Solving these and using $I_{Gzz}=m{R}^2/2$ :

$\begin{array}{l}N=mg\mu N=-m{a}_x\frac{1}{2}m{R}^2{\alpha }_z=-\mu NR\\ \Rightarrow a_x=-\mu g{\alpha }_z=-2\mu g/R\end{array}$

4.3   Find formulas for the velocity and angular velocity of the disk, during the period while the contact point is still slipping.

The acceleration and angular acceleration are constant, so we can use the constant acceleration formulas:

$v_x=v_0-\mu gt{\omega }_z={\omega }_0-2\mu gt/R$

 

4.4   Find a formula for the time at which the disk will reverse its direction of motion.

Velocity is reversed where v=0.   From the previous part, $v=v_0-\mu gt\Rightarrow t=v_0/\mu g$ at the reversal.

4.5   Find a formula for the time at which the disk begins to roll on the ground without slip.   Hence, show that the disk will reverse its direction only if $v_0<{\omega }_{0}R/2$

Rolling without sliding starts when $v_{xO}=-{\omega }_{z}R$.   We have that

$\begin{array}{l}{\omega }_z={\omega }_0-2\mu gt/R{v}_x=v_0-\mu gt\\ \Rightarrow v_0-\mu gt=-\left({\omega }_{0}R-2\mu gt\right)\text{when }{v}_x=-{\omega }_{z}R\\ \Rightarrow t=(v_0+{\omega }_{0}R)/3\mu g\end{array}$

The reversal will only occur if rolling without slip occurs after the reversal of velocity. This means

$(v_0+{\omega }_{0}R)/3\mu g>v_0/\mu g\Rightarrow v_0<{\omega }_{0}R/2$

 

Example 5:  The ‘Sweet Spot’ on a softball or baseball bat, or tennis or squash racket is a point that minimizes the reaction forces acting on the athlete’s hand when the ball is struck.   In fact, any rigid body has a sweet spot $–$ the magic point is called the ‘center of percussion’ of a rigid body.

For baseball and softball bats in particular, there is a standard ASTM test that can be used to measure the position of the sweet spot.   The test works like this: the bat is suspended from the knob on handle, so it swings like a pendulum.   The period of vibration $\tau$ of the swinging bat is then measured.  ASTM say that the center of percussion is then a distance

$d=\frac{{\tau }^{2}g}{4{\pi }^2}$

from the end of the handle.   Why does this work?   It seems that this test has nothing whatever to do with a ball hitting the bat!

We will solve this problem in two parts.  First, we will calculate a formula for the period of vibration in the ASTM test.   Then we will calculate the position of the center of percussion.  We will see that the ASTM test does indeed make the correct prediction.

We can calculate the period using the energy method.  The figure shows the ASTM pendulum test.  We assume that

·         The bat has a mass moment of inertia about its COM $I_{Gzz}$ 

·         The COM is a distance L from O

 

The bat pivots about O, so we can use the fixed axis rotation formula for the kinetic energy

$T=\frac{1}{2}{I}_{Ozz}{\left(\frac{d\theta }{dt}\right)}^2$

Here $I_{Ozz}=I_{Gzz}+M{L}^2$ (using the parallel axis theorem).

The potential energy is $V=-MgL\cos \theta$ .  

Energy conservation gives

$\begin{array}{l}T+V=\text{constant}\Rightarrow \frac{d}{dt}(T+V)=0\\ \Rightarrow I_{Ozz}\frac{d^2\theta }{d{t}^2}\frac{d\theta }{dt}+MgL\sin \theta \frac{d\theta }{dt}=0\end{array}$

If $\theta$ is small then $\sin \theta \approx \theta$ so the equation of motion reduces to

$\frac{I_{Ozz}}{MgL}\frac{d^2\theta }{d{t}^2}+\theta =0$

This is a standard ‘Case 1’ EOM.  The natural frequency is ${\omega }_n=\sqrt{\frac{MgL}{I_{Ozz}}}$ so the period is

 $\tau =\frac{2\pi }{{\omega }_n}=2\pi \sqrt{\frac{I_{Ozz}}{MgL}}$

Next, we find the position of the ‘sweet spot’. We can do this by calculating the reaction forces on the handle when the bat is struck, and finding the impact point that minimizes the reaction force. 

 

The figure shows an impact event.   We assume that:

• The bat rotates in the horizontal plane (so gravity acts out of the plane of the figure).  
• The bat rotates about the point O
• The ball impacts the bat a distance d from the handle.
• The ball exerts a (large) force $F_{impact}$ on the bat
• Reaction forces $R_x,R_y$ act on the handle during the impact.

 

This is a planar problem, so we can use the 2D equations of motion. The equation for translational motion gives

$(R_x-F_{impact})i+R_{y}j=M{a}_G$ 

 

For the rotational equation we can also use the short-cut for fixed axis rotation

 

$\begin{array}{l}{{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=I_{Ozz}{\alpha }_z}k\\ \Rightarrow -F_{impact}dk=I_{Ozz}{\alpha }_{z}k\Rightarrow {\alpha }_z=-\frac{F_{impact}d}{I_{Ozz}}\end{array}$

We can relate $a_G$ to ${\alpha }_z$ using the rigid body formula:

$a_G={\alpha }_{z}k\times r_G-{\omega }_z^2{r}_G=\alpha {}_{z}Li+{\omega }_z^{2}Lj$

 

We therefore see that

$\begin{array}{l}(R_x-F_{impact})i+R_{y}j=M\left(\alpha {}_{z}Li+{\omega }_z^{2}Lj\right)\\ \Rightarrow R_x=F_{impact}+M{\alpha }_{z}L\\ \Rightarrow R_x=F_{impact}\left(1-\frac{MdL}{I_{Ozz}}\right)\end{array}$

The sweet spot is at the position that makes $R_x=0$ , which shows that

$d=\frac{I_{Ozz}}{ML}$

For comparison, the ASTM formula gives

$d=\frac{{\tau }^{2}g}{4{\pi }^2}=\frac{g}{4{\pi }^2}{\left(2\pi \sqrt{\frac{I_{Ozz}}{MgL}}\right)}^2=\frac{I_{Ozz}}{ML}$

Example 6. The ‘Cubli’ is used to develop control algorithms used to stabilize aircraft and spacecraft.  It consists of a cube whose attitude can be controlled by spinning a set of reaction wheels inside the cube.

This simplified 1-D version of the device is used to test the algorithm that stands the cube up on one edge.   The goal of this problem is to do the preliminary design calculations needed to set up the system.

Idealize the rectangular frame as four rods with length L and combined mass M and the spinning wheel as a ring with radius R and mass m.  The corner at O is supported by a frictionless bearing.

Part 1: Find formulas for the mass moments of inertia of the frame and the wheel (about the center of the wheel).

 

The ring is easy $–$ we can use the formula $I_R=m{R}^2$ 

The frame is made up of four rods of mass M/4.   The moment of inertia of one rod about its center of mass is $\frac{1}{12}\frac{M}{4}{L}^2$ .  We need to shift the COM by a distance of L/2 to the center of the frame.  The total mass moment of inertia of the frame is therefore

$I_F=4\left(\frac{1}{12}\frac{M}{4}{L}^2+\frac{M}{4}{\left(\frac{L}{2}\right)}^2\right)=\frac{1}{3}M{L}^2$

 

Part 2: The frame is at rest and the wheel is spun up (clockwise) to an angular speed ${\omega }_0$ .   Find the total angular momentum of the system about the corner at O.

The formula for angular momentum is $h_O={\displaystyle \sum _{}^{}r\times m{v}_G}+{\displaystyle \sum _{}^{}I\omega }$ 

Since the frame is not moving only the second term contributes and we get $h=-m{R}^2{\omega }_{0}k$ 

 

Part 3: Thee wheel is then braked quickly, which causes the frame to rotate about the corner O at angular speed ${\omega }_f$ , while the motor driving the ring spins at (clockwise) angular speed ${\omega }_1$ (note that this is relative to the frame).  Write down the angular momentum of the system about O.

 

          Note that the frame rotates about O so the COM of the ring and frame are both in circular motion about O.   We know the speed of their COMs are therefore ${\omega }_{f}L/\sqrt{2}$ 

Use the formula again

$\begin{array}{l}h={\displaystyle \sum r_G\times m{v}_G+I_{Gzz}{\omega }_{z}k}\\ =-\left(\frac{1}{3}M{L}^2{\omega }_f+\frac{L}{\sqrt{2}}M\frac{L}{\sqrt{2}}{\omega }_f+m{R}^2({\omega }_1+{\omega }_f)+\frac{L}{\sqrt{2}}m\frac{L}{\sqrt{2}}{\omega }_f\right)k\\ =-\left(\frac{5}{6}M{L}^2+m(R^2+\frac{L^2}{2})\right){\omega }_{f}k-m{R}^2{\omega }_{1}k\end{array}$

We could also use the fixed axis rotation formula for the frame (using the mass moment of inertia about O) but this would not work for the ring, because O is not a stationary point on the ring.

Part 4: Explain why angular momentum is conserved about O during the braking.  Use momentum conservation to find an equation relating ${\omega }_f$ to $({\omega }_1-{\omega }_0)$

  The external forces acting on the frame and ring together are (1) gravity and (2) reaction forces at O.   We assume that the speed change of the rotor takes place over a very short time interval.  The force of gravity is constant and exerts a negligible impulse on the system during this time interval.  The reactions exert a finite impulse, but if we take moments about O the external angular impulse about O on the system vanishes.    This means angular momentum must be conserved. 

$\begin{array}{l}{h}_1-h_0=0\Rightarrow -\left(\frac{5}{6}M{L}^2+m(R^2+\frac{L^2}{2})\right){\omega }_{f}k-m{R}^2{\omega }_{1}k+m{R}^2{\omega }_{0}k=0\\ \Rightarrow {\omega }_f=\frac{m{R}^2({\omega }_0-{\omega }_1)}{\left(\frac{5}{6}M{L}^2+m(R^2+\frac{L^2}{2})\right)}\end{array}$

Part 5: For the special case ${\omega }_1=0$ show that the critical value of ${\omega }_0$ required to flip the frame (and ring) into the stationary vertical configuration is

${\omega }_0=\sqrt{\left(\frac{5}{6}M{L}^2+m(R^2+\frac{L^2}{2})\right)(m+M)}\frac{\sqrt{gL}\sqrt{\left(\sqrt{2}-1\right)}}{m{R}^2}$

Energy is conserved as the frame rotates up onto its edge. 

The formula for the kinetic energy of a system of rigid bodies is

$T={\displaystyle \sum _{}^{}\frac{1}{2}m{\left|v_G\right|}^2}+{\displaystyle \sum _{}^{}\frac{1}{2}\omega \cdot I_G\omega }$

For 2D problems we can replace the last term by $\frac{1}{2}{I}_{Gzz}{\omega }_z^2$ 

Assume that the frame is at rest in the upright state.  The total potential and kinetic energy in the upright state is therefore

$T_1+U_1=\frac{1}{2}m{R}^2{\omega }_1^2+(m+M)g\frac{L}{\sqrt{2}}$

In the initial state

$\begin{array}{l}{T}_0+U_0=\frac{1}{2}(m+M){\left(\frac{L}{\sqrt{2}}{\omega }_f\right)}^2+\frac{1}{2}m{R}^2{\left({\omega }_1+{\omega }_f\right)}^2+\frac{1}{2}\frac{1}{3}M{L}^2{\omega }_f^2+(m+M)g\frac{L}{2}\\ =\frac{1}{2}\left(\frac{5}{6}M{L}^2+m(R^2+\frac{L^2}{2})\right){\omega }_f^2+m{R}^2{\omega }_1{\omega }_f+\frac{1}{2}m{R}^2{\omega }_1^2+(m+M)g\frac{L}{2}\end{array}$

Energy conservation gives

$\begin{array}{l}\frac{1}{2}\left(\frac{5}{6}M{L}^2+m(R^2+\frac{L^2}{2})\right){\omega }_f^2+m{R}^2{\omega }_1{\omega }_f+\frac{1}{2}m{R}^2{\omega }_1^2+(m+M)g\frac{L}{2}=\frac{1}{2}m{R}^2{\omega }_1^2+(m+M)g\frac{L}{\sqrt{2}}\\ \Rightarrow \frac{1}{2}\left(\frac{5}{6}M{L}^2+m(R^2+\frac{L^2}{2})\right){\omega }_f^2+m{R}^2{\omega }_1{\omega }_f-(m+M)gL\left(\frac{1}{\sqrt{2}}-\frac{1}{2}\right)=0\end{array}$

For ${\omega }_1=0$ we get

${\omega }_f=\sqrt{\frac{(m+M)}{\left(\frac{5}{6}M{L}^2+m(R^2+\frac{L^2}{2})\right)}}\sqrt{gL}\sqrt{\left(\sqrt{2}-1\right)}$

From part 4 we get

${\omega }_0=\frac{\left(\frac{5}{6}M{L}^2+m(R^2+\frac{L^2}{2})\right)}{m{R}^2}{\omega }_f=\sqrt{\left(\frac{5}{6}M{L}^2+m(R^2+\frac{L^2}{2})\right)(m+M)}\frac{\sqrt{gL}\sqrt{\left(\sqrt{2}-1\right)}}{m{R}^2}$

 

 

 

 

6.8.2 Solutions to 3D problems

 

Example 1: The figure shows a wheel spinning on a frictionless axle.   The axle is supported on one side (at A) by a pivot that allows free rotation in any direction.   If the wheel were not spinning, it would simply swing about A like a pendulum.   But if the angular speed is high enough, the axle remains horizontal, and the wheel turns slowly about the vertical axis.   This behavior is called ‘precession’ and is a bit mysterious $–$ why does spin somehow hold the wheel up?   The goal of this example is to explain this, and to calculate a formula for the rotation rate of the axle.

 

We will do this by showing that steady precession satisfies all the equations of motion.  

 

1.1 Let n be a unit vector parallel to the axle.  Consider the disk at the instant when $n=i$ , and assume that

·         the disk spins at constant rate about the axle at $\nu$ radians per second, 

·         the disk rotates slowly at constant rate about k at $\Omega$ radians per second

Find the angular velocity and angular acceleration at the instant shown in the figure

 

The angular velocity is easy $–$ we just add the two vectors: $\omega =\nu n+\Omega k=\nu i+\Omega k$

 

The angular acceleration is harder.  Both $\omega$ and $\Omega$ are constant.   But this does not mean that the angular velocity vector is constant, because the axle is rotating about the k axis.  The direction of the angular velocity is changing, even though the magnitude is not.   We can calculate the rate of change of n by using the rigid body formula

$\frac{d}{dt}(r_B-r_A)=v_B-v_A=\omega \times (r_B-r_A)$

If we choose A and B to be a unit distance apart, then $(r_B-r_A)=n$ and therefore

$\frac{dn}{dt}=\omega \times n$

We can now calculate the angular acceleration

$\alpha =\frac{d\omega }{dt}=\frac{d}{dt}(\Omega k+\nu n)=\nu \frac{dn}{dt}=\nu (\Omega k+\nu n)\times n=\nu \Omega j$

 

1.2 Find a formula for the acceleration of the center of mass of the disk

 

We can use the rigid body formula

$\begin{array}{l}{a}_G-a_A=\alpha \times (r_G-r_A)+\omega \times \left[\omega \times (r_G-r_A)\right]\\ =\nu \Omega j\times (di)+\left(\nu i+\Omega k\right)\times \left[\left(\nu i+\Omega k\right)\times (di)\right]\\ =-{\Omega }^{2}di\end{array}$

A quicker way is to notice that the COM is in circular motion around A and use the circular motion formula, with the same result.

 

1.3 Draw a free body diagram showing the forces acting on the wheel

 

1.3 Write down the equations of translational and rotational motion for the disk

${{\displaystyle \sum _{i}F}}^{(i)}=M{a}_G\Rightarrow Ti+(N-Mg)k=-Md{\Omega }^{2}i$

$\begin{array}{l}{{\displaystyle \sum _i{r}_i\times F}}^{(i)}+{\displaystyle \sum _j{Q}^{(j)}=M{r}_G\times a_G+I_G\alpha }+\omega \times \left[I_G\omega \right]\\ \Rightarrow di\times (-Mgk)=M(di)\times (-{\Omega }^{2}di)+\left[\begin{array}{ccc}{I}_{Gxx}& 0& 0\\ 0& I_{Gyy}& 0\\ 0& 0& I_{Gzz}\end{array}\right]\left[\begin{array}{c}0\\ \nu \Omega \\ 0\end{array}\right]+\left[\begin{array}{c}\nu \\ 0\\ \Omega \end{array}\right]\times \left[\begin{array}{ccc}{I}_{Gxx}& 0& 0\\ 0& I_{Gyy}& 0\\ 0& 0& I_{Gzz}\end{array}\right]\left[\begin{array}{c}\nu \\ 0\\ \Omega \end{array}\right]\end{array}$

 

Working through the cross products and the matrix-vector products we get

$Mgdj=I_{Gyy}\nu \Omega j-(I_{Gzz}-I_{Gxx})\nu \Omega j$

We see that steady precession can indeed satisfy all the equations of motion.  Moreover, for a disk (or any solid of revolution) $I_{Gzz}=I_{Gxx}$ , so we can calculate the precession rate

$\begin{array}{l}Mgdj=I_{Gxx}\nu \Omega j\\ \Rightarrow \Omega =\frac{Mgd}{I_{Gxx}\nu }\end{array}$

 

 

Example 2:  The prism shown in the figure floats in space (no gravity).  At time t=0 its faces are perpendicular to the i,j,k axes as shown.   It is then given an initial angular velocity $\omega ={\omega }_{z0}k+{\omega }_{y0}j$ with ${\omega }_{y0}<<{\omega }_{z0}$  (i.e. we set the body spinning about the k axis, but give it a very small disturbance) .   Investigate the nature of the subsequent motion, with both hand calculations and by writing a MATLAB script that will animate the motion of the prism.

 

No forces or moments act on the prism.  We can use the equations of motion

$0=M{a}_{G}0=M{r}_G\times a_G+I_G\alpha +\omega \times \left[I_G\omega \right]$

The angular momentum equation can be written out explicitly

$\left[\begin{array}{ccc}{I}_{Gxx}& 0& 0\\ 0& I_{Gyy}& 0\\ 0& 0& I_{Gzz}\end{array}\right]\left[\begin{array}{c}d{\omega }_x/dt\\ d{\omega }_y/dt\\ d{\omega }_z/dt\end{array}\right]+\left[{\omega }_x,{\omega }_y,{\omega }_z\right]\times \left(\left[\begin{array}{ccc}{I}_{Gxx}& 0& 0\\ 0& I_{Gyy}& 0\\ 0& 0& I_{Gzz}\end{array}\right]\left[\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right]\right)$

(we could substitute values for $I_{Gxx},I_{Gyy},I_{Gzz}$ in terms of a,b,c and M but it is clearer to leave them) Expanding out the matrix products and cross product gives

$\begin{array}{l}{I}_{Gxx}\frac{d{\omega }_x}{dt}+\left(I_{Gzz}-I_{Gyy}\right){\omega }_y{\omega }_z=0\\ I_{Gyy}\frac{d{\omega }_y}{dt}-\left(I_{Gzz}-I_{Gxx}\right){\omega }_x{\omega }_z=0\\ I_{Gzz}\frac{d{\omega }_z}{dt}+\left(I_{Gyy}-I_{Gxx}\right){\omega }_x{\omega }_y=0\end{array}$

At time t=0 ${\omega }_x$ is zero and ${\omega }_y$ is small.  They might increase, but we will only consider behavior while they remain small.  In this case ${\omega }_x{\omega }_y$ is extremely small so we can assume $d{\omega }_z/dt\approx 0$ .   We can then decouple the first two equations like this:

1. Differentiate the second equation with respect to time  $I_{Gyy}\frac{d^2{\omega }_y}{d{t}^2}-\left(I_{Gzz}-I_{Gxx}\right)\frac{d{\omega }_x}{dt}{\omega }_{z0}=0$
2. Now we can substitute for $d{\omega }_x/dt$ using the first equation, and divide by $I_{Gyy}$ 

$\frac{d^2{\omega }_y}{d{t}^2}+\frac{\left(I_{Gzz}-I_{Gxx}\right)\left(I_{Gzz}-I_{Gyy}\right)}{I_{Gyy}{I}_{Gxx}}{\omega }_{z0}^2{\omega }_y=0$

This is an equation of the form

$\frac{d^2{\omega }_y}{d{t}^2}+\lambda {\omega }_y=0$

We recognize this as an undamped vibration equation (case I or case II from our table of solutions).  Its solution depends on the sign of $\lambda$ :

1. For $\lambda >0$ the solution is ${\omega }_y=A\sin \sqrt{\lambda }t+B\cos \sqrt{\lambda }t$ where A, B are constants.  This is stable motion - ${\omega }_y$ remains small.
2. For $\lambda <0$ the solution is ${\omega }_y=A\exp \sqrt{\lambda }t+B\exp (-\sqrt{\lambda }t)$.   This is unstable motion - ${\omega }_y$ will become very large.

 

 

The sign of $\lambda$ is determined by the product $\left(I_{Gzz}-I_{Gxx}\right)\left(I_{Gzz}-I_{Gyy}\right)$.  There are three possible cases:

1. $I_{Gzz}$ is greater than $I_{Gxx},I_{Gyy}$ (the k axis has the maximum inertia).  Motion is stable
2. $I_{Gzz}$ is less than $I_{Gxx},I_{Gyy}$ (the k axis has the minimum inertia).  Motion is stable
3. $I_{Gzz}$ is between $I_{Gxx},I_{Gyy}$.  Motion is unstable.

 

We can learn more about the motion by using MATLAB to solve the equations of motion for us.   Since there is no motion of the center of mass, we only need to consider rotational motion.  We know that we can describe the orientation of the prism by the rotation tensor R and its rate of change of orientation by the angular velocity $\omega$ .   The orientation and angular velocity are governed by the differential equations

$\begin{array}{l}\frac{dR}{dt}=WR\\ I_G\frac{d\omega }{dt}+\omega \times \left[I_G\omega \right]=0\end{array}$

where $I_G=R{I}_G^0{R}^T$ is the rotated inertia tensor for the block, and W is the spin tensor

$W=\left[\begin{array}{ccc}0& -{\omega }_z& {\omega }_y\\ {\omega }_z& 0& -{\omega }_x\\ -{\omega }_y& {\omega }_x& 0\end{array}\right]$

 

We need to set up the MATLAB ‘ode’ solver to calculate R and $\omega$ as functions of time by integrating these equations.

 

We can store the unknown rotation matrix and the angular velocity vector in a MATLAB vector:

$w=\left[R_{xx},R_{xy},R_{xz},R_{yx},R_{yy},R_{yz},R_{zx},R_{zy},R_{zz},{\omega }_x,{\omega }_y,{\omega }_z\right]$

We need to write a MATLAB function that will calculate the time derivatives of this vector, given its current value.  The calculation involves the following steps:

(1)    Assemble the vectors $\omega$ and the rotation tensor R from the Matlab solution vector w.  Matlab has a useful function that will automatically convert a matrix to a vector, and vice-versa.  For example, $R$ (a 3x3 matrix) can be converted to w (a 1x9 column vector) using

w = reshape(transpose(R),[9,1]))

To transform w (as a column vector) back to R, you can use

R = transpose(reshape(w,[3,3]))

(2)    Calculate the spin tensor W

(3)    Calculate the rotated inertia tensor  $I_G=R{I}_G^0{R}^T$ (Matlab will multiply the matrices for us)

(4)    Solve the equations for the angular acceleration $\alpha$

(5)    Calculate $dR/dt=WR$ 

(6)    Assemble the matlab vector $\frac{dw}{dt}=\left[{\dot{R}}_{xx},{\dot{R}}_{xy},{\dot{R}}_{xz},{\dot{R}}_{yx},{\dot{R}}_{yy},{\dot{R}}_{yz},{\dot{R}}_{zx},{\dot{R}}_{zy},{\dot{R}}_{zz},{\alpha }_x,{\alpha }_y,{\alpha }_z\right]$ 

 

This sounds complicated but actually MATLAB is great at doing this sort of calculation efficiently.  Here’s a function:

 

function  dwdt = rigid_body_eom(t,w)

   Rvec = w(1:9); % Rotation matrix, stored as a vector

   omega = w(10:12);  % Angular velocity

   R = transpose(reshape(Rvec,[3,3]));

   II = R*I0*transpose(R); %Current inertia tensor, in fixed coord system

   W = [0,-omega(3),omega(2);omega(3),0,-omega(1);-omega(2),omega(1),0];

   alpha = -II\(cross(omega,II*omega)); % Angular accel

   Rdot = W*R; % Rate of change of rotation matrix

   Rdotvec = reshape(transpose(Rdot),[9,1]);

   dwdt = [Rdotvec;alpha];

end

 

We just need to set up ode45 to integrate (numerically) the differential equation:

 

omega0 = [0.,0.01,1];   % Initial angular velocity

a = [4,1,2];            % Dimensions (a,b,c) of the prism

time = 20;

initial_w = [1;0;0;0;1;0;0;0;1;transpose(omega0(1:3))];

I0 = [a(2)^2+a(3)^2,0,0;0,a(1)^2+a(3)^2,0;0,0,a(1)^2+a(2)^2];

options = odeset('RelTol',0.00000001);

sol = ode45(@(t,w) rigid_body_eom(t,w,I0),[0,time],initial_w,options);

animate_rigid_body(sol,a,[0,time])

 

You can download the full script here.

 

The figures below show animations of the predicted behavior for the three possible types of behavior