Алгоритм расчета контактных границ с учетом эрозии конечных элементов при высокоскоростном взаимодействии тел Текст научной статьи по специальности «Физика»

УДК 539.3

DOI: 10.20310/1810-0198-2016-21-3-1275-1277

АЛГОРИТМ РАСЧЕТА КОНТАКТНЫХ ГРАНИЦ С УЧЕТОМ ЭРОЗИИ КОНЕЧНЫХ ЭЛЕМЕНТОВ ПРИ ВЫСОКОСКОРОСТНОМ ВЗАИМОДЕЙСТВИИ ТЕЛ

© П.А. Радченко

Томский государственный архитектурно--строительный университет, e-mail: radchenko@live.ru

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

Ключевые слова: математическое моделирование; разрушение; динамическое нагружение; метод конечных элементов; упругость; пластичность.

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

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

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

ков, определяется размерность и выделяется память под соответствующие динамические массивы для используемых в алгоритме контакта переменных. Далее определяются коэффициенты А, В, С (направляющие косинусы вектора нормали к поверхности), Б уравнения Ах + Ву + С2 + Б = 0 , описывающего в пространстве плоскость каждой поверхностной грани конечного элемента, образуемую узлами пь п2, п3 (рис. 1).

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

Рис. 1. Общая схема контактного взаимодействия «элемент-элемент»

ISSN 1810-0198. Вестник Тамбовского университета. Серия Естественные и технические науки

Vi =

(A ■ x4 + B ■ y4 + C ■ z4)

(A ■ (X4 - Xi) + B ■ (У4 - y ) + C ■ (Z4 - Zi))

, i = 1,2,3

Здесь х, у, 2 с индексом 4 - координаты внутреннего узла рассматриваемого тетраэдра, с индексом , - координаты поверхностных узлов ть т2, т3, а А, В, С, Б -коэффициенты уравнения плоскости, образуемой узлами пь п2, п3 соответственно. Если коэффициент 0 < ц < 1, то соответствующее ребро проникло сквозь поверхностную грань конечного элемента. Поскольку точка пересечения ребра тетраэдра с плоскостью может находится и вне грани щп2п3, происходит проверка принадлежности точки пересечения соответствующему треугольнику. Для этого определяется сумма углов, образуемых данной точкой с узлами пь п2, п3. При ее равенстве 2л точка лежит внутри грани данного тетраэдра. Так как через одну плоскость могут проникнуть несколько ребер тетраэдра, но само ребро может пересекаться лишь с одним треугольником в пространстве, то организовывая цикл перебора всех возможных взаимодействий по узлам т1, т2, т3 каждого тела, можно получить крайне высокое распараллеливание данного алгоритма. По окончании данного цикла все данные о полученных прониканиях элементов заносятся в соответствующий массив.

В соответствии с моделью скольжения поверхностей без трения [3] проникший узел т2 сносится на граничную поверхность по нормали в положение т2' (рис. 1), в дальнейшем обозначим эту точку как С. Изменение количества движения данного узла в нормальном направлении определяется соотношением:

Js = Ms AV"

где Ms - масса узла т2; AVS" - нормальное изменение скорости данного узла. Импульс 1 передается трем узлам элемента п1, п2, п3, обозначим их ,, ], к соответственно. Доли импульса, передаваемые узлам треугольного элемента, обозначим

Ri = , RJ = , ^=,

где

Ji = MAVn, Jj = Mj AVj, Jk = Mk av".

Здесь т,, ту, ть, АV" , АУ" , АУк - масса и изменение нормальных скоростей узлов ,,], к соответственно. Тогда для определения АУ,", АУ", АУк имеют место соотношения

AV" = RMs а v/m,, AV" = RjMs AV"jmj, AV" = RkMs AV"/mk.

Для определения Ау ", АУ", АУк необходимо определить доли момента ^, Rj, Rk и изменение ско-

рости узла т2 АУ" . Доли импульса Ri, Rj, Я к находятся из условий сохранения количества движения и момента количества движения. Изменение скорости узла т2 определяется из условия равенства нормальных скоростей данного узла и поверхности элемента в точке С. Для определения скорости поверхности элемента в положении т2 используется линейная аппроксимация. Тогда выражение для скорости поверхности элемента и в направлении X имеет вид:

ис = ЧМ, + Чуиу + Чкик,

где ч^ = А"СкIЛук . В данной формуле ЛуЛ является

площадью треугольника, образованного узлами с индексами ] и к, а также точкой С, на которую сносится узел т2. Данный треугольник является противолежащим по отношению к узлу с индексом Лук - площадь всей грани, образованной узлами с индексами ,, ], к. Аналогично определяются константы 4 и 4 . Другие компоненты скорости ( ус и ^ ) определяются из выражения, аналогичного выражению для и .

