Моделирование течений вязкой несжимаемой жидкости со свободными границами методом ISPH c использованием схемы расщепления

МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ СО СВОБОДНЫМИ ГРАНИЦАМИ МЕТОДОМ ISPH C ИСПОЛЬЗОВАНИЕМ СХЕМЫ РАСЩЕПЛЕНИЯ

Кемеровский государственный университет

*****@***ru

Введение

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

Наиболее распространенными на данный момент методами частиц являются метод сглаженных частиц (SPH) [1-4] и полунеявный метод движущихся частиц (MPS) [5-6]. Существует также множество модификаций метода SPH, в частности, методы RKPM [7] и MLSPH [8], призванные улучшить его аппроксимационные характеристики.

В данной работе рассматривается еще одна модификация метода SPH – метод сглаженных частиц для несжимаемых сред (ISPH) [9-10], позволяющий, в отличие от классического метода SPH, в точности удовлетворить условию несжимаемости. Наиболее значимые отличия двух методов приведены в таблице 1.

Таблица 1 – Сравнение методов SPH и ISPH

	SPH	ISPH
Схема интегрирования	Явная	Неявная
Давление	Уравнение состояния	Уравнение Пуассона
Условие Куранта	Скорость звука	Максимальная скорость частиц
Плотность	Переменная	Постоянная

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

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

; (1)

, (2)

где – числовые индексы координат, – компоненты вектора скорости, – компоненты вектора плотности объемных сил, – символы Кронекера, и – давление и плотность жидкости соответственно, – коэффициент динамической вязкости, а компоненты тензора напряжений вычисляются по формуле:

(3)

Следует отметить, что классический метод SPH применяется при моделировании движений сжимаемых сред, а для замыкания уравнений (1)-(2) используется некоторое уравнение состояния. Для расчета течений несжимаемых жидкостей чаще всего используется баротропное уравнение состояния в форме Тета [1]. Подбором коэффициента объемного расширения в уравнении состояния достигается эффект несжимаемой среды. Тем не менее, имеют место вариации плотностей частиц.

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

Аппроксимация функций

Первым этапом в процессе построения аппроксимационных формул для функций, входящих в уравнения динамики сплошной среды, является точное представление этих функций в виде интеграла:

, (4)

где – дельта-функция Дирака.

Далее следует замена -функции некоторой классической функцией с компактным носителем, что позволяет получить интегральную формулу аппроксимации функции по ограниченной области [3-4]:

, (5)

с весовой функцией , которая, по терминологии метода, называется функцией ядра.

На следующем этапе область расчета подвергается дискретизации набором лагранжевых частиц, а интеграл (5) заменяется квадратурой:

, (6)

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

Используя для аппроксимации плотности формулу (6) получаем:

(7)

Следует отметить, что в методе ISPH, в отличие от классического метода сглаженных частиц, в котором значение плотности в частицах чаще всего находится путем решения аппроксимированного аналога уравнения (2), плотность вычисляется именно по формуле (7).

Интегрирование уравнений движения

Для интегрирования по времени уравнений движения используется схема типа «предиктор-корректор». Первоначально рассматривается действие на частицы лишь массовых сил и сил вязкости, что позволяет получить предварительные значения плотностей, скоростей и координат частиц. На данном этапе давление в расчет не принимается, а в формуле (3), в силу несжимаемости среды, исчезает член с дивергенцией скорости. Указанные ограничения приводят к уравнениям:

; (8.1)

; (8.2)

. (8.3)

Здесь – изменение скорости за счет массовых сил и сил вязкости, и – предварительные значения вектора скорости и радиус-вектора центров частиц соответственно, и – значения вектора скорости и радиус-вектора центров частиц на -ом шаге по времени.

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

; (9.1)

; (9.2)

. (9.3)

На этом этапе, как видно из уравнения (9.1), возникает потребность в вычислении давления, что приводит к необходимости решения уравнения Пуассона, которое имеет вид:

. (10)

Вывод этого уравнения из уравнений (1)-(2) приведен в [10]. Данное уравнение приводится к системе линейных алгебраических уравнений, а условия на свободной границе, в качестве которых выступает условие Дирихле, внедряются в матрицу коэффициентов получившейся системы по тому же принципу, что и в методе конечных элементов [12].

Уравнение Пуассона на давление

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

(11)

Формула для аппроксимации сил вязкости получается аналогичным образом и принимает вид:

(12)

Постановка условий на твердой границе

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

Первый – виртуальные частицы Монагана [1]. Эти частицы располагаются вдоль твердой границы в один ряд, не меняют своих координат с течением времени и действуют на частицы подвижной среды посредством какого-либо потенциала взаимодействия. Наибольшей популярностью среди исследователей пользуется потенциал Леннарда-Джонса, хотя его выбор не является обязательным. В работе [3] Монаган предлагает новый потенциал взаимодействия, учитывающий особенности метода SPH.

