Моделирование динамики твердых тел и систем связанных тел с механическими соударениями

ДИПЛОМНАЯ РАБОТА НА ТЕМУ: «Моделирование динамики твердых тел и систем связанных тел с механическими соударениями»

Введение

Задачи моделирования динамики реальных объектов окружающего мира всегда были актуальны, и для их решения было разработано много методов и моделей. Однако когда нужно рассчитывать динамику в режиме реального времени, то существующие методы оказываются малопригодными. Например, для обучения операторов роботов-манипуляторов, предназначенных для работы с предметами (например, разгребание завалов) необходимы специальные тренажеры, в реальном времени рассчитывающие динамику виртуальных объектов. Для этого нужны специальные методы. Сейчас такие методы интенсивно развиваются в рамках компьютерных игровых программ.

Игровая индустрия интенсивно развивается в последние годы и все чаще становится двигателем прогресса технологий, например в области трехмерной графики. В игровом проекте любого жанра всегда создается модель игрового мира и, как правило, много внимания уделяется тому, как этот мир будет выглядеть визуально. Это стимулирует такие технологии, как динамический расчет освещения, теней, моделирование оптических свойств материалов и т. д. Визуализация виртуальных сцен происходит в реальном времени (25−30 кадров в секунду). Такие технологии прогрессируют с каждым днем одновременно с ростом производительности процессоров и видеокарт, и на сегодняшний день виртуальные миры выглядят весьма достоверно.

Но визуальная реалистичность не является единственным требованием к виртуальному миру. Этот мир интерактивен и на каждое действие игрока должен адекватно реагировать. Часто виртуальная сцена представляет собой не просто статические декорации, а наполнена физическими объектами, которые взаимодействуют между собой и с окружающей статической средой. Это могут быть, например, ящики, столы и стулья, которые игрок может подвинуть или опрокинуть, занавески, которые игрок может задеть, и которые реалистично колышутся от ветра. В игре могут присутствовать механизмы с подвижными частями и т. д.

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

Если несколько лет назад реалистичная физика требовалась только в играх таких жанров, как авиасимуляторы и автосимуляторы, где подходила достаточно частная физическая модель, отвечающая требованиям этого жанра, то теперь в физическом представлении нуждаются сложные сцены, с множеством взаимодействующих объектов разных типов — для чего требуется более общий подход к физической модели игрового мира.

Главная проблема таких физических моделей состоит в том, что они не должны быть вычислительно сложными, так как расчет динамики физических объектов должен происходить в реальном времени и только частично загружать процессор. Поэтому большинство методов, разработанных в других областях, не пригодны для решения таких задач.

1. Обзор предметной области

1.1 Моделирование твердых тел В существующих моделях твердых тел, обычно тело описывается следующими понятиями:

1) Геометрия тела

Это совокупность примитивов описывающих форму тела. Как правило, набор примитивов включает сферы, цилиндры, параллелепипеды и произвольные выпуклые многогранники. Геометрия тела используется для определения столкновений с другими телами. Если же тело в текущий момент ни с чем не соприкасается, то геометрия никак не влияет на его движение.

Важно отметить, что геометрия тела в физическом представлении не обязательно совпадает с той геометрией, которая будет отображаться на экране. Объект в визуальном представлении может быть очень детализированным, но в физической модели описываться набором примитивов приближенно.

2) Масса тела

3) Координаты центра масс Позиция тела задается позицией его центра масс, поэтому геометрия тела задается относительно этой точки.

4) Матрица ориентации тела (или кватернион) Это матрица поворотов 3×3, определяющая как будет повернута геометрия, вокруг центра масс.

5) Вектор линейной скорости центра масс

6) Вектор угловых скоростей тела

7) Тензор инерции Матрица 3×3, составленная из значений моментов инерции тела относительно 3-х осей, проходящих через центр масс и позволяющая вычислить момент инерции относительно произвольной оси.

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

Если к какой-то точке тела прикладывается сила, то эта сила воздействует на скорость центра масс и скорость вращения тела с учетом плеча и моментов инерции.

Основная же задача моделирования — обработка столкновений и контактного взаимодействия нескольких тел. Для этого существует несколько методов:

метод штрафов метод импульсов

LCP метод

Метод штрафов:

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

Преимущества:

— Простота реализации Именно из-за простоты реализации и малой вычислительной сложности штрафные методы были основными, когда от физической модели не требовалось высокой точности и очень правдоподобного поведения.

Недостатки:

— Тела проникают друг в друга.

Чтобы минимизировать глубину этого проникновения, необходимы очень большие силы выталкивания, что приводит к нестабильности системы.

— Низкая точность результата.

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

— Сложность настройки Коэффициенты, подобранные для одной пары тел, могут давать плохие результаты для другой пары. Большие сложности возникают, когда несколько тел одновременно взаимодействует друг с другом.

Метод импульсов:

Этот метод наиболее популярен и большинство физических движков основано именно на этом методе.

