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

Ссылка на статью:

// Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. №6. С. 329-345.

Б01: 10.7463/0615.0779350

Представлена в редакцию: 31.05.2015 © МГТУ им. Н.Э. Баумана

УДК 532.5-1

Использование метода конечных элементов с частицами для решения задач гидродинамики

Давыдова Е. В.1, Корчагова В. Н.1, * [email protected]

л л

Марчевский И. К.1'

1МГТУ им. Н.Э. Баумана, Москва, Россия

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

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

Введение

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

Эти трудности снимаются при использовании лагранжева описания среды. Суть такого подхода заключается в том, что узлы сетки движутся вместе со средой, что позволяет рассматривать их как частицы среды; при этом сама сетка деформируется либо перестраивается на каждом шаге решения. Одним из методов, использующих лагранжево описание среды, является РБЕМ — метод конечных элементов с частицами [1,2].

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

Наука и Образование

МГТУ им. Н.Э. Баумана

Сетевое научное издание

ISSN 1994-0408

указанных задач лагранжевы методы различных типов используются традиционно и весьма эффективно: в сопряженных задачах гидроупругости — методы вихревых элементов [3, 4], при моделировании течений со свободной поверхностью — метод сглаженных частиц БРИ [5, 6]. Преимущества и недостатки использования сеточных методов и методов частиц применительно к решению различных задач подробно рассмотрены в работе [7].

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

1. Общее описание алгоритма РЕЕМ

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

Будем считать, что расчетная область лежит в плоскости Оху. Пусть хк — положение частицы в момент времени Ьк, хк+1 — положение частицы в момент времени ¿к+1 = Ьк + ДЬ. Скорость частицы в момент Ьк считается известной: ик = и(хк,Ьк).

Тогда новое положение частицы можно определить по формуле

хк+1 = хк + ик д., (!)

представляющей собой приближенное решение уравнения X = и на одном шаге по времени явным методом Эйлера.

Приведем основные этапы решения задачи с помощью метода РББМ (блок-схема представлена на рис. 1).

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

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

3. Выполняется построение конечно-элементной сетки в каждой подобласти таким образом, что частицы образуют множество ее узлов.

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

Рис. 1. Блок-схема РББМ

5. Происходит перемещение частиц по формуле (1), все их характеристики (плотность, давление) остаются неизменными.

6. Выполняется переход к следующему шагу.

2. Метод возможных треугольников: построение треугольной сетки на заданном наборе узлов

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

Алгоритм, называемый методом возможных треугольников и основанный на алгоритме Форчуна [8], позволяет для заданного набора точек плоскости 5 = {вь в2,..., вп} построить триангуляцию, удовлетворяющую критерию Делоне (т. е. если вокруг каждой треугольной ячейки сетки с узлами в^, в^, вк описать окружность, то в нее не попадет ни одна точка из набора 5 \ {в^, в^, вк}).

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

Заметающая прямая — это вертикальная прямая, движущаяся от левой границы области к правой и «сканирующая» положение узлов. Узлы, лежащие по левую сторону от заметающей прямой, уже соединены ребрами с другими узлами.

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

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

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

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

Построение сетки на заданном наборе точек 5 = {в1, в2, ..., вп}, предварительно отсортированном по возрастанию их абсцисс, осуществляется по следующему алгоритму.

1. На г-м шаге алгоритма помещаем (перемещаем) заметающую прямую в точку вг.

2. Проверяем положение «сгшкроШ»' ов уже имеющихся возможных треугольников относительно заметающей прямой. Если необходимому условию выполнения критерия Делоне удовлетворяют сразу несколько возможных треугольников, то нужно «замкнуть» тот треугольник, «сгткро1М» которого находится левее всех остальных. Список возможных треугольников, таким образом, обновится: «замкнутый» треугольник перейдет в список построенных, а новое построенное ребро может породить новые возможные треугольники. Соответственно, такая процедура проверки повторяется до тех пор, пока все «сгткро1М»' ы имеющихся возможных треугольников не окажутся справа от заметающей прямой.

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

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

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

