CC BY
193
Вычислительные технологии
Том 14, № 5, 2009
Динамический метод построения треугольных сеток в многосвязных областях
Е.И. Краус, В. М. Фомин, И. И. ШАБАЛИН Учреждение Российской академии наук Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, Новосибирск, Россия e-mail: kraus@itam.nsc.ru, fomin@itam.nsc.ru, shabalin@sbras.nsc.ru
Динамический метод построения регулярных 2Б-треугольных сеток Делоне распространен на случай многосвязных областей. Модифицирован алгоритм расчета контактных поверхностей учетом количества элементарных актов взаимодействия для каждого граничного узла сетки. Это позволило успешно выполнить численные расчеты соударения модельных блоков реактора ядерной энергетической установки с деформируемой поверхностью.
Ключевые слова: треугольная сетка, динамический метод, многосвязная область.
Для автоматического построения разностной сетки, как правило, требуется односвязность области. Известны алгоритмы и программы, позволяющие получать равномерные сетки на односвязных областях [1], к достоинствам таких алгоритмов можно отнести их универсальность по отношению к форме границы области. Односвязность достигается разделением исходной многосвязной области на ряд односвязных подобластей путем введения дополнительных разрезов [2]. Очевидный недостаток подобных алгоритмов — большой объем ручной обработки: разделение области на подобласти, разбиение каждой границы, ввод информации и т. д. К тому же, согласно [3], полученное решение может оказаться неудовлетворительным, т. е. область останется многосвязной, так как разрезы могут касаться старых границ. Обзор существующих методов генерации сеток представлен в работах [4-6]. Так, в работе [4] приведена систематическая классификация, основанная на порядке, в котором узлы и элементы создаются. Этапы, реализованные в существующих кодах, включают: 1) размещение основных узлов, их порядок и связи; 2) декомпозицию области на крупные блоки; 3) установление связей между блоками; 4) декомпозицию блоков области до уровня элементов (триангуляция Делоне). Эффективные алгоритмы для триангуляции Делоне изложены в работах [7, 8], при этом в общем случае элементы имеют произвольные стороны.
В работе [9] показано, что “плохие” треугольники (сильно вытянутые) при лагран-жевом подходе значительно увеличивают погрешность решения и замедляют его сходимость. Подобные ситуации возникают либо из-за сложной геометрической формы области, либо из-за “упущений” при программной реализации алгоритмов. Поэтому требуется не только провести триангуляцию счетной области, но и добиться, чтобы элементы сетки имели форму, близкую к правильным треугольникам. Одним из методов создания треугольной сетки является динамический подход, который позволяет
© ИВТ СО РАН, 2009.
реализовать триангуляцию как последовательность двух шагов: 1) распределение частиц (узлов будущей сетки) по области; 2) установление связей центров сети Вороного — преобразование в триангуляцию Делоне.
Следуя [10], кратко изложим суть динамического метода. Основой для выбора сил межчастичного взаимодействия является близость частиц. Поэтому вид силы выбирается так, чтобы вся система частиц в равновесии обеспечивала плотную упаковку или касание (малое перекрытие) сфер влияния узлов. Размер ячейки определяется расстоянием между центрами смежных частиц (рис. 1, a).
Форма силы межчастичного взаимодействия f (г) представляется подобно силе Ван-Дер-Ваальса (рис. 1, б). Сила отталкивания (притяжения) действует, когда две частицы находятся на расстоянии ближе (дальше), чем равновесное расстояние го. Заметим, что при таком виде межчастичной силы все “плохие” элементы физически неустойчивы, поскольку частицы стремятся занять вершины равностороннего треугольника (места устойчивого равновесия), и потому шанс создания таких “плохих” элементов значительно снижен.
Для примера, математическое выражение для силы представим в следующем виде:
f (г) = аг3 + Ьг2 + сг + d при 0 < г < 1.5г0, f (г) = 0 при г > 1.5г0, (1)
f (го) = f (1.5го) = 0.
Запишем уравнение движения г-й частицы (узла)
Ш<^ + К'^гТ = f-(г)- (2)
- з=^ -
где шг — масса частицы; Кг — коэффициент демпфирования; f г = ^ f ^, N — число
3 = 1
соседей г-й частицы, fгз- = f (ггз-), ггз- = у^(хг — х3)2 + (У* — Уз)2 — расстояние между частицами г и ]; хг — положение г-й частицы в счетной области.
В [11] для систем с одной степенью свободы и линейной силой f (х) = ах(£)
Ш^ + К М> + ах(() = 0 (3)
d^2 d^
получена следующая связь:
К = 1.5^Ша.
Заметим, что линейная аппроксимация силы используется только для того, чтобы определить нужные значения параметра X, и не используется при вычислении сил при динамическом моделировании системы частиц. В нашем случае а = f;|Г=Г0 (штриховая линия на рис 1, б).
Без потери общности можно принять, что масса всех частиц постоянна и равна единице: шг = 1. Константы в уравнении (1) можно найти из решения системы уравнений
для ключевых точек, например:
d =1,
а + Ь + с + d = 0,
а ■ (1.5)3 + Ь ■ (1.5)2 + с ■ (1.5) + d = 0, а ■ (1.25)3 + Ь ■ (1.25)2 + с ■ (1.25) + d = —0.3.
Притяжение
Рис. 1. Установление связей между центрами частиц и вид силы взаимодействия частиц
Для решения системы дифференциальных уравнений второго порядка с правой нелинейной частью расщепим каждое уравнение системы на два уравнения первого порядка по времени. Тогда система (2) приведется к виду
т,
ЇХі (і) Їі
Їщ (і)
(4)
Їі
+К.щ (і) = /і (і).
Разностная форма уравнений (4) с первым порядком точности имеет вид
гп+1 п
х"+1 = хп + Аіп Щ, <+1 - ип+Аіп (т) ( иП+12+ иП
т.
А іп Ї п т.
(5)
Разрешив второе уравнение (5) относительно ип+ , найдем соотношение для вычисления вектора скорости на следующем шаге по времени.
По начальным данным расположения частиц проводится интегрирование дифференциального уравнения по времени, пока система не достигнет равновесия, т. е. состояния, в котором максимальное смещение всех узлов за один шаг по времени не станет меньше, чем заданное малое значение £сг:
тах
і
\! «+1 - хп)2 + (уп+1 - уп)2 <
(6)
Проиллюстрируем применение динамического подхода к построению треугольной сетки на примере. На рис. 2 показано построение треугольной сетки в многосвязной
£
сг •
области. Первичной задачей при построении сетки стало создание “хорошего” начального распределения частиц. Затем равномерное распределение частиц было “испорчено” введением дополнительной внутренней границы (триангуляция на рис. 2, а). Разбиение вблизи геометрических границ стало хаотическим, без использования иерархического разбиения области. На рис. 2, а видно, что внесение граничных узлов в конфигурацию приводит к случайному числу “плохих” треугольников. Их число сокращается по мере движения частиц под действием сил, т. е. происходит плотная переупаковка частиц (триангуляция на рис. 2, б). И уже после нескольких десятков шагов численного интегрирования уравнения движения (2) разностная сетка с топологической нерегулярностью перестраивается в более регулярное состояние как вблизи границ тела, так и внутри многосвязной области (см. выделенные зоны 1, 2 на рис. 2).
а б
Рис. 2. Пример применения динамического подхода
Рис. 3. Примеры построения треугольной сетки в многосвязной области
Полный алгоритм построения разностной сетки следующий:
— определение каждого элемента конструкции как счетной области;
— погружение регулярной счетной области в мажорантную область треугольной сетки;
— отбрасывание узлов сетки, не попадающих в счетную область (например, “вырезание” отверстий в теле);
— расчет движения внутренних узлов сетки по уравнениям (5) при фиксированном распределении узлов на границе счетной области до выполнения заданного критерия
(6). Если в процессе движения внутренние узлы выходят за пределы границ счетной области, то такие узлы отбрасываются и не учитываются в дальнейшем;
— повторение вышеописанной процедуры во всех счетных областях.
На рис. 3 показаны результаты построения треугольных сеток в многосвязных областях динамическим методом. Из рисунка видно, что метод успешно может быть применен на произвольных многосвязных геометрических объектах. Необходимо только задать узлы сетки на границах этих объектов. Внутренняя структура перестраивается автоматически во всей области.
Пример расчета соударения геометрически сложных тел о преграду
Рассмотрим боковой удар модельного реактора космической ядерной энергетической установки с поверхностью Земли в двумерной постановке. В аварийных ситуациях современные космические аппараты с термоэмиссионными преобразователями “отстреливают” ядерную энергетическую установку (ЯЭУ). Однако существует определенная вероятность того, что часть реактора, содержащая ядерное топливо, несмотря на значительные тепловые и механические нагрузки при прохождении плотных слоев атмосферы, может достичь поверхности Земли. Скорость соударения оставшейся части реакторного блока способна достигать 400 м/с.
20
10
0
-10
-20
-20 0 20
Рис. 4. Плоская геометрическая модель реактора
Несмотря на современный уровень развития вычислительной техники и наличие достаточно реалистических математических моделей поведения материалов, решение задач удара реальных технических объектов получить практически невозможно. Это обусловлено сложностью пространственного расположения деталей и наличием многих масштабов. В таких случаях используется процедура упрощения объекта, которая позволяет построить ряд моделей для изучения влияния ударных параметров на основные детали исследуемого объекта [12]. Упрощение заключается в том, что внутри реакторной зоны проведено осреднение материалов мелкомасштабных деталей по аддитивной теории смесей [13] (основные материалы: бериллий, диоксид урана и гидрид циркония составляют порядка 95-97% массы объекта). Далее считаем, что при входе в плотные слои атмосферы внешние элементы конструкции сгорают и от реактора остается объект со сложным внутренним строением, показанным на рис. 4 (плоская конструкция объекта включает 82 тела и 81 контактную границу). В этом случае модель блока реактора имеет вид кольцевой структуры, т. е. счетная область обладает многосвязностью и обилием контактных поверхностей.
Ш\ I I I I I I I
-0.001 0.001 0.003 0.005 0.007 0.009
Рис. 5. Соударение модели реактора с гранитной плитой
При численном моделировании процессов удара расчет контактного взаимодействия — очень важный момент, поскольку точность расчета контактных границ в значительной мере определяет точность решения всей задачи в целом. Сложность реализации граничных условий на контактной поверхности заключается в том, что необходимо удовлетворить как кинематические, так и динамические условия на поверхности, положение которой также подлежит определению.
В работе [14] на основе принципа локальности действия сил реакции построен симметричный численный алгоритм расчета контактных границ при динамическом взаимодействии деформаций с учетом трения. Принцип локальности состоит в том, что влияние граничного узла одного тела распространяется на узел или ребро ячейки другого тела. Алгоритм реализовывался в два этапа. На первом этапе вычислялись предварительные значения скоростей и координат граничных узлов с учетом граничных условий за исключением поверхности контакта, где предполагается отсутствие сил реакции. На втором — определялись номера граничных узлов, которые пересекли предварительную поверхность другого тела. И, применяя принцип локальности действия сил реакции, вы-
Рис. 6. Кинограмма расчета бокового удара модельного реактора о поверхность гранитной плиты со скоростью 100 м/с
числяли значения этих сил в узлах, расположенных на поверхности контакта тел. При этом сила реакции включала нормальную составляющую и силу трения. Применение этих сил при коррекции скоростей в граничных узлах позволяла находить истинное положение контактной поверхности между телами. Однако большинство граничных узлов участвует в нескольких элементарных актах взаимодействия, поэтому в алгоритм был введен учет количества элементарных актов взаимодействия для каждого узла. Силы реакции от каждого взаимодействия данного узла суммируются, и вычисляется средняя сила реакции, которая используется при коррекции векторов скорости и положения границ. Данная модификация алгоритма позволила проводить расчеты контактных границ даже с самопересечением, т. е. взаимодействием узлов одной границы.
Математическая постановка задачи, а также результаты расчета продольного соударения ЯЭУ с поверхностью Земли подробно рассмотрены в [12], а на рис. 5 показана кинограмма процесса бокового удара модели реактора по поверхности гранитной плиты для демонстрации успешной работы модифицированного алгоритма расчета контактных границ без учета разрушения в многосвязной области.
На рис. 6 приведена кинограмма процесса расчета бокового удара реактора о гранитную плиту при скорости встречи 100 м/с. Данный расчет показывает, что боковой удар реактора более критичен к скорости встречи, нежели торцевой [12], так как имеют место большие деформации в области контакта. Движение более тяжелого наполнения реактора давит на бериллиевый корпус, вызывая сдвиговые разрушения в области контакта с гранитной плитой и разрыв материала корпуса с внутренней стороны, что приводит к нарушению герметичности реактора. Заливка из ZrH значительно растрескивается, что приводит к разрушению твэлов и как следствие — к радиоактивному заражению места падения реактора.
Выводы
Предложен простой алгоритм применения динамического метода построения треугольных сеток в многосвязных областях.
Модифицирован симметричный алгоритм расчета поверхностей контакта и показана его работоспособность при моделировании сложных объектов с многочисленными контактами.
Список литературы
[1] Умлнский С.Э. Алгоритм и программа триангуляции двумерной области произвольной формы // Проблемы прочности. 1978. № 6. С. 83-87.
[2] Попов И.В., Поляков С.В. Построение адаптивных нерегулярных треугольных сеток для двумерных многосвязных невыпуклых областей // Матем. моделирование. 2002. Т. 14, № 6. С. 25-35.
[3] Shewchuk J.R. Adaptive precision floating-point arithmetic and fast robust geometric predicates // Discrete & Computational Geometry. 1997. Vol. 18, N 3. P. 305-363.
[4] Ho-Le K. Finite element mesh generation methods: a review and classification // Computer Aided Design. 1988. Vol. 20, N 1. P. 27-38.
[5] Shephard M.S. Trends in automatic three-dimensional mesh generation // Computers & Structures. 1988. Vol. 30, N 12. P. 421-429.
[6] Thacker W.C. A brief review of techniques for generating irregular computational grids // Intern. J. for Numerical Methods in Engineering. 1980. Vol. 15. P. 1335-1341.
[7] Fang T.P., Piegl L.A. Algorithm for constrained Delaunay triangulation // The Visual Computer. 1994. Vol. 10. P. 255-265.
[8] Скворцов А.В. Триангуляция Делоне и ее применение. Томск: Изд-во ТГУ, 2002. 128 с.
[9] Dey Т.К., Bajaj C.L. On good triangulations in three dimensions // Intern. J. of Computational Geometry & Applications. 1992. Vol. 2, N 1. P. 75-95.
[10] Shimada K., Gossard D. Bubble Mesh: Automated triangular meshing of non-manifold geometry by sphere packing // ACM Third Symposium on Solid Modeling and Applications. 1995. P. 409-419.
[11] Ogata K. Modern Control Engineering. N.Y.: Prentice-Hall, 1970. 836 p.
[12] Краус Е.И., Фомин В.М., Шабалин И.И. Моделирование процесса соударения сложных двумерных тел о деформируемую преграду // Вычисл. технологии. 2006. Т. 11. Спецвыпуск, посвященный 85-летию со дня рождения академика Н.Н. Яненко. C. 104-107.
[13] Дремин А.Н., Карпухин И.А. Метод определения ударных адиабат дисперсных веществ // ПМТФ. 1960. № 3. С. 184-188.
[14] Фомин В.М., Гулидов А.И., Сапожников Г.А. и др. Высокоскоростное взаимодействие тел. Новосибирск: Изд-во СО РАН, 1999. 600 с.
Поступила в редакцию 26 января 2009 г.
