Моделирование бессеточными методами взаимодействия потока жидкости с твердыми телами Текст научной статьи по специальности «Физика»

Известия Института математики и информатики УдГУ

2017. Том 50

УДК 519.63, 532.5

© А. С. Сармакеева, Л. Е. Тонкое, А. А. Чернова

МОДЕЛИРОВАНИЕ БЕССЕТОЧНЫМИ МЕТОДАМИ ВЗАИМОДЕЙСТВИЯ ПОТОКА ЖИДКОСТИ С ТВЕРДЫМИ ТЕЛАМИ1

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

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

DOI: 10.20537/2226-3594-2017-50-05 Введение

Необходимость решения сопряженных задач (FSI — Fluid-Structure Interaction) взаимодействия потоков жидкости с плавающими телами возникает при прогнозировании последствий прорыва и разрушения гидротехнических сооружений, что в свою очередь служит целям предотвращения и минимизации последствий чрезвычайных ситуаций подобного рода. Решение таких задач связано с моделированием катастрофических потоков в условиях сложной трехмерной геометрии, с учетом динамики множества вовлекаемых течением объектов и характера их взаимодействий.

Одним из устоявшихся подходов к решению сопряженных задач FSI является применение алгоритмов итерационного связывания решений физических подзадач, получаемых сеточными методами [1]. При этом подзадача динамики жидкости численно решается на основе эйлерова описания, методом контрольных объемов или его модификаций, а динамика твердых тел разрешается на основе лагранжева подхода методом конечных элементов. Однако согласование решений на границах твердого тела и жидкости является сложной и трудоемкой задачей, связанной с необходимостью консервативной интерполяции сеточных величин (давления, перемещений) как правило на несогласованных сетках. Дополнительные сложности возникают при моделировании течений со свободной поверхностью. Например, широко распространенным является метод объема жидкости в ячейке (VOF - Volume Of Fluid), при помощи которого задачу сводят к рассмотрению системы из двух несмешивающихся фаз, где объемная доля фазы в ячейке определяется скалярной функцией. Однако, такой подход связан с проблемой восстановления поверхности раздела фаз по значениям скалярной функции.

Альтернативой эйлеровым сеточным являются лагранжевы бессеточные методы. Для решения задач динамики жидкости получил распространение метод гидродинамики сглаженных частиц (SPH — Smoothed particle hydrodynamics [2,3]), который обеспечивает выполнение законов сохранения при моделировании течений с развитой свободной поверхностью без введения каких-либо специальных алгоритмов, и успешно применялся авторами для моделирования взаимодействия потока жидкости с неподвижной преградой [4].

Решение широкого круга задач динамики твердых тел может быть получено бессеточным методом дискретных элементов (МДЭ [5]). Такой подход, в отличие от конечно-элементных методов, не требует заранее заданных контактных поверхностей и вращение тел требует меньших

1 Работа выполнена при частичной финансовой поддержке гранта РФФИ № 16-41-180276_р_урал_а.

вычислительных затрат. В МДЭ твердое тело представляется как агломерат жестко связных сфер. Моменты инерции и масса агломератов вычисляется как сумма за вычетом пересечений составляющих сфер. Новое положение вычисляется с помощью системы уравнений движения и уравнения моментов.

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

§ 1. Постановка задачи

Одной из тестовых задач, имеющих подробное экспериментальное описание является взаимодействие потока с незакрепленными твердыми телами [6]. Рассмотрим задачу о движении потока, вызванного распадом начального уровня жидкости.

Пусть имеется область Qf П 0д € Я3 размерами 8 х 0.7 х 1м (рис. 1), где в подобласти Qf € [3.5 х 0.7 х 0.4] м находится жидкость плотности р = 1000 кг/м3. Рассматриваются два тестовых случая, различающиеся только количеством (три и шесть соответственно) и расположением в Од незакрепленных блоков.

Схематично расположение блоков представлено на рисунке 1. Все блоки одинаковые и имеют размер 0.15 х 0.15 х 0.15 м. Основные характеристики материалов (Е — модуль Юнга, V — коэффициент Пуассона и к — коэффициент жесткости при контактном взаимодействии, ^ — коэффициент трения): для блоков — Е = 3 ■ 109 Па, V = 0.3, к = 0.6 кг/с2, ^ = 0.45; стенки бокса — Е = 6.5 ■ 1010 Па, V = 0.23, к = 0.85 кг/с2, ц = 0.4; основание бокса — Е = 2.1 ■ 1011 Па, V = 0.3, к = 0.8 кг/с2, ^ = 0.35.

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

2. Метод БРН

В методе гидродинамики сглаженных частиц [2,3] за основу берутся уравнения сохранения д (ри)

дЬ

+ V- (рии) = —Ур + V- т + р% + FFsч, (2.1)