Идея заключается в том, чтобы представить любое взаимодействие тел как удар. Даже если это длительный контакт — он представляется как множество микроударов. Преимущество перед методом штрафов заключается в том, что результирующие скорости тел рассчитывается исходя не из глубины проникновения, а из скоростей, которыми обладали тела до удара. Импульс, полученный телами при ударе, рассчитывается физически корректно, учитываются линейные и угловые скорости центров масс тел и точек контакта, коэффициент упругости, нормали к поверхности в точках контакта — при этом учитывается трение. Самого пересечения тел так и не происходит, и в итоге тела не проникают друг в друга.

Длительный контакт в этой модели ничем не отличается — например, если объект покоится на какой-то поверхности, то он в каждый шаг по времени соударяется с этой поверхностью какой-либо точкой, но эффект этого соударения достаточно мал и фактически только поддерживает тело в нужном положении. При каждом ударе учитывается трение, и когда таких микроударов много — получается такой же эффект, как от постоянной силы трения.

Преимущества:

— Тела не проникают друг в друга Импульс мгновенно меняет скорости тел, и это позволяет избежать взаимопроникновения тел.

— Высокая точность Физически обоснованный расчет импульсов при ударах позволяет избежать неоднозначности результата, в отличие от метода штрафов. И получаемый результат больше соответствует тому, который был бы получен в реальном эксперименте. При обработке длительных контактов метод также ведет себя достаточно стабильно.

Недостатки:

— Нестабильность при большом количестве контактов Покоящееся на поверхности тело на самом деле не находится в полном покое, а дрожит. Эта дрожь достаточно мала и не заметна, она зависит от величины шага по времени. Однако когда одно на другом покоятся больше 2-х тел, то дрожь становится заметнее, и иногда тела вступают в «резонанс» и подпрыгивают — т. е. система становится очень нестабильной.

— Некорректная сила трения Тело, лежащее на наклонной поверхности, сползает вниз, даже если сила трения достаточно велика. Дело в том, что в результате микроударов с поверхностью тело отскакивает под углом, а не вертикально вверх и в итоге по пилообразной траектории движется в сторону наклона поверхности.

Разработано множество искусственных способов борьбы с этими недостатками — как правило, отслеживаются такие частные ситуации и принудительно корректируются. Например, тело, которое обрело три точки опоры и его колебания упали ниже какого-то порога, «выключается» и дальше считается как статический объект. В этом случае оно не будет ползти вниз по наклонной поверхности, и на него можно положить еще несколько тел одно за другим. Но в этом случае возникают сложности с определением момента, когда это тело будет выведено из равновесия и его надо «включить».

Подобные «переключатели» хорошо работают только в отработанных частных случаях, а в любой более-менее сложной ситуации взаимодействия нескольких тел либо не срабатывают, либо вызывают некорректное поведение системы.

Метод LCP

(Linear Complementarity Problem — линейная задача сопряженности) Это аналитический подход к решению задачи столкновения. Идея метода в том, чтобы представить ограничения непроникновения тел как линейную задачу сопряженности. Записанная система из уравнений и неравенств решается методом, предложенным Лотстедом (Lotstedt), который основан на методе вращения Котла (Cottle) и Данцига (Dantzig). Если при контактах отсутствует сила трения, то решение всегда существует и единственно.

Однако введение силы трения делает задачу нелинейной, так как касательные силы трения, зависят от нормальной составляющей силы взаимодействия. В итоге система может не иметь решения или решение может быть не единственным.

Чтобы избежать вырождения задачи и добиться единственности решения применяются разнообразные хитрости. Например, вводятся дополнительные параметры «мягкости», спасающие систему от вырождения, накладываются искусственные условия, обеспечивающие единственность и т. д. Все это отрицательно сказывается на точности решения.

Преимущества:

— Тела не проникают друг в друга

— Устойчивое поведение при нескольких контактах одновременно.

В отличие от метода импульсов, LCP метод рассчитан на обработку длительных контактов и при наличии нескольких одновременных контактов не возникает нестабильного поведения типа дрожания и «резонансов».

Недостатки (возникают при наличии силы трения):

— Нет гарантии, что решение будет найдено Чем больше одновременных контактов, тем больше вероятность, что решение не будет найдено. Также возникают «тупиковые» ситуации, когда единственное найденное решение дает нулевые скорости, и в итоге тела застревают в тех случаях, когда на самом деле застревать не должны. Например, катящийся шар может застрять на небольшом сгибе поверхности, когда касается поверхности сразу двумя точками, если сила трения достаточно большая.

— Невысокая точность Способы, которыми добиваются стабильной работы алгоритма, искусственно обеспечивающие существование и единственность решения, сильно понижают точность метода и часто вызывают странное поведение объектов.

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

1.2 Моделирование систем связанных твердых тел В системе связанных твердых тел помимо самих тел существуют суставы, связывающие тела попарно.

«Сустав» — это точка, с которой связана собственная система координат, в которой оба связываемых суставом тела имеют определенное число степеней свободы. Суставы бывают следующих типов:

1) Шарнир Этот сустав дает телам 3 степени свободы, позволяет им вращаться по 3-м углам.

