Эффект Джанибекова в SimInTech
Модель и визуализация движения в SimInTech гайки-барашка при её вращении вокруг оси с промежуточным моментом инерции, известного как “Эффект Джанибекова”.
Среда динамического моделирования SimInTech предлагает простой набор функций для визуализации движения моделируемых объектов. В этой статье рассматривается построение модели и визуализации движения гайки-барашка при ее свободном вращении вокруг собственной оси с “промежуточным” моментом инерции. Это движение известно как “Эффект Джанибекова”.
Эффект Джанибекова
Эффектом Джанибекова называют необычное движение вращающейся вокруг одной из своих осей гайки-барашка, выражающееся в её периодическом быстром перевороте на 180 градусов. Этой яркое явление наблюдал космонавт Владимир Александрович Джанибеков на борту орбитальной стации Салют-7 в 1985 году. На рисунке 1 показано движение гайки-барашка, построенное при помощи описываемой далее модели в среде динамического моделирования SimInTech.
Вокруг этого явления возникло множество различных безумных теорий, например, встречалось объяснение этого эффекта “на основании философии Дуализма Диалектики Абсолютного Парадокса” (все буквы, конечно, заглавные), однако этот эффект описывается в рамках классической механики. Вращение гайки-барашка вокруг оси с промежуточным или средним моментом инерции, например вокруг собственной оси Y, если значения главных моментов инерции удовлетворяют условиям \(J_z < J_y < J_x\), является неустойчивым, что и приводит к наблюдаемому перевороту гайки. Некоторой аналогией движения гайки является поведение физического маятника, выведенного небольшим возмущением из неустойчивого (верхнего) положения равновесия.
Мятник будет совершать размашистые колебания с амплитудой чуть меньше 180 градусов, подходя к точке неустойчивого равновесия то слева то справа, надолго “зависая” в этом положении. Нижнее устойчивое положение маятник будет проходить с большой скоростью. Подобным же образом ведет себя гайка, но гайка, в отличие от маятника, имеет устойчивые и неустойчивые движения, а не положения, поэтому у гайки длительные движения в окрестности неустойчивоой оси Y сменяются быстрыми переворотами на 180 градусов.
Анализ устойчивости
Рассмотрим движение твердого тела, которое имеет три различных главных момента инерции, при этом обозначения главных центральных осей тела выбраны так, чтобы выполнялись условия \(J_z \leq J_y \leq J_x\).
Вращение вокруг оси с максимальным значением момента инерции
Предположим, что тело вращающется вокруг главной оси с максимальным моментом инерции – $Cx$, с угловой скоростью \(\omega_0 = \text{const}\). Пусть проекции угловой скорости на все собственные оси получили малые возмущения, которые обозначим как \(\delta_x\), \(\delta_y\) и \(\delta_z\) соответственно: \[\omega_x = \omega_0 + \delta_x, \; \ \omega_y = 0+\delta_y, \; \omega_z = 0 + \delta_z.\]
Уравнения движения тела в возмущениях имеют вид: \[\left\{ \begin{aligned} & J_x \dot \delta_x - (J_y-J_z) \delta_y \delta_z = 0, \\ & J_y \dot \delta_y - (J_z-J_x) \delta_z (\omega_0 + \delta_x) = 0, \\ & J_z \dot \delta_z - (J_x-J_y) \delta_y (\omega_0 + \delta_x) = 0. \\ \end{aligned} \right.\]
Раскрывая скобки и пренебрегая произведениями малых отклонений (малыми второго порядка), получим: \[\left\{ \begin{aligned} & J_x \dot \delta_x = 0, \\ & J_y \dot \delta_y - (J_z-J_x) \delta_z \omega_0 = 0, \\ & J_z \dot \delta_z - (J_x-J_y) \delta_y \omega_0 = 0. \\ \end{aligned} \right.\]
Из 1-го уравнения следует, что \[\delta_x = \text{const}.\]
Из 2-го и 3-го уравнений можно получить два независимых линейных дифференциальных уравнения второго порядка: \[J_y \ddot \delta_y - \omega_0^2 \frac{(J_z-J_x) (J_x-J_y)}{J_z}\delta_y = 0.\] \[J_z \ddot \delta_z - \omega_0^2 \frac{(J_z-J_x) (J_x-J_y)}{J_y}\delta_z = 0.\]
Вид решений последних двух линейных дифференциальных уравнений зависит от знака множителя при \(\delta_y\) и \(\delta_z\). При \(J_z<J_y<J_x\): \[\omega_0^2 \frac{(J_x-J_z) (J_x-J_y)}{J_y J_z} > 0,\]
следовательно решения будут иметь вид: \[\delta_y = a_y \sin (\lambda t + \varepsilon_y), \quad \delta_z = a_z \sin (\lambda t + \varepsilon_z)\]
Т.е. отклонения \(\delta_y\) и \(\delta_z\) не будут возрастать и вращение вокруг оси \(Cx\) будет устойчивым. Подобным образом показывается устойчивость вращения тела вокруг главной оси с минимальным моментом инерции \(Cz\).
Вращение вокруг оси с промежуточным значением момента инерции
Теперь предположим, что тело вращается вокруг оси \(Cy\) с промежуточным значением момента инерции. Рассмотрим устойчивость такого движения при небольших возмущениях проекций угловой скорости: \[\omega_x = \delta_x, \; \ \omega_y = \omega_0+\delta_y, \; \omega_z = \delta_z.\]
Уравнения движения тела имеют вид: \[\left\{ \begin{aligned} & J_x \dot \delta_x - (J_y-J_z) (\omega_0 + \delta_y) \delta_z = 0, \\ & J_y \dot \delta_y - (J_z-J_x) \delta_z \delta_x = 0, \\ & J_z \dot \delta_z - (J_x-J_y) \delta_x (\omega_0 + \delta_y) = 0. \\ \end{aligned} \right.\]
Пренебрегая малыми второго порядка, получим: \[\left\{ \begin{aligned} & J_x \dot \delta_x - (J_z-J_x) \delta_z \omega_0 = 0, \\ & J_y \dot \delta_y = 0, \\ & J_z \dot \delta_z - (J_x-J_y) \delta_x \omega_0 = 0. \\ \end{aligned} \right.\]
Из 2-го уравнения следует: \(\delta_y = \text{const}\). Из 1-го и 3-го уравнений: \[\ddot \delta_x + \omega_0^2 \frac{(J_z-J_y) (J_x-J_y)}{J_x J_z} \delta_x = 0.\] \[\ddot \delta_z + \omega_0^2 \frac{(J_z-J_y) (J_x-J_y)}{J_x J_z} \delta_z = 0.\]
При \(J_z<J_y<J_x\): \[\omega_0^2 \frac{(J_z-J_y) (J_x-J_y)}{J_x J_z} < 0,\]
следовательно решения для возмущений будут иметь вид: \[\delta_x = C_1 e^{\lambda t} + C_2 e^{- \lambda t}, \quad \delta_y = C_3 e^{\lambda t} + C_4 e^{- \lambda t}\]
Мы видим, что возмущения неограниченно возрастают, поэтому вращение вокруг оси \(y_2\) будет неустойчиво.
Полученные уравнения показывают только устойчивость или неустойчивость движения твердого тела при вращении вокруг одной из собcтвенных осей и не описывают само движение при перевороте гайки. Среда моделирования SimInTech позволит нам не только убедиться в неустойчивости движения вокруг оси с промежуточным моментом инерции, построив графики изменения угловой скорости гайки, но и “увидеть” это движение при помощи модуля 3D-визуализации.
Модель
Общий вид модели в SimInTech показан на рисунке 3.
Главным блоком модели вляется новый блок “Гайка (XZ’X’’)”, созданный на основе субмодели. Этот блок может использоваться для моделирования движения твердого тела вокруг центра масс в углах Эйлера. На вход блока подаётся вектор внешнего момента, который в этой задаче равен нулю, и массив значений главных моментов инерции объекта-гайки. Выходными сигналами блока являются вектор проекций угловой скорости гайки на собственные оси [wx, wy, wz] и вектор углов поворота связанной с гайкой системы координат по отношению к неподвижной системе. Угловое положение системы координат, связанной с гайкой, определяется углами Эйлера последовательности поворотов XZ’X’’. Структура блока показа на следующем рисунке:
В блоке “Гайка (XZ’X’’)” используются блоки-интеграторы, начальные условия которых (начальная угловая скорость и начальные значения углов Эйлера) заданы в свойствах блока “Гайка (XZ’X’’)”:
Для того, чтобы добавить эти свойства к субмодели необходимо перевести SimInTech в режим разработчика (Меню “Вид” - “Режим разработчика”), в общих свойствах субмодели “Гайка” указать новое наименование типа блока, теперь это не “стандартная” субмодель, а блок типа “ТвердоеТело131”. Таким образом на основе субмодели создается новый блок со своими уникальными свойствами.
Блок динамических уравнений
В блоке “Динамические уравнения” в явном виде на встроенном языке программирования записываются динамические уравнения Эйлера в проекциях на главные центральные оси тела: \[\left\{ \begin{aligned} & J_x \dot \omega_x - (J_z-J_x) \omega_y \omega_z = M_x, \\ & J_y \dot \omega_y - (J_z-J_x) \omega_z \omega_x = M_y, \\ & J_z \dot \omega_z - (J_x-J_y) \omega_x \omega_y = M_z. \\ \end{aligned} \right.\]
В рассматриваемом случае движение тела происходит “по-инерции”, поэтому все проекции внешнего момента на собственные оси равны нулю \(M_x = M_y = M_z = 0\). Ниже приведен код блока “динамические уравнения” (блок “Язык программирования” из раздела “Динамические”):
// Динамические уравнения Эйлера
// в главных центральных осях
// Входы :
// Torque - Вектор момента в главных центральных осях [Нм]
// J - Главные моменты инерции [кгм2]
// w - угловые скорости в проекциях на главные оси [1/с]
input Torque[3], J[3], w[3];
// Выходы :
// dw - Вектор угловых ускорений [1/с2]
output dw[3];
dw[1] = (Torque[1]+(J[2]-J[3])*w[2]*w[3])/J[1];
dw[2] = (Torque[2]+(J[3]-J[1])*w[3]*w[1])/J[2];
dw[3] = (Torque[3]+(J[1]-J[2])*w[1]*w[2])/J[3];
На вход блока подаётся вектор проекций угловой скорости тела, на выходе получается вектор проекций угловых ускорений тела (правая часть дифференциальных уранвений).
Блок кинематических уравнений
Кинематические уравнения связывают производные параметров, определяющих ориентацию базиса, связанного с твердым телом по отношению к неподвижному базису, с проекциями угловых скоростей. В рассматриваемой задаче используется система углов Эйлера последовательности XZ’X’’, т.е. из исходного положения, когда связанный с телом базис совпадает с неподвижным базисом, угловое положение первого определяется тремя последовательными поворотами:
- вокруг оси \(Cx_2\) на угол прецессии \(\psi\);
- вокруг оси \(Cz_2\) на угол нутации \(\vartheta\);
- вокруг оси \(Cx_2\) на угол собственного вращения \(\varphi\).
Кинематические уравнения для рассматриваемой последовательности поворотов: \[\left\{ \begin{aligned} \dot \psi & = \frac{\omega_z \sin \varphi - \omega_y \cos \varphi}{\sin \vartheta}, \\ \dot \vartheta & = \omega_z cos \varphi + \omega_y \sin \varphi, \\ \dot \varphi & = \omega_x + \frac{\omega_y \cos \varphi - \omega_z \sin \varphi}{\tan \vartheta} \end{aligned} \right.\]
Блок кинематических уравнений написан также на встроенном языке программирования (блок “Язык программирования” из раздела “Динамические”):
// Кинематические уравнения для последовательности поворотов 131
input w[3], a[3];
output da[3];
da[1] = (w[3]*sin(a[3])-w[2]*cos(a[3]))/sin(a[2]);
da[2] = w[3]*cos(a[3])+w[2]*sin(a[3]);
da[3] = w[1]+w[2]*cos(a[3])/tg(a[2])-w[3]*sin(a[3])/tg(a[2]);
Переход от углов последовательности 131 к углам последовательности 232
В модуле 3D-визуализации SimInTech ориентация объекта в пространстве также задается тремя углами Эйлера, определяющими ориентацию связанного с телом базиса по отношению к неподвижному базису, но используется последовательность YZ’Y’‘, т.е. первый поворот выполняется вокруг оси Y на угол прецессии, второй поворот вокруг новой оси Z (поэтому она обозначается как Z’) на угол нутации и третий поворот на угол собственного вращения вокруг оси Y’’.
К сожалению, в справке описан только этот способ определения ориентации тела в модуле 3D-визуализации, поэтому если в разрабатываемой модели используется другой набор углов, как в рассматриваемом случае, или вместо углов используются катернионы или направляющие косинусы, то для визуализации движения необходимо выполнять пересчет используемых параметров в углы последовательности YZ’Y’’.
Для преобразования одного набора углов Эйлера в другой используется субмодель “Углы 131 в 232”, структура которой показана на следующем рисунке.
На вход субмодели поступает массив (вектор) углов последовательности XZ’X’’ (или 131). Внутри субмодели по значениям углов формируются матрицы поворота вокруг оси X на угол прецессии: \[A_x(\psi) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \psi & -\sin \psi \\ 0 & \sin \psi & \cos \psi \end{bmatrix}\]
вокруг оси Z на угол прецессии: \[A_z(\vartheta) = \begin{bmatrix} \cos \vartheta & -\sin \vartheta & 0 \\ \sin \vartheta & \cos \vartheta & 0 \\ 0 & 0 & 1 \end{bmatrix}\]
и вокруг оси X на угол собственного вращения \[A_x(\varphi) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \varphi & -\sin \varphi \\ 0 & \sin \varphi & \cos \varphi \end{bmatrix}\]
Матрицы перемножаются в прямом порядке: \[A = A_x(\psi) \cdot A_z(\vartheta) \cdot A_x(\varphi)\]
В результате перемножения трех матриц элементарных поворотов получается матрица направляющих косинусов базиса связанного с телом по отношению к неподвижному базису.
Матрица направляющих косинусов поступает на вход субмодели DCM-A232, которая написана на встроенном языке программирования. В субмодели определяются углы последовательности 232 по элементам полученной матрицы направляющих косинусов. Ниже приведён код субмодели DCM-A232:
// Вычисление углов последовательности 232
// по элементам матрицы направляющих косинусов
input u[3,3];
output a[3];
// Для каждого угла есть два решения
a21 = arccos(u[2,2]);
a22 = 2*pi-a21;
a31 = atan2(u[2][3]/sin(a21),u[2][1]/sin(a21));
a32 = atan2(u[2][3]/sin(a22),u[2][1]/sin(a22));
a11 = atan2(u[3][2]/sin(a21),-u[1][2]/sin(a21));
a12 = atan2(u[3][2]/sin(a22),-u[1][2]/sin(a22));
// Выбираем первое...
a[1] = a11;
a[2] = a21;
a[3] = a31;
Найденные углы передаются на выход субмодели “Углы 131 в 232”. С выхода значения углов записываются в базу сигналов: переменные A1, A2, A3 соответсенно, которые были предварительно объявленны в редакторе сигналов (меню “Сервис” → “Сигналы”):
Для записи сигналов используется блок “Запись в список сигналов”:
Скрипт модели
В скрипте модели в блоке инициализации задаются значения главных моментов инерции гайки, создается окно 3D-визуализации, загружается объект - гайка из файла screw.stl.
Вне блока инициализации записаны две строчки кода, которые будут вызываться на каждом шаге интегрирования. Функция Viewer3DSetLCSAngles поворачивает гайку на углы Эйлера, полученные на шаге интегрирования, функция Viewer3DSetLCSPosition перемещает центр масс гайки вдоль оси x имитируя её движение после схода с винта.
initialization
var Jx = 2000;
var Jy = 4000;
var Jz = 6000;
Id = Viewer3DCreate;
OId = Viewer3DPlotObject(Id, 0, 4, [0, 0, 0], [0, 0, 0], "screw.stl");
Viewer3DShowObjectAxes(Id, OId, true);
end
Viewer3DSetLCSAngles(Id, OId, A1*180/pi, A2*180/pi, A3*180/pi);
Viewer3DSetLCSPosition(ID, oID, -0.2+0.05*time, 0, 0);
Перед запуском модели необходимо в настройках “Параметры расчета” на закладке “Синхронизация” поставить галочку “Синхронизировать с реальным временем”.
После запуска модели появится окно “Модуль 3D Визуализации SimInTech”, в котором будет показано движение гайки. Графики изменения проекций вектора угловой скорости гайки на собственные оси показаны на рисунке 15.