V- и = 0, (2.2)

где р — плотность, и — вектор скор ости, р — давление, р = /л(Уи + Уит) — тензор вязких напряжений, у — динамическая вязкость, д — вектор ускорения свободного падения, FFSI — сила взаимодействия частиц жидкости и твердых тел. В методе частицы играют роль дискретных элементов, позволяющих свести уравнения динамики сплошной среды в частных производных к системе обыкновенных дифференциальных уравнений. Таким образом, БРН аппроксимация уравнений (2.1), (2.2) имеет вид:

- уравнение неразрывности

^ = " Е " из) ■ ~ гэ>Ь)> (2-3)

3

- уравнение движения

^ = - Е ^ (| + § + + / + &> (2-4)

где Ту — диссипативный член в форме Монагана [2], шг, рг — масса и плотность частицы жидкости, Ш(г — г',Н) — сглаживающее ядро (для решения рассматриваемой задачи использовалась функция ядра в виде кубического сплайна), г г ,гу — радиус-вектор частицы г и ], Н — радиус сглаживания. В этом случае новое положение частицы вычисляется по скоростям окружающих ее частиц

ж = ад

где е — константа.

Для расчета течений несжимаемой жидкости использовалась формулировка метода \¥СБРН, предполагающая замыкание системы уравнений (2.3)-(2.5) уравнением состояния Тейта

Pí > Ро,

Рг

рг ^ ро,

где В = росО— модуль объемного сжатия, ро — плотность жидкости на свободной поверхности (р = 0), со — константа, имеющая физический смысл скорости звука, 7 — константа (для рассматриваемых задач 7 = 7).

Для устранения нефизичного увеличения давления и плотности жидкости вблизи непроницаемой твердой поверхности и стабилизации свободной поверхности используется фильтр Шепарда.

р3

33

где Щ = (Шгзшз/рз)•

В методе БРН применяются различные типы граничных условий, в настоящей работе использовались динамические граничные частицы, способ учета и задание которых подробно рассмотрены в [4]. Динамические граничные частицы, как и частицы жидкости, удовлетворяют уравнениям сохранения, но при этом остаются неподвижными (либо движутся согласно функции движения границы). Особенностью данного вида задания границ является изменение плотности и давления в граничных частицах при приближении к ним частиц жидкости: когда жидкая частица приближается к границе, плотность граничных частиц увеличивается согласно

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

§ 3. Метод дискретных элементов

В методе дискретных элементов (МДЭ) тело, вместо сетки, представляется набором сфер разного диаметра (рис. 2), приближенно описывающих его поверхность [5]. Агломерат сфер P(dim, R) характеризуется двумя основными параметрами, это dim — размерность сетки, на которой строилась модель из объемных данных, и R — разреженность — значение, характеризующее расстояние между центрами сфер в агломерате, которое влияет на величину перекрытий в теле.

(а) (б)

Рис. 2. Агломерат сфер Р(64, 50): а) объемное представление: б) представление в разрезе

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

du v"^

mi~dt = ^ + SI + gmc'

I — - Т

' <11 "

где — масса агломе рата, — сила взаимодействия между частицам и агломерата, —

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

Выбор функции, описывающей силовое взаимодействие контактирующих сфер, является одним из важнейших в МДЭ. Сила взаимодействия между сферами играет такую же роль, что и определяющие уравнения в механике сплошной среды, однако, ее структура намного проще, чем у определяющих уравнений и представляет собой скалярную функцию расстояния. Силовое взаимодействие в точке контакта между сферами с номерами ^ ] представляется в векторной форме следующим образом:

е, = щ ■ Щ] + Е, ■ тч- = + ) ■+ + Ф) ■ т,

где Е, и Е, — нормальная и касательная компоненты контактной силы, каждая из которой состоит из вязкой Е и упруг ой Е *'е силы, их вычисление определяется в данном случае на основе контактной модели Герца-Миндлина [7]; п, — единичный вектор, определяющий плоскость контакта между двумя сферами, а именно вектор нормали к плоскости пересечения

сфер; Tij — единичный вектор, принадлежащий плоскости контакта. В работе рассматриваются случаи контакта с вязко-упругим взаимодействием контактирующих тел.

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

§ 4. Связывание МДЭ с SPH

Одной из проблем при вычислении взаимодействия между телами и жидкостью, в рамках бессеточных методов, является определение положения интерфейсной границы. В SPH расчетная частица определяется областью с радиусом 2h, h — длина сглаживания, а в МДЭ граница агломерата сфер — поверхностями элементов-сфер. Во избежание необходимости построения и использования трехмерной геометрии поверхности взаимодействия воспользуемся контактным потенциалом [8]. Взаимодействие между жидкостью и твердыми телами частицами определяется через введение дополнительного слагаемого, определяющую силу взаимодействия Ffs4, между жидкостью и твердыми телами, в уравнения движения для SPH и МДЭ. В алгоритме рассматривается условие выполнения взаимодействия для каждой пары расчетных частиц жидкости и твердого тела в пределах области поверхности тела. Для определения силы взаимодействия необходимо определить массу плотноеть pi и координаты Xi ,Xj расчетных частиц,