2) Петля Подобно дверной петле, оставляет телам одну степень свободы — они могут поворачиваться относительно одной оси.

3) Двойная петля Вариант сустава с 2-мя степенями свободы. Можно представить его как совмещение двух перпендикулярных петель, которые позволяют телам поворачиваться по двум углам.

4) Поршень Сустав, позволяющий телам скользить друг относительно друга вдоль выбранной оси.

Внешний вид суставов первых трех типов приведен на рисунке 4.

Более сложные суставы можно получить, комбинируя последовательно эти типы суставов.

Для моделирования систем связных тел разработан ряд методов, которые мы описываем ниже.

Метод ограничений (алгоритм Физерстоуна):

Это универсальный координатный метод, позволяющий ограничить степени свободы, связанных тел. Этот метод основан на вычислении относительных линейных и угловых скоростей и ускорений тел в системах координат, связанных с суставами с помощью уравнений Ньютона-Эйлера. При этом для каждого сустава вычисляется инерция, центробежные силы и силы Кориолиса (т.к. система отсчета связанная с суставом не инерционная), и эти силы передаются между звеньями. Для метода разработана специальная пространственная алгебра, без которой реализация алгоритма слишком сложна. Иерархия связанных тел может быть только древовидной (без циклов). Чтобы учесть влияние тел друг на друга, применяется рекурсивный алгоритм Физерстоуна, каждый шаг которого включает 4 прохода от корня дерева к листьям и обратно.

Преимущества:

— Высокая точность Так как алгоритм честный и физически обоснованный, то дает результаты, близкие к тем, которые могут быть получены в реальном эксперименте.

— Стабильность Алгоритм почти никогда не дает сбоев Недостатки:

— Нельзя образовывать циклы Это основное ограничение алгоритма и именно из-за этой особенности он не применим в некоторых случаях. Однако обычно таких случаев стараются избегать, либо применяют хитрости — например, к суставам, которые необходимо замкнуть, прикладывают силы навстречу друг другу пропорционально расстоянию между ними — иногда этого достаточно.

— Высокая вычислительная сложность

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

— Сложное совмещение с методами соударений тел Этот метод сам по себе не решает задачу столкновения систем связанных тел между собой. Обычно его комбинируют с методом импульсов. Но для расчета каждого соударения, приходится делать проход по всему дереву, чтобы рассчитать импульс в контактной точке, потом удар рассчитывается, как для 2-х твердых тел, а полученные результат пересчитывается на новые скорости всех звеньев дерева еще за два прохода. Это особенно уменьшает быстродействие метода в случаях, когда происходит длительный контакт в нескольких точках, который в методе импульсов представляет собой множество микроударов.

— Громоздкая реализация Алгоритм достаточно сложен в реализации и требует разработки специальной пространственной алгебры.

Метод LCP:

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

Преимущества:

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

— Хорошо совместим с соударением тел Контактная точка при соударениях и взаимодействии не связанных тел ничем принципиально не отличается от сустава, только имеет 5 степеней свободы. Это позволяет решать задачу столкновения систем связанных тел.

Недостатки:

— Высокая вычислительная сложность Каждая степень свободы системы связанных тел увеличивает размерность матрицы, которая раскладывается на множители. Например, каждый шарнирный сустав добавляет 3 степени свободы, а контактная точка — 5 степеней. В итоге при увеличении группы связанных тел матрица растет очень быстро, и большая часть времени уходит только на заполнение этой матрицы и на разложение на множители.

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

O (m1)+ O (m23)+ O (n),

где m1 — кол-во сочленений в группе, m2 — общее количество степеней свободы, n — количество тел в группе.

— Нет гарантии, что решение будет найдено Для борьбы с этим недостатком применяют разные хитрости, смягчающие ограничения задачи, в итоге суставы тел могут расходиться с правильных мест и для ликвидации этого приходится дополнительно прикладывать силы, стягивающие суставы на место.

— Невысокая точность Точность обычно сильно падает не только из-за борьбы с предыдущим недостатком, но и из-за попыток уменьшить вычислительную сложность различными способами — например, суставы разбиваются на подгруппы и решение для них ищется отдельно.

— Громоздкая реализация Со всеми средствами, обеспечивающими стабильность и высокую скорость алгоритма, его реализация получается достаточно громоздкой.

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

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

LCP метод слишком медленно работает при наличии силы трения и большом количестве контактов, и в этом случае высока вероятность вырождения задачи.

Интересное решение этой проблемы было предложено специалистами из Стэндфордского Университета. Это эвристический алгоритм, но он дает очень правдоподобные результаты. Идея заключается в том, чтобы на каждом шаге по времени, если несколько тел пересеклись, искать такое геометрическое положение тел, при котором ни одно из них не будет пересекаться с другими.

