Free Precession of a Rotating Rigid Body

PlotPrecession[h_,R_,u_,tx_,frame_,cones_,axislabel_,veclabel_]:=Module[{tspace,tbody,ωx,Ωb,ω0,ω3,r,x,y,z,e1,e2,e3,ωbody,ωspace,ωf,L,I1,I2,I3,α,θ,Ω0,x0,y0,z0},
I1 = 1+(1/4)(((h/R))^2); (* These are the principal moments of our boday, assuming M R^2/5 = 1 *)
I2 = 1+(1/4)(((h/R))^2); (* The body has radius R and height h *)
I3 = 2;
ω0 = Sin[u];
ω3 = Cos[u];
α = u; (* Angle between w and e3 *)
θ = ArcTan[I1 ω0/(I3 ω3)]; (* Angle between L = I.ω and e3 *)
Ω0 = (I1-I3) ω3/I1; (* precession rate *)
ωx := {-ω0,0,-ω3+Ω0}; (* relative angular velocity *)
Ωb := {0,0,-Ω0}; (* precession vector in body frame *)
ωbody[t_]:={ω0 Cos[Ω0 t],-ω0 Sin[Ω0 t],ω3}; (* angular velocity in body frame*)
L[t_]:={I1 ω0 Cos[Ω0 t],-I2 ω0 Sin[Ω0 t],I3 ω3}; (* angular momentum in body frame *)

r[t_]:=RotationMatrix[Norm[Ωb]t,If[Abs[Ω0]<1.0 (10^(-10)),{1,0,0},Ωb]].RotationMatrix[Norm[ωx]t,ωx].{x0,y0,z0}; (* This is a tricky command! *)

x[t_]:=r[t]/.{x0->1,y0->0,z0->0}; (* x-axis in body frame *)
y[t_]:=r[t]/.{x0->0,y0->1,z0->0}; (* y-axis in body frame *)
z[t_]:=r[t]/.{x0->0,y0->0,z0->1}; (* z-axis in body frame *)

e1[t_]:={1,0,0}.x[t]{1,0,0} + {1,0,0}.y[t]{0,1,0}+{1,0,0}.z[t]{0,0,1}; (* e1 axis in space frame *)
e2[t_]:={0,1,0}.x[t]{1,0,0} + {0,1,0}.y[t]{0,1,0}+{0,1,0}.z[t]{0,0,1}; (* e2 axis in space frame *)
e3[t_]:={0,0,1}.x[t]{1,0,0} + {0,0,1}.y[t]{0,1,0}+{0,0,1}.z[t]{0,0,1}; (* e3 axis in space frame *)

ωspace[t_]:=ωbody[t].x[t]{1,0,0}+ωbody[t].y[t]{0,1,0}+ωbody[t].z[t]{0,0,1}; 
(* angular velocity in space frame *)

(* In space frame, L is fixed. In body frame, e1,e2, and e3 are fixed. The If loop selects the moving vectors, and the form of the angular velocity vector. *)

If[frame=="body",

tspace=tx;
tbody = 0;
ωf[t_]:=ωbody[t],

tspace=0;
tbody=tx;
ωf[t_]:=ωspace[t]];

Show[
Graphics3D[{Red,Tube[{{0,0,0},e1[tbody]},0.015]},
PlotRange->{{-1.5,1.5},{-1.5,1.5},{-1.5,1.5}},ViewPoint->{1,1,2},ImageSize->{460,480}],
Graphics3D[{Red,Sphere[e1[tbody],0.1]}], (* e1 axis has a red sphere *)
Graphics3D[{Green,Tube[{{0,0,0},e2[tbody]},0.015]}],
Graphics3D[{Green,Sphere[e2[tbody],0.1]}], (* e2 axis has a green sphere *)
Graphics3D[{Blue,Tube[{{0,0,0},e3[tbody]},0.015]}],
Graphics3D[{Blue,Sphere[e3[tbody],0.1]}], (* e3 axis has a blue sphere *)
Graphics3D[{Yellow,Tube[Table[ e3[tbody]h zz ,{zz,-0.5,0.5,0.1}],Table[R Sqrt[1-4(zz^2)],{zz,-0.5,0.5,0.1}]]}],
(* Our body is an ellipsoid formed by a set of tubes of different radii, in yellow *)
Graphics3D[{Black,Arrow[Tube[{{0,0,0},1.5ωf[tx]}]]}], (* angular velocity in black *)
Graphics3D[{Gray,Arrow[Tube[{{0,0,0},1.5L[tspace]/Norm[L[tspace]]}]]}], (* angular momentum in gray *)

If[cones, Graphics3D[{Blue,Cone[{e3[tbody]*0.7,{0,0,0}},0.7Tan[α]],Gray,Cone[{L[tspace]/Norm[L[tspace]],{0,0,0}},Tan[Abs[θ-α]]]}],{}], 
(* body cone in blue, centered about e3, space cone is gray, centered about L *)

Graphics3D[{Red,Arrow[Tube[{{0,0,0},x[tspace]}]]}], (* x-axis is red *)
Graphics3D[{Green,Arrow[Tube[{{0,0,0},y[tspace]}]]}], (* y-axis is green *)
Graphics3D[{Blue,Arrow[Tube[{{0,0,0},z[tspace]}]]}], (* z-axis is blue *)
If[axislabel,
Graphics3D[{
Text[Style["x",Italic,FontSize->30],1.2x[tspace]],
Text[Style["y",Italic,FontSize->30],1.2 y[tspace]],
Text[Style["z",Italic,FontSize->30],1.2 z[tspace]],
Text[Style["(Subscript[e, 1])",FontSize->30],1.25e1[tbody]],
Text[Style["(Subscript[e, 2])",FontSize->30],1.25 e2[tbody]],
Text[Style["(Subscript[e, 3])",FontSize->30],1.25e3[tbody]]}],{}],
If[veclabel,Graphics3D[{Gray,
Text[Style["L",Italic,FontSize->30],1.7 L[tspace]/Norm[L[tspace]]],Black,Text[Style["ω",FontSize->30],1.7ωf[tx]]}],{}],
SphericalRegion-> True
]];

