Эффект Джанибекова

Численное интегрирование уравнений пространственного движения твёрдого тела вокруг неподвижной точки и его визуализация

Самарский университет
Кафедра теоретической механики.
Юдинцев В. В.

Постановка задачи

Рассмотрим вращение твёрдого тела (гайки-барашка) вокруг центра масс в условиях невесомости. Это движение эквивалентно сферическому движению твердого тела вокруг неподвижной точки, совпадающей с центром масс.

С твёрдым телом связана система координат Cxyz, начало которой совпадает с центром масс тела, а оси направлены так, как показано на рисунке.

Моменты инерции вокруг осей связанной системы координат удовлетворяют следующему условию:

(Subscript[J, x])<(Subscript[J, y])<(Subscript[J, z])

Кинематические соотношения

Ориентацию связанной системы координат Cxyz относительно неподвижной системы координат (Subscript[Cx, 0])(Subscript[y, 0])(Subscript[z, 0]) будет определяться тремя углами Эйлера последовательности X-Z’-X’’:

Первый поворот совершается вокруг оси (Subscript[x, 0]) на угол ψ

Второй поворот совершается вокруг новой оси z=z' на угол ϑ

Третий поворот  совершается вокруг новой оси x на угол φ

Эти три угла позволяют задать угловое положение гайки-барашка в пространстве относительно неподвижной системы координат.

Матрица преобразования координат

Каждый элементарный поворот можно описать матрицей преобразования координат из новой системы координат (повернутой) в исходную. Например для поворота вокруг оси x на угол ψ матрица преобразования координат (или матрица поворота) имеет вид:

(Subscript[A, x])=(({{1, 0, 0}, {0, cosψ, −sinψ}, {0, sinψ, cosψ}}))

Предположим, что система координат Cxyz развёрнута относительно неподвижной системы координат (Subscript[Cx, 0])(Subscript[y, 0])(Subscript[z, 0]) вокруг оси (Subscript[Cx, 0]) на угол ψ.

Матрица преобразования координат

Если известны координаты точки (вектора) в повернутой системе, то координаты этой же точки в исходной неподвижной системе координат  будут определяться выражением:

(({{x}, {y}, {z}}))=(({{1, 0, 0}, {0, cosψ, −sinψ}, {0, sinψ, cosψ}}))(({{x'}, {y'}, {z'}}))

Обратное преобразование:

(({{x'}, {y'}, {z'}}))=(((({{1, 0, 0}, {0, cosψ, −sinψ}, {0, sinψ, cosψ}})))^T)(({{x}, {y}, {z}}))=(({{1, 0, 0}, {0, cosψ, sinψ}, {0, -sinψ, cosψ}}))(({{x}, {y}, {z}}))

Матрицы элементарных поворотов

Поворот вокруг оси x:

RotationMatrix[ψ,{1,0,0}]
%//MatrixForm

Поворот вокруг оси y:

RotationMatrix[ϑ,{0,1,0}];
%//MatrixForm

Поворот вокруг оси z:

RotationMatrix[φ,{0,0,1}]
%//MatrixForm

Сложный поворот (последовательность X - Z’ - X”)

Матрицы элементарных поворотов

A1=RotationMatrix[ψ,{1,0,0}]

A2=RotationMatrix[ϑ,{0,0,1}]

A3=RotationMatrix[φ,{1,0,0}]

Сложный поворот

A1=RotationMatrix[ψ,{1,0,0}]
A2=RotationMatrix[ϑ,{0,0,1}]
A3=RotationMatrix[φ,{1,0,0}]

Матрица сложного поворота определяется произведением  матриц элементарных поворотов в прямом порядке.

Матрица преобразования координат из повернутой системы (после трех поворотов) в исходную:

A123=A1.A2.A3;
A123//MatrixForm

Динамические уравнения Эйлера

Уравнения свободного вращения твердого тела вокруг центра масс в отсутствии внешних моментов (движение “по инерции”):

eq={
Jx wx'[t]==(Jy-Jz) wy[t] wz[t],
Jy wy'[t]==(Jz-Jx) wx[t] wz[t],
Jz wz'[t]==(Jx-Jy) wy[t] wx[t]
};

Динамические уравнения Эйлера

Уравнения свободного вращения твердого тела вокруг центра масс в отсутствии внешних моментов (движение “по инерции”):

eq={
Jx wx'[t]==(Jy-Jz) wy[t] wz[t],
Jy wy'[t]==(Jz-Jx) wx[t] wz[t],
Jz wz'[t]==(Jx-Jy) wy[t] wx[t]
};

wx, wy, wz - это проекции угловой скорости тела на связанные оси.

Параметры системы и начальные угловые скорости:

params={Jx->1,Jy->2,Jz->3,wx0->0.001,wy0->1,wz0->0.001};

Вращение начинается вокруг оси Cy c угловой скоростью 1 рад/с. 
Также заданы небольшие возмущения угловой скорости вокруг других осей (0,001 рад/с)

Численное интегрирование

Подставляем в уравнения движения параметры:

params={Jx->1,Jy->2,Jz->3,wx0->0.001,wy0->1,wz0->0.001};
neq={eq,wx[0]==wx0,wy[0]==wy0,wz[0]==wz0}//.params//Flatten

tk=60;
sol=NDSolve[neq,{wx[t],wy[t],wz[t]},{t,0,tk}]//Flatten;

Plot[{wx[t],wy[t],wz[t]}//.sol//Evaluate,{t,0,tk},ImageSize->500,FrameLabel->{"t, c","ω, рад/с"},PlotLegends->Placed[{"(Subscript[ω, x])","(Subscript[ω, y])","(Subscript[ω, z])"},{Right,Bottom}],PlotTheme->{"Web","FrameGrid"}]

В результате интегрирования динамических уравнений Эйлера получаются проекции угловые скорости.

Как определить углы ориентации связанной системы координат относительно неподвижной?

Кинематические уравнения

Кинематические уравнения связывают проекции угловой скорости с производными углов.

Выразим проекции угловой скорости на связанные оси тела через производные углов

(Subscript[ω, x])=(OverDot[φ])+(OverDot[ψ])·cos ϑ 
(Subscript[ω, y])=(OverDot[ϑ])sin φ -(OverDot[ψ])sin ϑ cos φ  
(Subscript[ω, z])=(OverDot[ϑ])cos φ +(OverDot[ψ]) sin ϑ sin φ

Выражение производных углов через угловые скорости (кинематические уравнения):

({{(OverDot[ψ]), =((Subscript[ω, z])sinφ−(Subscript[ω, y])cosφ)/sinϑ}}) 
(OverDot[ϑ]) = (Subscript[ω, z])cos φ+(Subscript[ω, y])sinφ;    
(OverDot[φ])=(Subscript[ω, x])+((Subscript[ω, y])cosφ−(Subscript[ω, z])sinφ)/tanϑ

Кинематические уравнения

kinematicEq={
ψ'[t]==(wz[t] Sin[φ[t]]-wy[t] Cos[φ[t]])/Sin[ϑ[t]],
ϑ'[t]==wz[t] Cos[φ[t]]+wy[t] Sin[φ[t]],
φ'[t]==wx[t]+(wy[t] Cos[φ[t]]-wz[t] Sin[φ[t]])/Tan[ϑ[t]]
};

params={Jx->1,Jy->2,Jz->3,wx0->0.001,wy0->1,wz0->0.001,ψ0->0.0,ϑ0->π/2,φ0->0.0};

Уравнения движения: динамические + кинематические

neq={
eq,kinematicEq,
wx[0]==wx0,wy[0]==wy0,wz[0]==wz0,
ψ[0]==ψ0,ϑ[0]==ϑ0,φ[0]==φ0
}//.params//Flatten

Результат

neq={{eq,kinematicEq},wx[0]==wx0,wy[0]==wy0,wz[0]==wz0,ψ[0]==ψ0,ϑ[0]==ϑ0,φ[0]==φ0}//.params//Flatten

tk=60;
sol=NDSolve[neq,{wx[t],wy[t],wz[t],ψ[t],ϑ[t],φ[t]},{t,0,tk}]//Flatten;

Plot[{wx[t],wy[t],wz[t]}//.sol//Evaluate,{t,0,tk},ImageSize->500,FrameLabel->{"t, c","ω, рад/с"},PlotLegends->{"(Subscript[ω, x])","(Subscript[ω, y])","(Subscript[ω, z])"},PlotTheme->{"Web","FrameGrid"}]

Углы

Plot[{ϑ[t],φ[t]}*(180/π)//.sol//Evaluate,{t,0,tk},ImageSize->400,FrameLabel->{"t, c","Угол, градус"},PlotLegends->{"ϑ","φ"},PlotTheme->{"Web","FrameGrid"}]

Plot[ψ[t]*(180/π)//.sol//Evaluate,{t,0,tk},ImageSize->400,FrameLabel->{"t, c","Угол, градус"},PlotLegends->{"ψ"},PlotTheme->{"Web","FrameGrid"}]

Импорт 3D модели объекта

object=Import[NotebookDirectory[]<>"screw.stl",{"GraphicsComplex"}];

Graphics3D[{object},Axes->True]

Graphics3D[{EdgeForm[None],object},Axes->True,AxesLabel->{"x","y","z"}]

Поворот объекта в пространстве

Поворот объекта вокруг оси y

Graphics3D[
{
EdgeForm[None],
GeometricTransformation[object,RotationMatrix[30°,{0,1,0}]]
},
Axes->True,AxesLabel->{"x","y","z"}
]

Поворот объекта в пространстве

Матрица сложного поворота

A123=.;
A123[ψ_,ϑ_,φ_]:=RotationMatrix[ψ,{1,0,0}].RotationMatrix[ϑ,{0,0,1}].RotationMatrix[φ,{1,0,0}];

Manipulate[
Graphics3D[
{
{Green,Arrow[{{0,0,0},{0,1,0}}]},
{Red,Arrow[{{0,0,0},{1,0,0}}]},
EdgeForm[None],
GeometricTransformation[{object,Arrow[{{0,0,0},{0,1,0}}]},A123[ψ,ϑ,φ]]
},
Axes->False,AxesLabel->{"x","y","z"},PlotRange->{{-1,1},{-1,1},{-1,1}},Boxed->False
],
{ψ,0,2π},{ϑ,0,2π},{φ,0,2 π}]

Анимация

Animate[
Graphics3D[
{
EdgeForm[None],
GeometricTransformation[object,{A123[ψ[t],ϑ[t],φ[t]]/.sol,{0,0,0}}]
},
Axes->True,AxesLabel->{"x","y","z"},PlotRange->{{-1,1},{-1,1},{-1,1}}
]/.t->tt,
{tt,0,tk,0.2},DisplayAllSteps->True]

Экспорт анимации в файл (gif)

Анимированный GIF

Формируем список отдельных кадров

aniTable = Table[
   Graphics3D[
     {
      EdgeForm[None],
      GeometricTransformation[object, {A123[ψ[t], ϑ[t], φ[t]] /. sol, {0, 0, 0}}]
      },
     Axes -> True, AxesLabel -> {"x", "y", "z"}, PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}
     ] /. t -> tt,
   {tt, 0, tk, 0.5}];

Экспортируем в файл

Export[NotebookDirectory[]<>"res.gif",aniTable];

Можно и в avi (файл большого размера)

(*
Export[NotebookDirectory[]<>"res.avi",aniTable];
*)