Магія тензорною алгебри: Частина 18 — Математичне моделювання ефекту Джанібекова

Введення
Минула стаття повинна була бути про чисельному моделюванні ефекту Джанібекова, але мені раптом прийшла в голову думка, що цей ефект можна досліджувати якісно, нехай і досить наближеним першим методом Ляпунова. Однак, чисельне моделювання теж дуже цікаве питання, тим більше що лежить у площині моїх дослідницьких завдань. Тому, сьогодні ми
  1. Остаточно визначимося з тим, як використовувати параметри Родріга-Гамільтона для опис орієнтації тіла в просторі
  2. Розглянемо форми представлення рівнянь руху вільного тіла: покажемо як тензорні рівняння можна перетворити на матричні і компонентні.
  3. Виконаємо моделювання руху вільного твердого тіла при різних співвідношеннях між головними моментами інерції і покажемо, як проявляє себе ефект Джанібекова.

1. Диференціальні рівняння вільного руху в тензорною формі
Ми вже не раз розглядали ці рівняння у векторному вигляді

\begin{align*}
&m \, \vec a_{c} = \sum \vec F_k^{\,e} \\
&\mathbf I_c \, \vec\epsilon + \vec\omega \times \left(\mathbf I_c \, \vec\omega \right) = \sum \vec M_{c}(\vec F_k^{\,e}) 
\end{align*}
Векторна форма запису зручна для загального аналізу характеру залежностей, вона звична і в ній видно, що означає конкретне доданок. Однак, для подальшого перетворення рівнянь у форму, зручну для моделювання, перейдемо до тензорною запису

\begin{align*}
&m \, \left(\ddot x^{\,i} + \Gamma_{\,kl}^{\,i} \, \dot x^{\,k} \, \dot x^{\,l} \right) = X^{\,i} \\
&I_{\,j}^{\,i} \, \dot\omega^{j} + \varepsilon^{\,ijk} \, \omega_{\,j} \, g_{\,kl} \, I_{\,p}^{\,l} \, \omega^{\,p} = M^{i}
\end{align*}
де x^{\,i}— контравариантные координати центру мас тіла; X^{\,i}— контравариантные компоненти головного вектора зовнішніх сил, прикладених до тіла; M^{\,i}— контравариантные компоненти головного моменту зовнішніх сил, прикладених до тіла.

Система рівнянь (2) вже є замкнутою, інтегруючи її можна отримати закон руху центра мас і залежність кутової швидкості тіла від часу. Але, нас ще буде цікавити орієнтація тіла, тому доповнимо цю систему рівнянь

2 \, \varepsilon^{\,ijk} \, \lambda_{\,j} \, \dot\lambda_{\,k} + 2 \, \lambda_0 \, \dot\lambda^{\,i} - 2 \, \lambda^{\,i} \, \dot\lambda_0 = \omega^{\,i}
Рівняння (3) є ніщо інше як подання компонент кутової швидкості через параметри орієнтації Родріга-Гамільтона. Це вираження ми вже отримували у попередніх статтях. Тепер ми будемо розглядати його як диференціальне рівняння, що зв'язує параметри орієнтації з компонентами кутової швидкості.

Однак, параметри Родріга-Гамільтона є надлишковими — їх чотири, а для опису орієнтації тіла в просторі досить трьох координат. І число невідомих у системі (2), (3) перевищує число рівнянь на одиницю. Значить нам доведеться доповнити рівняння (2) і (3) рівнянням зв'язку між параметрами орієнтації. статті про параметри Родріга-Гамільтона ми показали, що поворот тіла зручно описувати одиничним кватернионом, що є

\lambda_0^2 + \lambda_1^2 + \lambda_2^2 + \lambda_3^2 = 1
або, в тензорному вигляді

\lambda_0^2 + \lambda_{\,i} \, \lambda^{\,i} = 1
Продиференціюємо (4) за часом

2\, \lambda_0 \, \dot\lambda_0 + \dot\lambda_{\,i} \, \lambda^{\,i} + \lambda_{\,i} \, \dot\lambda^{\,i} = 0
З урахуванням комутативності скалярного добутку вважаємо \dot\lambda_{\,i} \, \lambda^{\,i} = \lambda_{\,i}\,\dot\lambda^{\,i}, тоді

