Устойчивость одного конвективного течения во вращающейся трубе Текст научной статьи по специальности «Физика»

УДК 532.5

Устойчивость одного конвективного течения во вращающейся трубе

А.В. Проскурин1,2, А.М. Сагалаков2

1 Алтайский государственный технический университет им. И.И. Ползунова (Барнаул, Россия)

2 Алтайский государственный университет (Барнаул, Россия)

The Convective Flow Stability In Rotating Pipe

A.V. Proskurin1,2, A.M. Sagalakov2

1 Polzunov Altai State Technical University (Barnaul, Russia)

2 Altai State University (Barnaul, Russia)

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

Ключевые слова: конвекция, гидродинамическая устойчивость, метод дифференциальной прогонки.

БСТ 10.14258/izvasu(2014)1.2-38

Effective solution of hydrodynamic stability problems, as well as solution of the eigenvalue problem for non-selfadjoint operators in general cases can be obtained only numerically. However, in this case standard numerical methods are usually not applicable: the eigenfunctions of linear stability problem have poor properties. Thus, the differential sweep method is applied for solving this eigenvalue problem. This step reduces the problem to a sequence of Cauchy problems for nonlinear ordinary differential equations system which can be easily integrated numerically. At the same time, an iterative process for eigenvalues is carried out. It is found out that increasing of Prandtl number leads to stabilization of the flow with respect to small perturbations.

Key words: convection, flow stability, differential sweep method.

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

Рассмотрим течение в круглой вращающейся трубе. Ось г цилиндрической системы координат совпадает с осью трубы, г — радиальная координата, ^ — азимутальная. Система уравнений движения вязкой теплопроводящей жидкости в приближении Обербека-Буссинеска имеет вид

8V ~dt

дТ ~dt

+ (VV)V = -PrVp + PrAV - PrRa ■ rTer

+ (V V)T = AT, divV = 0,

где Ра = ^'х'1То — число Ре лея; Рг = ^ — число Прандтля; Д — радиус трубы; ш — угловая ско-

(1) (2) (3)

рость вращения трубы; йТ0 — разность температур между стенкой и осью трубы; V — кинематический коэффициент вязкости; х — коэффициент температуропроводности; в — коэффициент объемного расширения; ег — единичный радиальный вектор.

Решение (1)—(3) может быть представлено в виде

V = и + V, (4)

Т = Т0 + Т,

Р = Ро + Р,

(5)

(6)

96

На

и

гп _ _

0 ~ 9-128 др Да 2

(1 - г2)3,

(7)

Далее мы рассмотрим устойчивость этого течения.

Подставляя (4)-(6) в (1)-(3) и пренебрегая членами второго порядка малости, получим линейную систему

алгебраических преобразований получим систему уравнений

где и, Т0, Р0 — стационарные значения скорости, температуры, давления соответственно; V, Т, р — их малые возмущения.

В работе [1] было найдено решение задачи (1)-(3), которое в стационарном случае выражается при помощи элементарных функций

— iaCvr + 'аиvr = —Рг ■ — (тоГ)' ( т2 + 1

— Рг

+ а Vr —

2'т

г2 о<f

— РгНа ■ гТ, (12)

— ¡аСу^ + ¡аИу^ = —Рг • —</+

+ Рг

т +1 2\ 2'т

+ о. ] Юц, + —ррУг

— ''гаСгох + га.игих + и= —Рг ■ 'ац+ \r-vZ)' /т2

+ Рг

+ 1 2 + а

(13)

(14)

— 'аСТ + 'аЛТ + =

