Wheel on the plane

Construction of the Euler equations

(* moments of inertia of the wheel with respect to contact point *)
(Subscript[J, 1])=(J/2)+M (R^2);  (Subscript[J, 2])=(J/2);   (Subscript[J, 3])=J+M (R^2); 

(* torque components with respect to contact point *)
(Subscript[K, 1])=-R Fg Cos[θ]+R Fx Sin[θ]Sin[φ];
(Subscript[K, 2])=0;
(Subscript[K, 3])=-R Fx Cos[φ];

(* ω in the body frame *)
(Subscript[ω, 1])=(OverDot[θ]) ;
(Subscript[ω, 2])=(OverDot[φ]) Sin[θ];
(Subscript[ω, 3])=(OverDot[φ]) Cos[θ]+(OverDot[ψ]);

(* Rules for substitution of variables *)
anglest= {θ->θt[t],(OverDot[θ])->(θt)'[t],(Overscript[θ, ¨])->(θt)''[t],φ->φt[t],(OverDot[φ])->(φt)'[t],(Overscript[φ, ¨])->(φt)''[t],ψ->ψt[t],(OverDot[ψ])->(ψt)'[t],(Overscript[ψ, ¨])->(ψt)''[t]};
anglestinv= {θt[t]->θ,(θt)'[t]->(OverDot[θ]),(θt)''[t]->(Overscript[θ, ¨]),φt[t]->φ,(φt)'[t]->(OverDot[φ]),(φt)''[t]->(Overscript[φ, ¨]),ψt[t]->ψ,(ψt)'[t]->(OverDot[ψ]),(ψt)''[t]->(Overscript[ψ, ¨])};

(Subscript[EuEq, 1])=(Subscript[J, 1])(∂_t)((Subscript[ω, 1])/.anglest)==((Subscript[J, 2])-(Subscript[J, 3]))(Subscript[ω, 2])(Subscript[ω, 3])-(Subscript[J, 2])(Subscript[ω, 2])(OverDot[ψ])+(Subscript[K, 1]) /.anglest ;
(Subscript[EuEq, 2])=(Subscript[J, 2])(∂_t)((Subscript[ω, 2])/.anglest)==((Subscript[J, 3])-(Subscript[J, 1]))(Subscript[ω, 3])(Subscript[ω, 1])+(Subscript[J, 1])(Subscript[ω, 1])(OverDot[ψ])+(Subscript[K, 2])/.anglest ;
(Subscript[EuEq, 3])=(Subscript[J, 3])(∂_t)((Subscript[ω, 3])/.anglest)==((Subscript[J, 1])-(Subscript[J, 2]))(Subscript[ω, 1])(Subscript[ω, 2])+(Subscript[K, 3])/.anglest ;

(Subscript[EuEq, 1]) /.anglestinv 
(Subscript[EuEq, 2]) /.anglestinv 
(Subscript[EuEq, 3]) /.anglestinv

Parameters, initial conditions, and solution of the Euler equations

(Subscript[θ, 0])=π/2;   (Subscript[φ, 0])=π/2;   (Subscript[ψ, 0])=0π/4;
(OverDot[(Subscript[θ, 0])])=0.01;   (Subscript[(OverDot[φ]), 0])=1;    (Subscript[(OverDot[ψ]), 0])=-0.01;
tMax=20;
MM=1;  RR=1; II=MM (RR^2);  FFg=10;  FFx=-0.03;

Eqs={(Subscript[EuEq, 1]),(Subscript[EuEq, 2]),(Subscript[EuEq, 3])} /. {M->MM,R->RR,J->II,Fg->FFg,Fx->FFx};
IniConds={θt[0]==(Subscript[θ, 0]),(θt)'[0]==(OverDot[(Subscript[θ, 0])]),φt[0]==(Subscript[φ, 0]),(φt)'[0]==(Subscript[(OverDot[φ]), 0]),ψt[0]==(Subscript[ψ, 0]),(ψt)'[0]==(Subscript[(OverDot[ψ]), 0])};

Solution=NDSolve[Join[Eqs,IniConds],{θt,φt,ψt},{t,0,tMax},MaxSteps->3000000]
θtt[t_]:=θt[t] /. Solution[[1]]
φtt[t_]:=φt[t] /. Solution[[1]]
ψtt[t_]:=ψt[t] /. Solution[[1]]
Plot[θtt[t]/π,{t,0,tMax},PlotRange->{0,1},PlotPoints->200]
Plot[φtt[t]/(2π),{t,0,tMax},PlotRange->{-1,3},PlotPoints->200]
Plot[ψtt[t],{t,0,tMax},PlotRange->{-10,10},PlotPoints->200]

Linear motion of the wheel

Vxt[t_]:=R (Subscript[ω, 1])Sin[θ]Sin[φ]-R (Subscript[ω, 3])Cos[φ] /. anglest /. R->RR /. Solution[[1]]
Vyt[t_]:=-R (Subscript[ω, 1])Sin[θ]Cos[φ]-R (Subscript[ω, 3])Sin[φ] /. anglest /. R->RR  /. Solution[[1]]
Vzt[t_]:=R (Subscript[ω, 1])Cos[θ] /. anglest  /. R->RR /. Solution[[1]]

(*Plot[Vxt[t],{t,0,tMax},PlotPoints->200]
Plot[Vyt[t],{t,0,tMax},PlotPoints->200]*)
(*Plot[Vzt[t],{t,0,tMax},PlotPoints->200]*)
(* Linear displacements of the geometrical center *)
Solx=NDSolve[{Xt'[t]==Vxt[t],Xt[0]==0},{Xt},{t,0,tMax},MaxSteps->1000000]
Soly=NDSolve[{Yt'[t]==Vyt[t],Yt[0]==0},{Yt},{t,0,tMax},MaxSteps->1000000]

Xtt[t_]:=Xt[t] /. Solx[[1]];
Ytt[t_]:=Yt[t] /. Soly[[1]];
Ztt[t_]:=R Sin[θ]  /. anglest /. R->RR /. Solution[[1]]

OffSet=2.1;
XMax=Xtt[tMax];
XIni=If[XMax>0,-OffSet,OffSet];
XFin=If[XMax>0,XMax+OffSet,XMax-OffSet];

YMax=Ytt[tMax];
YIni=If[YMax>0,-OffSet,OffSet];
YFin=If[YMax>0,YMax+OffSet,YMax-OffSet];

(* coordinates of the contact point *)
ξx[t_]:=Xtt[t]+R Cos[θ]Sin[φ]  /. anglest /. R->RR /. Solution[[1]]
ξy[t_]:=Ytt[t]-R Cos[θ]Cos[φ]  /. anglest /. R->RR /. Solution[[1]]

Plot[Ztt[t],{t,0,tMax},PlotPoints->300,PlotRange->{0,1},PlotStyle->{Thick,Black}]

Plot[{Xtt[t],ξx[t]},{t,0,tMax},PlotPoints->300,PlotRange->All,PlotStyle->{{Thick,Black},{Thick,Red}}]

Plot[{Ytt[t],ξy[t]},{t,0,tMax},PlotPoints->300,PlotRange->All,PlotStyle->{{Thick,Black},{Thick,Red}}]

(* X-Y Trajectory of the geometrical center *)
TrajCenter:=ParametricPlot[{Xtt[t],Ytt[t]},{t,0,tMax},PlotRange->All,PlotPoints->1000,PlotStyle->{Thick,Black}]

(* X-Y Trajectory of the contact point *)
TrajContact:=ParametricPlot[{ξx[t],ξy[t]},{t,0,tMax},PlotRange->All,PlotPoints->1000,PlotStyle->{Thick,Red}]

Show[TrajCenter,TrajContact]

3d plotting

TrajectoryCM:=ParametricPlot3D[{Xtt[t],Ytt[t],Ztt[t]},{t,0,tMax},PlotRange->{{XIni,XFin},{YIni,YFin},{0,2}},PlotStyle->{Thick},ImageSize->{700,700}] (* Trajectory of the center of mass *)
 
(* ImageSize set here and used everyehere *)
TrajectoryCP:=ParametricPlot3D[{ξx[t],ξy[t],0},{t,0,tMax},PlotRange->{{XIni,XFin},{YIni,YFin},{0,2}},PlotStyle->{Thick,Red},PlotPoints->200] (* Trajectory of the contact point *)

rx[θ_,ϕ_,ψ_]:=Cos[ϕ]Cos[ψ]-Cos[θ]Sin[ϕ]Sin[ψ]
ry[θ_,ϕ_,ψ_]:=Sin[ϕ]Cos[ψ]+Cos[θ]Cos[ϕ]Sin[ψ]
rz[θ_,ϕ_,ψ_]:=Sin[θ](1+Sin[ψ])  (* Any point on the perifery of the wheel *)

anglesTime= {θ->θtt[t],ϕ->φtt[t]};
rxt[t_,ψ_]:=Xtt[t]+rx[θ,ϕ,ψ]/.anglesTime 
ryt[t_,ψ_]:=Ytt[t]+ry[θ,ϕ,ψ]/.anglesTime 
rzt[t_,ψ_]:=rz[θ,ϕ,ψ]/.anglesTime   (* The same with previously calculated time dependence of θ and ϕ *)

WheelCP[tt_]:=ListPointPlot3D[{{rxt[t,-π/2],ryt[t,-π/2],rzt[t,-π/2]}}/.t->tt,PlotRange->{{XIni,XFin},{YIni,YFin},{0,2}},PlotStyle->{{PointSize[0.015],Red}}] (* Position of the contact point *)

WheelCM[tt_]:=ListPointPlot3D[{{Xtt[t],Ytt[t],Ztt[t]}}/.t->tt,PlotRange->{{XIni,XFin},{YIni,YFin},{0,2}},PlotStyle->{{PointSize[0.015],Black}}] (* Position of the center of mass *)

WheelPerifery[tt_]:=ParametricPlot3D[{rxt[t,ψ],ryt[t,ψ],rzt[t,ψ]}/.t->tt,{ψ,0,2π},PlotRange->{{XIni,XFin},{YIni,YFin},{0,2}},PlotStyle->{Thick,Blue}] (* The wheel itself *)

Wheel[tt_]:=Show[WheelPerifery[tt],WheelCP[tt],WheelCM[tt]]  (* Wheel with CP and CM *)

WheelTrajectory[t_]:=Show[TrajectoryCM,TrajectoryCP,WheelPerifery[t],WheelCP[t],WheelCM[t]] (* Full trajectories of CP and CM and the wheel at time t *)

WheelTrajectory[tMax/4]

Manipulate[WheelTrajectory[t],{t,0,tMax}]

Export["Wheel.avi",Animate[WheelTrajectory[tt],{tt,0,tMax}]] (* Export of animated objects does not work *)

Show[TrajectoryCM,TrajectoryCP,Wheel[0],Wheel[2.25],Wheel[2.59],Wheel[3],Wheel[4],Wheel[7.4],Wheel[7.75],Wheel[8.5],Wheel[12.5],Wheel[12.9]]