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

Вычислительные технологии

Том 11, Специальный выпуск, 2006

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ ЖИДКОСТИ СО СВОБОДНЫМИ ГРАНИЦАМИ МЕТОДАМИ SPH И MPS

К. Е. Афанасьев, А. Е. Ильясов. Р. С. Макарчук. А. Ю. Попов

Кемеровский государственный университет, Россия e-mail: [email protected], makOkemsu.ru

This work is devoted to numerical solution of the free boundary problems in hydrodynamics using Smoothed Particle Hydrodynamics (SPH) and Moving Particles Semi-Implicit (MPS). These methods allow simulation of phenomena that incorporate large deformations of the domain boundaries, as well as alteration of its connectedness and intersection of its boundaries. Underlining ideas of the methods are described, and the results of numerical calculations of some free boundary problems are presented.

Введение

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

К одному из наиболее сложных для моделирования классов задач относятся задачи о течении жидкости со свободными границами, сопровождающиеся эффектами сильно нелинейных деформаций этих границ. В качестве примеров можно привести: взаимодействие волн с препятствиями, накат волны на наклонный берег, деформацию пузырей в жидкости вблизи свободных и твердых поверхностей и т. д. Ряд такого рода задач (в двумерной, осесимметричной и пространственной постановках) был решен с использованием разновидностей методов граничных элементов [1]. Привлекательность этих методов обусловлена тем, что дискретизации подвергается лишь граница области. С целью реализации такой возможности требуется переход от исходной краевой задачи для дифференциальных уравнений, описывающих некоторый процесс, к соотношениям, связывающим неизвестные функции на границе области. Эти соотношения, по существу, представляют собой граничные интегральные уравнения. Вместе с тем, если возникает необходимость отыскания

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

решения в любой внутренней точке, то это можно сделать, используя известные значения функций на границе области, В результате использования методов, построенных на основе граничных интегральных уравнений, размерность расчетной области понижается на единицу. Существенным недостатком этих методов является то, что невозможно продолжать моделировать процесс распространения волн после изменения связности расчетной области. Например, обрушение волны можно рассчитать только до момента соприкосновения гребня волны с ее подошвой, далее проводить расчет становится невозможным в силу изменения связности области течения и перехлеста границ. Указанными недостатками обладают также и хорошо разработанные сеточные методы: конечных разностей, конечных элементов и контрольного объема,

В последнее время для решения задач механики жидкости со свободными границами разрабатываются методы, позволяющие проводить расчеты течений с сильными деформациями границ расчетной области, допускающие изменение связности области расчета и перехлест ее границ. Эти методы называются бессеточными, для их реализации не требуется информация о связях между узлами, что позволяет избежать проблем, относящихся к построению сетки, а также вызывающих необходимость отслеживать межузловые связи, К чисто бессеточным относятся метод сглаженных частиц (Smoothed Particle Hydrodynamics [2-4]), полунеявный метод движущихся частиц (Moving Particles Semi-implicit [5, 6]) и метод лагранжево-эйлеровых частиц [7], К гибридным методам относится бессеточный метод конечных элементов (Meshless finite element method [8]), Подробный обзор и классификацию бессеточных методов можно найти в обзоре [9],

1. Основные уравнения динамики жидкости

Движение ньютоновской вязкой жидкости в некоторой области течения О описывается уравнениями Навье — Стокса, уравнением сохранения внутренней энергии, а также уравнением неразрывности и каким-либо уравнением состояния p = p(p):

где а,Ь = 1, 2, 3 — числовые индексы координат; уа — компоненты вектора скорости; ра — компоненты вектора плотности объемных сил; и — внутренняя энергия; ТаЬ = дуа/дхЬ + дуЬ/дха — 2 div у5аЬ/3 — компоненты тензора вязких напряжений; ёаЬ — символы Кронекера; р и р — давление и плотность жидкости соответственно.

Основная идея метода ЭРН состоит в дискретизации сплошной среды набором частиц, которые движутся в потоке жидкости, заполняющей область О, причем используется интегральное представление характеристик течения в виде

(1)

(2)

dP А— = —р div и, dt

(3)

2. Метод SPH

а б

Рис. 1. Взаимодействие частиц (а) и функция ядра (б).

