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

Рис. 9. Пример областей Вороного для плоского случая. Точки a и b ближайшие, если aV(Fb) и bV(Fa).

Установить, что точки a и b ближайшие, не проверяя остальных, можно в том случае, если выполняется следующее условие:

Точка a принадлежит области, образованной элементом, которому принадлежит точка b (в данном случае вершина) и точка b принадлежит области, образованной элементом, которому принадлежит точка a (в данном случае ребро). Т.е. aV(Fb) и bV(Fa).

В противном случае точки a и b не являются ближайшими и надо проверять соседние элементы.

В трехмерном пространстве областей Вороного будет три типа - для граней, для ребер и для вершин. Соответственно проверяться будут пары вершина-грань, ребро-ребро. В остальном все аналогично.

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

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

3.2.6 Обработка взаимодействия

Воздействие на произвольную точку

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

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

Рассмотрим пример с двумя частицами, соединенными связью.

Рис. 10. Перемещение произвольной точки p, лежащей на связи между частицами x1, x2 в произвольную точку q сдвигом частиц x1 и x2

Чтобы переместить частицу x1 в новую позицию x1, нам достаточно переместить только саму частицу, а x2 оставить на месте. После удовлетворения ограничению-равенству связи, обе частицы сдвинутся, но эффект будет таким, как если бы на частицу x1 действовала сила в том направлении, в котором мы ее сдвинули.

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

причем с1 + с2 = 1. Эти коэффициенты показывают степень близости точи p к крайним частицам. Например при с1=1, с2=0 точка p совпадет с частицей x1, а при с1 = с2 = 0.5 будет лежать ровно посередине.

Нам нужно переместить точку p в точку q. Для этого мы будем перемещать точки x1 и x2 вдоль направления

Сдвигаем частицы x1 и x2 пропорционально коэффициентам на некоторое расстояние вдоль направления ?.

(5)

где - некоторое неизвестное значение.

Но нам известно, что новое положение p после перемещения частиц будет

Вспомним, что p = q. Значит, мы должны выбрать значение таким, чтобы p точно совпало с q. Так как мы перемещаем обе частицы только в направлении ?, то и точка p переместится в направлении ?. Следовательно, решение уравнения p = q может быть найдено решением:

(6)

Распишем левую сторону:

Это вместе с правой частью (6) дает:

Подставляя в (5), получим новые позиции частиц, для которых p совпадает с q.

После удовлетворения ограничению-равенству связи, частицы и точка p сместятся, но эффект будет именно таким, как если бы к точке p была приложена сила в направлении ?.

Заметим, что если с1=1, с2=0, (т.е. точка p совпадает с частицей x1), то сдвигается только частица x1, а x2 остается на месте, т.е. приложение сил непосредственно к частице сдвигом самой частицы является лишь частным случаем этого способа. Важно, что этот способ можно использовать для объемной конструкции (например, тетраэдра), причем точка p может быть расположена произвольно относительно частиц.

Ее положение записывается как линейная комбинация частиц конструкции:

причем

Тогда для каждой частицы новое положение:

где

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

Сила трения

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

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

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

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

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

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

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

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

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

(7)

где - коэффициент трения.

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

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

Этот угол показан на Рис. 11.в.

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

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

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

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

Рис. 12

Эти кубы имеют 2 контакта «вершина-грань» и 2 контакта «ребро-ребро».

При контакте грань-грань, на самом деле рассматриваются несколько контактов типа вершина-грань или ребро-ребро. Например, куб лежащий на поверхности имеет 4 контактных точки в вершинах, а один куб, лежащий на другом со смещением (Рис. 12) имеет с ним две контактных точки типа ребро-ребро и еще две типа вершина-грань. Аналогично контакт ребро-грань обычно состоит из 2-х контактов.

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

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

3.2.7 Деформируемые тела

Из частиц и связей можно собирать конструкции с множеством степеней свободы. Например, веревки (Рис. 13.1) или ткани (Рис. 13.2). Также можно создавать объемные конструкции без лишних степеней свободы, но со связями небольшой жесткости и множеством частиц. (Рис. 13.3). При этом тело может состоять из множества тетраэдров и содержать множество вершин внутри себя. Но для оптимизации в большинстве случаев достаточно бывает соединить все частицы поверхности с одной внутренней «центральной» частицей.