Приравнивая нормальные скорости узла т2 и поверхности элемента в положении С, находим:

AV" =-

V" - V"

1 + qR,MjM, + qjRjMs/Mj + qkRtMjMt

тгП лтП

где У и У - нормальная скорость поверхности элемента и скорость узла т2, нормальная к поверхности элемента:

V" = Aus + Bvs + Cws,

V" = Au + Bv + Cw .

В данных уравнениях А, В и С - направляющие косинусы вектора нормали к поверхности элемента. Таким образом, сначала определяются А У" , а затем

однозначно вычисляются АУ" , АУ" и А У" .

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

Данный алгоритм позволяет адекватно описывать поведение тел различной геометрии при численном моделировании деформации и разрушения конструкций при импульсном нагружении [4]. Применение па-

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

СПИСОК ЛИТЕРАТУРЫ

1. Радченко П.А., Батуев С.П., Радченко А.В. Трехмерное моделирование деформации и разрушения гетерогенных материалов и конструкций при динамических нагрузках (EFES 1.0). Свидетельство о регистрации программы для ЭВМ № 2014614671.

2. Schoberl J. NETGEN - An advancing front 2D/3D-mesh generator based on abstract rules // Computing and Visualization in Science. 1997. V. 1 (1). Р. 41-52.

3. Johnson, G. R. Tree-dimensional computer code for dynamic response of solids to intense impulsive loads // Int. J. Numer. Methods Eng. 1979. V. 14. № 12. Р. 1865-1871.

4. Радченко П.А., Батуев С.П., Радченко А.В., Плевков В.С. Численное моделирование разрушения оболочки из бетона и фибробето-на при импульсном воздействии // Омский научный вестник. 2015. № 3 (143). С. 345-349.

БЛАГОДАРНОСТИ: Работа выполнена при финансовой поддержке РФФИ (грант № 16-31-00125 мол_а, грант № 16-38-00256 мол_а).

Поступила в редакцию 10 апреля 2016 г.

UDC 539.3

DOI: 10.20310/1810-0198-2016-21-3-1275-1277

ALGORITHM OF CALCULATION OF CONTACT BOUNDARIES TAKING INTO ACCOUNT THE EROSION OF FINITE ELEMENTS AT HIGH-VELOCITY INTERACTION OF BODIES

© P.A. Radchenko

Tomsk State University of Architecture and Building, Tomsk, Russian Federation, e-mail: radchenko@live.ru

In this paper we propose an algorithm of calculating contact boundaries of the interacting bodies. This algorithm maintains a high degree of parallelization which is important for high-quality finite element meshes. The erosion is to support high speed computations and also to maintain the stability of the solution and the time step at high loading rates. The adequacy of the proposed algorithm is confirmed by comparison with several experimental studies.

Key words: mathematical modeling; fracture; dynamic loading; finite element method; elasticity; plasticity.

REFERENCES

1. Radchenko P.A., Batuev S.P., Radchenko A.V. Trekhmernoe modelirovanie deformatsii i razrusheniya geterogennykh materialov i konstruktsiy pri dinamicheskikh nagruzkakh (EFES 1.0). Svidetel'stvo o registratsii programmy dlya EVM no. 2014614671.

2. Schöberl J. NETGEN - An advancing front 2D/3D-mesh generator based on abstract rules. Computing and Visualization in Science, 1997, vol. 1 (1), pp. 41-52.

3. Johnson, G. R. Tree-dimensional computer code for dynamic response of solids to intense impulsive loads. Int. J. Numer. Methods Eng., 1979, vol. 14, no. 12, pp. 1865-1871.

4. Radchenko P.A., Batuev S.P., Radchenko A.V., Plevkov V.S. Chislennoe modelirovanie razrusheniya obolochki iz betona i fibrobetona pri impul'snom vozdeystvii. Omskiy nauchnyy vestnik - Omsk Scientific Bulletin, 2015, no. 3 (143), pp. 345-349.

GRATITUDE: The work is fulfilled under financial support of Russian Fund of Fundamental Research (grant no. 16-3100125 Mon_a, grant no. 16-38-00256 Mon_a).

Received 10 April 2016

Радченко Павел Андреевич, Томский государственный архитектурно-строительный университет, г. Томск, Российская Федерация, кандидат физико-математических наук, доцент кафедры прикладной математики, e-mail: radchenko@live.ru

Radchenko Pavel Andreevich, Tomsk State University of Architecture and Building, Tomsk, Russian Federation, Candidate of Physics and Mathematics, Associate Professor of Applied Mathematics Department, e-mail: radchenko@live.ru