v^1 mj mi W (xi - Xj Г-1_ т„, ,

= Е77Ал" щь-гГ v"w(x'~

где Ffs4 — сила взаимодействия частиц г и j на поверхности Г^&ь W — ядро сглаживания, haw — осредненная длина сглаживания, K, n — параметры, определяющие вид ядра.

§ 5. Численное моделирование

В рамках верификации метода рассматривалась задача об обрушении столба жидкости с последующим его взаимодействием с несколькими твердыми телами, имеющая подробное экспериментальное описание [6]. Геометрические размеры расчетной области для двух рассмотренных вариантов эксперимента приведены на рисунке 1. Общее количество частиц, используемых в расчетах — 1 206 907, го них 240 907 для представления непроницаемых границ. Это обеспечило среднюю длину сглаживания havr = 0.021 м.

На рисунке 3 показано развитие течения жидкости после удаления перегородки. В момент времени t = 0.98 с поток достигает пирамиду блоков, в результате чего происходит смещение блоков и начинается их дальнейшее движение вместе с потоком. В момент времени t = 1.28 с становится хорошо заметно, что в численном эксперименте поток жидкости не смыкается, обтекая нижний блок, однако различие в динамике твердых тел незначительно: взаимное расположение блоков совпадает с наблюдаемым в эксперименте. Графики изменения положения центров масс блоков, приведенные на рисунке 4 также свидетельствуют о хорошем воспроизведении динамики тел в численном эксперименте, имеется лишь незначительное отставание в движении верхнего блока по оси x (рис. 4, б), для t > 2 с.

Далее рассмотрим эксперимент с шестью телами (рис. 6). Блоки расположены на расстоянии 0.005 м по оси у. Результаты вычислений показывают (рис. 6, a, t = 0.82 с), что при взаимодействии потока жидкости с твердыми телами, жидкость не проникает в узкие зазоры между блоками (рис. 6, a, t = 1.15с), что приводит уже к более существенным отличиям в динамике тел и потока жидкости. По-видимому, такое нефизичное поведение объясняется особенностями реализации граничных условий для непроницаемых поверхностей.

На рисунке 5 приведены графики изменения положения центра масс верхнего блока. В мо-t=1

0.075 мм по о си z, отсутствующий в физическом эксперименте, где блок не изменяет своего положения. Однако далее результаты расчета согласуются с физическим экспериментом.

г = 1.28с

а) б)

Рис. 3. Взаимодействие потока жидкости с тремя телами: а) расчет по представленной модели; б) физический эксперимент [6]

ж, м

1.5 г, с

ж, м 2.5

_ 1 2

5 ] 1.5

г, м 0.3 0.2 0.1

1.5 г, с

б)

с)

Рис. 4. Графики изменения положения центра масс в первом эксперименте: а) нижнего кубика по оси ж; б) и с) верхнего кубика по осям ж и г; 1 — физический эксперимент [6], 2 — численный расчет

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

г

а)

г

б)

ж

б) по оси г; 1 — физический эксперимент [6], 2 — численный расчет

t = 1.46c

a) 6)

Рис. 6. Взаимодействие потока жидкости с шестью телами: а) расчет по представленной модели; б) физический эксперимент [6]

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

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

1. Кузьмин И.М., Тонков Л.Е., Копысов С.П. Алгоритмическое и программное обеспечение решения задач взаимодействия конструкции с жидкостью/газом на гибридных вычислительных системах // Компьютерные исследования и моделирование. 2013. Т. 5. № 2. С. 153-164.

2. Monaghan J.J. Simulating free surface flows with SPH // Journal of Computational Physics. 1994. Vol. 110. P. 399-406. DOI: 10.1006/jcph.l994.1034

3. Gingold R.A., Monaghan J.J. Smoothed particle hydrodynamics: theory and application to non-spherical stars // Mon. Not. R. Astr. Soc. 1977. Vol. 181. P. 375-389. DOI: 10.1093/mnras/181.3.375

4. Копысов С.П., Тонков Л.Е., Чернова A.A., Сармакеева A.C. Моделирование взаимодействия с преградой потока несжимаемой жидкости методами VOF и SPH // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. 2015. Т. 25. Вып. 3. С. 405-420.

DOI: 10.20537/vml50311

5. Копысов С.П., Караваев A.C., Сармакеева A.C. Моделирование динамики произвольных тел методом дискретных элементов // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. 2015. Т. 25. Вып. 4. С. 473-482. DOI: 10.20537/vml50404