Рис. 13. Примеры деформируемых тел, составленных из частиц и связей.

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

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

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

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

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

Однако в нашем случае мы можем давать хорошее начальное приближение корня, которое позволит вычислить его быстрее при использовании метода Ньютона-Рафсона. Если нам известен квадрат текущего расстояния L2 между частицами, то, взяв некоторое начальное приближение L0, мы можем вычислить значение L с заданной точностью, используя следующую рекуррентную запись:

Процесс останавливается, если | Ln+1 - Ln | < ?

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

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

4. Исследование модели

4.1 Точность интегрирования Верле

Чтобы оценить порядок точности схемы Верле для траектории x(t) разложим приращения для момента времени t0 в ряд Тейлора:

Сумма этих выражений дает:

Мы получили выражение вида (2), т.е. схему Верле для постоянного шага - следовательно, эта схема имеет четвертый порядок точности для траектории при постоянном шаге по времени.

Если мы вместо суммы запишем разницу первых двух выражений, то получим выражение для скорости:

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

Найдем теперь глобальную ошибку схемы. Так как:

и

следовательно:

аналогично:

По индукции можно показать, что общая запись:

Если мы рассматриваем глобальную ошибку между позициями x(t0) и x(t0 + T), где T = n?T, то очевидно, что:



Следовательно, глобальная (совокупная) ошибка над интервалом времени T , разбитом на шаги величиной ?t:

Поскольку скорость определяется из отдельных позиций в схеме Верле, то ее глобальная ошибка также O(?t2).

4.2 Сохранение механической энергии системы

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

Рассмотрим систему из двух частиц, соединенных связью длины L, причем одна из частиц неподвижно закреплена в пространстве.

Рис. 14. Потеря скорости (а значит и кинетической энергии) при вращении одной частицы вокруг закрепленной неподвижной частицы.

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

Как видно из Рис. 14, до того как применяется ограничение связи, частица сдвигается из точки 2 в направлении, являющемся продолжением отрезка, соединяющей точки 1 и 2 на длину этого отрезка. После применения ограничения, расстояние между точками 2 и 3 получается меньше, чем между точками 1 и 2. В схеме интегрирования Верле это означает, что часть скорости потеряна. Величина этой потери зависит от угла между отрезком [1, 2] и касательной к окружности в точке 2. При малых углах поворота этот отрезок почти совпадает с касательной, поэтому потеря скорости становится незначительной.

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

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

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

Где E0 - начальное значение полной механической энергии, E - ее текущее значение. ЕК - текущее значение кинетической энергии.

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

где x - текущее положение частицы, xp - положение на предыдущем шаге, xp - скорректированное положение частицы на предыдущем шаге.

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

4.3 Сходимость

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

Проблема большого диаметра:

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

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

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

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

Эта проблема решается добавлением связей ограничений-неравенств типа «не больше чем» (Рис. 15). При этом L14 = L12 + L23 + L34; L13 = L12 + L23; L24 = L23 + L34

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

Рис. 15. Улучшение сходимости цепи на случай растяжения, созданием связей ограничений-неравенств типа «не больше чем» (показаны пунктирными линиями).

Проблема слабого звена:

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

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

Рис. 16

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

4.4 Зависимость свойств от параметров

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

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

Niter - кол-во итераций в алгоритме выравнивания связей

Nsqrt - кол-во итераций при вычислении корня

?t - величина шага по времени

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

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

Для анализа зависимости жесткости от параметров, проведем серию экспериментов с цепью из 15-и звеньев, один конец которой закреплен, а к другому концу прикладывается сила.

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

Рис. 17. Зависимость относительного удлинения цепи от количества итераций.

На Рис. 17 показан график зависимости относительного удлинения цепи L2/L, где L2 - установившаяся длина цепи, а L - ее реальная длина, от количества итераций Niter при сохранении остальных параметров постоянными. Эта зависимость соответствует функции f (N) = c/N + 1, т.е. является гиперболической.