д

Рис. 2. Графическая демонстрация выполнения основных операций алгоритма метода возможных треугольников

Об использовании обобщения триангуляции Делоне. В процессе тестирования алгоритма РБЕМ на треугольных сетках обнаруживаются трудности, связанные с классической триангуляцией: возможно появление неоднозначности описанного выше алгоритма в случае, когда несколько узлов лежат на одной окружности либо близки к этому (рис. 3). Для реализации РБЕМ в этом случае его разработчики (Е. Опа1е, Б.К Ыекоп) [1,9] рекомендуют использовать обобщение триангуляции Делоне, которое позволяет исключить неоднозначность построения сетки (см. рис. 3). При построении такой сетки появляются многоугольные ячейки, что требует использования специальных функций формы для соответствующих конечных элементов [10, 11].

а б в

Рис. 3. Варианты построения сетки на узлах, лежащих на одной окружности: а, б — равноправные варианты для классической триангуляции Делоне; в — четырехугольная ячейка, полученная при обобщении триангуляции Делоне

3. Аппроксимация уравнений гидродинамики

Общая постановка задачи. Рассмотрим плоское течение вязкой несжимаемой жидкости, занимающей область П. Плотность жидкости считается постоянной и равной р, коэффициент динамической вязкости равен При этом жидкость может находиться под воздействием внешних объемных сил, объемная плотность которых равна Ь = {61, 62}.

Движение такой жидкости описывается уравнением равновесия и уравнениями Навье — Стокса:

дХ=°; (2)

Ощ др д2щ .

рЖ = -дХ + ^дХ"дХ" + рк г = *,2 (3)

Бщ дщ дщ

Здесь —- = ——+ —--материальная (субстанциональная) производная.

иЬ дЬ дх^

Пусть п = {п1,п2} — единичный вектор нормали к границе области течения, направленный внутрь жидкости; т^ — компоненты тензора вязких напряжений в жидкости. Тогда граничные условия можно записать следующим образом:

1) т^щ — рщ = ~оП1, если на участке границы задано нормальное напряжение (давление); равенство аПг = 0 обеспечивает условие свободной поверхности;

2) иг = и — условие прилипания, если на участке границы задана скорость ее движения; аналогичное условие ставится на входной границе расчетной области, если набегающий поток на ней считается невозмущенным;

3) ^т^ = 0 на выходной границе области; дп

др

4) — = 0 на границе тела, обтекаемого жидкостью. дп

В начальный момент времени £ = 0 примем жидкость покоящейся. Таким образом, получаем начальное условие:

иг(0, х) = 0, X е П.

Дискретизация по времени: алгоритм дробных шагов. Уравнение (3) записано в лагранжевой постановке. Аппроксимируем полную производную по времени, для правой части уравнения используем полностью неявную схему:

£иг _ иг(хк+1,^+1) - иг(хк,£к) ик+1 - ик 1 дрк+1 , - д2ик+1

Д£ Д£ р дхг + р дxj дхj + ^ (4)

При интегрировании по времени уравнений (2), (3) возникают трудности, когда жидкость несжимаема или почти несжимаема: возникают нефизичные осцилляции давления, которые необходимо подавлять. Чтобы преодолеть эти трудности, удобно использовать алгоритм дробных шагов [12].

Суть метода заключается в расщеплении каждого шага по времени на две части. При этом вводится величина и* — прогноз скорости. Уравнения (2), (3) разбиваются на 3 уравнения так, чтобы процесс решения на каждом шаге обеспечивал бездивергентность поля скоростей, т. е. выполнялось уравнение (2):

1) находим прогноз скорости и*, решив уравнение

Т = + 1 =1' 2- (5)

Д£ р дxj дxj

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

2) находим значение давления, решив уравнение Пуассона с учетом граничных условий

д2рк+1 р ди*

дхгдхг Д£ дхг' 3) находим скорость ик+1, учитывая влияние давления

ик+ - и* 1 дрк+1

(6)

Д£ р дхг

(7)