где 8 — дельта-функция Дирака, Далее вышеприведенный интеграл аппроксимируется интегралом вида

А(т) = ! А(т')Ш(г — т',Ь)Лт' п

с весовой функцией Ш, которую часто называют функцией ядра, а сами интегралы заменяются конечной суммой [10]:

П

Шг

А3{г) = УА-^Ш{г-п,Ь), (5)

1=1 Рг

где тг, шг, рг — радиус-вектор, масса и плотность г-й частицы соответственно; п — количество соседних к г-й частиц. Соседними к данной г-й частице называют такие частицы, расстояния от центров которых до центра данной частицы не превосходят радиуса. Кг. На рис, 1, а, показано, что частица 1 является соседней к г-й, а частица 2 — нет. Величина — носитель функции Ш, и она называется сглаживающей длиной, В качестве веса Ш обычно используют полиномиальные сплайны. Из (5) следует, что градиент искомой функции представляется формулой

П

УАДг) = V Аг—У\У(г - п, К).

Л Рг

г=1 '

Используя вышеприведенные предположения, уравнения (1)-(3) можно привести к следующему виду:

^=ЁЙ + | + Па)ПЬУ^(г,-г„Л1)+

' П

Шз

+ ££— [(и,Т}‘ + ПТ*) УИ-'^г. - г„ Л,)] , (6)

4=1 3=1 РгРп

^иг _ (граЬ\2 _ 2;

(И ~ 2 рг^ 1 } 2

йрг

V п

ЕЕ Шз (у4 — 4(тг — Т,, Лг),

4=1 3=1

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

П

ТГ = У — \щ - «,^Г(Г4 - Ы) + Ы)+

3=1 Рз 1

+ (Vе - - Гз, К)

п

Т-Ь = У— - у?)Ч]¥ь(гг - г,, Ы) + - уьг)Ч]¥а(гг - г,, Ы)} ,

з=1 Рз

причем а = Ь = с.

Функция ядра. Для расчета задач методом ЭРН применяется много различных видов функции ядра (рис, 1,6), начиная с гауссовой функции и заканчивая сплайнами различного порядка. Помимо уже известных функций ядра можно разрабатывать собственные, но при этом необходимо руководствоваться предъявляемыми к ним требованиями:

\¥(г, К) = 0, г > /г,

