CC BY
17УДК 531.746.22:519.87
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ В КАНАЛЕ С ПОДВИЖНЫМИ ГРАНИЦАМИ
И. В. Бычин 2, Т. В. Гавриленко 1i2, В. А. Галкин ^ А. В. Гореликов ^ А. В. Ряховский 2
1 Федеральный научный центр Научно-исследовательский институт системных исследований
Российской академии наук, 2 Сургутский государственный университет, igor-bychin@yandex.ru, taras.gavrilenko@gmail.com, val-gal@yandex.ru, gorelikov_a@list.ru, echo47@rambler.ru
В статье рассматривается математическое моделирование течений вязкой несжимаемой жидкости в каналах с периодически изменяющимися во времени поперечными размерами: построение математической модели; алгоритм численного решения; проведение вычислительных экспериментов по исследованию зависимости характеристик течения от режима колебаний границ канала.
Ключевые слова: численное моделирование, течение в канале.
NUMERICAL SIMULATION OF VISCOUS INCOMPRESSIBLE FLUID FLOW IN CHANNEL WITH MOVING BOUNDARIES
I. V. Bychin 2, T. V. Gavrilenko 1 2, V. A. Galkin 1 2, A. V. Gorelikov 1 2, A. V. Ryakhovsky 2
1 System Research Institute, Russian Academy of Sciences, 2 Surgut State University,
igor-bychin@yandex.ru, taras.gavrilenko@gmail.com, val-gal@yandex.ru, gorelikov_a@list.ru, echo47@rambler.ru
The article considers a mathematical modeling of viscous incompressible fluid flows in channels with transverse dimensions alternating with time. The modeling consists of: construction of a mathematical model; numerical algorithm; carrying out computational experiments to study the dependence of flow characteristics on the vibration mode of edges of a channel.
Keywords: numerical simulation, channel flow.
Введение. Построение и обоснование математического описания динамики сложных систем, к которым относятся многофазные смеси в слабо сжимаемых жидкостях, представляют интерес с точки зрения моделирования технологических процессов. При этом одним из важных вопросов является управление разделением фракций за счет изменения параметров входящего и выходящего потока и деформацией области течения. В практически важных задачах, как правило, присутствует ряд усложняющих факторов: периодическое изменение внешнего силового поля, химические реакции, наличие поверхности раздела и т. д. В связи с этим представляется целесообразным исследование различных режимов течения жидкости в условиях сложной, изменяющейся во времени геометрии в упрощенной постановке для однофазных жидкостей.
Постановка задачи и математическая модель. Рассматривается ламинарное течение изотермической вязкой несжимаемой жидкости в прямоугольном канале D = {(.X, y) : X е (0, L), y е (0,0.25!)} . Верхняя и нижняя стенки канала на интервале x е (0.25L, 0.5L) синхронно совершают периодическое движение навстречу друг другу в направлении оси Oy с частотой щ, тем самым изменяя поперечные размеры канала на данном участке. На левой границе (x = 0) задается скорость входящего потока жидкости, которая также периодически изменяется с частотой щ.
Обезразмеривание производилось по следующим формулам:
т-2 2
* * v * L * pv v г
r = Lr, u = -u, t = — t, p = p, а>01 = — /од.
L V L L
г\ О. * *
Звездочкой отмечены размерные величины: r - радиус-вектор точки; t - время; u -
* Т
скорость жидкости; p - давление; р - плотность; v - кинематическая вязкость; L - размер канала вдоль оси х. В безразмерных переменных система уравнений гидродинамики вязкой несжимаемой изотермической жидкости имеет вид:
— + (u-V)u = — Vp + Au , (1) 8t
div u = 0. (2)
На верхней и нижней границах канала заданы условия прилипания. Скорости движения нижней (uby ) и верхней (ы'у ) границ канала на участке x е (0.25,0.5) определяются по формулам:
uby =sgn(sin 2^ и), uy = — sgn(sin 2^ /xt). (3)
о о
Из формулы (3) следует, что максимальное безразмерное расстояние, на которое вставки перемещаются внутрь канала, не зависит от частоты fl и составляет 1/16, т. е. 1/4 от высоты канала.
На входной границе задается следующее периодическое по времени распределение скорости:
UxL=o =—12Rey(4y — l) 3 + /ot)), (4)
где Re = — - число Рейнольдса; u - средняя по толщине жидкости скорость; /0 - частота v
изменения скорости входного потока жидкости. Коэффициент 12 в (4) получается из формулы среднего значения скорости в задаче о течении Пуазейля между двумя плоскостями [1].
Предполагается, что скорость выходящего потока (при х = 1) направлена вдоль оси Ох и рассчитывается таким образом, чтобы выполнялось уравнение неразрывности, т. е. чтобы расход на входе был равен расходу на выходе.
В начальный момент времени скорость жидкости во всем канале направлена вдоль оси Ох и определяется по формуле:
ux (x, y)\ t=0 =—12Re y (4 y — l) , (x, y) е D. (5)
Численные методы. Дискретные аналоги дифференциальных уравнений составлены при помощи метода контрольного объема [2] в прямоугольных декартовых координатах. Поперечное движение вставок в канале моделируется с помощью блокировки контрольных объемов [2], при этом используется аналог VOF (volume of fluid) метода [3]. Дискретные аналоги уравнений (1-2) получены с использованием неявной схемы и схемы со степенным законом для аппроксимации конвективных и диффузионных потоков на гранях контрольных объемов [2]. При получении дискретного аналога уравнения Навье - Стокса использовались разнесенные расчетные сетки, т. е. компоненты скорости рассчитывались на гранях основных контрольных объемов. Для решения уравнений гидродинамики используется алгоритм PISO (Pressure-Implicit with Splitting of Operators) [4-5]. Алгоритм PISO выбран как наиболее эффективный для решения нестационарных задач гидродинамики вязкой жидкости в рамках метода контрольного объема. Получающиеся в результате дискретизации СЛАУ решаются параллельной модификацией метода переменных направлений.
Рис. 1. Иллюстрация к алгоритму расчета скорости на выходе канала
Сеточное значение скорости на выходной границе и°ш рассчитывается по формуле:
и°хШ (У) = и? (]) + л, У = 2, М 2,
где и^2 - значение скорости на левой границе контрольного объема, прилегающего к выходной границе (рис. 1); L2, M2 - количество контрольных объемов вдоль осей Ox и Oy соответственно; и1^ - скорость жидкости на левой границе области D; H - размер области D вдоль оси Оу; Ауу - размер вдоль оси Оу у'-го контрольного объема. Константа Я
определяется из условия равенства объемного расхода жидкости на входной и выходной границах:
М 2 М 2 М 2 М 2
£ иХп(У)Ауу =£ и°°и (у)Ауу =£ ) + Я) Ау у =£ и^Ауу + ЯН,
у =2 у=2 у=2 у=2
M 2
х СО - "х С/))Avj
j=2
H
Вычислительные эксперименты по исследованию зависимости характеристик течения от режима колебаний боковых границ канала проводились на системе кластерного типа АПК-5 СурГУ. Для организации вычислительного процесса на кластере использована программная реализация в рамках разработанного ранее авторами программного комплекса для моделирования процессов динамики вязкой жидкости и магнитной гидродинамики на высокопроизводительных системах. Программный код комплекса написан на языках Фортран и Си. Комплекс спроектирован в рамках многоуровневого подхода, предполагающего сочетание различных технологий параллельного программирования для гибридных и традиционных вычислительных систем. Данный комплекс позволяет на базе технологий MPI [6-7], OpenMP [7-8], OpenCL [9-10] организовать отображение многоуровневого параллелизма вычислительного кода на многоуровневый параллелизм вычислительной системы. Распределение вычислений между узлами кластера выполнено методом геометрической декомпозиции расчетной области с использованием стандарта MPI, по схеме один узел - один процесс. В рамках одного узла, распараллеливание вычислительного процесса на многопоточное исполнение ядрами CPU проведено по технологии OpenMP. Рассматриваемый подход позволяет дополнительно на узле организовать гетерогенные вычисления на ускорителях с использованием стандарта OpenCL.
Таким образом, программный комплекс обеспечивает эффективный счет на всех доступных вычислительных ресурсах кластера. Программный комплекс прошел всестороннее тестирование на следующих тестах: нестационарные трехмерные задачи Коши для уравнений Навье - Стокса и уравнений магнитной гидродинамики с точными решениями [11-12]; сравнение с результатами натурных экспериментов [13]; сравнение с результатами вычислительных экспериментов, представленных в литературе [14]. Результаты тестов на численную сходимость на последовательности вложенных расчетных сеток и с разными шагами по времени демонстрируют второй порядок аппроксимации по координатам и первый порядок по времени [15]. Результаты тестирования позволяют сделать вывод о возможности использования данного программного комплекса для вычислительных экспериментов в рассматриваемой предметной области.
Результаты вычислительных экспериментов. Проведена серия вычислительных экспериментов для следующих значений частоты поперечных колебаний вставок: /1 = 50; 75; 100; 125; 150, при этом частота колебаний скорости входного потока и число Рейнольдса фиксировали^: /0 = 25, Re = 50. В каждом случае расчеты проводились до того
момента времени, на котором течение гарантировано переходило в установившийся режим, т. е. колебания всех величин становились периодическими (рис. 2).
Исследовалось влияние частоты поперечных колебаний вставок / на поля скорости и давления в установившемся режиме течения. На рис. 3 представлена эволюция структуры течения и распределение поперечной составляющей скорости и для частоты / = 75 за один период.
На рис. 4 показано типичное для всех проведенных расчетов распределение давления р и продольной составляющей скорости их при движении вставок навстречу друг другу и при движении вставок в противоположные стороны.
^ /1//0 = 2
Ь) /,//0 = 6
Рис. 2. Колебания продольной составляющей скорости в точках на оси канала: А(0.125, 0.125) (красная линия); В(0.375, 0.125) (зеленая линия); С(0.75, 0.125) (синяя линия)
в установившемся режиме течения.
I
-16 -14 -12 -10 - 8 - S -4 -2 0 2 4 6 10 12 14 16
a) t = 0.02T
■16-14 -12-10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16
b) t = 0.26T
Ш I I 1 ■
-25 -12.5 0 12 5 25
c) t = 0.5 T
-25 -12.5 0 12 5 25
d) X = 0.527;
0,25 0.5 0.75 1
I ■
-14 -7 0 7 14
e) X = 0.767
■9 -4.5 0 4.5 9
0 X = 0.987
Рис. 3. Эволюция структуры течения за период движения поперечных вставок 7 = 1 / /, т. е. показаны линии тока нестационарного поля скоростей в фиксированный момент времени цветом показано распределение поперечной составляющей скорости и (/ = 75)
а) г — 0.261;
в .я
Р: -8000 4000 0 4000 8000
Ь) г — 0.761
Рис. 4. Распределение давленияр (слева) и продольной составляющей скорости их (справа): а) при движении вставок навстречу друг другу г — 0.261 Ь) при движении вставок в противоположные стороны г — 0.761 (— 75)
Рис. 5. Зависимость приведенных амплитуд колебаний продольной составляющей скорости и давления (6) от отношения частоты колебаний границ к частоте пульсаций входящего потока — / /0 в точках А(0.125, 0.125),В(0.375, 0.125), С(0.75, 0.125) в установившемся режиме течения
На рис. 5 представлены зависимости приведенных амплитуд колебаний продольной составляющей скорости Ар и давления Ар от отношения частоты колебаний границ к
частоте пульсаций входящего потока (f — f / f) в точках A(0.125, 0.125), B(0.375, 0.125), C(0.75, 0.125) в установившемся режиме течения. Величины Д~ и Др определялись по формулам:
(max и - min и )| (max р - min р)|
Ди--f, Др--f, (6)
(max их - min их ^ f -0 (max р - min р)\f -0
где max мх , min мх , max р , min р - максимальное и минимальное значения продольной составляющей скорости и давления за период колебаний в точках (А, В, С); амплитуды мх и p нормируются на соответствующие амплитуды для течения в канале с неподвижными боковыми стенками f — 0 .
Заключение. Построена математическая модель течения вязкой несжимаемой жидкости в канале с периодически изменяющимися поперечными размерами. На основе метода контрольного объема разработан алгоритм численного решения задач течения в каналах с подвижными вставками. В рамках ранее созданного авторами программного комплекса для моделирования процессов динамики вязкой жидкости на высокопроизводительных вычислительных системах разработано программное обеспечение для решения рассматриваемого класса задач. Проведены вычислительные эксперименты по исследованию зависимости характеристик течения от режима колебаний боковых границ канала.
Работа выполнена при поддержке РФФИ, гранты: 16-29-15105 офи_м, 18-07-00504 А.
Литература
1. Ландау Л. Д., Лифшиц Е. М. Теоретическая физика : учеб. пособие ; в 10 т. Т. VI. Гидродинамика. 5-е изд., стереот. М. : Физматлит, 2001. 736 с.
2. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М. : Энергоатомиздат, 1984. 152 с.
3. Hirt C. W., Nicholls B. D. Volume of Fluid (VOF) method for dynamical free boundaries // Comput Phys. 1981. № 39. P. 201-225.
4. Issa R. I. Solution on the implicitly discretised fluid flow equations by operator-splitting // Journal of Computational Physics. 1985. № 61. P. 40-65.
5. Issa R. I., Gosman A. D., Watkins A. P. The Computation of compressible and incompressible recirculating flows by a non-iterative implicit scheme // Journal of Computational Physics. 1986. № 62. P. 66-82.
6. Gropp W. Using MPI: portable parallel programming with the message-passing interface. Cambridge : MIT Press, 2014. 336 p.
7. Антонов A. С. Технологии параллельного программирования MPI и OpenMP. М. : Изд-во Моск. ун-та, 2012. 344 с.
8. Chapman B. Using OpenMP: Portable Shared Memory Parallel Programming. Cambridge : MIT Press, 2008. 378 p.
9. Scarpino M. OpenCL in Action: How to Accelerate Graphics and Computations. Cambridge : Manning Publications, 2011. 456 p.
10. Kaeli D. R. Heterogeneous Computing with OpenCL 2.0. Cambridge : Morgan Kaufmann, 2015. 330 p.
11. Бетелин В. Б., Галкин В. А., Гореликов А. В. Алгоритм типа предиктор-корректор для численного решения уравнения индукции в задачах магнитной гидродинамики вязкой несжимаемой жидкости // ДАН. 2015. Т. 464. № 5. С. 525-528.
12. Галкин В. А. и др. Моделирование и управление разделением фаз в слабо сжимаемых вязких теплопроводящих жидкостях типа нефти в случае газообразных и твёрдых включений // Вестн. кибернетики. 2015. № 3. С. 27-45.
13. Ряховский А. В. Математическое моделирование естественной конвекции во вращающихся сферических слоях : дис. ... канд. физ.-мат. наук. Сургут. 2013. 143 с.
14. Бычин И. В., Гореликов А. В., Ряховский А. В. Исследование установившихся режимов естественной конвекции во вращающемся сферическом слое // ВАНТ. 2016. № 1. С. 48-58. Сер. Матем. моделирование физич. процессов.
15. Бычин И. В., Галкин В. А., Гавриленко Т. В., Гореликов А. В., Ряховский А. В. Программный комплекс численного моделирования конвекции в сферических слоях на гибридных вычислительных системах (CPU/GPU) // Матем. моделирование. 2014. Т. 26 (10). C. 95-108.
