CC BY
40
УДК 532.542
Метод исследования устойчивости течений в трубах А.В. Проскурин1,2, А.М. Сагалаков1
1 Алтайский государственный университет (Барнаул, Россия)
2 Алтайский государственный технический университет им. И.И. Ползунова (Барнаул, Россия)
Viscous Fluid Flow in a Coaxial Pipe
A.V. Proskurin1,2, A.M. Sagalakov1
1 Altai State University (Barnaul, Russia)
2 Polzunov Altai State Technical University (Barnaul, Russia)
Предложен метод для решения задач устойчивости течений вязкой несжимаемой жидкости в трубах. Метод основан на использовании функций Рвачева, при помощи которых составляется граничная функция. Ее роль аналогична роли сетки в методе конечных элементов, но она зависит только от формы области и ее не требуется изменять для увеличения точности. С помощью граничной функции и полиномов Чебы-шева составлены структуры, приближенно представляющие решение уравнений и удовлетворяющие граничным условиям. Метод использовался как для расчета ламинарного стационарного течения, которое возникает в трубе под действием постоянного продольного градиента давления, так и для решения линеаризованных уравнений эволюции возмущений. Давление исключалось при помощи уравнения Пуассона. Приведен пример исследования устойчивости течения в трубе прямоугольного сечения. Получены собственные числа, определяющие нарастание или затухание возмущений, и графики скорости. Предложенный алгоритм не зависит от формы сечения трубы и может использоваться для исследования устойчивости течений в произвольных трубах, в том числе с внутренними элементами. Предложенный алгоритм проще, чем метод коллокаций, тау-метод или спектрально-элементный метод, и поэтому должен работать быстрее и может быть быстро перенесен на новые компьютерные архитектуры, такие как CUDA, OpenCL, Xeon Phi.
Ключевые слова: течения в трубах, метод Рвачева, бессеточные методы, устойчивость течений.
DOI 10.14258/izvasu(2016)1-10
In this paper, we propose a method to solve problems of stability of viscous incompressible fluid in pipes. The method is based on the R-functions with which the boundary function is obtained. Boundary function role is similar to the mesh in the finite element method, but it depends on shape of the region only and does not require any changes to increase accuracy. With the help of the boundary function and Chebyshev polynomials we design and elaborate structures for an approximate solution of the equations that satisfy the boundary conditions. This method is used for calculation of laminar steady flow in a pipe due to the constant pressure gradient and for solution of the linearized eigenvalue problem for perturbations. The pressure is handled by the Poisson equation. An example of flow stability in a rectangular duct is investigated. We obtain the eigenvalues and plots of perturbation velocity. The proposed algorithm does not depend on the shape of duct cross-section and can be used to investigate the stability of flows in arbitrary pipes including internal elements. This algorithm is simpler than the collocation method, tau method or spectral element method. Therefore, it should work faster and can be quickly adapted to new computer architectures, such as CUDA, OpenCL, Xeon Phi.
Key words: pipe flow, Rvachev functions, mesh-free
methods, flows stability.
Введение. Проблема устойчивости течений является сложной и интересной задачей гидродинамики. Вопросы устойчивости плоских течений, сводящиеся к уравнениям типа Орра-Зоммерфельда, относительно хорошо исследова-
ны (обзор приведен, например, в работах [1,2]). Намного труднее обстоит дело с изучением устойчивости течений, в которых присутствует только одна однородная координата или таких направлений вообще нет, потому что они приводят к урав-
нениям в частных производных. Стандартные методы вычислительной гидродинамики, такие как, например, метод конечных объемов, непригодны для удовлетворительного нахождения собственных значений при больших числах Рейнольдса, так как требуется высокая точность приближенного представления решения, реализуемая только методами спектрального типа. Современное состояние проблемы изложено в обзоре [3].
В работе предложен спектральный метод на основе функций Рвачева [4-6] для исследования устойчивости течений вязкой несжимаемой жидкости в трубе прямоугольного сечения (см. рис. 1). Обзор применений метода функций Рвачева можно найти в работах [6-8]. Применение И-функций в компьютерной графике и техническом моделировании объектов широко обсуждается в научной литературе [9,10]. В гидродинамике данный метод почти не использовался, один из примеров приведен в исследовании авторов [11].
Ось г прямоугольной системы координат выбрана параллельно оси трубы, граница трубы обозначена П. Рассмотрим систему уравнений Навье-Стокса для несжимаемой жидкости
дУ 1 п
— + (ГУ) V = — Vр + ^АУ от р р
(%иУ = 0,
(1)
где У - вектор скорости, п - коэффициент вязкости, р - плотность жидкости.
У
-и
ЬУ <1
—Ьу
и
Рис. 1. Поперечное сечение трубы
1. О методе И-функций. Метод И-функ-ций позволяет построить аналитическое уравнение границы области ш(х, у)= 0, в которой ищется решение дифференциального уравнения. Здесь ш(х, у) — функция расстояния до границы, которую дальше мы будем называть граничной функцией.
В методе И-функций рассматриваемая область разбивается на несколько подобластей, границы которых описываются элементарными функциями. Из этих областей при помощи логических операций (конъюнкции, дизъюнкции и других) можно составить исходную область. При помощи И-операций, введенных в [5], аналогично можно составить уравнение границы области из уравнений границ подобластей. В данной работе использовалась И-конъюнкция, предложенная в работе [12]
Д<? = ш1А3ш2 = ~Ь+аШ1 +аШ2 - + акш*, (2)
где и шо - граничные функции двух пересекающихся областей, а > 1, 1 < Ь < 2 - }</2. Полагая прямоугольник О пересечением двух полос
Ш Эх —
2 т 2 ~ Ьх
2ЬХ
шву
2Ьу
, при помощи (2)
получим граничную функцию = швх Л в шву, график которой приведен на рисунке 2. Параметры а и Ь определяют радиус скругления углов прямоугольника. Большие значения соответствуют меньшим радиусам скругления. В данной работе мы приняли эти значения равными а = 50, к = 50, а Ь =1 на основе рекомендаций из статьи [12] и собственных оценок.
Рис. 2. График
2. Ламинарное течение в прямоугольной трубе. Найдем стационарное решение задачи (1) для течения в трубе прямоугольного сечения. Будем считать, что градиент давления к вдоль оси г постоянный, а скорость и направлена также вдоль оси г. Тогда система уравнений (1) в стационарном случае в примет вид
\дх2 ду2
и=к- = -1. п
(3)
Так как уравнение (3) линейно, его правую часть без уменьшения общности можно положена равной -1, а параметры течения будут зависеть только от геометрии трубы. Значения скоростей для
к
всех остальных величин - можно получить эле-
п 17
ментарным преобразованием уже вычисленных
скоростей. Далее, при исследовании устойчивости решение этого уравнения будет нормироваться так, что максимальная скорость итах будет
х
равна единице. Граничные условия имеют вид
и = 0 на Г. (4)
3. Постановка задачи устойчивости.
Рассмотрим устойчивость стационарного течения. Для этого решение уравнения (1) представим в виде
V = и + V,
р = Ро + р;
(5)
где и, Ро — не зависящее от времени решение задачи (3)—(4), V, р — малое возмущения скорости и давления. Линеаризованное уравнение Навье-Стокса примет вид
д 1
+ + = + - д., (б)
= 0,
итахХ^'
где Re =
число Рейнольдса. В качестве
характеристических параметров выбраны: полуширина канала d и максимальная скорость лами-
П
нарного течения Umax, v = —. Граничные условия
Р
имеют вид
v = 0 на Г. (7)
Давление можно исключить при помощи уравнения Пуассона (см. [13]). Для этого применим оператор V к обеим частям линеаризованного уравнения Навье-Стокса (первое в (6)). Используя уравнение непрерывности (второе в (6)), после простых преобразований получим
Ap = —2
dvx dU dvy dU dz dx dz dy
(8)
В качестве граничного условия на твердой стенке для давления используем выражение
dp dn
= 0 на Г.
(9)
Представим скорость и давление в виде
= Vx(x, y)e'
^i(z-Ct)
Vy = Vy (x,y)eia(z-Ct),
Vz = Vz (x, y)e p = p(x,y)eia(z-Ct),
ia(z-Ct)
(10)
ух,уу — проекции амплитуд возмущений скорости на соответствующие оси декартовой системы координат, р — возмущение давления, а — продольное волновое число, С = X + гУ — комплексная фазовая скорость, в которой аХ - собственно фазовая скорость, а аУ - декремент затухания возмущения (У < 0) или инкремент его нарастания (У > 0).
Тогда линеаризованное уравнение Навье-Стокса и уравнение Пуассона примут вид
э2 _ :
dx2 By2
. fTT dp 1 ( d2 д2 2
ia(U - C)vz + (vx— + Vy-^W
d_
' dx
д
= — iap +
1 ( d2 d2
Re
dx2 dy2
д2
dx2 dy2
pp — a p = —2ia
d_
' dx
vx— + vy
и (11)
4. Численный метод. Зададим равномерное разбиение прямоугольной расчетной области по горизонтали {хо,х\,... ,хпх} и по вертикали
{Уо, У1,...,Уиу} и обозначим это множество точек как X. Поле скорости ух(х,у) во всех этих
к
точках обозначим как ух. Введем также вектор Ф = {То(х)То (у),То(х)Т1(у),...,ТПх(х)ТПу (у)}, где Т(х), Т(у) — полиномы Чебышева первого рода степени г и ], приведенные к соответствующему интервалу, пх, пу — количество членов разложения по осям х и у. Вектор Ф, вычисленный во всех точках рассматриваемой области, обозначим к
как Ф. Тогда поле скоростей без учета граничных условий можно представить в виде
f
Ф © V x
(12)
где Vx — коэффициенты разложения, а знаком © обозначено скалярное произведение. Если значения скорости в точках X известны, можно вычислить значения Vx, решив линейную алгебраическую систему уравнений вида
AV x = b,
(13)
где
ff
A = (Ф, Ф),Лij ff
b = (Vx, Ф), bi
fL y ff
1 dx 1 dy Фj Фi,
-Lx -Ly
fLx Ly y ff dyVx Фi,
1 dx
-Lx J -Ly
(14)
к к
под знаком интеграла поля юх, Ф^ умножаются поэлементно, буквой Ь с соответствующими индексами обозначены границы области интегрирования. Интегралы вычислялись методом Симпсона на множестве точек X.
Для вычисления ламинарного течения, удовлетворяющего граничным условиям (4), представим скорость и в виде
f U
f
Sv © U
ff
ипФ © и.
V
x
V
y
V
z
x
V
x
Подставим выражение (15) в уравнение (3) и,
/
применяя операцию (, $*), получим систему вида (13), в которой
/
/
А — $ XX + $ уу — ($ хх + $
/ / 6 — {-1, $ .
уу
/
, $*),
(16)
Более подробно вычисление стационарного решения описано в работе [14], а график скорости и(х, у) приведен на рисунке 3.
Далее произведем дискретизацию уравнений (11). Для удовлетворения граничным условиям (7), (9), согласно работе [15], скорость и давление возмущений представим в виде
/
V х
/
— $ *
® V х
^у — ® Vу ,
/
V г
/
— $ *
® V г
(17)
'/ / /
/
р — $рр — Ф + Ф + Ф ® р,
Рис. 4. График первых собственных функций ух и уу, Яе = 100
Подставим (19) в (18). Получим
А — с В£,
(20)
где = + Здесь следует заметить,
что граничная функция является номализо-ваииой, то есть = 1. Чтобы это условие
выполнялось, взяты уравнения полос и , из которых составлен квадрат, удовлетворяющие этому условию. В работах [4,5,12] показано, что И-операция (2) сохраняет нормализованность граничной функции.
где
Рис. 3. График и(х,у)
Подставим (17) в (11). Применяя операцию
/
(, $ *), из первых трех уравнений получим алгебраическую задачу на собственные значения
Ац V х + Р х Р — СВ11V х,
А22 V у + Р у р — СВ 22 V у, (18)
А31V х + А32 V у + А33 V г + Шррр — С В 331) г,
Чтобы решить эту задачу, нужно выразить давление через скорость при помощи уравнения Пуассона (четвертого в (11)). Оно примет вид
Р — (Р хх + Р уу)
(#1V х + #2 V у ) — М1V х + М 2 V у.
А —
А11 + Р х М1 Р у М1
,А31 + шРМ 1
Р х М 2
А22 + Ру М 2 А32 + шРМ 2
0 0
А32 у
€ —
'V х
г> у
.V г,
В—
В11 0 0
0 В22 0
0 0 В33
Данная задача на собственные значения решалась итерационными методами: методом Арнольди, Крылова-Шура при помощи библиотеки ЯЬЕРС [16,17].
Рис. 5. График действительной части одной из собственных функций уг, Яе =100
/
5. Результаты вычислений. Рассмотрим течение в трубе при значениях параметров Ьх = Ьу = 1, то есть в квадратной трубе. Выберем равномерное разбиение по осям х и у по 200 точек. Для вычисления поля скорости стационарного течения возьмем количество членов разложения по 20 для каждой из осей. График вычисленного поля скорости и(х, у) приведен на рисунке 3. Данное поле скорости будет использовано при вычислении собственных чисел задачи устойчивости. Количество точек и количество членов разложения Чебышева выбрано аналогично, вычисления производились при Кв = 100.0 и а = 1.0. Наиболее опасное собственное значение: (0.86857, -0.14154), соответствующие собственные функции приведены на рисунках 4 и 5.
Заключение. В работе рассмотрена задача устойчивости течения вязкой несжимаемой жидкости в трубе прямоугольного сечения. Предложен метод на основе функций Рвачева, который может использоваться для исследования устойчивости течений в каналах произвольной формы и при наличии внутренних элементов.
В целом, данный метод может быть разделен на две части. В первой части происходит вычисление граничной функции ша, роль которой близка
к роли сетки в методе конечных элементов. Граничная функция составляется в начале и зависит только от формы рассматриваемой области, увеличение точности вычислений не требует изменения этой функции.
Во второй части мы производим дискретизацию задачи и решаем ее численно. Критически важная с точки зрения производительности часть предложенного алгоритма по сложности соответствует методу Галеркина с численным интегрированием. Широко используемый тау-метод или метод коллокаций, приводит к более сложной компьютерной реализации, так как необходима специальная часть кода, которая производит удовлетворение граничным условиям. Для областей непрямоугольной формы эти методы сильно усложняются или их применение становится невозможным, тогда как предложенный алгоритм не зависит от формы рассматриваемой области. Более гибкие сеточные методы, такие как спектрально-элементный, намного более громоздки. Таким образом, данный метод является более простым, чем существующие, и поэтому следует ожидать, что он может быть скорее перенесен на новые аппаратные платформы, такие как CUDA, OpenCL, Xeon Phi и при этом будет быстрее работать.
Библиографический список
1. Гольдштик М.А., Штерн В.Н. Гидродинамическая устойчивость и турбулентность. — Новосибирск, 1977.
2. Henningson D.S., Schmid P.J. Stability and transition in shear flows — New-York, 2001.
3. Theofilis V. Advances in global linear instability analysis of non-parallel and three-dimensional flows / / Progress in Aerospace Sciences. — 2003. — №39.
4. Кравченко В.Ф., Рвачев В.Л. Алгебра логики, атомарные функции и вейвлеты в физических приложениях. — М., 2006.
5. Рвачёв В.Л. Теория R-функций и некоторые ее приложения. — Киев, 1982.
6. Shapiro V. Semi-analytic geometry with R-functions // Acta Numerica. — 2007. — V. 16.
7. Tsukanov I., Shapiro V., Zhang S. A Mesh-free Method for Incompressible Fluid Dynamics Problems // Int. J. Numer. Meth. Engng. — 2003. — V. 58.
8. Shapiro V., Tsukanov I. The Architecture of SAGE — A Meshfree System Based on RFM // Engineering with Computers. — 2002. — V. 18., №4.
9. Fougerolle Y., Gribok A., Foufou S., Tru-chetet F. and others Boolean operations with implicit and parametric representation of primitives using R-functions / / Visualization and Computer Graphics, IEEE Transactions. — 2005. — V. 11., №5.
10. Freytag M., Shapiro V., Tsukanov I. Field modeling with sampled distances // Computer-Aided Design. — 2006. — V. 38., №2.
11. Proskurin A., Sagalakov A. The numerical investigation of the stability of the localized perturbation in Poiseuille flow // Computational technologies. — 2013. — V. 18., №3.
12. Слесаренко А.П. S-функции в обратных задачах аналитической геометрии и моделировании тепловых процессов // ВосточноЕвропейский журнал передовых технологий. — 2012. — Т. 1, №4(55).
13. Fletcher C. Computational techniques for fluid dymamics 2. — Springer, 1991.
14. Проскурин А.В., Сагалаков А.М. Математическое моделирование течений в трубах с помощью метода функций Рвачева //XI Всероссийский съезд по фундаментальным проблемам теоретической и прикладной механики : сборник
трудов (Казань, 20-24 августа 2015 г.). — Казань, 2015.
15. Proskurin A., Sagalakov A. A method for modelling MHD flows in pipes // Book of abstracts of Russian conference on Magnetohydrodynamics. June 22-25, 2015, Perm, Russia. — Perm, 2015.
16. Saad J. Numerical methods for large eigenvalue problems. — Manchester, 1992.
17. Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems // ACM Trans. Math. Software — 2005. — V. 31., №2.