\У(г, К)(1г = 1,

О

lim W(r, К) = 8(г). h^0

Здесь г = \r - п\, q = f/h. Кроме того, на функцию ядра можно накладывать дополнительные условия, которые обеспечивали бы лучшую устойчивость метода и более высокую степень аппроксимации функций, характеризующих течения. Такие дополнительные условия и вытекающие из них способы построения функций ядра привели, например, к появлению метода RKPM (Reproducing Kernel Particle Method) [11],

Для задач, рассматриваемых в данной работе, применяется кубический сплайн Монагана [12]:

Г 2/3 - q2 + q3/2, 0 < q < 1;

W(r,h) = —-гз { (2 — g)3/6, 1 < q < 2;

7nh2\ 0, q > 2.

3. Модификации метода сглаженных частиц

Искусственная вязкость. Для лучшей устойчивости метода SPH в правую часть формулы (6) добавляют дополнительный член, называемый искусственной вязкостью [13]:

п _ J ( — acij &ij + )/pij, vij Vij < 0;

ij — 1 0, Vij Vij < 0,

где aij = hvijrij/(rij + nf); cij = (ci + cj)/2 — средняя скорость звука в частицах і и j; П2 = 0.001h2; vij = vi - vj и rij = Гі - fj.

Модель несжимаемой жидкости. Оригинальный метод 81*11 пригоден только для моделирования сжимаемых сред, однако существуют некоторые модификации метода, позволяющие применить его и для случая слабосжимаемой среды, “Несжимаемость” в методе ЭРН моделируется при помощи выбора подходящего уравнения состояния, В рассматриваемых задачах скорости частиц много меньше скорости звука в воде. Для расчета течений жидкости чаще всего используется уравнение состояния в форме Тета [13, 14]:

р = В

Здесь модуль объемного сжатия:

В

-1

Ро

2

РоСо

1 '

где р0 — начальная плотность среды; с0 — начальная скорость звука; 7 = 7 — показатель адиабаты,

В

ле [13]

200родЯ

В

7

где д — ускорение свободного падения, аЯ — начальная высота столба жидкости.

Условия на твердой границе. В методе сглаженных частиц для постановки условий на твердой границе наиболее часто используется метод виртуальных частиц. Существует несколько способов задания виртуальных частиц. Виртуальные частицы Монагана располагаются вдоль твердой границы в один слой (рис, 2, о), в них не задаются никакие характеристики. Виртуальные частицы взаимодействуют с частицами среды посредством некоторого потенциала. Выбранный потенциал взаимодействия включается в уравнения движения в качестве дополнительной силы, В методе сглаженных частиц используется потенциал Леппарда — Джонса [13]:

и(г) = —

г

Го

12

— I + ( —

ГГ

Го

где В — глубина потенциальной ямы; г0 — радиус взаимодействия. Потенциал является двухпараметрическим, что позволяет независимо задавать два каких-либо параметра среды.

6

Рис. 2. Виртуальные частицы Монагана (а) и Морриса (б).

Виртуальные частицы Морриса [15] располагаются не только вдоль твердой границы, но и внутри нее (рис, 2,6), Они включаются в суммирование при нахождении характеристик свободных частиц и обладают всеми их характеристиками, однако не изменяют своих координат. Здесь возможны два варианта — виртуальные граничные частицы могут изменять свои характеристики (плотность, давление, скорость), а могут сохранять первоначальные значения.

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

Сглаживающая длина, переменная по времени. Сглаживающая длина играет значительную роль в методе SPH, Для устойчивости метода необходимо, чтобы у каждой частицы в каждый момент времени количество соседей было постоянным. При значительных деформациях такое условие может не выполняться, тогда для приближенного удовлетворения этого условия в алгоритм метода SPH включается уравнение для изменения сглаживающей длины со временем [17]:

dhi 1 hi dpi

dt v pi dt ’

где v — размерность задачи.

Для модифицированных методов частиц, предназначенных для расчета течений несжи-

h

dhi 1

—— = —hi div v. dt v

XSPH (extended SPH) модификация. При решении задач со свободными границами часто приходится прибегать к формуле коррекции скорости [13], которая позволяет упорядочивать расположение частиц и помогает избежать их столкновений:

pi + pj

4. Интегрирование системы уравнений

Для интегрирования полученных обыкновенных дифференциальных уравнений исполь-

зуется схема типа “предиктор-корректор” [10]

Предиктор

Корректор

pl-1/2 + (At/2)(dp“-1/dt); vl-1/0 + (At/2)(<fo"-7dt).

pl+1/2 = pl-1/2 + At(dp”/dt); vl+1/° = v1-1/° + At(dvn/dt); xT1 = x; + At(v“+1/0/dt).

n

p

Для первого шага по времени

Рг/2 = Р°г + (Л*/2)ММ); у]/2 = V0 + (Д^/2)(<^0 /(Щ; х1 = х0 + Д^(^1 /2/<И).

Шаг по времени выбирается исходя из условия Куранта — Фридрихса — Леви [15]:

шт(Кг)

Д£ < С------7“---гг—г^-, (7)

шах(е* + || V* ||)

г

где Кг, ег, V* — сглаживающая длина, скорость звука и скорость г-й частицы соответственно, Константа С € (0,1), устойчивый счет наблюдается при С = 0.3,

5. Тестирование метода

Для тестирования метода выбраны задача о ламинарном течении в плоском канале (течение Пуазейля) с твердыми границами и задача Овсянникова о деформации жидкого эллипса со свободной границей.

Для обеих задач использовались следующие коэффициенты: р0 = 1000 кг/м3 (начальная плотность частиц); е0 = 1440м/с (начальная скорость звука); а = 0.1 в = 0 (константы из искусственной вязкости).

5.1. Задача о ламинарном течении в плоском канале

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

■ результат расчета — аналитическое решение

В области О течение жидкости описывают уравнения движения и неразрывности:

dvl _ Pout ~ Рт /л<9У dt pL p dy2 ’

dP A— = —p cliv V,

dt H

где Pin, Pout — давление, заданное та входе и выходе канала; р L — плотность, коэффициент динамической вязкости и длина канала соответственно. Твердые границы моделируются виртуальными частицами Морриса, Бесконечность канала имитируется iiiik.iii-ческим возвращением частиц, вышедших за пределы правой открытой границы, на левую

О

собой квадрат со стороной 0.1м. Скачок давления: Pout — Pin = 3.75 ■ 10-7кг/(м-с2); коэффициент динамической вязкости ^ = 10-3кг/(м-с).

Тестирование задачи проводилось при разном количестве частиц (625, 900, 1600) и различных шагах по времени (At = 10-5, At = 5 ■ 10-5, At = 10-4, At = 5 ■ 10-4), Удовлетворительные результаты достигнуты при использовании 1600 расчетных частиц с шагом At = 10-4

линией — аналитическое решение. Максимальное отклонение численного и аналитического решений не превышало 0.8%.

5.2. Задача о деформации жидкого эллипса

Для тестирования метода в случае решения задач со свободными границами выбрана задача Овсянникова [17] о деформации жидкого эллипса, принадлежащая классу задач невесомой жидкости со свободными границами. Задача формулируется следующим образом: в начальный момент времени область расчета представляет собой круг радиуса 1, в который заключена несжимаемая жидкость. Деформация круга в эллипс начинается под действием начального распределения скоростей в отсутствие внешних сил. Для обеспечения несжимаемости необходимо, чтобы площадь эллипса оставалась постоянной, т, е, аЬ = 1 во все время расчета, где а, Ь — полуоси эллипса,

В любой момент времени £ уравнение эллипса имеет вид

x

a2 (t)

+ a2 (t)y2 = 1.

y, M

1.5 1

0.5

0

-0.5

-1

-1.5

y, m

y, m

........................................................................x, M

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

Полуось а(£) находится из решения системы обыкновенных дифференциальных уравнений:

а' = е, С = 2е2/(а5 + а) (8)

с начальными условиями

а(0) = 1, е(0) = 100.

Начальное распределение скоростей для решения задачи методом 81*11:

VI (0) = хе(0)/а(0), г>2(0) = —уе(0)/а(0).

Система уравнений (8), (9) решалась методом Рунге — Кутты четвертого порядка для получения свободной границы эллипса. Результаты сравнения расчетов методом 81*11 с решением по методу Рунге — Кутты приведены на рис, 4,

Расчеты задачи методом 81*11 проводились с различным начальным распределением частиц. На рис, 4 приведены результаты расчетов с использованием нерегулярного распределения, которое гарантирует для каждой частицы одинаковое количество соседних с ней частиц, распределенных в различных направлениях. При использовании регулярного распределения (частицы находятся в узлах прямоугольной сетки) для получения удовлетворительного результата необходимо использовать модификации метода для функций ядра с эллиптическим носителем.

Из приведенных рисунков видно, что результаты расчетов методом 81*11 хорошо согласуются с результатами, полученными методом Рунге — Кутты, Отклонение численного и аналитического решений не превышало 0,5%,

6. Решение модельных задач

Для расчета модельных задач использовались следующие значения параметров: ^ = 10-3кг/(м-с) (коэффициент динамической вязкости); р0 = 1000кг/м3 (начальная плотность частиц); е0 = 1440м/с (скорость звука); а = 0.1 в = 0 (константы для искусственной вязкости); г0 = <х (радиус взаимодействия для сил Леппарда — Джонса), где <х — начальное расстояние между частицами; В = 0.049 (глубина потенциальной ямы). Выбор шага по времени осуществляется согласно условию (7), Для постановки условий на твердых границах использовались виртуальные частицы Монагана,

6.1. Обрушение столба жидкости

В данной задаче столб вязкой жидкости с нулевым начальным вектором скорости в момент времени £ = 0 начинает обрушаться под действием силы тяжести. На рис, 5 приведены картины течений для различных моментов времени с количеством расчетных частиц, равным 900,

Для сравнения на рис, 6 и 7 показаны фрагменты течений для 625 и 1600 частиц соответственно для двух моментов времени. Из приведенных иллюстраций видно хорошее качественное совпадение результатов.

у,м £=0 у,м £ = 0.3с у,м

2.5 2.5 2.5:

2 2 2:

1.5 1.5 1.5:

1 1 1:

£ = 0.7 с

0.5

0

0 0.5 1 1.5 2 2.5 3

£ = 1.3с

у,м

2.5 2

1.5 1

0.5

0

0 0.5 1 1.5 2 2.5 3

0.5

х,м 0

у, и

2.5 2

1.5 1

0.5

х.м 0

ЯВЬ*.

-.......

0 0.5 1 1.5 2 2.5 3

£ = 1.4с

0 0.5 1 1.5 2 2.5 3

0.5 0

у,и

2.5 2

1.5 1

0.5

хм 0

0 0.5 1 1.5 2 2.5 3

£ = 1.6с

х,м

0 0.5 1 1.5 2 2.5 3

Рис. 5. Обрушение столба жидкости (Ж = 900).

у,м

2.5 2

1.5 1

0.5

0

£ = 1.3с

у,м

2.5

2

1.5

1

.. 5И 0.5

1«Д—«я хм 0

£ = 1.6с

0 0.5 1 1.5 2 2.5 3

0 0.5 1 1.5 2 2.5 3

х,м

Рис. 6. Обрушение столба жидкости (Ж = 625).

На рис, 8 представлены картины течений для 900 частиц с постоянным шагом по времени Д£ = 10-4 с. Хорошо видно также качественное совпадение результатов и их устойчивость.

£=0

£ = 0.3с

£ = 0.7 с

у,м

2.5 2

1.5 1

0.5

0

0 0.5 1 1.5 2 2.5 3

у,м

2.5 2

1.5 1

0.5

0

0 0.5 1 1.5 2 2.5 3

у,м

2.5 2

1.5 1

0.5

0

0 0.5 1 1.5 2 2.5 3

£ = 1.3с

£ = 1.4с

£ = 1.6с

у,м

2.5 2

1.5 1

0.5

0

0 0.5 1 1.5 2 2.5 3

у,м

2.5 2

1.5 1

0.5

0

0 0.5 1 1.5 2 2.5 3

у,м

2.5 2

1.5 1

0.5

0

0 0.5 1 1.5 2 2.5 3

Рис. 8. Обрушение столба жидкости (Д£ = 10 4 с).

у, м 0.1

0.08

0.06

0.04

0.02

0

£ = 0.0018 с

0 0.02 0.04 0.06 0.08 0.1 0.12

^ х,м

у, м 0.1 Ь

0.08

0.06

0.04

0.02

0

£ = 0.0099 с

0 0.02 0.04 0.06 0.08 0.1 0.12

1^. х,м

у, м 0.1

0.08

0.06

0.04

0.02

0

£ = 0.0297 с

0 0.02 0.04 0.06 0.08 0.1 0.12

£ = 0.0477 с

£ = 0.0837 с

£ = 0.1017с

у, м 0.1

0.08

0.06

0.04

0.02

0

у, м 0.1

0.08

0.06

0.04

0.02

0

0 0.02 0.04 0.06 0.08 0.1 0.12

. х,м

у, м 0.1

0.08

0.06

0.04

0.02

0

х.м

х.м

Рис. 10. Падение капли: а результат авторов; б результат работы [18].

6.2. Падение капли в резервуар с жидкостью

Задача формулируется следующим образом: круглая капля диаметром 0.01 мм падает с начальной скоростью 2 м/с в резервуар, наполненный жидкостью. Начальная плотность капли и жидкости в резервуаре р0 = 1000кг/м3, коэффициент динамической вязкости ^ = 10-3кг/(м-с). Проведена серия расчетов для различного количества частиц и разных временных шагов. На рис. 9 приведены картины течений для различных моментов времени. Жидкость в резервуаре содержала 5400 частиц, капля состояла из 139 частиц. На рис. 10 приведены фрагменты течений для момента времени і = 0.0297 с и результаты работы [ 181. Видно качественное сходство результатов, однако наличие количественного различия говорит о необходимости использования дополнительных модификаций метода, в частности применения метода СЭР [19|.

7. Метод М Р8

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

Далее следует стандартная для метода ЭРН процедура замены 5-функции функцией ядра Ш(г, Ь). а также замены интеграла конечной суммой [5|:

П

(г - Г, К)

Мг) = Ч------------------, (10)

^2^(г — гг, К)

г=1

где г — радиус-вектор; п — число соседних частиц; К — носитель функции ядра (размер ядра).

Функция ядра должна удовлетворять следующим требованиям:

IV(г, Н) = 0, г > /і,

Иш IV(г, Н) = 8(г).

Ь.^0

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

П

по формуле

П* = У W(г* - Г, , Л). (11)

3=1

В качестве функции ядра обычно используется следующая функция [5, 6]:

л*т<- 1Л Г Ь/г—!, О <г<1г; ,1ГЛ

И''(г,Л) = (0, г>Л. <12>

Выражение для аппроксимированной плотности в методе XII’8 записывается как [6]

рг = тЫ (г* ), (13)

где N — среднее количество частиц в единице объема; масса частицы т в рамках метода

считается постоянной. Приближенное значение N можно получить, используя формулу

П

W(г* - Vу , Л*)

= 1}У(г)(1г ' ^

п

Таким образом, из (11), (13) и (14) получается следующее выражение для аппроксимации плотностей:

рг = т Г . (15)

] W(г)йг п

Из (15) вытекает линейная зависимость плотности от количественной плотности частицы, а, так как в силу несжимаемости жидкости плотность в частице остается постоянной, значит, и количественная плотность должна сохраняться. Эта константа обозначается через П0,

7.1. Аппроксимация градиента функции

В методе МРЭ вектор градиента функции Р в частице г находится как взвешенное среднее величины (Р, — Р)(Г3 — гг)/|гз — г*|2 по всем соседним к г-й частицам, с использованием формулы аппроксимации (10):

По

3=г

Р — Р

3 \(г, - ri)W{ri - Гз,Ы)

|г3 — Гг\

где й — размерность пространства.

Отличие способа аппроксимации градиента функции в методе X11 * 8 от метода 81*11 заключается в том, что он основывается на разности между значениями функции в частицах и не зависит от градиента функции ядра, К тому же в методе 81*11 возможна группировка частиц, так как, начиная с некоторого расстояния между частицами, силы отталкивания, возникающие из-за градиента давления, уменьшаются, В X11 *8 этот эффект невозможен.

7.2. Аппроксимация оператора Лапласа

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

2 (] 1

= -т----------------------------------------------------------------------------------^(п - Гз, /г*)/*, (16)

А Ні

где

/ Ш(т)т2Лт А п

Перенос от частицы і к і:

2с1 1

Ь

А= -г—{гі - Гі, ^-)Л- (17)

А По

Согласно (16), (17), полное изменение физической характеристики в частице г за время Аъ представляется следующим выражением:

А/г = У (А/л^г - А/г^') = Г■?’ ~~ (18)

Если в качестве функции ядра выбрана функция (12), то значение параметра А следует вычислять в зависимости от размера ядра по следующим формулам: для “одномерного случая”

А = -/г2,

6 ’

для “•двумерного случая”

А = 0.3Л2.

А

нет необходимости, так как функция Гаусса является весовой в теоретической модели диффузии.

7.3. Модель несжимаемости

Условие несжимаемости (р/(Ъ = 0 эквивалентно сохранению количественной плотности частиц щ = п0. Это достигается путем вычисления давлений при помощи уравнения Пуассона, Модифицированное для метода МРЭ уравнение выглядит следующим образом:

= (19) Аъ2 п0

Его левая часть расписывается согласно модели аппроксимации оператора Лапласа, В итоге получается система линейных уравнений, решая которую, можно найти давления в частицах. Затем полученные значения используются для нахождения градиента давления и вычисления ускорения в частицах.

7.4. Граничные условия

Если для частицы выполняется неравенство щ < Рщ, где Р — эмпирический параметр, значения которого лежат в интервале от 0,95 до 0,99, то эта частица принадлежит свободной поверхности. Для давления в этой частице выставляется значение 0 в системе линейных алгебраических уравнений, полученной из уравнения (19), в качестве условия Дирихле, Твердая стенка моделируется тремя слоями неподвижных частиц Морриса, Если раскрыть левую часть уравнения Пуассона (19) согласно формуле аппроксимации оператора Лапласа (18), то получится система линейных алгебраических уравнений:

Матрица системы (20) является самосопряженной (симметричной), а также разреженной, в силу того что функция ядра имеет конечный носитель.

8. Решение модельных задач

В расчетах использовались следующие параметры: шаг по времени Дг = 10-3 с; коэффициент динамической вязкости ^ = 10-3кг/(м-с).

г = 0 г 0. 3 О г = 0.7с

У,Ы у,м у,м

3 3 3

2.5 2.5 2.5

2 и 2 2

I

1.5 1.5 1.5

1 1 1

0.5 0.5 0.5

0 X, М 0 X, М 0 х, м

0 0.5 1 1.5 2 2.5 3 3.5 0 0.5 1 1.5 2 2.5 3 3.5 0 0.5 1 1.5 2 2.5 3 3.5

г = 1.3с г = 1.4с г = 1.6 с

У,ы у,м у,м

3 3 3

2.5 2.5 2.5

2 ... 2 ! 2

1.5 1.5 1.5

1 1 1

0.5 0.5 0.5

0 \,х, М 0 х, М 0 х, м

0 0.5 1 1.5 2 2.5 3.5 0 0.5 1 1.5 2 2.5 3 3.5 0 0.5 1 1.5 2 2.5 3 3.5

Рис. 11. Обрушение столба жидкости (Ж = 900).

у,м

3

2.5

г = 1.3с

у,м

3

2.5

2 11 = = 2

1.5 1 ' . 1 | 1.5

1 Н | 1

0.5 | 0.5 1

0 1 Х М 0

г = 1.6с

0 0.5 1 1.5 2 2.5 3 3.5

0 0.5 1 1.5 2 2.5 3 3.5

Рис. 12. Обрушение столба жидкости (Ж = 625).

у,м

3

2.5 2

1.5 1

0.5

0

г = 1.3с

у,м

3

2.5 2

1.5 1

0.5

0

г = 1.6с

0 0.5 1 1.5 2 2.5 3 3.5

0 0.5 1 1.5 2 2.5 3 3.5

Рис. 13. Обрушение столба жидкости (Ж = 1600).

г = 0

г = 1.3с

у,м

3

2.5 2

1.5 1

0.5

г = 0.3с

у,м

3

2.5 2

1.5 1

0.5

0

0 0.5

Ж

1.5 2 2.5

1.4

у,м

3

2.5 2

1.5 1

0.5

0

г = 0.7с

3 3.5

у,м

3

2.5 2

1.5 I

0.5

01

0 0.5

1.5 2

1.6

з 3.5

0 0.5 1 1.5 2 2.5 3 3.5

у,м

3

2.5 2

1.5 1

0.5

0

0 0.5 1 1.5

8.1. Обрушение столба жидкости

Формулировка задачи аналогична формулировке, использованной для решения задачи методом ЭРН, На рис, 11 приведены результаты расчетов для разных моментов времени с количеством частиц, равным 900, Для сравнения па рис, 12 и 13 показаны фрагменты течений для моментов времени г = 1.3 и г = 1.6 с с использованием 625 и 1600 расчетных частиц соответственно. Также для сравнения па рис, 14 представлены результаты расчетов для 900 частиц с шагом по времени Дг = 10-4 с.

Приведенные картины течений хорошо согласуются с результатами расчетов методом ЭРН, что позволяет говорить о высоких возможностях обоих методов, а сходство результатов расчетов с различным количеством частиц (удовлетворительные результаты получаются при использовании 900 частиц) говорит о сходимости данных методов.

8.2. Падение капли в резервуар с жидкостью

Проведена серия расчетов для разного количества частиц и различных шагов по времени. Использовались те же значения величин, что и для решения данной задачи методом ЭРН, результаты которого приведены выше. На рис, 15 показаны картины течений в различные моменты времени. На раппих этапах счета результаты хорошо согласуются с результатами, полученными методом ЭРН, а также приведенными в работе [ 191. Однако далее виден неустойчивый характер счета методом МРЭ, что оставляет необходимость дальнейшего глубокого анализа алгоритмов и их возможной модификации.

у,м

0.1

0.08

0.06

0.04

0.02

0

г = 0.0018с

г = 0.0099 с

у, м У, м ;

0.1 0.1:

0.08 0.08:

0.06 0.06:

1 А 1 0.04 0.04.

0.02 1 " I 0.021

х, м 0 х, м 01

г = 0.0297с

0 0.02 0.04 0.06 0.08 0.1 0.12

г = 0.0477с

О 0.02 0.04 0.06 0.08 0.1 0.12

г = 0.0837с

О 0.02 0.04 0.06 0.08 0.1

г = 0.1017с

0.12

у,м 0.1 у, М 0.1 у, м 0.1

0.08 0.08 0.08

0.06 0.06 0.06

0.04 0.02 У I 0.04 0.02 | ' 42 О О 0. 0.

0 X, М 0 1, ,щ ,р, р, м 0

'

0 0.02 0.04 0.06 0.08 0.1 0.12

х, м

Заключение

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

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

[1] Afanasiev К.Е., Stukolov S.V. Simulation of problems with free surfaces by a boundary element method // Computational Science and High Performance Computing. Series: Notes on Numerical Fluid Mechanics and Multidisciplinary Design (NNFM). 2005. Vol. 88. P. 307-338.

[2] Lucy L.B. A numerical approach to the testing of fusion process // Astronom. J. 1977. Vol. 82, N 12. P. 1013-1024.

[3] Gingold R.A., Monaghan J.J. Smoothed particle hydrodynamics: theory and application to non-spherical stars // Monthly Notices of the Royal Astronom. Society. 1977. Vol. 181. P. 375-389.

[4] Vila J.P. On particle weighted methods and smooth particle hydrodynamics // Math. Models and Methods in Appl. Sci. 1999. Vol. 9, N 2. P. 161-209.

[5] Koshizuka S., Tamako H., Oka Y. A particle method for incompressible viscous flow with fluid fragmentation // Comp. Fluid Dynamics J. 1995. P. 29-46.

[6] Koshizuka S., Nobe A., Oka Y. Numerical analysis of breaking waves using the moving particle semi-implicit method // Intern. J. for Numer. Methods in Fluid. 1998. P. 751-769.

[7] Франк A.M. Дискретные модели несжимаемой жидкости. М.: Физмат.пп. 2001.

[8] Idelsohn S., Onate Е., Del Pin F. A lagrangian meshless finite element method applied to fluid structure interaction problems // J. Computer and Structures. 2003. Vol. 81. P. 655-671.

[9] Liu G.R. Mesh Free Methods: Moving Beyond the Finite Element Method. CRC Press, 2003.

[10] Liu G.R., Liu M.B. Smoothed particle hydrodynamics: a meshfree particle method // World Sci., 2003.

[11] Liu W.K., Jun S., SlHLlNG D.T. ET al. Multiresolution reproducing kernel particle method for computational fluid dynamics // Intern. J. of Numer. Method in Fluids. 1997. Vol. 24. P. 1-25.

[12] Monaghan J. J. Smoothed particle hydrodynamics // Ann. Rev. Astron. and Astrophysics. 1992. Vol. 30. P. 543-574.

[13] Monaghan J.J., Thompson M.C., Hourigan K. Simulation of free surface flows with SPH // J. of Comput. Phys. 1994. Vol. 110. P. 399-406.

[14] КОУЛ. P. Подводные взрывы. М.: Изд-во иностр. лит., 1950.

[15] MORRIS J.P., Fox P.J., Zhu Y. Modeling low Reynolds number incompressible flows using SPH // J. of Comp. Phys. 1997. Vol. 136. P. 214-226.

[16] Benz W., Bowers R.L., Cameron A.G.W., Press W.H. Dynamic mass exchange in doubly degenerate binaries // Astrop. J. 1990. Vol. 348. P. 647.

[17] Овсянников Л.В. Общие уравнения и примеры: Задачи о неустановившемся движении жидкости со свободной границей. Новосибирск: Наука, 1967.

[18] Cueto-Felgueroso L., Colominas I., Mosqueira G. et al. Fluid flow simulation by Lagrangian particle methods.

https://dspace.udc.es/bitstream/2183/410/l/2004_-ECCOMAS_-SPH.pdf

[19] MORRIS J.P. Simulating surface tension with smoothed particle hydrodynamics // Intern. J. for Numerical Methods in Fluids. 2000. Vol. 33, N 3. P. 333-353.

[20] SlJEYOSHl М., Naito S. A numerical study of violent free surface problems with particle method for marine engineering // 8th Intern. Conf. on Numer. Ship Hydrodynamics. Busan, Korea, 2003.

Поступила в редакцию 20 февраля 2006 г., в переработанном виде — 4 апреля 2006 г.