Дискретизация по пространству: метод конечных элементов. Неизвестные функции скорости и давления представляются в виде линейных комбинаций функций формы и

значений скоростей и давления в частицах следующим образом:

п п

пг = £ N(х)(иг)1 = [Щ{щ}, р = £ N1 (х)(р) = [М][р], 1=1 1=1

где п — количество частиц.

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

Для получения систем линейных алгебраических уравнений, аппроксимирующих задачу (5)-(7), используем метод Бубнова — Галеркина: системы формируются из условия ортогональности невязки Я соответствующего уравнения проекционным функциям N(х1, х2), в качестве которых выступают функции формы.

Дискретные аналоги уравнений (5)-(7) можно записать в виде:

([М + ^Д и) {<} = [М]К} + Д^}, (8)

трк+1} = Д ([в]т {«*}), (9)

[Ы]{пк+1} = [М]{м*} - Д [В]{рк+1}, (10)

Р

где

[М ] =

[ММ] 0

о [М

[В] = {[В1], [В2]}т , {^} = {[Л], [^2]}т, г = 1, 2.

Матрицы [М], [5], [В1], [В2] и векторы {и}, [^2] получают ансамблированием элементных матриц по всем конечным элементам:

КеУ 1п(б)]

М] = Е^Т (е)]Т[N(б)] [а( е Че) '

П(е)

[Вг] = Е^ ( / [N(е)] ^^(в))[«(е)];

6 П(е)

= Е[а(е)]^ / Ь ¿Ц + В«(6)]Т( /[^Т^1 , г = 1, 2.

6 П(е) 6 Г(е)

Здесь [N(б)]

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

Отметим, что граничные условия, задающие компоненты скорости в граничных узлах, вводятся в матрицы систем (8) и (10) после ансамблирования напрямую (уравнение, соответ-

ствующее граничному узлу, заменяется на нужное соотношение). Число обусловленности матриц при этом невелико (сопоставимо с их размером). Решение систем проводится методом BiCGStab с iLU-предобуславливанием [13, 14], поэтому несимметричность матриц не является критичной.

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

Выбор функций формы: <<non-Sibsonian interpolation». На каждом шаге расчета по времени происходит изменение положений частиц и перестроение сетки, поэтому функции формы, используемые в PFEM, также должны соотвествующим образом изменяться. Ячейки сетки, получаемой в результате работы обобщенного алгоритма триангуляции Делоне, не обязательно будут треугольными. По этой причине необходимо применять функции формы, которые можно построить для многоугольных ячеек сетки с произвольным количеством узлов в ячейке по единому алгоритму.

Удобно применять функции формы бессеточного метода конечных элементов (MFEM), получаемые с помощью так называемой «поп^Ьвошатьинтерполяции [10]. Их построение основано на определении соседства узлов конечного элемента и его внутренних точек путем составления диаграммы Вороного.

Пусть P = [п\,... ,nm} — множество узлов, образующих вершины многоугольника (конечного элемента), x — внутренняя точка этого многоугольника (рис. 4, а). Вычисление значений функции формы Ni(x) узла п в точке x происходит в несколько этапов:

• строится ячейка Вороного для точки x на множестве P U {x};

• вычисляются значения функций

sj(x)

Уз =

hj(x)

, j = 1,

, m,

где — длина стороны ячейки Вороного, разделяющей точки х и щ; Н^ — расстояние от х до щ;

а б

Рис. 4. К построению функций формы: а — ячейка Вороного; б — вид одной из функций формы на четырехугольных конечных элементах

вычисляются значения функции формы:

Уг(х)

Щх)

£ Уз(х) 3=1

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

При сборке матриц систем линейных алгебраических уравнений (8)-(10) необходимо вычислять интегралы от произведения значений функций формы и их производных. Целесообразно считать эти интегралы численно, используя квадратурную формулу Гаусса для расчета интеграла функции по треугольнику [15]:

I/(х) ¿Б « Бд ]г ир /(хр).

р=1