Второй тип – виртуальные частицы Морриса [2]. Эти частицы располагаются вдоль твердой границы в несколько слоев. Количество слоев зависит от сглаживающей длины частиц подвижной среды. Это позволяет решить одну из главных проблем метода SPH – несимметричность функции ядра вблизи границ. Действие частиц Морриса на частицы подвижной среды отличается от действия частиц Монагана тем, что нет необходимости использовать какой-либо потенциал взаимодействия. Вместо этого, в частицах Морриса вычисляются значения характеристик исходя из их значений в частицах подвижной среды.

В методе ISPH для постановки условий на твердых границах используются виртуальные частицы Морриса. В работе [10] применяется подход, использованный при решении задач методом MPS [5-6], с 5 слоями частиц. На твердых границах для решения уравнения Пуассона на давление ставится условие Неймана, означающее равенство давлений в частицах, расположенных вдоль нормали к твердой границе. В данной работе используется 3 слоя частиц. Решение уравнения Пуассона, а также и другие вычисления проводятся только с использованием внутреннего слоя граничных частиц, остальные слои нужны лишь для вычисления плотностей как граничных так и внутренних частиц подвижной среды по формуле (7).

Постановка условий на свободной границе

Для идентификации частиц, принадлежащих свободной границе, можно применять различные способы. Одним из таких способов является использование алгоритма Дильтса [13], основанного на том, что каждая частица имеет размер, который в методе ISPH определяется сглаживающей длиной. Однако чаще всего, как в методе MPS, так и в методе ISPH, частицы, принадлежащие свободной границе, идентифицируются по их плотности. Данный подход основан на том, что плотность частиц, принадлежащих свободной границе, меньше плотности внутренних частиц среды в силу недостатка соседей с одной стороны границы.

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

Решение модельной задачи. Задача об обрушении столба жидкости

Решаются уравнения (1)-(2). В начальный момент времени столб вязкой жидкости начинает обрушаться под действием силы тяжести. Для задачи использовались следующие значения физических характеристик: – плотность жидкости, – коэффициент динамической вязкости. На рис. 1 приведены картины течений в различные моменты времени с расчетным количеством частиц – 1357. Численные результаты, полученные при решении данной задачи методами SPH и MPS, приведены в работе [14].

а)б)в)

г)д)е)

Рисунок 1 – Задача об обрушении столба жидкости: а) ; б) ; в) ; г) ; д) ; е)

Заключение

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

1.  Monaghan, J. J. Simulation of free surface flows with SPH [Текст] / J. J. Monaghan, M. C. Thompson, K. Hourigan // Journal of computational physics. – 1994. – № 000. – P. 399-406.

2.  Morris, J. P. Modeling low Reynolds number incompressible flows using SPH [Текст] / J. P. Morris, P. J. Fox, Y. Zhu// Journal of computational physics. – 1997. – № 000. – P. 214-226.

3.  Monaghan, J. J. Smoothed Particle Hydrodynamics [Текст] // Reports on Progress in Physics. – 2005. – № 68, P. .

4.  Liu, G. R. Smoothed particle hydrodynamics: a meshfree particle method [Текст] / G. R. Liu, M. B. Liu // World Scientific, 2003.

5.  Koshizuka, S. A particle method for incompressible viscous flow with fluid fragmentation [Текст] / S. Koshizuka, H. Tamako, Y. Oka // Computational Fluid Dynamics. – 1995. – P. 29-46.

6.  Koshizuka, S. Numerical analysis of breaking waves using the moving particle semi-implicit method [Текст] / S. Koshizuka, A. Nobe, Y. Oka // International Journal for Numerical Methods in Fluids. – 1998. – P. 751-769.

7.  Liu, W. K. Multiresolution Reproducing Kernel Particle Method for Computational Fluid Dynamics [Текст] / W. K. Liu, S. Jun, D. T. Sihling, Y. Chen, W. Hao // International Journal for Numerical Methods in Fluids. – 1997. – № 24. P.

8.  Dilts, G. A. Moving-least-squares-particle hydrodynamics i: consistency and stability [Текст] // International Journal for Numerical Methods in Engineering. – 1999. – 44(8). – P. 1115 – 1155.

9.  Cummins, S. J. An SPH projection method [Текст] / S. J. Cummins, M. Rudman // Journal of computational physics. – 1999. – № 000. – P. 584-607.

10.  Shao, S. Incompressible SPH method for simulating Newtonian and non-Newtonian flow with a free surface [Текст] / S. Shao, E. Y.M. Lo // Advances in Water Resources. – 2003. – № 26. – P. 787-800.

11.  Brookshaw, L. A Method of Calculating Radiative Heat Diffusion in Particle Simulation [Текст] // Proc. ASA. – 1985. – 6 (2). – P. 207-210.

12.  Коннор Дж., Метод конечных элементов в механике жидкости. / Пер. с англ. – Л.: Судостроение, 1979. – 264 с.

13.  Dilts, G. A. Moving-least-squares-particle hydrodynamics ii: conservation and boundaries [Текст] // International Journal for Numerical Methods in Engineering. – 2000. – 48(10). – P. 1503 – 1524.

14.  Афанасьев, моделирование течений жидкости со свободными границами методами SPH и MPS. [Текст] / , , // Вычислительные технологии. – 2006. – Т. 11. – Спец. Выпуск. – С. 26-44.