Рассматривается каждая пересекающаяся пара тел и вычисляется минимальный поворот и сдвиг, необходимый для того, чтобы ликвидировать пересечение. Но после этих операций перемещенные тела могут пересечься с другими телами, и тогда процесс повторяется. Новые позиции и повороты тел, полученные после этого итерационного процесса, сравниваются с позициями и поворотами на предыдущем шаге по времени — исходя из этого, вычисляются новые линейные и угловые скорости тел. В итоге мы неявным образом приложили к телам силы, которые не дали им пересечься на текущем шаге по времени.

Неявно выраженная скорость и силы в совокупности с итерационным процессом удивительным образом приводят к тому, что одни тела давят на другие, их давление передается третьим и далее по сети контактов. В результате даже большие кучи тел ведут себя очень естественно и стабильно.

Однако алгоритм достаточно неоднозначен — многое зависит от порядка обработки пар. Чтобы не появлялось заметной ассиметричности в поведении системы, порядок обработки пар выбирается случайно. Подобная неоднозначность результата, несомненно, является основным недостатком метода.

Также не было предложено, как расширить этот метод на задачу обработки взаимодействия систем связанных тел.

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

1.3 Моделирование деформируемых тел Динамика деформируемого тела не может быть описана одной материальной точкой, как это делалось при моделировании твердых тел. Поэтому деформируемые тела состоят из множества точек, обладающих собственными массами и скоростями. При этом геометрия тела должна быть привязана к этим точкам, которые являются узлами сетки.

К деформируемым телам относятся такие объекты, как ткани (представляющие собой плоскую сетку), веревки (линейная цепь), а так же более сложные структуры — например растения или резиновый мячик, составленный из множества треугольников.

Метод пружин Материальные точки, из которых составлены деформируемые тела, соединяются пружинами. На каждом шаге времени к точкам, связанным пружиной, прикладываются силы равной величины и направленные в разные стороны. Величина силы — функция от разницы текущего расстояния между точками и длины пружины в расслабленном состоянии, а также иногда от относительной скорости сближения точек. Выбор функции зависит от того, какой результат требуется. Это может быть и обычная модель пружины, сила действия которой прямо-пропорциональна растяжению (сжатию) с коэффициентом жесткости k, но обычно нужны пружины, быстро гасящие колебания, и при отсутствии внешнего воздействия, быстро возвращающие тело в расслабленное состояние.

Если этим методом моделируются объемные тела (а не ткани или веревки), то структуру материальных точек и связывающих их пружин надо выбирать таким образом, чтобы эта конструкция была устойчивой и сохраняла форму. При неудачной структуре, тела могут складываться или выворачиваться, что обычно нежелательно.

Особо отметим, что основной проблемой является взаимодействие деформируемых и твердых тел — на стыке двух физических моделей возникают дополнительные трудности. Обычно их решают грубо — деформируемые тела, ткани и веревки считаются невесомыми по отношению к твердым телам и никак на них не влияют. В итоге дизайнеры игровых сцен ограничены в средствах и не могут объединить в одну более-менее сложную конструкцию одновременно и деформируемые и твердые тела. Когда же, например, твердое тело надо подвесить на веревках, то веревки моделируются как система связанных твердых тел, а это не очень экономно, так как рассмотренные выше методы моделирования систем связанных тел медленно работают с длинными цепями тел.

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

1) динамика твердых тел

2) динамика систем связанных тел

3) динамика деформируемых тел.

Основные сложности возникают на стыках этих методов — причем часто разработчики вынуждены грубо обходить многие проблемы. Например, из-за медленного расчета столкновений систем связанных тел, такие столкновения полностью исключаются и системы тел сталкиваются только со статической геометрией сцены. Про невесомость тканей и веревок упоминалось выше.

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

2. Постановка задачи Существующие системы расчета динамики систем тел в реальном времени используют описанные выше методы и обладают множеством недостатков. Свободно распространяемые системы, такие как Newton, Tokamak и Open Dynamics Engine (ODE), обладают слишком низким качеством моделирования и множеством ошибок, причем открытые исходные коды есть только у ODE, однако ODE основан на методе LCP и имеет принципиальные ограничения с силой трения и это делает его неприемлемым для целого ряда задач.

Коммерческие системы, такие как Havok, Meqon, Novodex тоже не избавлены от недостатков, но большинство проблем тех методов, на которых они основаны, обходятся разными хитростями и специальными надстройками. Наиболее качественной системой является Havok, но лицензия на его использование дорого стоит. Его исходные тексты закрыты, и поэтому не известно, как он устроен и насколько «честно» моделируется динамика.

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

Задачей данной работы является разработка собственной системы расчета динамики для использования в нескольких проектах. За основу решено было взять нестандартный подход к моделированию динамики твердых тел, недавно предложенный Томасом Якобсеном.

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

Исследовать метод Якобсена и реализовать его для расчета динамики систем связанных тел.

Модифицировать метод Якобсена для обработки механических столкновений и длительных взаимодействий систем связанных тел.

Выбрать и реализовать метод обнаружения столкновений Разработать оптимизацию обнаружения столкновений, учитывая специфику метода Якобсена.

Провести аналитический и экспериментальный анализ полученного метода, оценить его производительность и точность.

Основные требования к методу:

