Назад, на численное моделирование
Корректность численного решения задачи зависит от правильности задания начальных и граничных условий и, как правило, расчет кинематических параметров на контактных поверхностях является наиболее сложной и трудоемкой операцией, так как для контактирующих узлов требуется производить дополнительные действия, согласующие скорости и перемещения.
Желающих более подробно ознакомиться с алгоритмами и подходами к численному моделированию контактных взаимодействий, рекомендую недавно вышедшую, достаточно полную обзорную статью Бураго Н.Г., Кукуджанов В.Н. Обзор контактных алгоритмов// МТТ, 2005, №1, с.45-87. Помимо неплохого обзора, в ней описана некоторая классификация подходов и алгоритмов, которая позволит читателю лучше оценить текущее состояние этого вопроса и направления, в которых ведутся дальнейшие поиски.
При численном моделировании задач интенсивного взрывного и ударного нагружения, для определения скоростей и координат узлов, находящихся на границе раздела двух сред, используются различные подходы.
Метод сквозного счета является наиболее простым и, в то же время, корректным методом. При данном подходе узлы расчетных сеток на контактной границе совпадают и определение скоростей проводится, как для обычных внутренних узлов. Однако, область применения метода сквозного счета сильно ограничена по ряду причин: необходимость совпадения узлов требует совпадения размеров ячеек контактирующих тел и не позволяет рассчитывать скольжение одного материала по другому. Фактически, метод применим только при моделировании образца со статической трещиной. Кроме того, для уменьшения погрешности требуется дополнительное согласование размеров ячеек (в зависимости от плотности контактирующих областей).
Метод зеркальных ячеек не требует совпадения узлов контактирующих сред и отличается от метода сквозного счета введением слоя фиктивных «зеркальных» ячеек. При этом возможно рассчитывать скольжение расчетных сеток друг по другу, отскок при ударе, однако трудоемкость определения напряжений в «зеркальных» ячейках и неприменимость метода для расчета узлов острых осколков не позволяет использовать его для задач фрагментации.
Подход Master- Slave основан на том, что для сильно разноплотных сред, массой менее плотной среды (например, газа) можно пренебречь, учитывая только ее давление на более плотную среду, снося при этом по нормали граничные узлы газа ( Slave- линия) на поверхность более плотного материала ( Master - линия). Если плотности контактирующих сред близки, то можно на каждом шаге по времени менять Master на Slave и наоборот, однако в задачах фрагментации, при сильно развитой поверхности разрушения, невозможно разделить контактные границы на Master и Slave, так как они зачастую являются продолжением друг друга.
Также при моделировании взаимодействия тел широко используются комбинации различных методов и подходов:
Особенности реализации контактных условий для тел с высокой криволинейностью поверхности. Коррекция перемещений.
Классический подход к реализации граничных условий на контактных поверхностях в Лагранжевой системе координат состоит в том, что скорости граничных узлов устанавливаются в соответствии с законами сохранения импульса и момента импульса, а их координаты корректируются путем сноса узлов по нормали к границе контакта.
Рассмотрим недостатки классического подхода на примере задач взрывного или высокоскоростного ударного дробления, когда образуется множество осколков с очень развитой поверхностью разрушения. Если взаимодействующие поверхности равноправны (например, два фрагмента одного и того же тела), то к ним невозможно применить поход Master- Slave, то есть выделить Master-поверхность, на которую производится снос точки. Таким образом, требуется взаимная коррекция координат узлов. Однако, если поверхность контакта криволинейна, невозможно обеспечить выполнение требования расположения всех узлов на контактной поверхности: либо для некоторых узлов (узлы 3,5 на рис.1) условие контакта не будет выполняться, либо некоторые узлы будут заходить за поверхность противолежащего тела (узел 7 на рис.1).
рис.1
В случае, если поверхность контакта прямолинейна (соударение пластин) или узлы на противолежащих поверхностях совпадают (вложенные соосные оболочки), можно расположить узлы обоих поверхностей вдоль линии контакта, соблюдая условие непроникновения узлов в противолежащий материал, однако эти случаи являются весьма редкими и при моделировании задач дробления на них ориентироваться не стоит. Даже когда поверхность контакта первоначально плоская (например, соударение пластин), за счет погрешности, обусловленной выбором начальной точки (узла на контактной поверхности, с которого начинается процедура коррекции), поверхность контакта может слегка отойти от плоскости.
Если считать допустимым, что при реализации итерационной процедуры сноса узлов на контактную поверхность некоторые из них могут отойти от линии контакта, то вывод точки из «зацепления» приводит к резкому изменению ее скорости (пропорциональному давлению в точке контакта) и, как следствие, к нефизическим осцилляциям, для устранения которых требуется значительно увеличивать искусственную вязкость.
Возможна ситуация, когда два остроугольных осколка, первоначально совпадающих вершинами (рис.2,а), «зашли друг за друга» (рис.2,б,в). С классической точки зрения вершины этих осколков не взаимодействуют ни с какой поверхностью, они являются свободными узлами, так как не находятся внутри материала. Если в дальнейшем осколки начнут расходиться и один из узлов окажется внутри материала (рис.2,г), попытка реализовать его снос на поверхность может оказаться весьма плачевной (такого резкого изменения формы и плотности элемента программа может не выдержать). Данная ситуация встречается редко, однако реализация условия проверки подобного взаимодействия значительно снижает скорость численного расчета.
Рис.2
Исходя из вышеизложенного, можно сделать вывод о том, что при реализации граничных условий на контактных поверхностях коррекция перемещений не является жизненно необходимой процедурой, более того, в некоторых случаях она приводит с существенным погрешностям, не говоря уже о машинном времени, затрачиваемом на данную операцию. Для обеспечения в задачах численного моделирования условия непроникания по нормали необходимо, чтобы при возникновении контакта не происходило дальнейшее проникновение узла в материал. Для выполнения данного условия вполне достаточно обеспечить равенство нормальных составляющих скоростей в точке контакта. Захождение узла в материал на глубину порядка нескольких процентов от размера ячейки не только не должно запрещаться, но наоборот, должно приветствоваться, ибо в данном случае контакт получается более надежным, узел не сможет «случайно» выйти из контакта, тем самым исключаются нефизические осцилляции, вызванные резким изменением скорости узла.
При выполнении условия контакта (т.е. при захождении узла в ячейку противолежащей области) требуется дополнительно определить нормаль к поверхности контакта и набор узлов, с которыми должен осуществляться обмен импульсами. Эта, на первый взгляд, элементарная задача, при определенных обстоятельствах (контактные границы весьма сложной формы, при развитой поверхности разрушения, при взаимодействии острых осколков) вызывает значительные затруднения, так как в этих случаях понятие нормали к контактной поверхности становится несколько неопределенным.
Неопределенность вызвана следующими причинами:
- в ряде случаев точку контакта, а, следовательно, и нормаль, нельзя определить однозначно;
- нормаль к поверхности в точке для каждой из контактирующих поверхностей может иметь различное значение, особенно если одна из контактирующих поверхностей представляет из себя осколок, а их сближение происходит под углом;
- поскольку граница дискретна, «нормаль в точке» может не совпадать с реальной нормалью (например, в случае с окружностью), а нормаль в узле вообще не определена.
Граничный узел считается контактирующим с противолежащей поверхностью, если он находится внутри одной из ее ячеек. Математически выразить условие контакта можно разными способами. Например, если считать, что нормали всех плоскостей тетраэдра ориентированы наружу (см. рисунок), то необходимым и достаточным условием того, что узел U находится внутри (или на границе тетраэдра), будет отрицательное (или равное нулю) значение скалярного произведения для каждой плоскости 
, где Ri - радиус-вектор от произвольной точки i-ой плоскости до узла U.
Проверять контакт граничного узла с каждой ячейкой, естественно, нереально. Для каждого граничного узла определяется его ближайшее окружение – набор ячеек, центры которых попадают в сферу радиусом Rs (центр сферы находится в узле). Из списка привязки стоит исключить ячейки, в которые входит сам узел.
Для достаточно подробных сеток (сотни тысяч ячеек) определение набора ячеек для каждого узла является достаточно долгой процедурой, поэтому, она проводится через определенный промежуток времени. Число шагов по времени, через которое переопределяются списки ячеек, находится из условия сохранения актуальности списка – узел не должен за указанный промежуток времени приблизиться к другим ячейкам, поэтому радиус сферы привязки я беру как Rs=0.4*dx+0.75*Rmax (dx - средняя сторона ячейки, Rmax - максимальное ребро ячейки, соответственно 0.75*Rmax - максимальный радиус окружности описанной вокруг ячейки, с центром, совпадающим с центром ячейки). Допустим, что максимальная (из всех узлов) скорость равна Vmax, тогда за момент времени 
, определяемый из условия 
-0.75*Rmax, узел гарантированно не зайдет ни в одну ячейку, не попавшую в список привязки. Учитывая, что вероятность того, что узел и ячейки из списка привязки движутся в противоположных направлениях, причем с максимальными скоростями, ничтожно мала, промежуток времени можно увеличить, смягчив условие: 
-0.75*Rmax.
Как показывает практика, величину Rs-0.75*Rmax (я беру ее равной 0.4*dx) необходимо выбирать из условия пустого списка привязки на свободной поверхности, т.е. меньше средней высоты тетраэдра. При больших значениях резко увеличивается число ячеек в списке привязки, а при меньшем – слишком маленьким становится период Ts (при числе ячеек порядка нескольких сотен тысяч слишком частое переопределение списка привязки слишком замедляет счет).
На каждом шаге по времени, для узлов, находившихся в контакте, проверяется условие сохранения контакта (достаточно экономичная процедура). Для узлов, не контактировавших, проверяется факт захождения в ячейки из списка привязки.
В качестве нормали к поверхности контакта берется средняя нормаль поверхностных плоскостей, окружающих граничный узел.
Применение индексов при определении списка привязки.
Прямой перебор ячеек, даже при таком редком действии, как определение списка привязки, является слишком трудоемкой операцией. Имеет смысл рассматривать только ячейки, лежащие в определенной близости от узла.
Данная идеология широко применяется в бессеточных (MeshFree) методах, для обозначения используемых при этом подходов употребляются термины "метод динамической сетки" и "Grid-алгоритм".
Как же определить понятие близости для неструктурированных сеток (да и для структурированных в большинстве случаев)? R- tree индексы применять смысла нет, система динамическая, переопределение индексов будет слишком дорогой операцией, да и позиционирование в R-дереве тоже операция не из дешевых. Напрашивается самый простой выход – поделить пространство плоскостями, перпендикулярными осям X, Y, Z и на каждом участке иметь список указателей на попавшие в него ячейки. Для двумерных задач это еще прокатывает, но для трехмерных уже нет – несмотря на то, что относительная доля заполненных участков мала (большинство участков пусты), приходится хранить их все в виде трехмерного массива и объем этого массива превышает все разумные пределы (в некоторых из моих задач более 2 Gb). Я применяю следующую схему индексации, которая позволяет снизить время составления списка привязки на порядок (по сравнению с прямым перебором, даже если при нем отсеиваются ячейки, лежащие относительно далеко), тем самым это процедура перестает быть самой критичной с точки зрения машинного времени:
1) Берется обычный пустой ListYchNoGroup: TStringList (свойство Sorted по умолчанию имеет значение false и делать его true ни в коем случае нельзя, иначе каждое добавление нового элемента будет перестраивать индекс – даже для 10тыс. элементов это будет ужасно).
2) Область считается виртуально поделенной на кубические участки, имеющие индексы i, j, k с шагом, равным Rs по каждой координате.
3) Для каждой ячейки определяются индексы i, j, k участка, в котором находится ее центр и добавляется в ListYchNoGroup указатель на ячейку, используя AddObject. Строка, соответствующая ячейке, представляет из себя IntToStr(i)+';'+IntToStr(j)+';'+IntToStr(k).
4) Сортируем ListYchNoGroup, используя метод Sort. Данная процедура весьма незначительна по трудоемкости. Кстати, это правило справедливо практически для всех баз данных и используемых в них индексов: при добавлении (а зачастую, как ни странно и при удалении) большого количества элементов, время добавления элементов с отключенным индексом и последующей переиндексацией (пересозданием индекса) может отличаться от времени добавления со включенным индексом на порядки!
5) Берем пустой ListIndex:THashedStringList (свойство Sorted по умолчанию имеет значение false и менять его не в коем случае нельзя, по тем же причинам). THashedStringList – это TStringList, имеющий внутренний хранимый индекс для обеспечения быстрого поиска, именно это нам и требуется.
6) Идем по отсортированному списку ListYchNoGroup и для каждого участка i, j, k, имеющего ячейки, в ListIndex добавляется массив этих ячеек, с соответствующей строкой. Если проводить аналогию с базами данных, накопление элементов выполнялось с отключенным индексом в одной базе, далее она была отсортирована и в другую базу перенесены уже сгруппированные значения.
7) По завершении, ListYchNoGroup очищается – больше он нам не понадобится. ListIndex сортируется, используя метод Sort.
8) В дальнейшем, при определении списка привязки для конкретного узла, рассматриваем только ячейки, расположенные в участках с индексами, попадающие в диапазон [ i0-1.. i0+1, j0-1.. j0+1, k0-1.. k0+1] от индекса участка ( i0, j0, k0), в котором находится сам узел. Собственно, именно для этого нам и нужен был THashedStringList, чтобы IndexOf работал очень быстро.
Получается, что время определения списка привязки зависит от числа узлов (или ячеек) практически линейно, в то время как при прямом переборе – квадратично. Учитывая, что определение списка привязки производится еще и не на каждом шаге по времени, получается достаточно щадящая процедура…
Коррекция скоростей по принципу Master- Slave мало приемлема потому, что:
Рассмотрим более общий алгоритм. На первом этапе все узлы на контактных поверхностях считаются не контактирующими и их скорость определяется согласно стандартному алгоритму (как для узлов на свободной поверхности). Затем проверяется условие контакта, в качестве которого принимается факт захождения узла в противолежащую область.
рис.3
Рассмотрим узел U0 (рис.3), который пересек граничную плоскость и находится внутри материала. Если скорость данного узла превышает скорость поверхности противолежащего материала, то есть вектор относительной скорости направлен внутрь материала, включается процедура коррекции скоростей. Для трехмерной численной схемы с тетраэдрической сеткой граничная плоскость (с которой контактирует узел) будет представлена в виде треугольника и контактирующий узел обменяется импульсами с тремя узлами этой плоскости. В ряде случаев, более удобным является подход, когда считается, что узел обменивается импульсами с ячейкой, за поверхность которой он зашел.
При взаимодействии контактирующих поверхностей (обмене импульсами) выполняются следующие соотношения:
1) Закон сохранения импульса 
(3 независимых уравнения) - данный закон выполняется во всех случаях;
2) Условие распределение приращений скоростей между взаимодействующими узлами. При выборе данного условия необходимо, чтобы соблюдался закон сохранения момента импульса 
. Поскольку закон сохранения момента импульса представляет всего 3 уравнения, а всего неизвестных компонент 12 (в случае с плоскостью), или 15 (в случае с ячейкой), его приходится дополнять другими соображениями (например, для абсолютно твердого тела тангенциальные приращения скоростей совпадают, а для деформируемого зависят от положения U0). Даже в 2 D расчетах, при определении тангенциальных составляющих скоростей, приходится учитывать, что они находятся на одной линии и закон сохранения момента импульса выполняется тождественно. Наиболее естественно было бы использовать в качестве этих соображений условие пропорционального (в зависимости от расстояния до точки контакта) распределения приращений скоростей 
, или условие пропорционального распределения по массе 
(здесь в качестве xi поочередно выступают проекции на все 3 оси).
В случае, если масса узлов не учитывается, невозможно обеспечить соблюдение закона сохранения момента импульса, поэтому мы делаем выбор на условии пропорционального распределения по массе (в этом случае закон сохранения момента импульса является частным случаем данного распределения)
6 независимых уравнений (в случае контакта с плоскостью)
В случае контакта с ячейкой добавляется еще 3 уравнения 
3) Условие равенства скоростей в точке контакта (при выполнении условия прилипания):
3 независимых уравнения, здесь pi – относительный вклад каждого узла в значение скорости в точке контакта ( 
).
рис.4
pi= h/ H отношение перпендикуляра из точки контакта на противолежащую сторону к высоте, опущенной из i-ого узла (рис.4).
При скольжении данная система вырождается в одно уравнение, для нормальной составляющей скорости.
4) Тормозящее действие силы трения (при скольжении)
2 независимых уравнения, здесь 
- единичный вектор, совпадающий с направлением 
, 
- нормальная составляющая силы, действующей на узел U0, Fn - сила нормального давления, k-коэффициент трения.
Стоит обратить внимание, что для данной системы уравнений не выполняется закон сохранения энергии. Обмен импульсами (т.е. торможение узлов), вызывает рост давления в приграничных областях и изменение потенциальной энергии. Если не учитывать это обстоятельство, в определенных задачах дисбаланс энергии может достигать значительных величин. Это возможно в ситуациях, когда менее плотная среда с большим градиентом скоростей тормозится на более плотной среде. Например, при подрыве безоболочечного заряда ВВ (на расстоянии от плотной поверхности) и последующем расширении продуктов детонации, практически вся энергия системы заключена в кинетической энергии, причем максимальной скоростью, а, значит, и максимальной энергией, обладают именно поверхностные узлы. При их торможении на плотной поверхности, величина дисбаланса энергии сравнима с кинетической энергией поверхностных узлов в продуктах детонации. Если сетка в газе слишком грубая, это вызывает большой градиент скоростей у поверхности и величина дисбаланса может достигать десятков процентов.
При попытке учесть закон сохранения энергии dЕп1+dЕк1+dЕп2+dЕк2=0, появляется дополнительно только одно уравнение, а неизвестных два - изменение потенциальной энергии в 1-ом материале и во 2-ом. Поэтому требуется еще одно уравнение, связь между изменениями потенциальной энергии в обоих телах (материалах)
Предположив dЕп1= dЕп2=0, dЕк1=- dЕк2, мы получим абсолютно упругий удар, без переходных процессов, что, конечно, не соответствует постановке задачи.
На самом же деле, закон сохранения энергии при соударении формулируется в виде
dЕп1=- dЕк1;
dЕп2=- dЕк2.
То есть вся разность энергии перешла в потенциальную энергию, причем обмена потенциальной энергией между
взаимодействующими телами нет. При таком подходе закон сохранения энергии выполняется автоматически и в явном виде не требуется
при определении изменения скоростей
Рассмотренные выше соотношения удобно разделить на содержащие нормальные компоненты скорости и на содержащие тангенциальные соотношения скорости. В этом случае, система для определения нормальной составляющей не меняется при смене условия идеального скольжения на скольжение с трением.
При взаимодействии с плоскостью, система линейных уравнений четвертого порядка для определения скоростей 4-х узлов будет состоять из 12 уравнений:
Компоненты скорости, нормальные к границе раздела, находятся из системы, состоящей из 4-х уравнений:
(2.1)
здесь 
(проекция векторов 
на 2 взаимно перпендикулярных направления 
, лежащих в плоскости контакта).
В случае, если используется условие идеального скольжения, тангенциальные компоненты скоростей остаются без изменения.
В случае условия прилипания для тангенциальных компонент скоростей решается система уравнений, аналогичная системе для нормальных компонент:
(2.2)
Если используется скольжение с трением, вначале, предполагая прилипание, определяется 
. Если 
, то скольжение отсутствует и условие прилипания остается в силе. В случае, если 
(силы трения не хватает, чтобы полностью затормозить узел), для определения тангенциальных скоростей используется система (2.3), в которой, по сравнению с системой (2.2), условие равенства скоростей в точке контакта заменено уравнением, определяющим тормозящее действие силы трения.
(2.3)
В дальнейшем узел двигается совместно с данной плоскостью до момента возникновения в зоне контакта растягивающих напряжений. При этом вектор относительной скорости узла будет направлен от контактной поверхности и узел постепенно отойдет от нее.
Погрешность, связанная с выбором начальной точки.
При решении вышеописанной системы уравнений 2.1-2.3 для каждого узла (например, в случае контакта двух пластин), изменение скорости для первого узла, к которому применили процедуру коррекции скоростей, и для последующих, будет различаться. Это связано с тем, что для первого узла решается система с невозмущенными значениями скоростей, а для последующих – с уже скорректированными скоростями.
Данная проблема проявляется при шахматном расположении узлов контактирующих поверхностей (когда узлы одной поверхности соответствуют серединам сторон другой), при торможении менее плотной среды возможен даже отскок. Рассмотрим этот вопрос на примере двухмерной осесимметричной задачи: коррекция скоростей граничных узлов внутреннего (индекс А) и внешнего (индекс В) тел, соприкасающихся по окружности. Допустим, что узлы расположены в шахматном порядке, что тела находятся в равновесии и определим погрешность вышеописанного подхода к коррекции скоростей.
Пусть граница соприкосновения каждого тела состоит из 160 узлов. Это достаточно мелкая сетка для того, чтобы не принимать во внимание погрешности дискретизации. Поскольку тела находятся в равновесии, импульсы, которыми они должны обменяться, равны maVa= mbVb. При отсутствии погрешности, в результате обмена импульсами мы должны получить равенство скоростей нулю на обеих поверхностях.
При однократном решении для каждого узла обоих поверхностей системы уравнений 2.1-2.3, наблюдается следующая закономерность:
1) выделяется зона начальной точки – несколько узлов, где остаточная скорость составляет до 12%;
2) во всех остальных узлах, которые назовем регулярной зоной, остаточная скорость практически одинакова и составляет около 2%.
Повторная итерация снижает погрешность в зоне начальной точки примерно в 2 раза, в регулярной зоне – примерно в 50 раз.
При 9 итерациях погрешность в зоне начальной точки снижается до 1%, в регулярной области – до 0.01%.
В случае, если массы узлов областей А и В различны, остаточная погрешность меньше в том случае, если первой система уравнений решается для области с меньшей массой узлов.
В трехмерном случае остаточная погрешность меньше, чем в двухмерном, т.к. узел обменивается импульсами не три раза, а четыре – это аналогично дополнительной итерации. Рассмотрим варианты подходов к уменьшению погрешности процедуры коррекции скоростей и, в частности, погрешности выбора начальной точки.
1) Данная проблема практически полностью (если нормаль к поверхности контакта определена) снимается, если узлы совпадают – в этом случае каждая система уравнений охватывает только 2 взаимодействующих узла и не влияет на остальные. Таким образом, можно рекомендовать зеркальное (или близкое к зеркальному) расположение узлов на контактирующих поверхностях.
2) Для корректности определения скоростей и устранения погрешности, связанной с выбором начальной точки, необходимо решать полную систему уравнений, включающую в себя изменения скоростей всех узлов данного пятна контакта.
3) Как показано выше, повторная коррекция скоростей значительно уменьшает погрешность как в регулярной зоне, так и в зоне начальной точки. Поэтому, если нет возможности реализовать предыдущие два варианта, можно остановиться на дублировании процедуры коррекции скоростей.
Стоит заметить, что для трехмерных расчетов при хорошем контакте погрешность, связанная с выбором начальной точки незначительна, а при взаимодействии осколков гораздо большую погрешность вносит неопределенность нормали, поэтому в задачах фрагментации я дублирование не использую. Однако в особо точных расчетах, где погрешность нежелательна даже в одном узле, итерации не помешают.
Назад, на численное моделирование