6. Canelas R.B., Crespo A.J.С., Dominguez J.M., Ferreira R.M.L., Gomez-Gesteira M. SPH-DCDEM model for arbitrary geometries in free surface solid-fluid flows // Сотр. Phys. Comm. 2016. Vol. 202. P. 131-140. DOI: 10.1016/j.cpc. 2016.01.006

7. Mindlin R.D. Compliance of elastic bodies in contact // Journal of Applied Mechanics. 1949. No. 16. P. 259-268.

8. Vignjevic R., De Vuyst Т., Campbell J.C. A frictionless contact algorithm for meshless methods // ICCES. 2007. Vol. 3. No. 2. P. 107-112.

Поступила в редакцию 01.10.2017

Сармакеева Анастасия Семеновна, младший научный сотрудник, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: asarmakeeva@gmail.com

Тонков Леонид Евгеньевич, к. ф.-м. н., заведующий лабораторией, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: tnk@udman.ru

Чернова Алена Алексеевна, к. т. н., старший научный сотрудник, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: alicaaa@gmail.com

A.S. Sarmakeeva, L.E. Tonkov, A. A. Chemova Meshfree methods for simulation fluid-structure interactions

Citation: Izv. Inst. Mat. Inform. Udmurt. Cos. Univ., 2017, vol. 50, pp. 36-44 (in Russian). Keywords: smooth particle hydrodynamic, discrete element method, solid body dynamic, fluid-structure interaction.

MSC2010: 76D27, 76M25

DOI: 10.20537/2226-3594-2017-50-05

The fluid-structure interaction algorithm (FSI) by meshless methods is presented. Dynamics of fluid is simulated by Smoothed Particle Hydrodynamic Method (SPH) and, for moving solid bodies, Discrete Element Method (DEM) is used. The interaction between the fluid and solid bodies is determined by additional contact forces in motion equations for SPH and DEM, respectively. The special condition of contact algorithm is considered for each pair of liquid and solid particles within the interface, the area. The features and calculation procedures of these numerical methods are described. During simulations, the dam break experiment was carried out. The fluid-body interactions are considered according to the boundary conditions introduced by dynamic boundary particles. Obtained results show that the algorithm gives good agreement with the physical experiment.

REFERENCES

1. Kuzmin I.M., Tonkov L.E., Kopysov S.P. Algorithms and software for solving coupled fluid-structure interaction problems on hybrid HPC platform, Komp'yuternye Issledovaniya % Modelirovanie, 2013, vol. 5, no. 2, pp. 153-164 (in Russian).

2. Monaghan J.J. Simulating free surface flows with SPH, Journal of Computational Physics, 1994, vol. 110, pp. 399-406. DOI: 10.1006/jcph.l994.1034

3. Gingold R.A., Monaghan J.J. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon. Not. R. Astr. Soc., 1977, vol. 181, pp. 375-389. DOI: 10.1093/mnras/181.3.375

4. Kopysov S.P., Tonkov L.E., Chernova A.A., Sarmakeeva A.S. Modeling of the incompressible liquid flow interaction with barriers using VOF and SPH methods, Vestnik Udmurtskogo Universiteta. Matematika. Mekhanika. Komp'yuternye Nauki, 2015, vol. 25, issue 3, pp. 405-420 (in Russian).

DOI: 10.20537/vml50311

5. Kopysov S.P., Karavaev A.S., Sarmakeeva A.S. Simulation of arbitrary solid body dynamics by discrete element method, Vestnik Udmurtskogo Universiteta. Matematika. Mekhanika. Komp'yuternye Nauki, 2015, vol. 25, issue 4, pp. 473-482 (in Russian). DOI: 10.20537/vml50404

6. Canelas R.B., Crespo A. J.C., Dominguez J.M., Ferreira R.M.L., Gomez-Gesteira M. SPH-DCDEM model for arbitrary geometries in free surface solid-fluid flows, Comp. Phys. Comm., 2016, vol. 202, pp. 131-140. DOI: 10.1016/j.cpc. 2016.01.006

7. Mindlin R.D. Compliance of elastic bodies in contact, Journal of Applied Mechanics, 1949, no. 16, pp. 259268.

8. Vignjevic R., De Vuyst T., Campbell J.C. A frictionless contact algorithm for meshless methods, ICCES, 2007, vol. 3, no. 2, pp. 107-112.

Received 01.10.2017

Sarmakeeva Anastasiya Semenovna, Junior Researcher, Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia. E-mail: asarmakeeva@gmail.com

Tonkov Leonid Evgen'evich, Candidate of Physics and Mathematics, Head of Laboratory, Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia. E-mail: tnk@udman.ru

Chernova Alena Alekseevna, Candidate of Engineering, Senior Researcher, Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia. E-mail: alicaaa@gmail.com