Работа в реальном времени с частотой не менее 30 Гц Выполнение законов механики Ньютона.

Допустимы ограничения на высокие линейные и угловые скорости.

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

3. Построение математической модели

3.1 Метод Якобсена Для начала опишем метод Якобсена в его изначальном виде, как он был предложен автором.

Идея метода достаточно необычна и предполагает модель твердого тела, отличающуюся от рассмотренной выше. Ее преимущества не очевидны для случая отдельных твердых тел, но в моделировании систем связанных твердых тел преимущества значительные.

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

Этот факт не удивителен, так как попарное взаимодействие частиц в этом методе удовлетворяет всем трем законам Ньютона, а значит и система, составленная из таких взаимодействующих частиц, будет вести себя в соответствии законам Ньютона. В реальности твердые тела состоят из молекул и силы взаимодействия любых двух молекул направлены вдоль одной линии, равны по величине и противоположны по направлению. И из этой структуры возникают все макро-понятия твердого тела, такие как моменты инерции, линейные и угловые скорости движения.

Основная задача в данном случае добиться большой жесткости пружин, чтобы они не растягивались — тогда кристаллическую решетку твердого тела можно упростить вплоть до тетраэдра с четырьмя частицами в вершинах — и эта фигура сохранит динамические свойства твердого тела.

Томас Якобсен предложил в модели твердого тела материальную точку с моментами инерции заменить на тетраэдр с материальными точками в вершинах, каждая из которых имеет массу и линейную скорость. Тетраэдр выбран как самая простая из объемных фигур, имеющих только 6 степеней свободы, как и требуется для твердого тела, однако твердое тело может быть представлено и любой другой фигурой с таким же количеством степеней свободы — более сложные фигуры понадобятся при моделировании суставов.

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

Интегрирование Верле Это метод численного интегрирования уравнений движения, часто использующийся в задачах молекулярной динамики и относится к методам конечных разностей. Существует три его разновидности — позиционный, чехарда (leapfrog) и скоростной. Нас интересует позиционный вариант, его запись выглядит так:

(1)

если шаг по времени постоянный, то выражение примет следующий вид:

(2)

Преимущества этой схемы в довольно высокой стабильности, особенно при наложении граничных условий. В ней нет явно выраженной скорости, и это основная причина, по которой она была выбрана. Кроме того, она быстро вычисляется, примерно так же как интегрирование по Эйлеру.

Рассмотрим простой случай двух частиц, которые мы хотим связать жесткой пружиной длиной L. Допустим в текущий момент частицы находятся на правильном расстоянии L друг от друга и координаты частиц на предыдущем шаге по времени не совпадают с текущими — т. е. у частиц есть какая-то скорость (Рис. 2). Также на частицы действует сила гравитации.

Рис. 2. Шаг интегрирования Верле на примере 2-х частиц со связью

При интегрировании Верле происходит следующее:

Скорость, которой частица обладала на предыдущем шаге и которая может быть записана через разницу текущих и предыдущих координат (v = (xi — xi-1) /i -1) изменяется на величину ai, где a — это ускорение от силы действовавшей на протяжении текущего шага i. Новая позиция получается сдвигом частицы на расстояние, пройденное частицей за время i, т. е.

Таким образом, мы получили выражение (1).

После того как мы вычислили новые позиции частиц летящих пока независимо, мы хотим добиться того, чтобы частицы на этом шаге не разлетелись, а остались на расстоянии L друг от друга. Удовлетворить этому ограничению мы можем простым способом — сдвинем полученные позиции частиц относительно их общего центра масс таким образом, чтобы расстояние между ними вновь было равно L. Этим действием мы неявно изменили скорость частиц на текущем шаге, которая вычисляется из разницы текущих координат и предыдущих, причем изменение этой скорости равносильно неявному применению силы в точности такой по величине и направлению, чтобы удержать между частицами необходимое расстояние L — это и есть сила нашей бесконечно жесткой пружины. С этого момента пружина перестает быть пружиной по своей сути, поэтому будем называть ее связью.

Если мы продолжим этот процесс, то увидим, что эта система из 2-х частиц движется подобно твердому телу — центр масс системы будет двигаться по баллистической кривой, а частицы будут вращаться вокруг центра масс с постоянной скоростью. На самом деле из-за ошибки метода скорость вращения постепенно затухает, и чем больше частицы поворачиваются за каждый шаг по времени, тем затухание сильнее. Однако при достаточно малом шаге по времени или медленном вращении скорость затухания достаточно мала. О том, откуда берется ошибка метода и как бороться с затуханием будет описано ниже.

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

Теперь рассмотрим более сложный случай с 3-мя частицами и 2-мя связями (Рис. 3). Тут проблема заключается в том, что частица 2 испытывает на себе влияние сразу двух связей. Каким-то образом мы должны удовлетворить одновременно двум ограничениям длин.