{гГУ ^ + (15)

(гVr )

гт

= ——г^ - юмх.

(16)

После некоторых упрощений получим

Dvr = —Ргд'—

— Рг 'аго'х + 7

(г^У

— РгНа ■ гТ, (17)

д^и

т

+ (и V)v + =

= —РгУр + PгAv — РгНа ■ гТ, (8)

дТ

+ (и V)T + (vV)Tо = ДТ,

(9)

= о, (10)

решение которой можно представить в виде

{V,р, Т} =

= {vr(г), (г), vz(г), ц(г),1(г)}в1а(г-С1)+1т^,

где vr, vv, vz, д, Т — амплитуды скорости, давления и температуры соответственно; а — осевое волновое число; т = 0,1,2, 3,... — азимутальное волновое число; С = X + 'У — комплексная фазовая скорость, в которой X — собственно фазовая скорость, а аУ - декремент затухания возмущения (У < 0) или инкремент его нарастания (У > 0). После подстановки и простых

Dvv = —^Ргц + Рг ^ ^

Dvz + и V =

= —'аРгц + Рг

DlT + T0vr =

'(г^ У

27 Н--уг

(18)

' (гVz )'

(гТ)

+

+

Т

(гvr)' 'т - =--Ую - ю:их

(19)

(20)

(21)

(И) где В = га(и - С) + Рг{+ а2), В1 = га(11 -

С) + ^т + о?, 7 = —. Введем обозначения

В = гог , К = го^,

Ф

К' (п^)'

И =

И' (гvz)' Ф = — = Б = гТ,

2

г

г

г

г

г

г

г

г

V

z

2

г

г

2

г

г

г

г

z

г

г

г

г

г

г

Система уравнений (17)-(21) примет вид

Р га

д'=--— — г-^Ф-ПаБ, (23)

Ф'

2 + —р^В + гад,

В' = —7К — га.2,

О' =

Е±

г

Б+^В г

Я' = Фг, 2' = Фг, 5' = Ог.

(24)

(25)

(26)

(27)

(28)

Запишем систему (23)-(28) в матричном виде, для этого введем векторы Ш = {В, К, Г} и V = [д, Ф, Ф, О}.

Ш' = М1Ш + M2V, V' = М3Ш + М4У,

(29)

где

М1 =

0 —1 —га 0 0 0 0 0

0 0 0 0 , М2 = 0 г 0 0

0 0 0 0 0 0 г 0

0 0 0 0 0 0 0 г

М3 =

/__Щ-

I гРг

щ

и'

г Рг

Г

0

р

г Рг 0

0

га

Р гРг

Вх

г

—Яа \ 0 0

— гт/

/0 —7 —га 0\

М4 =

7 0 0 га 0 0 0 0 0 0 0

Полагая возмущения скорости и температуры равными нулю на стенке трубы и ограниченными на ее оси, запишем граничные условия в виде

Ш = [В, К, Г} = 0 при г = 0,1. (30)

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

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

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

Введем уравнения для подпространства решений

Ш = Л(г^, (31)

где Л - матрица размера 4 х 4.

Выражение (31) обычно называют схемой прогонки. В процессе вычислений схему прогонки можно изменять, добиваясь наиболее «благоприятного» поведения коэффициентов матрицы Л.

На границе и оси трубы заданы условия Ш = 0, из чего немедленно следуют граничные условия Л = 0 при г = 0,1. Нередко применяют непосредственное интегрирование системы (29) на небольших отрезках от границ с последующим использованием той или иной схемы прогонки.

Численные эксперименты позволили установить, что схема (31), определенная непосредственно граничными условиями, неоптимальна при интегрировании в средней части отрезка [0,1]. Наиболее экономичным с учетом затрат машинного времени и простоты алгоритма оказался следующий вариант прогонки: на небольших расстояниях от границ прогонка велась по схеме (31), определенной граничными условиями, а в средней части использовалась «обращенная» схема прогонки

V = СШ,

(32)

причем в точке стыковки в силу непрерывности собственных функций и их производных С = Л-1. Так как уравнения содержат коэффициенты вида интегрирование велось не точно от оси трубы, а с некоторым малым отступом, эта техника подробно описана в [2]. Системы дифференциальных уравнений для матриц Л и С можно получить, продифференцировав (31), (32) и подставив в них выражения (29). Находим

Л' =М1Л + М2 — Л (М3Л + М4), (33) С' =М4С + М3 — С (М2С + М1). (34)

3

г

1

3

г

0

1

3

г

0

0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1.2

0.01 0.02

0.03 0.04 X

0.05 0.06

0

Рис. 2. Нейтральные кривые а(Яа) при Рг = 0.001, сплошная линия соответствует т = 0, пунктирная - т = 1,

точечная - т = 2

Рис. 3. Критические зависимости Яа*(Рг), сплошная линия соответствует т = 0, пунктирная - т = 1,

точечная - т = 2

Интегрирование уравнений (33), (34) велось численно от границ к некоторой средней точке, в которой, так же как и в точках инверсии схемы прогонки, векторы V и Ш должны быть непрерывны, что позволяет записать систему алгебраических уравнений для величин Ш

(С+ — С-)Ш = 0, (35)

где знаками «+» и « —» обозначены прогоночные коэффициенты С, полученные встречным интегрированием от оси и стенки трубы. Поскольку Ш = 0, можно записать дисперсионное соотношение

det (С+ (С) — С- (С)) = 0. (36)

Численные эксперименты показали, что при такой схеме организации прогонки дисперсионное соотношение (36) обладает «хорошими» свойствами. Для его решения использовались итерацио-ные методы.

После определения собственных значений легко найти Ш и V на всем отрезке интегрирования, что составляет содержание обратной прогонки.

На рисунке 1 представлены собственные значения при Рг = 0.001, Яа = 10, 1д а = 0.5.

При У < —0.2 собственные числа имеют одну фазовую скорость и поэтому выстраиваются в вертикальную линию. При У > —0.2 собственные значения расположены в некотором диапазоне по X, крайнее левое соответствует неустойчивому возмущению. На рисунке 2 изображены нейтральные кривые а(Яа) при Рг = 0.001 , сплошная линия соответствует т = 0, пунктирная - т = 1, точечная - т = 2. Верхние ветви кривых практически совпадают, но критические значения достигаются при разных значениях а, с увеличением т также изменяется область генерации неустойчивости за счет положения нижней ветви нейтральных зависимостей. На рисунке 3 приведены критические зависимости Яа*(Рг), сплошная линия соответствует т = 0, пунктирная - т = 1 , точечная - т = 2. С увеличением числа Рг увеличивается соответствующее критическое число Релея ЯаФ, зависимость близка к линейной. Возмущения с т = 1 наиболее опасны. Таким образом, увеличение числа Прандтля приводит к стабилизации течения по отношению к малым возмущениям.

Авторы благодарят В.Н. Штерна за постановку задачи.

Библиографический список

1. Бирих Р.В., Пухначев В.В. Осевое конвективное течение во вращающейся трубе с продольным градиентом температуры // Доклады Академии наук. - 2011. - Т. 436, №3.

2. Гольдштик М.А., Штерн В.Н. Гидродинамическая устойчивость и турбулентность. — Новосибирск, 1977.

3. Проскурин А.В., Сагалаков А.М. Устойчивость течения Пуазейля при наличии продольного магнитного поля // Прикладная механика и техническая физика. — 2008. — Т. 49, №3.

4. Проскурин А.В., Сагалаков А.М. Новая ветвь неустойчивости магнитогидродинамическо-го течения Пуазейля в продольном магнитном поле // Письма в ЖТФ. — 2008. — Т. 34, вып. 5.