\lambda_0 \, \dot\lambda_0 + \lambda_{\,i} \, \dot\lambda^{\,i} = 0
і є шукане рівняння зв'язку. Повна система рівнянь руху вільного твердого тіла в тензорною формі буде мати вигляд

\begin{align*}
&\dot x^{\,i} - v^{\,i} = 0 \\
&m \, \left(\dot v^{\,i} + \Gamma_{\,kl}^{\,i} \, v^{\,k} \, v^{\,l} \right) = X^{\,i} \\
&2 \, \varepsilon^{\,ijk} \, \lambda_{\,j} \, \dot\lambda_{\,k} + 2 \, \lambda_0 \, \dot\lambda^{\,i} - 2 \, \lambda^{\,i} \, \dot\lambda_0 - \omega^{\,i} = 0 \\
&\lambda_0 \, \dot\lambda_0 + \lambda_{\,i} \, \dot\lambda^{\,i} = 0 \\
&I_{\,j}^{\,i} \, \dot\omega^{j} + \varepsilon^{\,ijk} \, \omega_{\,j} \, g_{\,kl} \, I_{\,p}^{\,l} \, \omega^{\,p} = M^{i}
\end{align*}
Досить страшнувато — (6) містить 13 нелінійних диференціальних рівнянь першого порядку з 13 невідомими величинами. Страшно виглядає з-за загальною тензорною запису, але при переході до конкретних координат, в нашому випадку декартовым, система (6) значно спроститься.

2. Матрична форма диференціальних рівнянь руху твердого тіла в декартовом базисі
Введемо вектор-стовпець фазових координат тіла

\mathbf y = 
\begin{bmatrix}
\mathbf x^T && \mathbf q^T && \mathbf v^T && \mathbf \omega^T
\end{bmatrix}^T
де \mathbf x = \begin{bmatrix} x_c && y_c && z_c \end{bmatrix}^T\mathbf v = \begin{bmatrix} v_{cx} && v_{cy} && v_{cz} \end{bmatrix}^T— положення і швидкість центру мас тіла; \mathbf q = \begin{bmatrix} \lambda_0 && \lambda_1 && \lambda_2 && \lambda_3 \end{bmatrix}^T\mathbf \omega = \begin{bmatrix} \omega_x && \omega_y && \omega_z \end{bmatrix}^T— орієнтація та кутова швидкість тіла.

У декартовом базисі метричний тензор представлений одиничною матрицею а символи Кристоффеля дорівнюють нулю, тому система рівнянь (6) в матричній формі запишеться так:

\begin{align*}
&\mathbf{\dot x} - \mathbf v = \mathbf 0 \\
& 2\, \mathbf B \, \mathbf{\dot q} - \mathbf D \, \mathbf \omega = \mathbf 0 \\ 
&m \mathbf{\dot v} = \mathbf X \\
&\mathbf I_c \, \mathbf{\dot\omega} + \mathbf \Omega \, \mathbf I_c \, \mathbf \omega = \mathbf M_c
\end{align*}
де введені матриці

\mathbf B = 
\begin{bmatrix}
\lambda_0 && \lambda_1 && \lambda_2 && \lambda_3 \\
-\lambda_1 && \lambda_0 && -\lambda_3 && \lambda_2 \\
-\lambda_2 && \lambda_3 && \lambda_0 && -\lambda_1 \\
-\lambda_3 && -\lambda_2 && \lambda_1 && \lambda_0 \\
\end{bmatrix}, \quad \mathbf D = 
\begin{bmatrix}
\mathbf 0^T \\
\mathbf E
\end{bmatrix}, \quad \mathbf \Omega = 
\begin{bmatrix}
0 && -\omega_z && \omega_y \\
\omega_z && 0 && -\omega_x \\
-\omega_y && \omega_z && 0
\end{bmatrix}
Вирішуючи систему (7) стосовно перших похідних, отримуємо

\begin{align*}
&\mathbf{\dot x} = \mathbf v \\
& \mathbf{\dot q} = \frac{1}{2} \, \mathbf B^{-1} \, \mathbf D \, \mathbf \omega \\ 
&\mathbf{\dot v} = \frac{1}{m} \, \mathbf X \\
&\mathbf{\dot\omega} = \mathbf I_c^{-1} \, \left(\mathbf M_c - \mathbf \Omega \, \mathbf I_c \, \mathbf \omega \right)
\end{align*}
систему рівнянь руху у формі Коші.

