среда, 24 августа 2016 г.

Моделирование вращения твёрдого тела: С++, SFML, OpenGL.

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

Довольно давно я уже сталкивался с библиотекой OpenGL и даже участвовал в проекте, где в десктопном приложении была реализована 3D визуализация с использованием этой библиотеки, но это было очень давно. Тогда последняя версия OpenGL была 2.0, а в нашем коде, если мне не изменяет память, использовалась версия 1.1. С тех пор утекло много воды и я неожиданно обнаружил, что мои знания об OpenGL сильно устарели: в уроках OpenGL так и написано "Забудьте все, что вы знали об OpenGL ранее, если ваши знания касаются glBegin() и подобных функций". Что ж, программировать трёхмерную графике по крайней мере весело! Почему бы и не обновить знания?

При создании анимации, компьютерных игр или решении инженерных задач может потребоваться промоделировать и визуализировать движение каких-то твёрдых предметов. Хотя сейчас для компьютерного моделирование механики существуют разнообразные физические движки, которые с успехом решают эти задачи, всё равно интересно и полезно представлять как движутся твёрдые тела. Здесь мы рассмотрим свободное движение твёрдого тела, на которое не действуют никакие силы и моменты. Этакий брошенный и свободно летящий "кирпич". В этом случае движение можно разделить на два компонента: поступательное движение и вращательное.

Свободное вращательное движение в общем случае более сложное чем вращение вокруг фиксированной оси: объект крутится вокруг какой-то оси, при этом сама ось вращения тоже совершает цикличекие движения.

Я воспользуюсь следующими известными из курса механики фактами:
1. Момент импульса объекта не меняется, так как на него не действуют моменты внешних сил:
2. Момент импульса связан с угловой скоростью формулой $\vec{L} = \hat{J} \vec\omega $, где $J$ - тензор инерции. Следовательно, угловую скорость в каждый момент времени можно выразить формулой $\vec\omega = \hat{J}^{-1}\left(t\right)\vec{L}$.
3. Ориентацию тела в пространстве можно описать ортогональной матрицей $\vec{Q}$ с положительным определителем.

При расчётах будем использовать неподвижную систему координат. В ней все компоненты вектора $\vec{L}$ будут постоянными, а компоненты тензора инерции $\hat{J}_{ij}$ будут меняться, так как будет меняться ориентация тела. Угловую скорость в каждый момент времени можно определить по формуле $\vec\omega = \hat{J}^{-1}\left(t\right)\vec{L}$. Зная угловую скорость можно примерно определить ориентацию тела через короткий промежуток времени $\Delta t$.


 Соединив всё это воедино полуаем такой алгоритм:
1. Инициализируем тензор инерции тела, момент имульса, матрицу ориентации тела:
$$\hat{J^{-1}_0}=\left(\begin{array}{ccc}1/J_x & 0 & 0 \\
0 & 1/J_y & 0 \\
0 & 0 & 1/J_z \end{array}\right) $$

$$Q_0=\left(\begin{array}{ccc}1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \end{array}\right) $$

2. Делаем следующие итерации:
$$\vec{\omega_{\tau}} = \hat{J^{-1}_\tau} \vec{L}$$
$$dQ_\tau = R\left(\vec{n}_\omega, |\vec\omega|\Delta t\right)$$
$$\hat{J^{-1}_{\tau+1}} = dQ_\tau \hat{J^{-1}_\tau} dQ_\tau^{-1}$$
$$Q_{\tau+1} = dQ_\tau Q_\tau $$


Где $\vec{n}_\omega$ единичный вектор, совпадающий с направлением $\vec{\omega}$, $|\vec \omega|$ - абсолютная величина вектора $\vec{\omega}$. Угол, на который поворачивается тело за время $\Delta t$ примерно равен $|\vec\omega|\Delta t$. Матрица трёхмерного поворота вокруг заданной оси $\vec n$ на заданный угол $\alpha$ обозначена как $R\left(\vec n, \alpha\right)$. Выражение для неё несложно найти, например, в википедии. В библиотеке glm есть готовая функция для создания матрицы поворота по заданной оси и углу: glm::rotate

Реализация. 

Я пошёл по пути наименьшего сопротивления и просто объединил два примера кода из сети: пример из туториала по OpenGL (есть русская версия), в котором создаётся трёхмерное изображение куба и код из туториала SFML, в котором показывается как инициализировать окно OpenGL.

Для работы с матрицами и векторами используем библиотеку glm. Она тесно интегрирована с OpenGL и уже и так используется в примере с кубом. В то же время вполне подходит и для любых других других применений если вам не нужны матрицы размером больше чем 4x4 и вектора размером больше чем 4. 

Хотя в расчётах используется только лишь обратный тензор инетрции , я пересчитываю на каждом щаге и прямой тензор. Может, потом пригодится.

Матрица ориентации $Q$ это и есть матрица модели. В приемере со статическим кубом из OpenGL туториала матрицы модели всегда единичная, а матрицы модель-вид-проекция (MVP) вычисляется так:
В нашем случае матрица модели должна пересчитываться при каждой перерисовке кадра:
Операция glm::scale тоже часть вычисления матрицы модели. Она растягивает наш куб в прямоугольный параллелепипед в соответствии с главными моментами инерции.

Результат. 

Исходники доступны на GitHub. Хотя я создавал этот проект в Microsoft Visual Studio, код должен быть переносим на Linux и OS X, хотя я не пробовал собирать его в других системах.

 

В результате видим, что параллелепипед вращается, и движение зрительно похоже на то, что ожидаешь увидеть. Так же воспроизводится эффект гайки Джанибекова: закрученное вокруг оси близкой к оси со средним моментом инерции тело периодически "переворачивается". Однако если понаблюдать подольше, видно, что симуляция неустойчива: через какое-то время вращение уходит от оси со средним моментом инерции и стабилизируется вокруг оси с максимальным или минимальным моментом инерции. По моему убеждению это артефакт используемого вычислительного метода. В следующих постах я попробую бороться с этой неустойчивостью.