Manipulate[
PlotPrecession[h,R,u,tx,frame,cones,axislabel,veclabel],
"height",
{{h,0.1,""},0.02,2,.01, Appearance-> "Labeled", ImageSize-> Tiny},
"radius",
{{R,0.3,""},0.02,1,.01, Appearance-> "Labeled", ImageSize-> Tiny},
"time",
{{tx,0,""},0,1000,.01,Appearance-> "Labeled", ImageSize-> Tiny},
{{u, .1, Style["u", Italic]},0,π,.01, Appearance-> "Labeled", ImageSize-> Tiny},
Delimiter,
"frame",
Control[{{frame, "body", ""},{"body","space"}}],
Delimiter,
"show cones",
Control[{{cones,True, ""},{True,False}}],
"axis labels",
Control[{{axislabel,True,""},{True,False}}],
"vector labels",
Control[{{veclabel,True,""},{True,False}}],
SaveDefinitions->True, ControlPlacement-> Left]

The free precession of a rotating rigid body is a classic application of Euler's equations. Since the precession is free, there is no torque, and the angular momentum is constant in the space frame. In the body frame, however, the body axes (Subscript[e, 1]), (Subscript[e, 2]), (Subscript[e, 3]) (indicated by the red, green, and blue spheres) are fixed, and the angular momentum L and space axes (x,y,z) (indicated by the red, green, and blue arrows) will change with time. For an axisymmetric ellipsoid, the angular velocity ω and the angular momentum L precess about the symmetry axis (Subscript[e, 3]). This precession can be visualized by the rotation of a space cone (centered about L) around the body cone (centered about (Subscript[e, 3])), with the angular velocity ω at their point of contact. This Demonstration shows the free precession as a function of time in both the body and space frames. The body's height and radius can be varied, as can the initial angular velocity ω={sin u, 0, cos u}.

Euler's equations for a rotating rigid body can be exactly solved for the free motion (with zero torque) of an axisymmetric body, with principal moments of inertia (Subscript[I, 1])=(Subscript[I, 2])≠ (Subscript[I, 3]). These give the angular velocity ω and angular momentum L in the body frame, and will exhibit free precession if the angular velocity is not initially along one of the symmetry axes. The body frame and space frame can be related by a product of two rotation matrices. This is used to transform ω, L, and the space and body axes (x,y,z) and ((Subscript[e, 1]), (Subscript[e, 2]), (Subscript[e, 3])), and to visualize their dynamics as a function of time. In the body frame, both ω and L precess about (Subscript[e, 3]), while in the space frame, ω and (Subscript[e, 3]) precess about L.

The solution to Euler's equations can be found in the following references. For more information, see the Wikipedia entries for "Euler's equations (rigid body dynamics)" and "Precession".

References

[1] J. R. Taylor, "Rotational Motion of Rigid Bodies," Classical Mechanics, Herndon, VA: University Science Books, 2005 pp. 367–416.

[2] S. T. Thornton and J. B. Marion, "Dynamics of Rigid Bodies," Classical Dynamics of Particles and Systems, Pacific Grove, CA: Brooks/Cole, 2004 pp. 411–467.

Contributed by: Frederick W. Strauch