На рисунке показан график зависимости удлинения цепи от величины шага по времени, из которого видно, что зависимость является функцией f (?t) = c?t2 + 1.

Рис. 18. Зависимость относительного удлинения цепи от величины шага по времени.

Рис. 19. Зависимость относительного удлинения цепи от коэффициента упругости.

Коэффициент упругости k ухудшает сходимость, и поэтому, при фиксированном количестве итераций, уменьшает жесткость цепи (Рис. 19). Функция зависимости f (k) = c/(1 - k) + 1. Это гипербола и при k = 1 она уходит в бесконечность. Это значит, что при абсолютной упругости связей длина цепи под воздействием силы не может окончательно установиться на конкретном значении независимо от числа итераций. Поэтому для сходимости коэффициент упругости всегда должен быть меньше единицы. Однако во время расчета ударов несколько первых итераций можно сделать с k = 1 для того, что бы импульс распространился по связям, но после этого необходимо еще несколько итераций с меньшим либо вообще нулевым k, чтобы обеспечить жесткость связей.

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

Рис. 20 Зависимость относительного удлинения цепи от количества шагов приближения корня методом Ньютона-Рафсона.

Исследование зависимости относительного удлинения цепи от количества звеньев дало неожиданную зависимость - скорость роста относительного удлинения уменьшается с ростом количества звеньев.

Рис. 21. Зависимость относительного удлинения цепи от количества звеньев.

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

Когда звено одно - оно всегда абсолютно жесткое и этим отличается от общего случая нескольких звеньев. Этим объясняется особенность поведения графика при N = 1. В нашем эксперименте случай был особым и при 2-х звеньях, т.к. оба являются крайними - одно закреплено, а ко второму прикладывалась сила. Удлинение этих звеньев отличается от удлинения средних, поэтому график при N = 2 отклоняется от линейного.

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

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

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

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

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

Например, если количество итераций Niter = 1 и при приближении квадратного корня Nsqrt = 1, то мы получаем самую быструю с вычислительной точки зрения систему. При этом жесткость связей такой системы вполне достаточная, чтобы моделировать, например, ткань или какое-нибудь эластичное деформируемое тело. Поскольку визуальное представление ткани треугольниками обычно совпадает с ее физическим представлением сеткой из частиц и связей, то требуется достаточно большое количество частиц и связей. Но при единственной итерации за шаг по времени и замене корня лишь первым шагом приближения сводит количество операций до 4N делений за один шаг по времени, где N - количество связей - это настолько быстро, что размер такой конструкции ограничен в основном тем, сколько треугольников может нарисовать визуализирующая часть программы.

4.5 Тесты на производительность

Все тесты выполнялись на машине Pentium 4 с тактовой частотой 1600 МГц. Производительность измерялась без визуализации объектов. Таким образом, для расчета динамики было высвобождено 90% времени центрального процессора.

Тест 1:

В статическую квадратную коробку 500х500х500 м. бросаются отдельные кубики 20х20х20 м. Каждый кубик содержит 4 частицы и 6 связей (тетраэдр). Количество итераций Niter = 15, количество приближений корня Nsqrt = 1.

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

Рис. 23. Графики зависимости среднего времени на расчет одного шага (сек.) от количества кубиков в сцене. Серый график - кубики связаны в цепь, черный - отдельные кубики с проверкой столкновений.

Из графика видно, что с частотой 30 раз в секунду может обрабатываться сцена, содержащая до 220 кубиков.

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

Тест 2:

Кубики последовательно связаны в одну цепь шарнирными суставами. Каждый кубик содержит в себе 4 частицы и 10 связей (3 дополнительные связи для сустава и 1 связь-ограничитель). Столкновения между кубиками не проверяются. Как и в предыдущем тесте Niter = 15, Nsqrt = 1.

График теста показан на рисунке 23 серым цветом. Как видно, график линейный. Время расчета линейно зависит от количества связей и числа итераций и может быть вычислено теоретически по формуле ?t = cNlinkNiterNsqrt; При этом для случая с проверкой столкновений формула примет вид ?t = c1Niter (NlinkNsqrt + c2Ncoll), где Ncoll - количество взаимодействующих пар.

