CC BY
40
Моделирование осесимметричных течений вязкой несжимаемой жидкости методом конечных элементов с частицами PFEM-2 в программном комплексе Kratos с открытым кодом
1 Е.В. Смирнова <alena.davidova@ispras.ru> 1,2 И.К. Марчевский <iliamarchevsky@mail.ru> 2 В.О. Бондарчук <viktoriia.bondarchuk@gmail.com>
1 Институт системного программирования им. В.П. Иванникова РАН, 109004, Россия, г. Москва, ул. А. Солженицына, д. 25. 2Московский государственный технический университет им. Н.Э. Баумана, 105005, Россия, Москва, ул. 2-я Бауманская, д. 5
Аннотация. Получены необходимые соотношения для реализации алгоритма моделирования осесимметричных течений вязкой несжимаемой жидкости в методе конечных элементов с частицами PFEM-2. Осесимметричная модель реализована в программном комплексе с открытым исходным кодом Kratos на основе существующего модуля для решения плоских задач. Была проведена валидация осесимметричной модели на тестовых задачах. В качестве модельных рассмотрены задачи о течении в трубе (задача Пуазейля) и задача о моделировании падения капли в глубокий слой жидкости. Численное решение, полученное методом конечных элементов с частицами, для задачи Пуазейля показало удовлетворительное согласие с аналитическим; решение задачи о падении капли в слой глубокий жидкости в комплексе Kratos сравнивалось с результатом моделирования в программном пакете Gerris с открытым исходным кодом. Результаты расчетов также удовлетворительно согласуются.
Ключевые слова: метод конечных элементов с частицами; вязкая несжимаемая среда; осесимметричное течение; задача Пуазейля; падение капли в слой жидкости.
Б01: 10Л5514/ISPRAS-2018-30(2)-13
Для цитирования: Смирнова Е.В., Марчевский И.К., Бондарчук В.О. Моделирование осесимметричных течений вязкой несжимаемой жидкости методом конечных элементов с частицами PFEM-2 в программном комплексе Kratos с открытым кодом. Труды ИСП РАН, том 30, вып. 2, 2018 г., стр. 263-284. DOI: 10Л5514ДSPRAS-2018-30(2)-13
1. Введение
Метод конечных элементов с частицами (Particle Finite Element Method, PFEM [1]) представляет собой сеточный эйлерово-лагранжев метод вычислительной гидродинамики несжимаемой среды. Наиболее эффективным является его модификация PFEM-2 [2], в рамках которой расчет ведется на неподвижной конечноэлементной сетке, через которую перемещаются частицы - «носители» скорости среды. Движение частиц позволяет моделировать движение среды, обусловленное конвективным переносом, а влияние остальных процессов (вязкость, градиент давления) учитывается путем решения соответствующих дифференциальных уравнений методом конечных элементов.
На сегодняшний день единственным программным комплексом, в котором реализованы PFEM/PFEM-2, является пакет Kratos c открытым исходным кодом (http://kratos-wiki.cimne.upc.edu). Kratos - конечноэлементный пакет, позволяющий решать задачи теории упругости, механики разрушения, гидродинамики и т.д. В текущей версии Kratos можно моделировать плоские и трехмерные течения жидкостей, однако возможность решения задач в осесимметричной постановке до настоящего времени не реализована. Целью данной работы является разработка алгоритма моделирования несжимаемых осесимметричных течений методом конечных элементов с частицами и его реализация в рамках программного комплекса Kratos.
2. Постановка задачи
2.1 Постановка и принцип решения задачи методом PFEM-2
Движение вязкой несжимаемой жидкости описывается уравнениями Навье -
Стокса и уравнением неразрывности [3]
+^'=v ■(мУ*) ~Vp+р;
V-v = 0,
где р - плотность, v = (u,v,w) - скорость среды; ц - коэффициент динамической вязкости; p - давление; g = (gj,g2,S3) - вектор массовой плотности объемных сил. Система уравнений дополняется необходимыми начальными и граничными условиями, чаще всего в их роли выступают:
• условие непротекания: v и|р = 0,
• условие прилипания: v Г|г2 = 0,
• заданные величины для скорости или давления: v |р = v 0, p |р = p 0-
Здесь индексы п и т указывают на нормальные и касательные компоненты вектора скорости соответственно.
Смысл использования лагранжевых методов заключается в возможности учета конвективного слагаемого в уравнении (1) путем моделирования движения частиц - «носителей» скорости - по расчетной области. Вводя понятие
а д _
материальной (субстанциональной) производной — =--+ (у -V),
б! д!
учитывающей изменение соответствующей величины в лагранжевой частице, систему уравнений (1)-(2) можно записать в виде
Р— = V-(^) ^р + р;
а/
V-V = 0.
Полную производную по времени заменим разностным аналогом: /) и /(! + Д!) - /) + 0(Д/
а/ л
С учетом (4) уравнения (3) переходят в систему уравнений (в частных произвоных по пространственным координатам), которая должна решаться на каждом п -м шаге по времени:
V(пъ1) , /14 , V(п)
Р « V - (п+в)) - Vp(n+1)+р+Р у -,
Д Д
V-V (п+1)=0,
Параметр в определяет способ учета «вязкого» слагаемого: в = 0 и в = 1 для явной и неявной аппроксимации соответственно.
2.2 Постановка осесимметричной задачи
Примем, что вязкая несжимаемая жидкость с плотностью р и вязкостью V занимает область цилиндрической формы и находится в поле массовой силы £ , действующей вдоль оси области. Введем систему координат {г, ф, г], направив ось Ог по оси симметрии области, ось Or - перпендикулярно ей. Будем расматривать осесимметричные незакрученные течения среды (Уф =0), когда окружное перемещение частиц равно нулю. Тогда можно перейти к двумерной задаче гидродинамики в области О = {{г,г]: г е {г^Г2],г е {г^г2]] (рис. 1). Течение жидкости описывается уравнениями (3), при этом все дифференциальные операторы подразумеваются действующими в цилиндрической системе координат {г,г]. Граничные условия (ГУ) соответствуют указанным выше, однако ГУ на оси области (при г = 0) является условие ограниченности радиальной скорости.
Рис. 1. Расчетная область Fig. 1. The computational domain
3. Краткое описание метода PFEM-2
Эйлерово-лагранжев метод PFEM-2 [2, 4] предполагает рассмотрение в области течения неподвижной конечноэлементной сетки и набора частиц, движущихся сквозь ячейки сетки. Частицы перемещаются по полю скоростей, одновременно являясь его «носителями», что обеспечивает учет конвективного слагаемого из уравнения (2). При этом система уравнений (5) решается на неподвижной сетке методом конечных элементов [5]. Отметим, что при моделировании многофазных течений частицы также являются маркерами соответствующей фазы. Расчет течения методом PFEM-2 производится по следующему алгоритму (операции 2-5 выполняются в цикле по времени).
1) В расчетной области строится конечноэлементная сетка. Задается начальное распределение частиц, при котором на каждую ячейку сетки приходится, как правило, несколько десятков частиц.
2) По известному полю скоростей с n -го шага по времени частицы переносятся вдоль линий тока в новые положения.
3) Скорости частиц проецируются (интерполируются) на узлы конечноэлементной сетки.
4) Решается система уравнений (5); найденные ранее скорости являются начальным приближением для значений скоростей на следующем шаге.
5) По найденному полю скоростей в узлах сетки корректируются скорости частиц. При необходимости некоторые частицы исключаются из расчета, а в ячейки с малым числом частиц вводятся дополнительные частицы.
Алгоритмы перечисленных операций с особенностями их реализации описаны в публикациях, посвященных PFEM-2. В данной работе приведем лишь их краткое описание, необходимое для понимания общей идеологии метода PFEM-2. При этом остановимся более подробно на решении методом конечных элементов системы уравнений (3). Для простоты будем рассматривать задачу в двумерной (плоской) постановке и считать, что массовые силы направлены вдоль второй координаты, т.е. вектор g имеет координаты (0, g).
3.1 Перемещение частиц
Рассмотрим частицу жидкости с номером р. Зная координаты частицы в некоторый момент времени Уп, ее положение в момент времени ¿и+1 > Уп можно определить выражением
. . Уп+1
г(рй+1)= хрп)+ | Vt(xtp)dt,
Уп
где V - скорость, а индекс У означает, что вычисление должно производиться непрерывно по времени. Скорость частицы в момент времени Уп+1 определяется аналогичным образом:
^ v(pn) + | (
где а - ускорение. Однако уравнения (6)-(7) не могут быть решены напрямую, так как перемещение частицы зависит от ее скорости в каждый момент времени, а скорость - от координаты частицы. К тому же для решения уравнений производится дискретизация по времени, поэтому значения величин известны только в определенные моменты времени Уп_1, Уп, Уп+1 и т.д. Для приближенного расчета положения частицы р в момент времени Уп+1 из формулы (6) получаем:
хрп+1)« Х(рп)+ | Vп(хУр^.
Уп
Такой подход фактически означает интегрирование движения частицы вдоль линии тока (т.е. в постоянном поле скоростей). Скорость частицы может быть найдена как
Уп+1 Уп+1
v(pn+1)«vpn)+ (1 _в) | а(п)(хр^у+в | а^(хр^, о<в< 1.
Уп Уп
Перемещение частиц может осуществляться двумя способами. Первый, более простой способ заключается в аппроксимации интегралов в правой части выражения (9) следующим образом (принимаем в простейшем случае в = 0 и
дополнительно полагаем, что для ускорения справедливо а(п) « а (п+1)):
v(pn+1) и ^+а (п)( х(рп))ду,
где Д = Уп+1 _ Уп. Затем согласно (8) находится положение частицы выражению. Однако при использовании данного способа положение частицы
в момент времени tn+i может существенно отличаться от истинного. Поэтому для повышения точности на этапе перемещения частиц применяют второй способ - интегрирование вдоль линий тока.
Данный способ предполагает перемещение частицы по известному на шаге tn полю скорости. При этом отрезок \tn ,tn+i ] дробится на M более мелких:
L о л л M— M I
|f n , tn , tn ,•■■, tn , tn }
где tm = tn + mx, x = At/M. Перемещение частицы вдоль линий тока происходит согласно выражению
x Pm+1}= x Pm)+ v Pn)(X (pm))x, 0 < m <M-1.
Иллюстрация интегрирования вдоль линий тока приведена на рис. 2. Стоит отметить, что скорость частиц, найденная таким образом, позже будет скорректирована.
Рис. 2. Интегрирование вдоль линий тока Fig. 2. Integration along the streamlines
3.2 Проецирование скоростей с частиц на сетку
Для описания механизма проецирования величин с частиц на узлы заданной в расчетной области сетки рассмотрим скалярное поле Л. С каждой г -й частицей ассоциировано значение ^ . Для того, чтобы вычислить значение Л у для каждого у -го узла сетки, необходимо использовать все частицы в соседних к узлу у элементах (рис. 3). Полагается, что частицы, расположенные ближе к узлу у , дают больший вклад в значение Л у. Для
выполнения данного требования вводится весовая функция, с помощью которой будет проводиться вычисление значения поля в узлах. В качестве
весовой функции выбрана стандартная функция формы элемента Л1 (х); при этом вес вклада в значение Л у от частицы г равен N (х^ ) = Л/ .
Рис. 3. Узел с индексом J и соседние к нему конечные элементы с частицами Fig. 3. Node with the index J and the neighbouring finite elements with the particles
Таким образом, значение Я в узле j определяется выражением
n
TAN/ л j = "-,
J n
TNJ
i
где верхний индекс у функции формы указывает на принадлежность функции конкретному узлу конечного элемента. Спроецированная на узлы сетки с помощью (10) скорость определяется аналогичнно выражением
n
TVnNj
V п+1= Л-
у ] п •
ЕЫ/
I
Обозначение (") указывает на то, что вычисленная скорость является первым приближением к скорости на новом шаге по времени. Данная величина уточняется в итерационном процессе решения уравнений движения.
3.3 Решение системы уравнений
При использовании метода конечных элементов (МКЭ) значение некоторой величины Я в точке х вычисляется как сумма произведений значений этой величины в узлах на значения соответствующих узлам функций формы [5]:
к
Я( х,Г ) = (О Ы-1 (х).
j=l
Здесь K - количество узлов в элементе (рассматриваем локальную нумерацию узлов в конечном), Ij - значение величины Я в узле j , Nj (x) -функция формы j -го узла.
Ограничимся рассмотрением симплексных (треугольных в плоском и тетраэдральных в пространственном случае) конечных элементов 1-го порядка с линеными функциями формы, частные производные которых постоянны. Отметим, что именно в таком виде PFEM-2 реализован в KRATOS. Представляя u,v,p в виде (12), уравнения (5), полагая в = 1, принимают вид:
Р fuf+N - V • (^vfuf+N) + VZpn+1Ni = p fu/N,
At i i i At i
P^N-V • (MV^N1) + v|pn+1Ni = pg + ^b>lNl,
At i i i i i i At i i
K , . V^Yyrl+1N1 = 0. i
Выполняя стандартную процедуру МКЭ и умножая уравнения на функции формы NJ и далее интегрируя их по конечному элементу Q , получаем:
\PYyT1NjN1 dQ - JNjV • (M^f+1VNi)dQ + J^pn+1NJVNi dQ = QAt i Q i Q i
P K . . = \PTv"NJNidQ, At QAt i
JPivn+1NjNi dQ - JNJV • (M^vf+1VNi)dQ + J^Pr^VN dQ = QAt i Q i Q i
P K . . ( )
= JpgN] dQ + J^YyfN^1 dQ,
Q QAt i
K
jKVf1 • VN1 )dQ = 0.
Q i
С учетом допущений о постоянстве градиентов функций для интегралов в (14) можно ввести следующие обозначения:
M ч = J NJNidQ;
Q
GiJ = J NJVNidQ = 1 VN'V ;
Q 3
11] = | Ы] У<УЫ* )ёО = |У{Ы7 УЫ1 )ёО-|УЫ7 УЫ1 ёО =
= {Ы^УЫ1 • пёБ- |УЫ7УЫ1 ёО « -[УЫ-'УЫ1 ёО = -УЫ7УЫ1У. 5 О О
(в последнем соотношении пренебрегли поверхностным интегралом.) С учетом введенных обозначений систему (14) удобно записать в матричной (блочной) форме:
[О] {г}"+11 Г^да+-£■ {¥}"1 ДРГ1 Г оМ
[М ] + М[Ь]
[О]
о
здесь {V}" = {м" у" ,ы" у" }т - вектор-столбец из компонент скорости в узлах конечного элемента, {Р}" = {р",р",р"}Т - вектор-столбец из значений
1 т
давления в узлах КЭ, {Я} = — V{0,1,0,1,0,1} - вектор-столбец, компонентами
которого являются нули и третья часть объема КЭ. Все матрицы, входящие в (16), состоят из соответствующих им элементов (15), упорядоченных необходимым образом. Матрицы [М ] и [Ь] - квадратные размером (2"х2"), где " = 3 - количество узлов в конечном элементе, а двойка указывает на количество компонент вектора скорости. Матрицы [О ] и [О ] имеют размер
(2"х") и ("х2") соответственно и состоят из элементов О^ . По сложившейся
в МКЭ традиции будем называть [М ] матрицей масс, [Ь ] - матрицей вязких напряжений, [О ] - матрицей дивергенции, [О ] - матрицей градиента. Для составления глобальной системы линейных алгебраических уравнений необходимо для каждого конечного элемента записать локальную систему уравнений (16). Затем локальные матрицы с помощью специальной процедуры аггрегации «встраиваются» в глобальную матрицу [5]. Полученная в итоге система уравнений решается итерационными методами. Для того, чтобы реализовать возможность расчета осесимметричных течений, на базе готовых модулей приложения р£ет2_аррИса'Ыоп для решения задач гидродинамики методом РРЕМ-2 необходимо модифицировать соответствующие программные модули, отвечающие за расчет коэффициентов соответствующих матриц и интерпретацию результатов решения.
О
О
О
3.4 Коррекция положения частиц
После решения системы уравнений скорости частиц, определенные изначально по значениям старого поля скоростей, должны быть обновлены. Для этого используется поправка, вносимая в скорости частиц [4]:
я"+1= vn+1-v n+1.
Утверждается, что использование для коррекции скоростей частиц величины поправки позволяет уменьшить численную диффузию, особенно на
сетках с разноразмерными ячейками.
3.5 Об алгоритме поиска межфазной границы
Для моделирования течений со свободной поверхностью, а также многофазных течений требуется определение положения границы фаз. Для этого частицам приписывается параметр - маркер фазы, к которой они принадлежат. Положение межфазной границы при этом не должно определяться параметрами отдельных частиц. Например, на рис. 4 показано, что некоторые частицы могут располагаться не с нужной стороны границы раздела фаз.
Рис. 4. Конечные элементы, частицы и положение границы раздела фаз Fig. 4. Finite elements, particles and the interface between phases position
Предполагая, что в расчетной области присутствует две фазы, можно промаркеровать частицы легкой фазы параметром +1, а частицы тяжелой -параметром -1. Тогда для вычисления положения границы раздела фаз можно ввести «псевдо»-функцию уровня ф , изоповерхности ф = 0 которой будут соответствовать межфазная граница. Значения функции ф в узле j конечноэлентной сектки находится по формуле
Zsign/^/
фj =
_ I
т
где суммирование ведется по частицам в соседних к узлу у элементах, а sign1■ определяет знак маркера 1 -й частицы.
4. Решение осесимметричной задачи
4.1. Представление функций формы конечных элементов и дифференциальных операторов в осесимметричной задаче
При рассмотрении плоской задачи функция формы Ы' ' -го узла конечного элемента {',],к} в точке Р(Р^,Р2), лежащей внутри конечного элемента, определялась как отношение площадей треугольников {Р,],к} и {',],к}. В осесимметричном случае такой функции формы соответствует функция, определяемая как отношение объемов фигур, полученных из треугольников {Р,],к} и {',],к} их вращением вокруг оси симметрии 02 (рис. 2). Объем подобной фигуры можно найти точно, но в первом приближении он равен ¥']к = 2лТ'Г1]'£1 Б']к , где S']k - площадь треугольника {',],к}, а гт-
радиальная координата его средней точки. Таким образом, функция формы осесимметричного конечного элемента определяется выражением
N (г,2) = Ур]к/У']к • Если приближенно принять гттк ~ гРт]'к, то функции формы формально остаются прежними, N (г,г) = Sp]kIS']k •
J
Ж
b
Рис. 2. К определению функций формы: a - двумерные функции формы, b - осесимметричные функции формы Fig. 2. To the shape functions definition: a - two-dimensional shape functions, b - axisymmetric shape functions
Система линейных алгебраических уравнений, полученная после ансамблирования локальных систем (16), сохраняет свою структуру для любой двумерной постановки. Однако чтобы перейти к цилиндрическим координатам, необходимо учесть изменение правила действия дифференциальных операторов и якобиан преобразования при вычислении интегралов.
В цилиндрической системе координат [г,ф,z} с условием независимости всех рассматриваемых величин от окружной координаты (X(r$,z) = Ä(r,z)) дифференциальные операторы градиента скалярного поля и дивергенции векторного поля записываются следующим образом:
dF ^ dF
dar
daz
VF(r,z) = — e — e 2, V- a(r,z) =-+ — +-
dr
dz
dr
dz
a
r
r
где е I и е 2 - орты системы координат {г,г}, а обозначения аг и а1 отвечают соответствующим компонентам векторного поля а. Якобиан перехода к цилиндрическим координатам J = г.
4.2. Интегрирование уравнений движения в конечно-элементном представлении
Примем следующие допущения и обозначения:
• как и ранее, считаем функции формы Ы' линейными функциями координат г и 1;
обозначим за r элемента Q
V
mid
радиальную координату средней точки конечного
• |гёО. = — = V , где V - объем конечного элемента, а V -нормированный объем;
• |Ы' (х)гйО« ; □
• в первом приближении считаем координату г постоянной на конечном элементе и равной гт.
Перейдем к интегрированию компонент матриц, входящих в локальную систему уравнений (16), и определим «добавочные члены», возникающие в осесимметричных задачах. Впоследствии эти «добавочные слагаемые» учтены в соответствующем модуле для программного комплекса Kratos. Действие оператора градиента в цилиндрических координатах {г, г} формально не отличается от действия градиента в «плоской» системе координат {х,у}. Таким образом, матрица [О], входящая в систему (16), не изменится при переходе к цилиндрическим координатам, а при ее интегрировании производные функции формы необходимо умножить на объем конечного элемента V.
Рассмотрим уравнение неразрывности в цилиндрических координатах и запишем его с учетом конечноэлементного представления скорости:
„ _ dv' v
V-v =-+ —
dr r
r dvz + -
dz
к j=1
f
\
dNj Nj -+ —
dr r v /
+v
dNj
к
(
= z
j=1
Проинтегрируем получившееся выражение: 274
dNj z dNj
dr
+v
dz
K Nj + Zvi—.
j=1
r
m
r
v
r
K
(
JE
qj=i
dN} z dNJ
dr
- + Vi
dz
K NJ
rdQ+jE vrJ-rdQ =
QJ=1 J r
K
= E
J=1
(
dNJ 7 dNJ
Л
dr
+V
dz
K
V + E (vr JNJ dQ).
J=1
Интеграл от функции формы можно посчитать аналитически или численно интегрированием по гауссовым точкам; с учетом выбранных допущений имеем
К 1 к
Е V] I Ы^П)«1е ],
]=1 О 3 ]=1
где S - площадь треугольного конечного элемента.
Выражение (18) содержит два слагаемых. Первое из них после формальной подстановки ]—х, г—у будет являться конечноэлементным представлением
дивергенции в координатах {х,у}. Выражение (19) представляет собой «добавочное слагаемое» к матрице дивергенции которая учитывает
переход к цилиндрической системе координат.
Рассмотрим «вязкое» слагаемое. Проинтегрируем его по объему конечного элемента с учетом якобиана цилиндрической системы координат:
IМ->У{УЫ')]dО. = {V • (]Ы]ш' )йЮ - |У(Ы->]=
ООО
= |гЫ]VN' •ndS- |V(Ы->]у^ЫЧО. « -|V(Ы-> г)VN'dО.
S О О
Для полученного интеграла можно записать приближенное выражение
JV(NJr)VNl dQ = Q
dNl dNJ dNl dNJ --+--
dr dr dz dz
T. dNl ,NJ
V +-J-rdQ г
dr Q r
dNl dNJ dNl dNJ +
dr dr dz dz
V +1S N 3 dr
Первое слагаемое в (20) - такое же, как и в случае плоской задачи с точностью до замены обозначений координат; второе - «добавочное слагаемое», которое необходимо ввести в матрицу вязкости \Ь ] при решении задачи в осесимметричной постановке.
К интегрированию компонент матрицы масс вернемся ниже.
Локальную систему линейных алгебраических уравнений (16) можно записать
в сокращенной форме
[A\x} = {f },
r
V
r
V
Q
где [л] - матрица левой части; {x} - вектор-столбец неизвестных; {f} -вектор-столбец правой части. При переходе к осесимметричной постановке матрицу левой части можно представить в виде суммы двух матриц, одна из которых равна матрице, получающейся при решении двумерной задачи [Л^д , а вторая учитывает описанные выше «добавочные слагаемые», возникающие из-за осесимметричности [л]г:
[Л] = [A]2D .
4.3. Обеспечение устойчивости схемы
Несмотря на простоту изложенной процедуры приближенного интегрирования, решение системы уравнений, полученной ансамблированием матриц, записанных вышеописанным образом, не всегда будет устойчивым. Для этого часть слагаемых записывают в модифицированном виде.
Рассмотрим следующее слагаемое из уравнения (14):
□ '
При составлении локальной системы данное слагаемое распадется на три компоненты. Однако при непосредственном его интегрировании в составленной матрице не будет диагонального преобладания, что приведет к ее плохой обусловленности. Однако (22) можно преобразовать таким образом, чтобы не нарушить условие диагонального преобладания:
iы] n\+ы2ы2+ы3ы3 « iы]ы] n2+ы3 = и] |ы] сш и - ы]¥. □ □ □ 3
При этом было принято допущение о малом различии величин скоростей в узлах конечного элемента.
Запись компонент матрицы вязкости в виде (20) также приводит к плохой обусловленности матрицы [А] Авторы работы [6] предлагают использовать следующий способ вычисления компонент матрицы вязкости:
Ч, ] = м№ ]т [С ][В] ]ап, □
где матрицы [В' ] и [С ] имеют следующий вид:
[B ]=
Nir 0 "2 0 0 0"
0 N iz , [c]= 0 2 0 0
N!z N]r 0 0 1 0
N1 /r 0 0 0 0 2
Здесь запись или означает дифференцирование I -й функции формы по координате г или 2 соответственно. При такой записи компонент матрицы вязкости в [Ь ] появляются компоненты вида
| Ы^Ы^а = | ЫЫ^а.
□ г г а г
При точном интегрировании данного выражения возникают слагаемые вида
г
1п — , где г и г! - радиальные координаты узлов I и 1 конечного элемента.
О 1
При удалении элемента на достаточно большое расстояние от оси величина логарифма будет стремиться к нулю и данные слагаемые не будут приводить к появлению каких-либо сложностей. Однако при приближении элемента к оси величина логарифма неограниченно растет, что вносит ошибку в расчеты. При реализации осесимметричной модели РРБМ-2 при численном интегрировании выражений использовались формулы приближенного интегрирования по одной, трем и шести гауссовым точкам [5]. Стоит заметить, что при добавлении в локальную матрицу аппроксимации уравнения неразрывности в форме (18) с учетом (19) соответствующие уравнению строки получаются приблизительно одинаковыми, что также приводит к плохой обусловленности системы. Чтобы в какой-то степени решить эту проблему, в уравнение неразрывности можно добавить дополнительное слагаемое, «уравновесив» его слагаемым в правой части. К примеру, таковым может быть лапласиан давления. Запишем его в конечно-
элементной формулировке, домножим на функцию формы N и
проинтегрируем по объему конечного элемента:
( \
| ы1 др"+хгйа = рги+11 Ы1 жг^а = Ли+1 ^Ы^уЫ-МБ - /уы^уы^а а а Чб а у
г \
« - РП+1 /уы^УЫ1 ¿а = - РП+1 | гУЫ1 уы1 ¿а +1 ЫУГУЫ ¿а
а ча а у
«-р^УИ1 УИгУ +1 ^б! (23)
Ч 3 дг У
Полученное выражение умножается на стабилизирующий параметр т и добавляется в строки локальной матрицы [Л], соответствующие уравнению неразрывности. В соответствующие компоненты вектора правой части войдет похожее выражение. Полагая значение Ур" известным и в простейшем случае равным £ , получаем:
JV- grdQ«-VNJ • gV. Q
С учетом (23)-(24) локальную систему уравнений можно записать в виде
At
-T[D]g
О [M ]+m At [G]
[D] t[L]
{P}
n+1
( О ^
pg{R} + О {V }n At
В рассматриваемой реализации РРБМ-2 стабилизирующий параметр т был выбран равным [20]
Г Л-1
= 1 А! И2 И '
V И У
где И - длины наименьшей грани конечного элемента; уп - модуль средней на элементе скорости в момент времени tn . Зависимость параметра от А! необходима для случаев невязких течений в ячейках с нулевой скоростью.
4.4 Обеспечение баланса массы в осесимметричном случае
Так как схема метода конечных элементов не является консервативной, необходимо использовать инструмент, позволяющий корректировать объем жидкости с целью обеспечения его постоянства. Это актуально при моделировании двухфазных течений и течений со свободной границей. В программном комплексе Kratos для того, чтобы обеспечить сохранение объема жидкости, достаточно заменить в соответствующей процедуре объемы плоских конечных элементов на объемы осесимметричных.
5. Решение модельных задач
5.1 Течение жидкости в цилиндрической трубе
Рассмотрим установившееся течение вязкой несжимаемой жидкости в цилиндрической трубе радиусом Я и длиной Ь (рис. 3). Течение происходит под действием постоянной разности давлений р1 и р2, заданных на концах трубы. Коэффициент динамической вязкости равен ц . На стенках трубы
заданы граничные условия прилипания V |]=к = 0 . На оси - равенство нулю радиальной скорости у] |]=о= 0.
Рис. 3. Моделирование течения Пуазейля в круглой трубе Fig. 3. Simulation of the Poiseuille flow in a pipe
Такое течение (при малых скоростях, когда обеспечиваеся ламинарный режим течения) называется течением Пуазейля и имеет аналитическое описание (профиль Пуазейля). Распределение скорости в трубе в зависимости от расстояния до оси симметрии r дается выражением
v(r) = ^(^2 -r2). 4/lIL
Для решения задачи выберем следующие параметры L = 1,0 м, R = 0,15 м, р = 1000,0 кг/м3, Ap = p2 -Р1 = 0,013 Па, ^ = 0,001 Па с. Задача решалась при шаге пространственной дискретизации Ah = 0.01 м.
На рис. 4 приведены полученные результаты. Для построения графиков использовались данные о скоростях, взятые в центральном сечении трубы (при z = L/2). Из графиков видно удовлетворительное совпадение результатов расчета с аналитическим решением.
Рис. 4. График зависимости скорости от радиальной координаты (синяя кривая - профиль Пуазейля, красная кривая - результат расчета) Fig. 4. Flow velocity dependence on the radial coordinate (blue curve - Poiseuille velocity profile; red curve - result of the flow simulation)
5.2 Падение капли в слой жидкости
Рассмотрим каплю жидкости с параметрами pf =1179.0 кг/м3 и
/df = 1.86 -10 5 Па с, падающую в бассейн такой же жидкости под прямым углом (рис. 5). Скорость капли в момент соударения vimp = 2,56 м/с.
Параметры окружающего жидкость газа pa =1.204 кг/м и /а = 1.90 -10 Па с. На поведение жидкости в общем случае сильно влияет значение поверхностного натяжения [7], однако, так как в программном комплексе Kratos отсутствуют модели поверхностного натяжения для PFEM-2, примем коэффициент поверхностного натяжения равным нулю. На боковой и нижней стенках области заданы граничные условия прилипания, но оси -
ограниченность радиальной скорости v , на верхней границе равенство нулю давления p .
Рис. 5. Моделирование падения капли в глубокий слой жидкости
Fig. 5. Simulation of the droplet impact onto a deep pool
Задача о падении капли в слой жидкости аналитического решения не имеет, и из-за допущения о равенстве нулю коэффициента поверхностного натяжения становится невозможным сравнение полученных результатов с данными экспериментов. Поэтому решение задачи о падении капли с помощью PFEM-2 будем сравнивать с результатами численного эксперимента, поставленного в программном пакете с открытым исходным кодом Gems [8, 9], который позволяет решать задачи гидродинамики со свободной поверхностью интегро-интерполяционным методом. Его отличием от других программ этого класса является использование структурированных динамически перестраиваемых сеток. Gems успешно зарекомендовал себя в решении задач о падении капель как в слой жидкости, так и на твердую поверхность [10, 11]. Решения, полученные в задаче о падении капли в моменты времени t = 0,0; 0,00117; 0,00234; 0,00351; 0,00468; 0,00585 c, приведены на рис.6.
\
t=0,0 с
t=0,00234 с
1=0,00351 с
t=0,00468 с
/
t=0,00585 с
Рис. 6. Моделирование падения капли в глубокий слой жидкости Fig. 6. Simulation of the droplet impact onto a deep pool
Для того, чтобы количественно сравнить полученные результаты, построим графики зависимости глубины кратера от времени (рис. 7). Сравнение графиков для глубины кратера показывает также удовлетворительное соответствие между результатами численных экспериментов, проведенных в двух программных пакетах, Kratos и Gerris, в основе которых лежат разные численные методы решения задач гидродинамики.
h/d
Рис. 7. Схема для нахождения глубины кратера (a) и график зависимости приведенной глубины кратера от безразмерного времени (b) Fig. 7. The scheme for the crater depth determination (a) and the normalized crater depth dependency on dimensionless time (b)
Благодарности
Работа выполнена при частичной финансовой поддержке Российского фонда
фундаментальных исследований (проект 17-08-01468)
Список литературы
[1]. Idelsohn S.R., Onate E., Pin F.D. 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б vol. 61, No. 7, 2004, pp. 964-989. doi: 10.1002/nme.1096
[2]. Idelsohn S.R., Nigro N.M., Gimenez J.M., Rossi R., Marti J.M. A fast and accurate method to solve the incompressible Navier-Stokes equations. Engineering Computations, vol.30, No. 2, 2013, pp.197-222. doi: 10.1108/02644401311304854
[3]. Милн-Томсон Л.М. Теоретическая гидродинамика. М.: Мир, 1964. 660 стр.
[4]. Becker P. An enhanced Particle Finite Element Method with special emphasis on landslides and debris flows. PhD thesis, 2015, Universitat Politecnica de Catalunya.
[5]. Зенкевич О.С. Метод конечных элементов в технике. М.: Мир, 1975. 543 стр.
[6]. Hughes T.J.R., Liu W.K., Brooks A. Finite Element Analysis of Incompressible Viscous Flows by the Penalty Function Formulation. Journal of Computational Physics, vol. 30, No. 1, 1979, pp. 1-60. doi: 10.1016/0021-9991(79)90086-X
[7]. Yarin A.L. Drop Impact Dynamics: Splashing, Spreading, Receding, Bouncing. Annual Review of Fluid Mechanics, vol. 38, 2006, pp. 159-192. doi: 10.1146/annurev.fluid.38.050304.092144
[8]. Popinet S. Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries. Journal of Computational Physics, vol. 190, No. 2, . 2003, pp. 572600. doi: 10.1016/S0021-9991(03)00298-5
[9]. Korchagova V.N., Kraposhin M.V., Marchevsky I.K., Smirnova E.V. Simulation of droplet impact onto a deep pool for large Froude numbers in different open-source codes. Journal of Physics: Conference Series, vol. 918. art. 012037, 2017. doi: 10.1088/1742-6596/918/1/012037
[10]. Agbaglaha G., Thoravala M.-J., Thoroddsen S.T., Zhanga L.V., Fezzaaa K., Deegana R.D. Drop impact into a deep pool: vortex shedding and jet formation. Journal of Fluid Mechanics, vol. 764, 2015m pp. 1-12. doi: 10.1017/jfm.2014.723
[11]. Renardy Y., Popinet S., Duchemin L., Renardy M., Zaleski S., Josserand C., Drumright-Clarke M.A., Richard D., Clanet C., Quere D. Pyramidal and toroidal water drops after impact on solid surface. Journal of Fluid Mechanics, vol. 484, . 2003, pp. 69-83. doi: 10.1017/S0022112003004142
Axisymmetric viscous incompressible flow simulation by using the Particle finite element PFEM-2 method in the open source Kratos code
1 E.V. Smirnova <alena.davidova@ispras.ru> 1,21.K. Marchevsky <iliamarchevsky@mail.ru>
2 V.O. Bondarchuk <viktoriia.bondarchuk@gmail.com> 1 Ivannikov Institute for System Programming of the Russian Academy of Sciences, 25, Alexander Solzhenitsyn st., Moscow, 109004, Russia
2 Bauman Moscow State Technical University,
5, 2-ndBaumanskaya st., Moscow, 105005, Russia
Abstract. In this paper, the particle finite element method (PFEM-2) for the simulation of the axisymmetric flows of a viscous incompressible fluid is considered. The necessary equations for the description of the axisymmetric flows by using the particle finite element method with some assumptions were obtained. The numerical model for solving axisymmetric flows of a viscous incompressible fluid was implemented in the Kratos open-source code on the basis of the existing application for solving two-dimensional problems. The numerical model for the simulation of the axisymmetric flows of a viscous incompressible fluid was validated on two problems. First of them is the Poiseuille problem about viscous flow in the pipe. The numerical solution obtained by particle finite element method are in the satisfactory agreement with the analytical solution for this problem. The second test is a problem of a droplet impact onto a deep liquid pool with similar fluid. As the comparison with analytical results is impossible, the results of particle finite element method simulation were compared with results of the numerical simulation obtained by using the open-source package Gerris (Volume of Fluid method). The results of comparing of numerical simulations of droplet impact onto a liquid pool obtained by two different codes also in satisfactory agreement.
Keywords: particle finite element method; viscous incompressible flow, axisymmetric flow;; Poiseuille flow; droplet impact onto a deep pool.
DOI: 10.15514/ISPRAS-2018-30(2)-13
For citation: Smirnova E.V., Marchevsky I.K., Bondarchuk V.O. Axisymmetric viscous incompressible flow simulation by using the Particle finite element PFEM-2 method in the open source Kratos code. Trudy ISP RAN/Proc. ISP RAS, vol. 30, issue 2, 2018, pp. 263-284 (in Russian). DOI: 10.15514/ISPRAS-2018-30(2)-13
References
[1]. Idelsohn S.R., Onate E., Pin F.D. The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves. International journal for numerical methods in engineering6 vol. 61, No. 7, 2004, pp. 964-989. doi: 10.1002/nme.1096
[2]. Idelsohn S.R., Nigro N.M., Gimenez J.M., Rossi R., Marti J.M. A fast and accurate method to solve the incompressible Navier-Stokes equations. Engineering Computations, vol.30, No. 2, 2013, pp.197-222. doi: 10.1108/02644401311304854
[3]. Milne-Thomson L.M. Theoretical Hydrodynamics. London: Macmillan & Co, 1968.
[4]. Becker P. An enhanced Particle Finite Element Method with special emphasis on landslides and debris flows. PhD thesis, 2015, Universitat Politecnica de Catalunya.
[5]. Zienkiewicz O.C. The finite element method in engineering science. London: McGraw-Hill, 1971. 521 p.
[6]. Hughes T.J.R., Liu W.K., Brooks A. Finite Element Analysis of Incompressible Viscous Flows by the Penalty Function Formulation. Journal of Computational Physics, vol. 30, No. 1, 1979, pp. 1-60. doi: 10.1016/0021-9991(79)90086-X
[7]. Yarin A.L. Drop Impact Dynamics: Splashing, Spreading, Receding, Bouncing. Annual Review of Fluid Mechanics, vol. 38, 2006, pp. 159-192. doi: 10.1146/annurev.fluid.38.050304.092144
[8]. Popinet S. Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries. Journal of Computational Physics, vol. 190, No. 2, . 2003, pp. 572600. doi: 10.1016/S0021-9991(03)00298-5
[9]. Korchagova V.N., Kraposhin M.V., Marchevsky I.K., Smirnova E.V. Simulation of droplet impact onto a deep pool for large Froude numbers in different open-source codes. Journal of Physics: Conference Series, vol. 918. art. 012037, 2017. doi: 10.1088/1742-6596/918/1/012037
[10]. Agbaglaha G., Thoravala M.-J., Thoroddsen S.T., Zhanga L.V., Fezzaaa K., Deegana R.D. Drop impact into a deep pool: vortex shedding and jet formation. Journal of Fluid Mechanics, vol. 764, 2015m pp. 1-12. doi: 10.1017/jfm.2014.723
[11]. Renardy Y., Popinet S., Duchemin L., Renardy M., Zaleski S., Josserand C., Drumright-Clarke M.A., Richard D., Clanet C., Quere D. Pyramidal and toroidal water drops after impact on solid surface. Journal of Fluid Mechanics, vol. 484, . 2003, pp. 69-83. doi: 10.1017/S0022112003004142