В методе Якобсена это решалось итерационным процессом, на каждом шаге которого длины всех связей поочередно выравнивались. В итоге процесс сходился к некоторому результату, удовлетворяющему сразу всем длинам. Недостатком этого способа является зависимость результата от порядка обхода связей, в результате чего в поведении системы могут появиться нежелательные закономерности. Чтобы избежать этого, можно каждый раз случайным образом менять порядок обхода, но это тоже не является хорошим решением, так как в этом случае даже при одинаковых начальных условиях результат моделирования может различаться.

Для решения этой проблемы был разработан специальный алгоритм с расщеплением частиц.

3.2 Модификация метода Якобсена

3.2.1 Алгоритм с расщеплением частиц Для предыдущего примера проведем итерационный процесс, на каждом шаге которого выполняются следующие действия:

1) Расщепляем частицу 2 на две частицы той же массы и выравниваем длины обеих связей независимо. Поскольку сдвиг позиций каждой пары выполняется относительно общего центра масс пары, то общий центр масс системы остается на месте.

2) Частица 2 теперь находится сразу в двух разных местах с координатами x1 и x2. Назначим ей новую позицию x = (x1 + x2) / 2. Как видно, положение общего центра масс снова не меняется.

3) Поскольку на шаге 2 длины связей снова нарушились, переходим к шагу 1 и повторяем процесс.

Рис. 3. Итерационный процесс с расщеплением частицы 2

Этот итерационный процесс быстро сходится к требуемому результату, с заданной точностью удовлетворяющему ограничениям длин. Центр масс этой системы тоже движется по баллистической кривой, но сама система из 3-х частиц теперь имеет 4 степени свободы (в 2х-мерном пространстве), несмотря на это, если мы посчитаем общий момент вращения системы, то заметим, что он тоже сохраняется постоянным, если не учитывать постепенное затухание из-за ошибки метода.

Если мы добавим еще одну связь между частицами 1 и 3, то получим треугольник, который уже имеет всего 3 степени свободы, т. е. он ведет себя, как твердое тело в 2-х мерном пространстве, сохраняя свою форму.

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

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

3.2.2 Представление суставов Если у двух фигур, например тетраэдров, выберем общую частицу, то получим полную аналогию сустава «шарнир», который дает телам 3 степени свободы движения относительно друг друга. (Рис. 4. а) Рис. 4. Представление суставов совокупностью частиц и связей

Две общие точки ограничивают степень свободы до одной — это соединения типа «петля» (Рис. 4. б), тела могут поворачиваться только вдоль одной оси. «Двойная петля» (Рис. 4. в) — сустав с 2-мя степенями свободы. Мы пока не будем затрагивать сустав типа «поршень», так как для этого необходимо расширение модели связей, которое будет сделано ниже.

Таким образом, мы можем создавать сложные системы связанных тел, при этом выбор структуры ничем не ограничен, можно образовывать циклы и никакое сочетание суставов не приведет к конфликту, которые бывают в классических методах. При этом расчет динамики такой системы значительно более быстрый, чем в методе ограничений и LCP методе — время вычисления линейно зависит от количества связей в системе (при фиксированном кол-ве итераций выравнивания).

Для того чтобы ограничить сгиб каждого сустава, нам нужно добавить еще один тип связи. Сейчас мы использовали только связь, восстанавливающую свою длину, как в случае ее превышения, так и в случае уменьшения — т. е. эта связь выполняет ограничение-равенство. Также можно ввести связи с ограничениями-неравенствами. Например, связь типа «не меньше чем» восстанавливает свою длину L, только в случае, если ее текущая длина стала меньше L, но никак не реагирует, если длина превышается. Аналогично для связей «не больше чем»

Рис. 5. Пример связи ограничения-неравенства, ограничивающей максимальный угол сгиба шарнирного сустава

Связывая определенные частицы суставов этими связями, мы можем добиться угловых ограничений для любых типов суставов.

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

Для связи с произвольной жесткостью новая длина L в итерационном процессе выбирается по следующему закону:

где

L — текущая длина связи

L — новая длина связи

Lrelax — требуемая длина связи

N — кол-во итераций

q — жесткость связи, q [ 0, 1 ]

Важно, что в отличие от классического коэффициента жесткости пружины, область значений которого находится в пределах от нуля до бесконечности, в нашей модели коэффициент жесткости q меняется от нуля до единицы и, когда он равен единице — это равносильно бесконечной жесткости пружины.

3.2.3 Связь с классической моделью От совокупности частиц, представляющих твердое тело, можно перейти к классическому представлению твердого тела. Для этого достаточно вычислить координаты его центра масс и получить матрицу ориентации, взяв в качестве базиса два любых вектора — например v1 = x2 — x1 и v2 = x3 — x1. Это стандартная функция в большинстве библиотек для работы с векторами и матрицами.

Сложнее поставить в соответствие конкретному телу с произвольными моментами инерции, несколько частиц, расположенных таким образом, что моменты инерции полученной системы соответствуют требуемым.

Допустим у нас тело массы M и моментами инерции Ix, Iy, Iz. Возьмем начальную конфигурацию тела в форме тетраэдра, либо более сложной конструкции (если есть суставы) с частицами массой m = M/n, где n — кол-во частиц. Вычислим координаты центра масс. Из этих частиц выделим те, позиция которых относительно центра масс не является принципиальной (например, частицы не относятся к суставам и с ними не связаны ограничители). Таких частиц должно быть не меньше 4-х.

