CC BY
170
Свободное программное обеспечение для моделирования жидкости со свободной поверхностью
Е.В. Давыдова < [email protected]> В.Н. Корчагова <[email protected]> МГТУ им. Н. Э. Баумана, 105005, Россия, г. Москва, ул. 2-я Бауманская, дом 5, стр. 1
Аннотация. Задачи течения вязкой несжимаемой жидкости со свободной поверхностью представляют собой отдельный класс задач механики сплошной среды и имеют большую практическую значимость. При моделировании таких процессов исследователь сталкивается с достаточно большим числом особенностей и ограничений, критически важных для получения корректного решения. Целью данной работы является обзор существующих численных методов, которые можно применить к решению задач течения жидкости со свободной поверхностью, и программных комплексов с открытым исходным кодом, в которых эти методы реализованы, а также выявление границ применимости рассмотренных программных комплексов. Были рассмотрены несколько численных методов (Volume of Fluid, Smoothed Particle Hydrodynamics, Particle Finite Element Method v.2), реализованные в пяти разных пакетах программ с открытым исходным кодом: OpenFOAM, Gerris, pySPH, DualSPHysics, Kratos. В качестве тестовых задач были выбраны следующие примеры: обрушение колонны жидкости и падение капли в слой жидкости. Решения, полученные в выбранных пакетах, были сопоставлены с известными результатами экспериментов. Проводилось сравнение времени, затраченного на решение задач в различных пакетах. Наилучшие результаты для выбранных тестов показали пакеты OpenFOAM и Gerris: в них содержатся необходимые для решения задач инструменты - возможность расчета задач в двумерной, трехмерной либо осесимметричной постановке, а также корректный учет поверхностного натяжения жидкости. При этом пакет Gerris дает существенное ускорение решения "крупных" задач за счет использования динамически перестраиваемых сеток. Пакет DualSPHysics направлен на решение задач прибрежной инфраструктуры, где рассматривают, как правило, трехмерные задачи и нет необходимости учета поверхностного натяжения. Пакет pySPH разрабатывался с целью наглядной демонстрации работы численного метода SPH, поэтому исходный код записан на языке Python без оптимизации. Пакет Kratos, в частности, модуль PFEM-2, в настоящий момент находится в разработке, поэтому некоторые возможности в нем еще не реализованы.
Ключевые слова: гидродинамика; переменная область течения; жидкость со свободной поверхностью; VOF; SPH; PFEM; открытый исходный код.
Davydova E., Korchagova V. Open-source software for modelling of free surface flows. Trudy ISP RAN /Proc. ISP RAS, 2014, vol. 28, issue 1, pp. 243-258
DOI: 10.15514/ISPRAS-2016-28(1 )-14
Для цитирования: Давыдова Е.В., Корчагова В.Н. Свободное программное обеспечение для моделирования жидкости со свободной поверхностью. Труды ИСП РАН, том 28, вып. 1, 2016 г., с. 243-258. DOI: 10.15514/ISPRAS-2016-28(1)-14
1. Введение
Исследование течений жидкости со свободной поверхностью представляет собой отдельный класс задач механики сплошной среды. Такие задачи имеют большую практическую значимость и возникают, к примеру, при моделировании плотин, дамб, решении экологических задач (загрязнение поверхностей, распространение примесей и т. д.). Как правило, при решении подобных задач рассматривают процессы брызгообразования, волнообразования, возникновения струйных течений и т. п. С математической точки зрения специфика таких задач состоит в том, что область течения является переменной, форма ее поверхности заранее неизвестна, и на этой поверхности задаются граничные условия специального вида. Во многих случаях в различные моменты времени область течения может состоять из разного числа связных компонент.
Даже в сравнительно простом случае моделирования ламинарного течения несжимаемой изотермической среды разработчик численного алгоритма или инженер-исследователь, выполняющий расчеты, может столкнуться с рядом условий, соблюдение которых может быть критически важным для успешного решения задачи. Этими условиями являются обеспечение консервативности и устойчивости разностной схемы, а также как можно более низкой ее диссипативности и монотонности в области разрыва (межфазной границы). Не менее важными для практики являются требования высокой вычислительной производительности и масштабируемости модели.
К дополнительным факторам, повышающим трудоемкость решения подобных задач, могут относиться необходимость учета движения границ расчетной области; высокие значения отношения плотностей жидкой и газообразной фаз; требования тщательного разрешения мелкомасштабных структур (пленок, пузырьков, капель), характерный размер которых может быть намного меньше, чем характерный масштаб длины в основной задаче и т. п. Весьма нетривиальным может быть обеспечение возможности сопряжения решения или вычислительного алгоритма в целом с математическими моделями из других разделов механики сплошной среды в рамках решения многодисциплинарных задач.
К сожалению, на практике удовлетворить всем перечисленным требованиям в полной мере невозможно. Поэтому численный метод решения поставленной задачи выбирается на основе некоторого компромисса в зависимости от целей исследования.
Давыдова Е.В., Корчагова В.Н. Свободное программное обеспечение для моделирования жидкости со свободной поверхностью. Труды ИСП РАН, том 28, вып. 1, 2016 г., с. 243-258
Целью данной работы является обзор существующих численных методов, которые можно применить к решению задач моделирования течения жидкости со свободной поверхностью, и программных комплексов с открытым исходным кодом, в которых эти методы реализованы, а также выявление границ применимости рассмотренных программных комплексов. В данной работе рассмотрены три численных метода:
• метод контрольных объемов с расчетом переноса объёмной доли жидкой фазы (FVM + VOF);
• метод сглаженных частиц (SPH);
• метод конечных элементов с частицами (PFEM).
Рассмотрены 5 программных средств, реализующих три вышеупомянутых численных метода решения уравнений Навье - Стокса для моделирования несжимаемых течений со свободной поверхностью: OpenFOAM, Gerris, pySPH, DualSPHysics, Kratos.
Для сравнения программ было выбрано два расчетных случая: задача моделирования обрушения двумерного столба жидкости и задача моделирования падения капли в слой жидкости.
2. Определяющие соотношения
Рассмотрим задачу о расчете течения вязкой несжимаемой жидкости со свободной поверхностью. Для описания течений жидкости и газа с учетом выбранных допущений о несжимаемости, изотермичности, наличии вязких сил и сил поверхностного натяжения используются следующие уравнения:
- условие несжимаемости:
V-U = 0; (1)
- уравнение сохранения импульса (уравнение Навье - Стокса), учитывающее наличие массовых сил:
p(jTt + (°-V)I}) = -VV + ^^ + P9- (2)
Здесь U — поле скоростей среды; р — поле давления; т = ß(VU + VUT) — тензор вязких напряжений в соответствующей фазе; р — плотность; ß — коэффициент динамической вязкости; g — ускорение свободного падения. Начальное распределение скоростей считается заданным. Граничные условия на бесконечности записываются следующим образом:
U ^ 0, Vp ^ —рд при Irl ^ ю. На поверхностях раздела фаз ставятся следующие условия [1]: Ür = Um, (Рг — Рж)п + (Гг — = —2акп,
где к — средняя кривизна поверхности раздела фаз; о — коэффициент поверхностного натяжения; индекс "ж" обозначает значение физической величины в жидкости, индекс "г" — значение физической величины в газе. Предположение о несжимаемости газа справедливо в случае, когда U « сг, где сг - скорость звука в газовой фазе.
3. Численные методы и программные комплексы
3.1. Метод контрольных объемов с расчетом переноса объёмной доли жидкой фазы
Суть подхода, основанного на расчете переноса объемной доли жидкости — Volume of Fluid (VOF), заключается в следующем: для аппроксимации скачков физических величин в рамках метода контрольных объемов вблизи поверхности раздела фаз используется «приближение смеси», «сглаживающее» границу раздела фаз в пределах некоторой области и позволяющее искать решение методом «сквозного счета». Для этого вводят маркерную функцию у (рис. 1), физический смысл которой — объемная доля жидкости [2]:
• у = 0 в областях, заполненных газом;
• у = 1 в областях, заполненных жидкостью;
• 0 < у < 1 в области перехода между фазами.
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
ТГ8 0.9 0 0 0
1 1 1 У2 0 0
1 1 1 0 0
Рис. 1. Значения маркерной функции в ячейках сетки Fig.1. Values oof marker function in the cells oof grid
Уравнения (1) и (2) рассматриваются для смеси, скорость, плотность и динамическая вязкость которой определяются следующим образом:
Л = уЛж +
Р = УРж+ (1-у)рГ; v = w ж+ (1-к)^г.
Давыдова Е.В., Корчагова В.Н. Свободное программное обеспечение для моделирования жидкости со свободной поверхностью. Труды ИСП РАН, том 28, вып. 1, 2016 г., с. 243-258
Для маркерной функции решают уравнение ее переноса жидкой фазой
ду
которое можно записать в виде:
d-£+V-(U Y) + V-(V Y(1-Y)) = 0, (3)
где 7 — относительная скорость движения тяжелой фазы (жидкости) относительно легкой (газа). Эта величина является неизвестной, поэтому при аппроксимации соотношения (3) ее заменяют на некоторое поле скоростей, вводимое искусственным образом. Введенное поле скоростей «сжимает» границу раздела фаз, не допуская слишком сильного ее «размазывания». При этом оно мало влияет на решение в целом, т. к. его действие проявляется в тонком слое вокруг поверхности раздела.
Для численного интегрирования определяющих соотношений используется метод контрольного объема. Подробно эта процедура описана в [2].
3.2. Метод сглаженных частиц
Метод сглаженных частиц — Smoothed Particle Hydrodynamics (SPH) принадлежит к классу бессеточных методов. Суть метода заключается в том, что сплошная среда, находящаяся в расчетной области, аппроксимируется множеством частиц [3]. Каждая частица ассоциирована с плотностью, скоростью, давлением среды в точке, в которой она расположена в текущий момент времени. Для восстановления непрерывных полей физических величин между частицами производится интерполирование на основе выкалывающего свойства дельта-функции:
А(х) = J А(х')5(х'-х) dnin VxERd.
Rd
Здесь А(х) — произвольное непрерывное распределение; 5(х) — d--мерная дельта-функция Дирака с присущими ей свойствами; d — размерность пространства.
В методе SPH вместо дельта-функции 8(х) используется некоторая гладкая функция W(x,h), называемая функцией ядра, которая зависит от пространственной координаты и параметра h, отлична от нуля лишь в малой окрестности начала координат (либо экспоненциально убывает на бесконечности) и удовлетворяет условию нормировки
J W(x,h) dn$ = 1.
Параметр h, называемый длиной сглаживания, определяет размер окрестности центра функции ядра, вне которой она обращается в ноль либо становится пренебрежимо малой. В качестве функции ядра может выступать функция Гаусса, полиномы и сплайны разных порядков и др.
Таким образом, интерполирование неизвестной физической величины А (х) по ее значениям в точках, в которых располагаются частицы, производится следующим образом (рис. 2):
А(х) « ^ A(xj)w(xj - х, h) AVj = ^ A(xj)w(xj - х, h)
i i Pi Здесь m.j — масса частицы с координатами Xj, p¡ — плотность среды в
Рис. 2. Круги радиусами h и характерный профиль функции ядра Fig.2. Circles of radius h and the characteristic profile of the core function
точке Xj.
3.3. Метод конечных элементов с частицами
Метод конечных элементов с частицами — Particle Finite Element Method (PFEM) — это гибридный эйлерово-лагранжевый метод, объединивший в себе достоинства бессеточных методов и метода конечных элементов [4]. Первая версия этого метода появилась в 2001 году. В области, занятой жидкостью, вводится конечно-элементная сетка, узлы которой рассматриваются как лагранжевы частицы. Частицы движутся вместе со средой, что позволяет "автоматически" учесть конвективное слагаемое в уравнениях Навье — Стокса. Это приводит к необходимости перестроения (деформирования) сетки на каждом шаге расчета по времени. На вновь построенной сетке решаются уравнения гидродинамики (которые в отсутствие
Давыдова Е.В., Корчагова В.Н. Свободное программное обеспечение для моделирования жидкости со свободной поверхностью. Труды ИСП РАН, том 28, вып. 1, 2016 г., с. 243-258
конвективного члена становятся линейными) методом конечных элементов. Сетка может разделяться на несвязанные между собой подобласти (например, при моделировании брызг), в этом случае уравнения решаются независимо в каждой подобласти (рис. 3).
Рис. 3. Подобласти, занятые жидкой фазой, в методе PFEM
Fig.3. Subregions employed by liquid phase in the PFEM method
Такой подход позволяет упростить процесс решения уравнений (так как нет необходимости аппроксимации нелинейных слагаемых), но процедура перестроения сетки является нетривиальной и трудоемкой задачей: время перестроения сетки сопоставимо со временем решения уравнений и даже может превышать его. Также частицы в процессе решения могут сильно сближаться, что приводит к уменьшению размера ячеек сетки и, как следствие, шага расчета по времени.
Для ускорения счета в 2012 году была разработана модификация метода — PFEM-2 [5]. В ней конечно-элементная сетка является неподвижной и достаточно грубой, а частицы, которые вводятся для учета конвективного члена, перемещаются из одной ячейки сетки в другую вместе со средой (рис. 4). Движение частиц рассчитывается явным методом с малым шагом по времени, соответствующим числу Куранта CFL ~ 0.1. Учет остальных сил (вязкие силы, градиент давления, внешние массовые силы) осуществляется путем решения соответствующего линейного уравнения на сетке с большим числом Куранта (CFL ~ 1...10), что позволяет существенно снизить время выполнения расчета. Разработаны специальные процедуры проецирования характеристик частиц на узлы сетки и наоборот.
О о с/ oo ¡So <yo 0 ( / 3 о г/ ?Уо /О О S п ^ Уп ° у о / А0
J о / °o¥ JO О С о о/ < °° у° О щ
°QJr> Ус / о; гг о ( иоПоУ ' "Уо Уос УЬ о 0 • /У •у3 о К'.9
/ О fjr О ¡/о ° 0 '-V о о/ Ро о О oj ш Щ
.•Л v., % »
Рис. 4. Частицы в ячейках неподвижной сетки в методе PFEM-2 Fig. 4. Particles in cells of fixed grid in the method PFEM-2
Важной особенностью метода PFEM-2 является наличие "монолитных" решателей, не предполагающих использования процедур расщепления по физическим процессам.
3.4. Программные комплексы
Рассмотрены пять программных комплексов с открытым исходным кодом, в которых реализованы указанные выше численные методы:
• OpenFOAM (http://www.openfoam.org/) - классический метод Volume of Fluid, реализованный на C/C++ и успешно себя зарекомендовавший во многих практических приложениях;
• Gerris (http://gfs.sourceforge.net/wiki/) - метод Volume of Fluid с динамическими адаптивными сетками, реализованный на C/C++, позиционируемый как пакет, позволяющий максимально точно учитывать капиллярные явления;
• pySPH (https://pysph.readthedocs.org/) - метод SPH, реализованный с использованием языка программирования Python и технологии компиляции при исполнении Cython;
• DualSPHysics (http://www.dual.sphysics.org/) - метод SPH, реализация на C/C++ с возможностью проведения расчетов как на CPU, так и на графических ускорителях, поддерживающих технологию CUDA;
• Kratos (http://www.cimne.com/kratos/) - метод PFEM-2, реализация — на C/C++ для дискретизации и решения уравнений в частных производных и на языке Python для подготовки исходных данных и управления процессом интегрирования.
4. Тестовые задачи
Для сравнения программ было выбрано два расчетных случая. Первым тестовым примером является задача моделирования обрушения двумерного столба жидкости (воды), в которой влияние поверхностного натяжения пренебрежимо мало. Данная задача является одним из наиболее известных бенчмарк-тестов численных моделей механики жидких сред с разделом фаз, для которого в большом количестве доступны экспериментальные данные. Во втором тестовом примере - задаче о моделировании падения капли в слой жидкости - поверхностное натяжение играет существенную роль. Во всех вышеперечисленных программных комплексах вопрос учёта поверхностного натяжения и связанных с ним эффектов стоит особняком. В некоторых пакетах не реализованы численные модели поверхностного натяжения либо их корректность вызывает сомнения, поэтому этот тестовый пример был рассмотрен только в двух из пяти представленных пакетов: OpenFOAM и Gerris.
4.1. Обрушение колонны жидкости
Рассмотрим двумерную задачу об обрушении колонны жидкости (воды) в прямоугольном резервуаре с гладкими стенками [6]. Геометрия расчетной области указана на рис. 5. В качестве характерного размера выбрана ширина колонны Ь = 14.6 • 10-2 м. Плотность жидкости р равна 999 кг/м3, коэффициент динамической вязкости ^ равен 10-3 Па-с, коэффициент поверхностного натяжения а равен 0.0727 Н/м.
А
<м
Рис. 5. Расчетная схема задачи об обрушении колонны жидкости Fig. 5. Computational scheme of the problem of liquid column collapse
В [6] приведены экспериментальные данные для Т < 1 с. Шаг по времени выбирался так, чтобы число Куранта не превышало значения 0.2. В расчетной области строится сетка либо задается распределение частиц так, чтобы в начальный момент времени в области, занятой жидкостью, было около 3000 частиц (либо узлов сетки).
Полученные результаты хорошо соотносятся с экспериментальными данными (рис. 6, 7). На рис. 6 представлены качественные результаты, позволяющие визуально сравнить полученные решения. На рис. 7 указаны соответствующие значения длины растекания жидкости в разные моменты времени для трех серий экспериментов и для полученных в пакетах численных решений. Характерные оценки времени счета в разных пакетах следующие.
• OpenFOAM: 1 мин;
• Gerris: 4 мин;
• pySPH: 9 мин;
• DualSPHysics: 3 мин;
• Kratos: 3 мин.
Пакет pySPH решил задачу за наибольшее время. Это может быть вызвано тем, что пакет написан на языке программирования высокого уровня Python, а
также тем, что код не подвергался оптимизации, так как пакет разрабатывался для учебных целей (наглядной демонстрации работы метода SPH). Задача была решена в пакете Gerris в четыре раза медленнее, чем в OpenFOAM, несмотря на то, что в этих пакетах используется один и тот же численный метод (VOF). Это связано с тем, что на "маленькой" задаче время, затрачиваемое на адаптивное перестроение сетки в Gerris, существенно превышает время решения уравнений гидродинамики.
experiment DualSPHysics Gerris
Kratos pySPH OpenFOAM
а
experiment DuaiSPHysics Gerris
Kratos PySPH OpenFOAM
б
Рис. 6. Качественное сравнение результатов расчетов в различных пакетах с экспериментом для задачи обрушения водяной колонны: а -1=0.2 с; б -1=0.8 с Fig. 6. A qualitative comparison of results of calculation in a variety ofpackages for the experiment with the problem of collapse of the water column: a -1=0.2 s; б -1=0.8 s
Рис. 7. Зависимость приведенной длины растекания жидкости от безразмерного времени для полученных численных решений и сравнение их с экспериментальными
данными
Fig. 7. Dependence of the reduced length of liquid spreading on the dimensionless time for obtaining numerical solutions and their comparison with experimental data
4.2. Падение капли в слой жидкости
Капля диаметром й = 2.67 • 10-3 м из вязкой несжимаемой жидкости (смеси глицерина и воды) падает в слой такой же жидкости толщиной Н = 2й вертикально вниз (рис. 8). Ширина слоя жидкости считается бесконечно большой по сравнению с размерами падающей капли. Скорость капли и в момент соударения с жидкостью равна 2.56 м/с. Плотность жидкости р равна 1179 кг/м3, коэффициент динамической вязкости ^ равен 18.57 • 10-3 Па-с, коэффициент поверхностного натяжения а равен 0.0668 Н/м.
% Э\
a ц р Н
Рис. 8. Расчетная схема для задачи о падении капли в слой жидкости Fig. 8. Computational scheme for the problem offalling drops into a layer of liquid
В численной постановке задача считается осесимметричной. Размеры расчетной области подбираются так, чтобы возмущения решения на краях области были пренебрежимо малы. В соответствии с [7] размеры расчетной области составляют 45*45 мм. Время расчетов Т составляет 50 мс. Число Куранта, по которому производится автоматический выбор шага по времени, равно 0.2. Сетка, используемая в расчетной области, подбирается так, чтобы в радиусе капли содержалось около 30 ячеек.
В пакетах DualSPHysics и Kratos нет моделей поверхностного натяжения и возможности решения задач в осесимметричной постановке. Модели поверхностного натяжения, реализованные в пакете pySPH, требуют дополнительного анализа, так как результаты, показанные им, вызывают сомнения. В связи с этим задача о падении капли в жидкость была решена в пакетах OpenFOAM и Gerris.
На рис. 9, 10 показаны количественные результаты расчетов: измерены глубина кратера, появляющегося после падения капли в жидкость, и диаметр этого кратера, измеренный на половины высоты слоя жидкости. Наблюдается хорошее согласие между экспериментальными данными и расчетами.
Рис. 9. Зависимость приведенной глубины кратера от безразмерного времени
Рис. 10. Зависимость приведенного диаметра кратера, измеренного на половине высоты слоя жидкости, от безразмерного времени Fig. 10. The dependence of the reduced diameter of a crater measured at half the height of the liquid layer on dimensionless time
Задача о падении капли в жидкость является достаточно "большой": сетка содержит около 250 000 ячеек (в OpenFOAM). Расчеты проводились в параллельном режиме, поэтому целесообразно указать характерные оценки времени решения в расчете на одно ядро.
Характерное время решения задачи в OpenFOAM составило около 12 часов, а в Gerris - 0.5 часа. Высокая эффективность пакета Gerris обеспечивается адаптивным перестроением сетки, которое позволяет существенно снизить ее размер (в 30-100 раз) и, соответственно, сильно снизить время расчетов.
5. Заключение
Рассмотрены три основных метода решения задач расчета течения вязкой несжимаемой жидкости со свободной поверхностью — Volume of Fluid, Smoothed Particle Hydrodynamics, Particle Finite Element Method, и разные их реализации в пяти программных комплексах с открытым исходным кодом. В этих пакетах программ были решены тестовые задачи. На основе полученных решений и при рассмотрении исходных кодов и документации к проектам сделаны следующие выводы о каждом из рассмотренных программных комплексов.
• OpenFOAM: универсальный пакет — реализованы математические модели многих физических эффектов; хорошо масштабируется.
• Gerris: узкоспециализированный пакет для решения задач расчета течений жидкости со свободной поверхностью; на "больших" задачах работает существенно быстрее, чем OpenFOAM.
• DualSPHysics: инженерный пакет для исследования прибрежной инфраструктуры; не реализован учет поверхностного натяжения, отсутствует возможность решения осесимметричных задач; есть возможность использования CUDA.
• pySPH: позволяет решать широкий класс задач; прост в использовании; имеются несколько моделей поверхностного натяжения; производительность — низкая.
• Kratos (PFEM-2): возможность решения ограниченного класса задач, в том числе сопряженных; наличие «монолитных» решателей; отсутствует возможность решения осесимметричных задач, нет учета поверхностного натяжения; высокая производительность.
Список литературы
[1]. J.U. Brackbill, D.B. Kothe, C. Zemach. A Continuum Method for Mode. Journal of
Computation Physics, 100(2), 1992. P. 335-354. DOI: 10.1016/0021-9991(92)90240-Y
[2]. J.H. Ferziger, M. Peric. Computational Methods for Fluid Dynamics. Berlin: SpringerVerlag, 2002. 426 p. DOI:10.1007/978-3-642-56026-2
[3]. Shaofan Li, Wing Kam Liu. Meshfree and Particle Methods. Berlin: Springer-Verlag,
2004. 502 p. DOI:10.1007/978-3-540-71471-2
[4]. S.R. Idelsohn, E. Onate, F. Del Pin. 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, 61(7), 2004. P. 964-989. D01:10.1002/nme. 1096
[5]. J.M. Gimenez, N.M. Nigro, S.R. Idelsohn. Evaluating the performance of the particle finite element method in parallel architectures. Comp. Part. Mech., 1(1), 2014. P. 103116. DOI: 10.1007/s40571-014-0009-4
[6]. S. Koshizuka, Y. Oka. Moving-Particle Semi-Implicit Method for Fragmentation of Incompressible Fluid. Nuclear Science and Engineering, 123, 1996. P. 421-434.
[7]. E. Berberovic, N.P. van Hinsberg, S. Jakirlic, I.V. Roisman, C. Tropea. Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution. Physical Review, 79, 2009. P. 1-15. DOI: 10.1103/PhysRevE.79.036306
Open-source software for modelling of free surface flows
E. Davydova <[email protected]> V. Korchagova <[email protected]> BMSTU, ul. Baumanskaya 2-ya, 5/1, Moscow, 105005, Russian Federation
Abstract. Problems of free-surface flow of viscous incompressible fluid are very useful in different practical cases. There are many specifies and limitations in these problems which are critically important for correct solving. The main goal is the review of existing numerical methods which can apply for modeling of free-surface flows and open-source programs where these methods are realized. Three methods for solving problems of free-surface flow were considered: Volume of Fluid, Smoothed Particle Hydrodynamics, Particle Finite Element Method v.2. They are realized in five open-source packages: OpenFOAM, Gerris, pySPH, DualSPHysics, Kratos. These packages were compared by modeling of two chosen cases: breaking of a dam and droplet impact to the liquid layer. Results of computations were compared with experimental results. There are good coincidence between them. The best results were obtained in OpenFOAM and Gerris. All main tools for modeling of free-surface flow are realized in these packages - the possibility of computations in 2D, 3D and axisymmetric model setup and also correct modeling of surface tension. Gerris can significantly accelerate computations in "big cases" due to dynamically adaptive remeshing. Further, DualSPHysics is the package for modeling of problems of coastal infrastructure where the most number of cases is 3D and the surface tension effect is negligible. The package pySPH was designed for clear demonstration of SPH working. The pySPH source code is on the Python language and not optimized. Kratos is the new package, which is in development now, therefore some tools are not developed in this moment.
Keywords: hydrodynamics; variable flow region; free-surface flow, VOF; SPH; PFEM; open-source packages.
DOI: 10.15514/ISPRAS-2016-28(1 )-14
For citation: Davydova E., Korchagova V. Open-source software for modelling of free surface flows. Trudy ISP RAN /Proc. ISP RAS, 2014, vol. 28, issue 1, pp. 243-258 (in Russian). DOI: 10.15514/ISPRAS-2016-28(1)-14
References
[1]. J.U. Brackbill, D.B. Kothe, C. Zemach. A Continuum Method for Mode. Journal of Computation Physics, 100(2), 1992. P. 335-354. DOI: 10.1016/0021-9991(92)90240-Y
[2]. J.H. Ferziger, M. Peric. Computational Methods for Fluid Dynamics. Berlin: SpringerVerlag, 2002. 426 p. D0I:10.1007/978-3-642-56026-2
[3]. Shaofan Li, Wing Kam Liu. Meshfree and Particle Methods. Berlin: Springer-Verlag, 2004. 502 p. D0I:10.1007/978-3-540-71471-2
[4]. S.R. Idelsohn, E. Onate, F. Del Pin. 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, 61(7), 2004. P. 964-989. D0I:10.1002/nme. 1096
[5]. J.M. Gimenez, N.M. Nigro, S.R. Idelsohn. Evaluating the performance of the particle finite element method in parallel architectures. Comp. Part. Mech., 1(1), 2014. P. 103116. DOI: 10.1007/s40571-014-0009-4
[6]. S. Koshizuka, Y. Oka. Moving-Particle Semi-Implicit Method for Fragmentation of Incompressible Fluid. Nuclear Science and Engineering, 123, 1996. P. 421-434.
[7]. E. Berberovic, N.P. van Hinsberg, S. Jakirlic, I.V. Roisman, C. Tropea. Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution. Physical Review, 79, 2009. P. 1-15. DOI: 10.1103/PhysRevE.79.036306