Здесь БД — площадь треугольника, ир = 1/3 — весовой множитель, хр — гауссовы точки, ¿-координаты которых получены из последовательности {2/3,1/6,1/6} циклической перестановкой ее элементов.

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

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

4. Численный эксперимент: задача о течении жидкости в квадратной каверне

Рассмотрим модельную плоскую задачу о течении вязкой несжимаемой жидкости в квадратной каверне П, длина стороны которой равна Ь (рис. 5). Верхняя граница каверны Г1 движется с постоянной скоростью и0. Плотность жидкости считается постоянной и равной р, коэффициент динамической вязкости равен р.

При построении математической модели течения жидкости использовались следующие безразмерные величины (знаком «~» обозначены размерные величины):

• и = й/и0 — скорость течения жидкости;

Р

• р = ——к — давление;

ри.о

• хз = Х3 /Ь, ] = 1, 2, — безразмерные координаты;

• Ь = ¿/¿о — безразмерное время;

• Яе = — число Рейнольдса.

Рис. 5. Расчетная схема задачи о каверне

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

V- и = 0;

—и „ 1 Л

— = + - Аи;

(11) (12)

и

Г1

(1;0)

т

и

Г2

(0;0)т, ^ д п

Г1УГ2

где = + (и ' V)u — полная (материальная) производная по времени.

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

Расчеты были проведены при АЬ = 0.01, К = 0.02, Яе = 400, где АЬ и К — шаги по времени и пространству соответственно. Диаграмма распределения скоростей на первом шаге по времени, полученная при первоначальной регулярной четырехугольной сетке, приведена на рис. 6.

Расчеты показали, что имеет место линейная сходимость решения. На рис. 7 приведены графики, демонстрирующие зависимость равномерной нормы разности решений е от количества частиц п. Значения компонент скорости щ и и2, полученные в процессе решения, сравнивались со значениями, полученными на сетке, состоящей из 3136 частиц, в одних и тех же точках. Из приведенных графиков видно, что при достаточно большом количестве частиц (больше 1000) норма разности решений невелика — имеет место численная сходимость.

0

Рис. 6. Диаграмма распределения скоростей при Ь = 0.01

б

Рис. 7. Разность решений е в зависимости от количества частиц п: а — для регулярной четырехугольной сетки; б — для треугольной сетки

а

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

Заключение

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

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

Список литературы

1. Idelsohn S.R., Onate E., Del Pin F. A Lagrangian meshless finite element method applied to fluid-structure interaction problems // Computer and Structures. 2003. Vol.81, no. 8-11. P. 655-671. DOI: 10.1016/S0045-7949(02)00477-7

2. Idelsohn S.R., Onate E., Del Pin F. The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves // International Journal for Numerical Methods in Engineering. 2004. Vol. 61, no. 7. P. 964-989. DOI: 10.1002/nme.1096

3. Ермаков A.B., Щеглов Г.А. Моделирование методом вихревых элементов динамики цилиндрической оболочки в пространственном потоке жидкости // Известия ВУЗов. Машиностроение. 2014. №3. C. 35-41.

4. Кузьмина К.С., Марчевский И.К. Об ускорении вычислений при решении двумерных сопряженных задач гидроупругости вихревыми методами // Вестник Пермского национального исследовательского политехнического университета. Сер. Аэрокосмическая техника. 2014. №39. С. 145-163.

5. Monaghan J.J. Smoothed Particle Hydrodynamics // Reports on Progress and Physics. 2005. Vol. 68. P. 1703-1759. DOI: 10.1088/0034-4885/68/8/R01

6. Li Shaofan, Kam Wing Liu. Meshfree Particle Methods. Springer Berlin Heidelberg, 2004. 502 p. DOI: 10.1007/978-3-540-71471-2

7. Idelsohn S.R., Onate E. To mesh or not to mesh. That is the question... // Computer Methods in Applied Mechanics and Engineering. 2006. Vol.195, iss. 37-40. P. 4681-4696. DOI: 10.1016/j.cma.2005.11.006

8. deBergM., vanKreveldM., OvermarsM., Schwarzkopf O. Computational Geometry. Springer Berlin Heidelberg, 2000. 379 p. DOI: 10.1007/978-3-662-04245-8

9. Onate E., Idelsohn S.R., Celigueta M.A., Rossi R. Advances in the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows // Computer Methods in Applied Mechanics and Engineering. 2008. Vol. 197, iss. 19-20. P. 1777-1800. DOI: 10.1016/j.cma.2007.06.005

10. Idelsohn S.R., Onate E., Calvo N., Del Pin F. The meshless finite element method // International Journal for Numerical Methods in Engineering. 2003. Vol.58, no. 6. P. 893912. DOI: 10.1002/nme.798

11. Kraus M., Rajagopal A., Steinmann P. Investigations on the polygonal finite element method: Constrained adaptive Delaunay tessellation and conformal interpolants // Computers and Structures. 2013. Vol. 120. P. 33-46. DOI: 10.1016/j.compstruc.2013.01.017

12. ЯненкоН.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967. 196 с.

13. Марчевский И.К., Пузикова В.В. Анализ эффективности итерационных методов решения систем линейных алгебраических уравнений // Математическое моделирование и численные методы. 2014. №4. С. 37-52.

14. Saad Y. Iterative Methods for Sparse Linear Systems. 2nd edition. Philadelphia, Pennsylvania: SIAM, 2003. 547 p.

15. Галанин М.П., Савенков Е.Б. Методы численного анализа математических моделей. М.: Изд-во МГТУ им. Н.Э. Баумана, 2010. 590 с.

Science ¿¿Education

of the Bauman MSTU

Science and Education of the Bauman MSTU, 2015, no. 6, pp. 329-345.

DOI: 10.7463/0615.0779350

Received:

31.05.2015

Electronic journal

© Bauman Moscow State Technical University

ISSN 1994-0408

On Using Particle Finite Element for Hydrodynamics Problems Solving

Davidova E.V.1, Korchagova V.N.1, Marchevskii I.K.1'* * [email protected]

1 Bauman Moscow State Technical University, Russia

Keywords: hydrodynamics, Lagrangian description, particle finite element method, mesh reconstrution

The aim of the present research is to develop software for the Particle Finite Element Method (PFEM) and its verification on the model problem of viscous incompressible flow simulation in a square cavity. The Lagrangian description of the medium motion is used: the nodes of the finite element mesh move together with the fluid that allows to consider them as particles of the medium. Mesh cells deform when in time-stepping procedure, so it is necessary to reconstruct the mesh to provide stability of the finite element numerical procedure.

Meshing algorithm allows us to obtain the mesh, which satisfies the Delaunay criteria: it is called "the possible triangles method". This algorithm is based on the well-known Fortune method of Voronoi diagram constructing for a certain set of points in the plane. The graphical representation of the possible triangles method is shown. It is suitable to use generalization of Delaunay triangulation in order to construct meshes with polygonal cells in case of multiple nodes close to be lying on the same circle.

The viscous incompressible fluid flow is described by the Navier — Stokes equations and the mass conservation equation with certain initial and boundary conditions. A fractional steps method, which allows us to avoid non-physical oscillations of the pressure, provides the time-stepping procedure. Using the finite element discretization and the Bubnov — Galerkin method allows us to carry out spatial discretization.

For form functions calculation of finite element mesh with polygonal cells, "non-Sibsonian interpolation" is used, which allows to construct and reconstruct the functions with desired properties for finite elements with arbitrary number of nodes.

The model problem of the flow simulation of viscous incompressible fluid in a square cavity is considered. The diagram of the velocity distribution after the first time step is shown and the numerical convergence of solutions for different initial meshes (regular quadrilateral and triangular) is discussed.

As the conclusion, the considered algorithm is useful for solving a number of problems in the domains with complex geometry, which can be either broken into several domains or merged. A Lagrangian approach us allows to take into account the convective terms of the governing equations automatically. However, the procedure of the mesh reconstruction has large computational cost and time step have to be reduced significantly if particles are located close one to each other. So, there are some possibilities of algorithm improvement.

References

1. Idelsohn S.R., Onate E., Del Pin F. A Lagrangian meshless finite element method applied to fluidstructure interaction problems. Computer and Structures, 2003, vol.81, no. 8-11, pp. 655-671. DOI: 10.1016/S0045-7949(02)00477-7

2. Idelsohn S.R., Onate E., Del Pin F. The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves. International Journal for Numerical Methods in Engineering, 2004, vol. 61, no. 7, pp. 964-989. DOI: 10.1002/nme.1096

3. Ermakov A.V., Shcheglov G.A. The application of the three-dimensional vortex element method to the fluid dynamic analysis of cylindrical shell elements. Izvestiia vysshikh ucheb-nykh zavedenii. Mashinostroenie = Proceedings of Higher Educational Institutions. Machine Building, 2014, no. 3, pp. 35-41. (in Russian).

4. Kuz'mina K.S., Marchevskii I.K. On computation speeding up when solving two-dimensional hydroelastic coupled problems by using vortex methods. Vestnik Permskogo natsional'nogo issledovatel'skogo politekhnicheskogo universiteta. Ser. Aerokosmicheskaya tekhnika = PN-RPU Aerospace Engineering Bulletin, 2014, no. 39, pp. 145-163. (in Russian).

5. Monaghan J.J. Smoothed Particle Hydrodynamics. Reports on Progress and Physics, 2005, vol. 68, pp. 1703-1759. DOI: 10.1088/0034-4885/68/8/R01

6. Li Shaofan, Kam Wing Liu. Meshfree Particle Methods. Springer Berlin Heidelberg, 2004. 502 p. DOI: 10.1007/978-3-540-71471-2

7. Idelsohn S.R., Onate E. To mesh or not to mesh. That is the question. .. Computer Methods in Applied Mechanics and Engineering, 2006, vol. 195, iss. 37-40, pp. 4681-4696. DOI: 10.1016/j.cma.2005.11.006

8. deBergM., vanKreveldM., OvermarsM., Schwarzkopf O. Computational Geometry. Springer Berlin Heidelberg, 2000. 379 p. DOI: 10.1007/978-3-662-04245-8

9. Onate E., Idelsohn S.R., Celigueta M.A., Rossi R. Advances in the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows. Computer Methods in Applied Mechanics and Engineering, 2008, vol. 197, iss. 19-20, pp. 1777-1800. DOI: 10.1016/j.cma.2007.06.005

10. Idelsohn S.R., Onate E., CalvoN., DelPinF. The meshless finite element method. International Journal for Numerical Methods in Engineering, 2003, vol.58, no. 6, pp. 893-912. DOI: 10.1002/nme.798

11. Kraus M., Rajagopal A., Steinmann P. Investigations on the polygonal finite element method: Constrained adaptive Delaunay tessellation and conformal interpolants. Computers and Structures, 2013, vol. 120, pp. 33-46. DOI: 10.1016/j.compstruc.2013.01.017

12. Yanenko N.N. The Method of Fractional Steps. The Solution of Problems of Mathematical Physics in Several Variables. Springer Berlin Heidelberg, 1971. 160 p. DOI: 10.1007/978-3-642-65108-3 (Russ. ed.: Yanenko N.N. Metod drobnykh shagov resheniya mnogomernykh zadach matematicheskoi fiziki. Novosibirsk, Nauka Publ., 1967. 196 p.).

13. Marchevskii I.K., Puzikova V.V. Performance analysis of iterative methods of combined linear algebraic equations solution. Matematicheskoe modelirovanie i chislennye metody, 2014, no. 4, pp. 37-52. (in Russian).

14. Saad Y. Iterative Methods for Sparse Linear Systems. 2nd edition. Philadelphia, Pennsylvania: SIAM, 2003. 547 p.

15. Galanin M.P., Savenkov E.B. Metody chislennogo analiza matematicheskikh modelei [Methods of numerical analysis of mathematical models]. Moscow, Bauman MSTU Publ., 2010. 590 p. (in Russian).