Теперь проведем следующий итерационный процесс:

1) Вычислим момент инерции относительно оси x.

где ryz — расстояние от i-й частицы до центра масс в плоскости yz.

2) Назначим частицам новую позицию xi = xc + (xc — xi)(Ix /Ix)2

3) Повторим пункты 1 и 2 для осей y и z.

4) Перейдем к пункту 1 с осью x и повторим процесс.

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

3.2.4 Механические соударения Способ обработки столкновений, предложенный Якобсеном, позволял моделировать только неупругие удары, хотя и неплохо справлялся с обработкой длительных контактов. Неупругость взаимодействия до сих пор считалась принципиальным ограничением метода. Как выяснилось, так оно и есть — в том виде, в каком алгоритм был предложен автором, упругость ударов не возможна без специальных алгоритмов обработки столкновений, поэтому зачастую используют подобие метода штрафов (основанном на приложении сил), а он дает низкую точность результата.

Теперь покажем, как можно расширить метод до обработки столкновений произвольной упругости.

Сначала построим упрощенную модель столкновений шаров без трения:

Допустим, у нас есть две частицы и между ними связь типа «не меньше чем», которая не допускает сближение частиц ближе, чем на расстояние L. Пока частицы находятся на большем расстоянии, то они никак не взаимодействуют. Но стоит им сблизиться на расстояние меньшее L, то связь восстановит это расстояние между ними. Промоделировав этот процесс, можно заметить, что в этом случае частицы ведут себя как при неупругом столкновении, как если бы они представляли собой шарики, причем их радиусы R1 + R2 = L.

Во-первых, восстановление длины происходит относительно общего центра масс частиц, таким образом, из-за использования интегрирования Верле скорость центра масс не изменяется, а значит, выполняется закон сохранения импульса.

Во-вторых, восстановление длины связи, при интегрировании Верле аннулирует относительную скорость сближения, т. е. относительная скорость разлета равна нулю, а это в совокупности с законом сохранения импульса является полным аналогом неупругого удара.

Рис. 6. Связь ограничение-неравенство между двумя шариками препятствует их взаимопроникновению, что в итоге приводит к модели неупругого удара Легко убедиться, что в нашей модели при центральном ударе выполняется известный закон для неупругого удара где V1 и V2 — исходные скорости шаров, V1 и V2 — новые скорости, а m1, m2 — массы шаров.

Чтобы получить модель абсолютно упругого удара, нам нужно добиться того, чтобы относительная скорость разлета шаров после удара по модулю была равна относительной скорости сближения до удара. Чтобы получить такой результат, запишем L = L + 2(Lp — L) для новой длины связи и получим выражение:

(3)

где L — текущая длина связи, Lp — длина связи на предыдущем шаге, L — новая длина связи.

Таким образом, мы получили модель абсолютно упругого удара. Чтобы получить модель удара произвольной упругости с коэффициентом упругости k, запишем следующее выражение для длины связи:

(4)

где Lrelax — требуемая длина связи.

Заметим, что если в выражении (4) подставить k = 1, то оно принимает вид (3), а при k = 0 превращается в L = Lrelax, т. е. в обычное восстановление длины связи.

Итак, мы получили модель механических соударений шариков произвольной упругости, в которой, как показали эксперименты, выполняются все основные свойства соударений.

Однако эксперименты с шариками соединенными в пары обычными связями по методу Якобсена показали, что даже при коэффициенте упругости удара k = 1 пары при столкновении часто ведут себя как неупругие, особенно когда связи лежат на линии, вдоль которой произошло столкновение. Объясняется это тем, что связь «ограничение-равенство» в методе Якобсена, которая просто восстанавливает длину связи, является моделью неупругой связи. В этом и кроется основная проблема метода — такие связи поглощают энергию удара при передаче частицам, которые они связывают и в итоге превращают удар в неупругий.

Возможно, именно в этом кроется причина неудач попыток применения метода Якобсена в чистом виде для моделирования упругих ударов.

Если же в методе Якобсена заменить способ восстановления длины с простого L = Lrelax на выражение (4), то связи можно сделать упругими, и пример с парами шариков в этом случае работает корректно.

Но эксперименты с шариками, соединенными в треугольники (по 3 шарика), вновь показали поглощение энергии удара даже при k = 1. В этот раз причиной стало расщепление частиц, введенное нами при модификации метода Якобсена. Энергия теряется при стягивании расщепленной частицы в точку:

где xi координаты i-й копии расщепленной частицы, N — кол-во копий частицы (равно количеству входящих в частицу связей) Тут по аналогии со связью нулевой длины, можно заменить простое перемещение копий частицы в точку, отражением их относительно этой точки на изначальное расстояние, умноженное на коэффициент упругости k.