3. Моделювання ефекту Джанібекова
У відсутність зовнішніх силових факторів права частина системи (8) дорівнює нулю, то рівняння руху центру мас інтегрується легко, з урахуванням початкових умов

\mathbf x(t) = \mathbf x_0 + \mathbf v_0 \, t
Обертання гайки описується системою семи рівнянь першого порядку, які отримуємо з (8), вводячи безрозмірні моменти інерції i_y = \frac{I_y}{I_x}і i_z = \frac{I_z}{I_x}

\begin{align*}
&\dot\lambda_0 = -\frac{1}{2} \, \left( \lambda_1 \, \omega_x + \lambda_2 \, \omega_y + \lambda_3 \, \omega_z \right) \\
&\dot\lambda_1 = \frac{1}{2} \, \left( \lambda_0 \, \omega_x + \lambda_3 \, \omega_y - \lambda_2 \, \omega_z \right) \\ 
&\dot\lambda_2 = -\frac{1}{2} \, \left( \lambda_3 \, \omega_x - \lambda_0 \, \omega_y - \lambda_1 \, \omega_z \right) \\ 
&\dot\lambda_3 = \frac{1}{2} \, \left( \lambda_2 \, \omega_x - \lambda_1 \, \omega_y + \lambda_0 \, \omega_z \right) \\
&\dot\omega_x = \left(i_y - i_z) \, \omega_y \, \omega_z \\
&\dot\omega_y = \frac{i_z - 1}{i_y} \, \omega_x \, \omega_z \\
&\dot\omega_z = \frac{1 - i_y}{i_z} \, \omega_x \, \omega_y
\end{align*}
Для чисельного інтегрування системи (9) задамо початкові умови

\begin{align*}
&\lambda_0(0) = 1, \quad \lambda_1(0) = \lambda_2(0) = \lambda_3(0) = 0 \\
&\omega_x(0) = \omega_0, \quad \omega_y(0) = \Delta\omega_y, \quad \omega_z(0) = 0
\end{align*}
де \omega_0— кутова швидкість гайки після сходження з різьби; \Delta\omega_y— початкова обурення кутової швидкості

При значеннях параметрів i_y = 2.0, \, i_z = 0.5, \, \omega_0 = 1.0, рад/с, \Delta\omega_y = 1 \cdot 10^{-10}, рад/с, рух гайки відбувається наступним чином

Параметри орієнтації Родріга-Гамільтона








Проекції кутової швидкості на власні осі



З графіків видно, що при I_y > I_x > I_z, досить незначне обурення вектора кутової швидкості призводить до періодичного лавиноподібного зміни орієнтації гайки в просторі.

Порівняємо отриманий результат з рухом тіла, закручених навколо осі з максимальним моментом інерції, тобто покладемо I_x > I_y = I_z, задавши наступні значення параметрів i_y = 0.5, \, i_z = 0.5, \, \omega_0 = 1.0, рад/с, \Delta\omega_y = 1 \cdot 10^{-2}, радий/з

Параметри орієнтації Родріга-Гамільтона









Проекції кутової швидкості на власні осі



Видно, що при досить значному обуренні кутової швидкості рух залишається стійким обертанням навколо осі xз невеликою прецесією.

Схожа картина спостерігається для тіла, закручиваемого навколо осі з мінімальним моментом інерції (I_x < I_y = I_z) i_y = 1.5, \, i_z = 1.5, \, \omega_0 = 1.0, рад/с, \Delta\omega_y = 1 \cdot 10^{-2}, радий/з

Параметри орієнтації Родріга-Гамільтона









Проекції кутової швидкості на власні осі



Частота прецесії істотно менше, ніж при закрутки навколо осі з максимальним моментом інерції, що логічно, так як коливання відбуваються навколо осі з великим моментом інерції, ніж у випадку I_x > I_y = I_z.

Висновок
Всі розрахунки виконані автором в СКА Maple 18. Графіки побудовані з лода розрахунку зв'язкою Kile + LaTeX + gnuplot.

Хотілося б ще зробити анімацію, однак досвід автора у цьому питанні вкрай малий. Тому хотів би задати питання читачам — існує для Linux/Windows), за допомогою якого маючи набір значень параметрів кватерниона орієнтації в залежності від часу зробити анімаційний ролик, що ілюструє рух тіла? Підозрюю, що подібне можна провернути з Blender 3D, але не впевнений.

А поки що, дякую за увагу!

далі буде...

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.