Содержание к диссертации
Введение
1 Постановка задачи 10
1.1 Недеформируемое приближение 10
1.2 Деформируемое приближение 10
Математическая модель
2.1 Недеформируемое приближение 12
2.1.1 Вывод якобиана на примере контактного взаимодействия 18
2.1.2 Трение скольжения, качения и верчения 19
2.2 Упругопластическое приближение 22
2.2.1 Определяющие уравнения 22
2.2.2 Граничные условия 24
2.2.3 Контактные условия 24
2.2.4 Спектральные свойства системы 25
Численный метод
3.1 Недеформируемое приближение 27
3.1.1 Эффективная формулировка LCP 29
3.1.2 Решение LCP 30
3.1.3 Решение ограничений, не входящих в LCP 31
3.1.4 Временная когерентность, горячий старт 32
3.1.5 Разрешение проникновений, метод псевдоскоростей 34
3.1.6 Ускорение метода Гаусса-Зейделя решения LCP 35
3.1.7 Альтернативные алгоритмы решения LCP 36
3.2 Упругопластическое приближение 37
3.2.1 Разрывный метод Галёркина 40
3.2.2 Метод характеристик 48
4 CLASS Поиск контактов в динамических задачах CLASS 52
4.1 Недеформируемое приближение 52
4.1.1 Представление физической геометрии 52
4.1.2 Экстремальное отображение 53
4.1.3 Геометрия Минковского 54
4.1.4 Использование экстремального отображения для неявного задания геометрии 56
4.1.5 Расстояние пересечения, установление факта пересечения 58
4.1.6 Алгоритм Гильберта-Джонсона-Кёрти(GJK) 63
4.1.7 Поиск расстояния пересечения, существующие подходы 66
4.1.8 Расширение GJK на случай пересекающихся геометрий, SGJK 67
4.1.9 Аппроксимация контактной площадки 70
5 Программные комплексы 72
5.1 Программный комплекс для моделирования твёрдых недеформируемых тел 72
5.2 Программный комплекс для моделирования твёрдых деформируемых тел 79
6 Расчёты 86
6.1 Экспериментальная проверка порядка сходимости по сетке 86
6.2 Расчёт задач распада разрыва 90
6.3 Задача сейсмического мониторинга 94
6.4 Расчёт сейсмоотклика от метановой бомб 104
6.5 Расчёт каверны под железнодорожным подстилающим покрытием 105
6.6 Расчёт процесса формирования отклика в процессе дефектоскопии 106
6.7 Расчёты конечных деформаций 107
6.8 Расчёт трения качения 108
6.9 Серия верификационных тестов Destruction 109
6.10 Программный комплекс RobSim 110
6.11 Расчёт взаимодействия консолидирующихся ледоваых образований с опорой нефтедобывающего сооружения 112
6.12 Расчёты сеточно-характеристическим методом 115
Заключение 118
Список литературы
- Деформируемое приближение
- Трение скольжения, качения и верчения
- Временная когерентность, горячий старт
- Использование экстремального отображения для неявного задания геометрии
Деформируемое приближение
Несмотря на то, что локальные уравнения эластики, лежащие в основе первого программного комплекса, достаточно точно описывают очень широкий круг явлений, включая динамику систем большого количества твёрдых деформируемых тел, численные методы, пригодные для этих моделей, как правило, требуют больших вычислительных ресурсов и малопригодны моделирования систем большого количества тел в реальном времени. Обычно количество тел, взаимодействие которых адекватно описывать такими точными моделями, исчисляется единицами или десяткам. Однако, второй программный комплекс ориентирован именно на моделирование сотен и тысяч тел на одном персональном компьютере в интерактивном режиме реального времени, когда пользователь может влиять на моделируемые процессы прямо во время расчёта.
Более того, макроскопические параметры, такие, как коэффициент упругости или коэффициент трения скольжения, полученные в результате численных экспериментов, проведённых первым вычислительным комлексом, могут прямым образом использоваться во втором комплексе. Так, например, полноволновой расчёт одной секунды задачи трения качения колеса с помощью первого программного комплекса, потребовал более 5 часов реального времени. Однако, данные, полученные в нём, могут использоваться как входные значения параметров во втором программном комплексе, способном рассчитать эту же задачу за доли миллисекунды. Также несмотря на принципальную разницу подходов, лежащих в основе программных комплексов, их области применения на некоторых задачах пересекаются и при достаточном количестве входных данных их возможно верифицировать при помощи друг друга.
Существует ряд сторонних библиотек, решающих те же уравнения динамики системы твёрдых тел, что и реализованный программный комплекс Destruciton. Наиболее перспектив-4 ными среди них являются являются: Havoc, PhysX, Bullet, ODE, Box2D, Chipmunk. Однако Destruction выгодно отличает от них ряд принципиальных возможностей: Ни одна из существующих библиотек не реализует один программный интерфейс как для 2д, так и для 3д физики. Существуют подходы, где третья координата просто игнорируется, у трёхмерного тензора инерции тел зануляются z-координаты, ставится квазидвумерная постановка, однако, в таком режиме работы библиотека потребляет столько же ресурсов, как и при трёхмерном расчёте. Destruction позволяет работать как в полноценном трёхмерном режиме, так и в двумерном, потребляя ровно столько вычислительных ресурсов, сколько было бы необходимо специальным реализациям.
Среди существующих библиотек лишь Bullet поддерживает работу с геометриями, заданными экстремальным отображением. Однако, реализованный в Bullet алгоритм обнаружения столкновений работает эффективно лишь в случаях малых взаимопроникновений твёрдых тел и вносит характерного рода артефакты в определение столкновений. Разработанный в данной работе алгоритм позволяет одинаково эффективно работать при любой величине взаимопроникновения и не создаваёт артефактов.
В прикладных задачах, для которых разрабатывался Destruction, требуется большое количество экзотических нестандартных соединений(торсионные качающиеся подвески, трапециидальные демпферы и другие). В то время как их с некоторой точностью можно аппроксимировать системой тел, в Destruction реализованы специальные типы соединений, выполняющие это гораздо более эффективно без введения дополнительных тел.
Наиболее мощные из приведённых библиотек, а именно Havoc и PhysX – коммерческие разработки с закрытым исходным кодом, следовательно, их область применения сильно ограничена, и условия лицензирования не позволяют применять их в ряде проектов.
Проект Volume, реализующий метод Галёркина и метод характеристик, разрабатывался изначально с целью сравнительного анализа этих двух методов. Большинство существующих программных комплексов, занимающихся решением уравнений эластики, либо давно не поддерживаются, так как написаны на устаревших языках программирования, либо специализируются на структурированных сетках, либо являются коммерческими закрытыми решениями, либо реализуют схемы низкого порядка точности, либо не полностью поддерживают параллельность с распределённой памятью, либо решают статическую задачу. Одним словом, найти готовый проект, удовлетворявший огромному количеству выдвигаемых требований, не представлялось возможным, поэтому было решено продолжать дорабатывать Volume.
Трение скольжения, качения и верчения
В выражении (2.44), наконец, видно, какой физический смысл несёт величина А; - при безразмерной нормировке J, Л имеет смысл величины обобщённой силы которой тела обмениваются посредством связи г.
Для дальнейших выкладок удобнее проинтегрировать (2.44) по времени и работать с величинами обменных импульсов, а не сил: AT = J XAt = J Л (2.45) в полученном выражении вектор множителей лагранжа Л имеет смысл величины обменных импульсов, а вектор AT - суммарный вектор изменения импульсов тел системы, которыми они обмениваются в результате воздействия через связи.
При выборе соответствующей нормировки функций ограничения CІ, матрица J будет безразмерной, а величину ЛІ можно трактовать как импульс, которым обмениваются тела через соединение г за малый промежуток времени dt, поэтому такой импульс принято называть обменным. Чтобы найти вектор обменных импульсов, подставим (2.35), (2.34) и (2.45) в (2.29) в момент времени t + dt:
Второе уравнение полученной системы означает совместность уравнений связи и начальных условий. Если оно выполняется для t = 0, то далее его можно не рассчитывать. Первое уравнение системы перепишем в виде: JM. J X = с — JV (2.47) Окончательно оформив последнее уравнение, получаем: АХ = Ь, (2.48) где А = JM.JT, Ъ = с— JV. Последняя система уравнений играет ключевую роль в моделировании динамики системы твёрдых тел, потому что, найдя из неё вектор обменных импульсов Л, используя уравнения (2.35) и (2.45), возможно рассчитать, как изменятся за промежуток времени dt скорости всех тел системы.
До данного момента мы рассматривали лишь т.н. двунаправленные(англ. Bilateral) соединения, обменный импульс в которых может принимать как положительные, так и отрицательные значения. Многие же реально существующие связи называются однонаправленными, так как действуют только в одну сторону. Например, сила реакции опоры при взаимодействии пары тел не может быть отрицательной. Гибкая верёвка, связывающая два тела, не может их расталкивать. Для описания системы однонаправленных связей используется система, называемая системой линейных дополнений или LCP(Linear Complementarity Problem)( [29]), которую для нашего случая можно сформулировать так:
Исходя из (2.23) и (2.48), первое уравнение этой системы требует выполнения условия неуменьшения C, то есть C 0, второе условие требует, чтобы обменные импульсы не at были отрицательными и, наконец, последнее условие можно трактовать как
Физический смысл этих трёх уравнений для случая моделирования контактного взаимодействия заключается в том, что в каждой точке контакта взаимодействующие тела либо отдаляются и сила реакции опоры равна нулю, либо они не сближаются и сила реакции опоры неотрицательна. Другими словами, контактные точки либо отдаляются, либо обмениваются неотрицательным нормальным импульсом, но не одновременно. Систему (2.49) можно записать и в более симметричном виде: О Л _1_ АХ — 6 0, (2.51) где под знаком _1_ понимается ортогональность в смысле равного нулю произведения строки на вектор.
Таким образом, окончательную систему уравнений, определяющую эволюцию во времени системы твёрдых недеформируемых тел, движущихся под действием ограничений, пригодную для численного решения, объединив (2.1), (2.2), (2.45) и (2.49), можно сформулировать следующим образом:
Явный вид якобиана jT, соответствующего контакту пары тел можно получить из (2.18). Физический смысл контактного взаимодействия - непроникновение контактирующих точек. Другими словами, относительная скорость точек контакта в проекции на нормаль контакта должны быть неотрицательной. Пусть контакт тел с номерами кх и к2 произошёл в точке q и нормаль к поверхности в этой точке - п. Тогда полная скорость скорость точки q тела с номером п\ складывается из линейной и угловой по закону: k = к! + (wid, q — Ркі) j (2.53) где () - векторное произведение. Аналогичным образом выражается скорость этой же точки q на втором теле: здесь сверху индексами к\ и к2 показаны номера элементов в строке jT, а дополнительными скобками показаны пары векторов, относящиеся к линейной и угловой скоростям одного твёрдого тела системы соответственно. Якобиан такого вида будет встречаться в большом количестве и других соединений.
В случае неупругого соударения контактное ограничение на скорость записывается следующим образом:
Упругое же соударение моделируется, используя коэффициент отскока а Є (0..1). Данный коэффициент определяется геометрией и материалом, либо устанавливается на основе экспериментальных данных. Коэффициент а = 0 соответствует абсолютно неупругому соударению, а а = 1 - абсолютно упругому. Абсолютно неупругий удар характеризуется тем, что сразу после столкновения относительная скорость тел в проекции на нормаль взаимодействия обнуляются. Абсолютно упругий удар означает, что в результате столкновения проекции скорость на нормаль контакта меняет знак. Общее уравнение, объединяющие эти два случая, можно записать следующим образом:
Выведенный в главе (2.1.1) якобиан ограничивает лишь нормальное движение тел, то есть учитывает лишь силу реакции опоры. В практических задачах часто используется более сложное контактное условие, в котором присутствует трение. Существуют различные типы трения, в данной главе будет рассмотрен подход к моделированию трения скольжения, трения качения [30] и трения верчения [31].
Временная когерентность, горячий старт
Для составления разностной схемы метода Галёркина [37] мы будем рассматривать гиперболическую систему в общем виде (2.93). Область интегрирования П разбивается на конформные тетраэдральные элементы Т , где m - уникальный номер каждого тетраэдра. Матрицы Apq, Bpq и Срд считаются постоянными внутри каждого отдельного тетраэдра.
Каждому тетраэдру ставится в соответствие локальная система косоугольных координат (, rj, (), где базисными векторами являются рёбра, выходящие из одной из его вершин. Численное решение Q в каждом тетраэдре представляется в координатах (, rj, () как Qp ( V) С) t) = QjT W i ( V) С) (3.32) где Ф[ (, rj, () - функциональный базис, зависит только от пространственных координат, Qpf (t) - коэффициенты при базисных функциях для тетраэдра га, зависят только от времени. Выбор конкретного функционального базиса определяется характером расчёта. В данной работе использовались как ортогональные полиномиальные базисы, так и неортогональные. Здесь и далее, если не обозначается обратное, то суммирование подразумевается только по повторяющимся нижним индексам. Домножая (2.93) на одну из базисных функций Фк и интегрируя по тетраэдру Т(т), получаем:
Решение на границе тетраэдров разрывно, поэтому вводится численный поток через границу с соседним тетраэдром под номером щ, j є [1..4]: Н = TLAt)Qi (3.35) где – матрица перехода в локальную систему координат, в которой осью является нормаль общей грани тетраэдров и , – предел решения задачи Римана распада разрыва на границе этих тетраэдров со стороны тетраэдра , в локальной системе координат. В рассматриваемом нами случае многомерную задачу Римана распада разрыва можно приближённо свести к одномерной задаче распада разрыва по оси .
Обратим внимание, что для краткости здесь и далее мы будем рассматривать один тетраэдр под номером и вместо громоздкой полной записи, в которой бы фигурирвали индексы обоих контактирующих тетраэдров и , мы будем указывать только индекс . Например: Tpq Ti)qn mj (3.36) Tlj - Щт,т-) (3.37) Матрица перехода T , появившаяся в (3.35), при домножении справа на вектор неизвестных Qi переводит его из локальной системы координат, где базисом являются п = Hj, s = Sj, t = tj, в мировую. При этом тройка п, s, t является ортонормированной, а векторы s и t определяются лишь с точностью до поворота относительно
Элементы матрицы Лрд, соответствующие отрицательным собственным числам, берутся из матрицы Ард ), а положительные – из Apj). Матрица Rpq состоит из столбцов ДЙ), соответствующих отрицательным собственным числам, и из столбцов щ , соответствующих положительным.
Рассмотрим структуру решения в области разрыва в системе координат, где осью х является внешняя нормаль Hj общей грани тетраэдров Т -"1) и Т т).
Благодаря гиперболичности системы, каждому её собственному числу будет соответствовать скачок в кусочно-постоянном решении, двигающийся со скоростью . Обозначим области постоянного решения между разрывами как - , , , , , + . При некоторые из обозначенных величин тождественны ранее введённым:
Обратим внимание, что для расчёта потока нам понадобится лишь величины существенны лишь для промежуточных расчётов.
Согласно соотношениям Ранкина-Гюгонио, для ненулевых собственных чисел выполня ется: область постоянного решения слева от границы разрыва, соответствующей собственному числу , а (+) – справа. Под матрицей sgn() понимается матрица той среды, в которой распространяется волна , в зависимости от её знака.
Для установления соотношения между областями решения и , разделяемыми переходом, соответствующим нулевому собственному числу, используются граничные условия вроде (2.99) (2.100).
Дополняя последнюю систему уравнений контактными условиями, например, условием слипания (2.99), получаем: Q — Qp = a1Rp1 + Q!2 »2 + a3Rp3) Qp Qp = a7Rp7 + a8Rp8 + 9 „9, (3.48) С Uj = (JC Uj Решая полученную систему, получаем: Q = ( 4( "I as (-Мі - st ) (R3)tr ) Qr "I as ( l it \) (R3)tr Q+ = II 2y qr МаГ 5; + С V, (3.49) где было введено обозначение: M3 n = A() -\—Ris (A3st — \A3st\) (R3)tr , 4 4 4 I (3.50) (3.51) иде: (3.52) j — 1 j\/]_ = — _/xTV + — TV + (it ) 4 4 \ I 2 y Наконец, подставляя в определение численного потока (3.35), поток представляется в виде: F3 = contactJ(Qc , QJnj ) = Т3М3Сп (Т3) О с + T3M3fu (Т3) О Р Р v s pq qr rs s pq qr rs s Подставляя получившееся значение численного потока (3.52) вместе с дискретизацией решения (3.33) в (3.34) и разбивая интеграл через по границе тетраэдра как 4 интеграла по его граням, получаем:
Уравнение (3.53) записано в мировой системе координат, однако если перейти к локальной системе координат, связанной с единичным тетраэдром, то все фигурирующие в нём интегралы можно предрассчитать заранее, что позволяет существенно сократить время расчёта. Для определённости пронумеруем грани тетраэдра следующим образом: Грань Вершины
Использование экстремального отображения для неявного задания геометрии
Важно отметить, что данный алгоритм очень быстро сходится на реальных задачах. При создании иллюстрации оказалось достаточно сложно вбрать такие начальные условия, чтобы было возможно проиллюстрировать хотя бы 5 итераций алгоритма, так как он сходится настолько быстро, что точки симплекса сливаются при уменьшении ошибки уже на 2-3 итерации. Другое важное свойство GJK - возможность использования временной когерентности системы. Приняв гипотезу, что за малое время At геометрии сдвинутся на малое расстояние, алгоритм GJK в момент времени t + At может стартовать с того вектора, на котором сошёлся в момент времени t. На практике это означает, что в случае медленного движения алгоритм практически всегда сходится за 1 итерацию
Как уже замечалось в (4.1.5), GJK – очень мощный алгоритм, с прекрасными показателями сходимости, но имеющий принципиальное ограничение. Ограничение заключается в том, что он разрабатывался для поиска расстояния непересекающихся геометрий. Основные подходы, решающие задачу поиска расстояния пересечения для геометрий, заданных с помощью экстремального отображения - это GJK-EPA и алгоритм фильтрации порталов Минковского. В ходе данной работы автором был разработан альтернативный подход, названный SGJK. Более подробно о каждом из этих алгоритмов пойдёт речь ниже. GJK-EPA Для расширения GJK на случай пересекающихся геометрий используется ставший уже стандартным алгоритм GJK-EPA [52], который можно сформулировать следущим образом:
Algorithm 14 Алгоритм GJK-EPA поиска расстояния между пересекающимися геометриями данный алгоритм носит название алгортма расширяющегося многогранника(англ. GJK-EPA, GJK-Expanding Polytope Algorithm). Основное отличие GJK-EPA от GJK заключается в том, что в случае непересекающихся геометрий достаточно хранить лишь 4 точки симплекса, на каждом шаге отсеивая лишние с помощью функции ClosestSubsimplex(), в EPA же приходится хранить весь симплекс, в результате чего он обладает куда худшими показателями производительности, нежели GJK, так как на каждом шаге необходимо находить величину -- (Simplex), сложность этой операции возрастает с увеличением числа граней симплекса. Данная проблема, судя по всему, носит принципиальный характер, так как автору диссертации на момент её написания не известно о существовании эффективных подходов к оптимизации этого алгортма. Самый распространённый способ обхода проблем, связанных с низкой производительностью GJK-EPA – ввод так называемого пограничного слоя (англ. Skin Width), который констатирует контакт, если расстояние между телами стало меньше определённой
Разумеется, данный подход не является полноценным решением, так как требует достаточно больших значений для устойчивой работы, что порождает существенные артефакты и не исправляет принципиальную проблему: при пересечении более чем на , такой алгоритм по-прежнему применять нельзя и приходится использовать медленный
Существует также особый алгоритм, называемый фильтрацией порталов Минковско-го(англ. Minkowsky Portal Refinement [54]), ищущий вместо глобального минимума расстояния от начала координат до поверхности конфигурационного пространства, некоторый локальный минимум, используя итерационную схему вида:
Другими словами, алгоритм в качестве приближения vi+\ выбирает такое направление v, для которого SBQA{V) является точкой пересечения CSO и луча щ. Ключевой недостаток алгоритма, из-за которого он не получил широкого распространения, кроется именно в поиске точки vi+1, а точнее итерационной схемы, использованной для её аппроксимации. Схема заключается в последовательном разбиении треугольника, пересекаемого лучом щ, точки которого лежат на конфигурационном пространстве, называемого порталом, пока его площадь не будет достаточно малой. Треугольник разбивается на три путём добавления вершины !Аев(щ), где щ - нормаль к плоскости треугольника на г-й итерации. При последовательном разбиении треугольник часто вырождается практически в отрезок с неопределённой нормалью, в результате чего из-за численной ошибки алгоритм может плохо сходиться или зациклиться.
Идея поиска локального минимума вместо глобального, используемая в алгоритме фильтрации порталов Минковского, показала себя на практике исключительно хорошо. Её успех объясняется тем, что в практических приложениях тела редко пересекаются на значительную глубину, поэтому локальный минимум при практически любом начальном приближении совпадает с глобальным. Если же пересечение настолько существенно, что локальных минимумов несколько, то это на практике соответствует случаю, что тела очень сильно пересеклись в результате нештатной работы программы и при этом по большому счёту нет особой разницы, в каком именно направлении следует их расталкивать.
Автором был разработан алгоритм, названный SGJK, эффективно использующий сильные стороны GJK, но в отличие от него способный искать расстояние пересечения. Как и алгоритм фильтрации порталов Минковского, в целях существенной оптимизации скорости работы, он использует поиск локального минимума вместо глобального.
Суть алгоритма заключается в последовательном движении некоторой точки относительно CSO таким образом, чтобы она гарантированно находилась вне CSO, но с каждой итерацией гарантированно приближалась к началу координат. При этом на каждой итерации необходимо отыскивать расстояние от текущей точки до CSO, для чего можно использовать GJK в чистом виде со всеми его достоинствами, так как точка гарантированно находится вне CSO.