В этом случае эксперименты с треугольниками из шариков дают правильные результаты, энергия при абсолютно упругих ударах не поглощается. Теперь можно соединять шарики и в более сложные структуры с множеством степеней свободы, и столкновения таких систем тоже происходит правильно, причем импульсы таких систем в контактных точках сложно посчитать аналитически (а это приходится делать в классических методах), а в нашем случае импульс учитывается автоматически в течение итерационного процесса.

Таким образом, мы построили новую модель связей, которые теперь обладают упругостью и при абсолютно упругом ударе не поглощают энергию. Это позволит нам расширить метод, чтобы с его помощью, не прибегая к другим методам, моделировать механические соударения и длительные контакты.

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

3.2.5 Проверка пересечений Задачей блока проверки пересечений является регистрация самого факта пересечения примитивов и определение точек наибольшего проникновения. Для решения этой задачи существуют уже готовые библиотеки, но в нашем случае выбор метода очень важен.

Основная проблема нашего метода обработки соударений и контактов заключается в том, что проверка на столкновение происходит несколько раз за один шаг по времени. Нагрузка на блок проверки пересечений возрастает пропорционально количеству итераций. Но можно учесть, что в рамках одного шага по времени объекты движутся мало, на каждой итерации фактически лишь уточняется положение тел в пространстве. С учетом этой особенности, был разработан специальный алгоритм оптимизации проверки столкновений.

Рис. 7. Ограничивающие объемы для оптимизации проверки столкновений На Рис. 7 изображено тело, форма которого задана неким примитивом (или их набором). Вокруг этого тела описывается сфера наименьшего возможного радиуса. Тело вместе со сферой заключено в прямоугольный параллелепипед, стороны которого ориентированы по осям мировых координат. Назовем этот параллелепипед ограничивающим ящиком, а сферу — ограничивающей сферой.

Предварительное выделение кандидатов на столкновение В отличие от ограничивающей сферы, ограничивающий ящик не движется вместе с телом. Он находится на месте до тех пор, пока сфера не выйдет за его пределы. Тогда для ящика рассчитывается новое положение и размеры с учетом величины и направления скорости движения тела (Рис. 8 справа). При этом гарантируется, что тело не вылетит из ящика в течение времени T, если продолжит лететь свободно и не столкнется с другими телами. В момент расчета нового положения ящика проверяется его пересечение с ограничивающими ящиками других тел. Те тела, с чьими ограничивающими ящиками было обнаружено пересечение, образуют пару с рассматриваемым телом и такие пары заносятся в список кандидатов на столкновение.

Рис. 8. Использование ограничивающих ящиков для выделения кандидатов на столкновение. Справа показано, что размеры и смещение ящика выбираются исходя из скорости объекта, и ящик остается неподвижным, пока объект не выйдет за его пределы Если тело по тем или иным причинам уменьшило свою скорость на значение больше некоторой величины, не достигнув края ящика, то строится вложенный ящик меньших размеров — при этом проводить проверку на пересечение имеет смысл только с ящиками тех тел, которые были в списке кандидатов на столкновение с этим телом, т. е. их ящики пересекались с ящиком этого тела до момента расчета нового размера.

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

Таким образом, вместо проверки пересечений между всеми парами объектов (N2 проверок) мы получаем метод, выделяющий пары, выполняя в среднем лишь N проверок за один шаг по времени.

Причем если объекты занимают достаточно устойчивую позицию или колеблются в рамках ограниченной области, то для них ограничивающие ящики не пересчитываются, а пары взаимодействия остаются постоянными.

Например, на Рис. 8. пирамидка из шариков достаточно стабилизировалась, и их ящики не пересчитываются, пунктирными линиями показываются установившиеся пары взаимодействия, и только они и проверяются на столкновения.

Но самое важное в этом методе, что внутри итерационного процесса, когда уточняются положения тел, пересчет ящиков практически никогда не происходит и поэтому проверка выполняется только для выбранных пар тел.

Точное определение столкновений На каждом шаге по времени у кандидатов на столкновение проверяется пересечение их ограничивающих сфер, сравнением расстояния между центрами сфер и суммой их радиусов.

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

Но основная проблема возникает, когда нужно проверить пересечения произвольных многогранников. Обычно это связано с большей вычислительной сложностью, чем в случае с обычными примитивами. Как правило, эти многогранники должны быть выпуклыми, чтобы не могло быть несколько глубочайших точек при взаимопроникновении.

Существует множество методов проверки пересечений, но в нашем случае нас интересуют именно такие, оптимизация которых учитывает, что изменения положений тел могут быть малы. С учетом этого был выбран метод Лина-Кенни.

Метод Лина-Кенни (Lin-Canny):

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

Идея в том, чтобы для данной пары многогранников проверять сначала те элементы, который при прошлой проверке были ближайшими. И если они теперь не являются ближайшими, то проверять их соседние элементы.

Проверка того факта, что рассматриваемые элементы ближайшие, основана на применении областей Вороного.

На Рис. 9 для наглядности показан плоский пример, в котором области Вороного делятся на два типа — прямоугольные для ребер и треугольные для вершин.