В этом тесте с частотой 30 раз в секунду может обрабатываться сцена, содержащая до 340 кубиков.

Таким образом, при 15 итерациях за один шаг по времени с частотой 30 раз в секунду может обрабатываться около 3'400 связей. При этом если делать только по 1-й итерации, то количество связей возрастает до 51'000. На самом деле это далеко не предел, если использовать не универсальный алгоритм, а специальную его модификацию без коэффициентов упругости и жесткости и без расщепления частиц, количество операций сокращается в 5 раз и количество связей можно поднять до 250'000.

Приведенные результаты не учитывают затраты времени на визуализацию сцены. В реальных игровых проектах визуализация занимает много времени центрального процессора, несмотря на то, что большинство расчетов происходит на процессоре видеокарты. Таким образом, если под расчет динамики будет выделено 10% процессорного времени, то вместо 200 объектов в реальном времени возможно рассчитывать только 20 объектов, что тоже не так уж и мало.

Однако с ростом мощности процессоров и видеокарт под физические расчеты в современных игровых проектах уже все чаще выделяется от 30% до 50% времени центрального процессора.

4.6 Общие тесты

В результате различных тестов были проверены следующие свойства и особенности модели.

Передача импульса:

На рисунке 24 показан тест - ряд из 4-х шариков лежит неподвижно и с этим рядом сталкивается шарик, движущийся со скоростью V вдоль оси ряда. Если удар абсолютно упругий, то импульс за один шаг по времени передается через весь ряд последнему шарику, который в итоге и приобретает скорость V, а остальные шарики остаются неподвижны (все шарики одинаковой массы).

Рис. 24. Импульс правого шарика до соударения (a) полностью передается левому шарику после соударения (б) если удар абсолютно упругий.

Распределение давления и импульса:

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

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

Рис. 25. Двухмерный пример установившейся структуры шариков (слева) и распространение импульса вглубь этой структуры и в стороны (справа).

Распространение волн:

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

Стабильность при нагромождениях:

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

Рис. 26. Пример небольшой стопки кубиков, в этой сцене отсутствует дрожание и дрейф объектов.

Корректная сила трения:

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

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

Качение объектов:

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

Реалистичное поведение ткани:

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

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

Выводы

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

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

1) Универсальность, единый алгоритм для расчета:

- динамики абсолютно твердых тел

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

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

- динамики деформируемых тел

2) Возможность настройки баланса между точностью и скоростью расчета

3) Высокая стабильность

4) Простота реализации

Недостатки:

1) Многократная проверка пересечений

2) Сложный этап подготовки к использованию

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

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

моделирование якобсен математический тест

Заключение

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

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

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

Выбран и реализован метод обнаружения столкновений.

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

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

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

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

Список использованной литературы

1. Baraff, David, Dynamic Simulation of Non-Penetrating Rigid Bodies, Ph.D. thesis, Dept. of Computer Science, Cornell University, 1992. http://www.cs.cmu.edu/~baraff/papers/index.html

2. Mirtich, Brian V., Impulse-base Dynamic Simulation of Rigid Body Systems, Ph.D. thesis, University of California at Berkeley, 1996. http://www.merl.com/people/mirtich/papers/thesis/thesis.html

3. Jacobsen, Thomas Advanced Character Physics, Copenhagen K Denmark http://www.teknikus.dk/tj/gdc2001.htm

4. Verlet, L. Computer experiments on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules, Phys. Rev., 159, 98-103 (1967).

5. Witkin, Andrew and David Baraff, Physically Based Modeling: Principles and Practice, Siggraph '97 course notes, 1997. http://www.cs.cmu.edu/~baraff/sigcourse/index.html

6. Кривцов А.М., Кривцова Н.В., Метод частиц и его использование в механике деформируемого твердого тела, Дальневосточный математический журнал ДВО РАН, 2002, Т.3, №2 с. 254-276.

7. Самарский А.А., Гулин А.В., Численные методы математической физики, Научный мир, 2003.

Размещено на Allbest.ru