CC BY
18
Вычислительные технологии
Том 18, № 6, 2013
О возникновении движения в конечном цилиндре*
Е. П. Магденко
Институт вычислительного моделирования СО РАН, Красноярск, Россия e-mail: magdenko_evgeniy@icm.krasn.ru
Рассмотрена задача о возникновении конвекции в цилиндрическом контейнере. Для её решения применён метод разделения переменных. В результате получено однородное дифференциальное уравнение шестого порядка с постоянными коэффициентами со сложными граничными условиями. Для случая монотонного возмущения получено аналитическое выражение для критических чисел Марангони. Рассмотрен случай, когда система находится в состоянии невесомости, для которого доказано, что с увеличением радиуса цилиндра критические числа Марангони стремятся к известным числам Марангони для бесконечного слоя.
Ключевые слова: критические числа Марангони, метод разделения переменных, свободная граница, конвекция.
1. Постановка задачи
Пусть цилиндрический контейнер заполнен покоящейся жидкостью с верхней свободной деформируемой границей, на которой задан теплообмен с окружающей средой.
Обозначим через П = (0, а) х (0, 2п) х (-к, 0) область, занимаемую жидкостью (рис. 1).
Рис. 1. Схема области конвекции
* Работа выполнена при поддержке Интеграционного проекта СО РАН № 38 и проекта РФФИ № 1101-00283.
Равновесное состояние системы в рамках модели Обербека — Буссинеска описывается уравнениями [1]
u = 0, (1)
Pz = pogße, (2)
0 = Az + B, (3)
где u — вектор скорости; p — отклонение давления от гидростатического; p0 — плотность, ß — коэффициент теплового расширения, 0 — температура жидкости; A, B — температурные коэффициенты, которые находятся из условия на нижнем основании контейнера и условия теплового контакта на свободной границе:
A = Bi (0oi - 002) B = 001 + Bi0o2 (4)
(1 + Bi)h , 1 + Bi , ()
здесь Bi = Yh/k — число Био, y — коэффициент межфазного теплообмена, 0o1 и 0o2 — температура на нижнем и верхнем основании соответственно. Поверхностное натяжение линейно зависит от температуры, а свободная граница является плоской (круг). Давление квадратично зависит от z и имеет вид
p = p0gß (^z2 + Bz^ + c, c = const. 2. Возмущённое решение
При некоторой критической разности температур на основаниях цилиндра 0o1 - 0o2 возникает движение — конвекция. С целью определения этой разности рассматривается линеаризованная на равновесном состоянии задача (1)-(4) о малых возмущениях системы в рамках модели Обербека — Буссинеска, решение которой находится в виде нормальных волн
(U, P, T, R) = (U (r, z), P (r,z) ,T (r, z), N) exp [i (s<p - Ct)], (5)
где U, P, T — возмущения основного решения u, p и 0; R — нормальная составляющая вектора возмущений на свободной границе; N — отклонение амплитуды возмущений свободной границы по нормали; s — азимутальное волновое число; C — комплексный декремент. Тогда для осесимметричного случая (s = 0) монотонных возмущений (C = 0) в безразмерных переменных (в качестве масштаба длины, времени, скорости, давления и температуры выбраны соответственно h, h2/v, v/h, p0v2/h2, Ah) задача (1)-(4) описывается уравнениями (U = (U, 0, W))
Pr = (Urr + 1 Ur + Uzz - rQ , (6)
Pz - GT = ( Wrr + 1 Wr + Wzz ) , (7)
r
Ur + U + Wz = 0, r
W = ^ f Trr + 1 Tr + Tzz Pr V r
где О = двАк4/у2 — число Грасгофа, Рг = у/х — число Прандтля. На свободной границе выполняются следующие условия [1]:
dW
-P + 2
dz (G' + Ga) + We ^Nrr + 1 N^j , (10)
dW dU М/л7
w + nzz = -pr+T-) ■ (11)
Tz + Bi (T + N) = 0, (12)
W = 0, (13)
здесь G' = g@Bh3/v2 — число Грасгофа; Ga = gh3/v3 — число Галилея; We = ah/pov2 — число Вебера; M = &Ah2/povx — число Марангони (заметим, что в силу (4) оно прямо пропорционально искомой разности температур на нижнем и верхнем основаниях цилиндра). Условия на нижнем основании дают
U (r,-1) = W (r,-1) = 0, T (r,-1) = 0. (14)
На боковой поверхности выполняются следующие условия
и(1 ,z) =0, W[ 1 ,z) =0, T[ 1 ,z) =0, \a J \a J J
т.е. жидкость может просачиваться по нормали к стенке, при этом её общий поток через всю боковую поверхность равен нулю.
Задача (6)-(14) допускает разделение переменных:
U =1 R (r) Fz (z), (15)
r
W = -1 Rr (r) F (z) , (16)
r
T = 1 Rr (r) D (z), (17)
r
где
R = Rn (r) = rJ\ (mr) , (18)
здесь m = aSn = hSn/a, Sn, n =1, 2..., — решение уравнения
Jo(S) = 0, (19)
первые корни которого равны [2] Si = 2.4048255577, S2 = 5.5200781103, S3 = 8.6537279129, S4 = 11.7915344391, S5 = 14.9309177086. В общем случае ¿П° ~ nn + 3п/4 при n ^ ж. Из (18), (19) выводим равенство Rnr (a) = 0. Таким образом, условия на боковой поверхности для возмущения температуры и касательной скорости заведомо выполнены. Заметим также, что согласно (15)-(17) величина N пропорциональна r-1Rr(r), т.е.
N = Nor-1Rr (r) = mNoJo (m), No = const.
Подстановка выражений (15)-(17) в уравнения (6)-(9) приводит к обыкновенному дифференциальному уравнению шестого порядка
-1 L3D - m2GD = 0, (20)
Pr
где Ь = — т2. Функция Р(г) вычисляется из уравнения
Р (г) = — ^(21)
В результате функция Д определяется с точностью до шести постоянных, которые находятся из семи граничных условий (N0 входит в число неизвестных постоянных) (10)-(14).
3. Зависимость числа Марангони от геометрии контейнера и физических параметров жидкости
Решение уравнения (20) выглядит следующим образом:
^ H 2m(cosh Aiz — cos A3z cosh A2z) — A3 sin A3z sinh A2z^
8m3
A23
H2 2m(sinh Aiz — cos A3z sinh A2z) — A3 sin A3z cosh A2z
8m3
A23
H3 sin A3z sinh A2z H4 sin A3z cosh A2z
2m
A3
+
2m
A3
+ H5 cosh Aiz + H6 sinh Aiz,
где
A2 = -7- m
A3
Ai = m(1 + b)i/2, 1 — 2 + ((1 — b)2 + b)i/2
i/2
m
b
1 — 2 +((1 — b)2 + b)i/2
i/2
b
PrG
m4
а Hi, i = 1,..., 6, — неизвестные постоянные. Формула для функции F имеет более громоздкий вид и поэтому в статье не приводится. Подставляя найденное решение в условия (10)-(14), получим систему уравнений, которая будет однородной относительно постоянных Hi, i = 1,..., 6. Нетривиальное решение существует тогда и только тогда, когда её определитель равен нулю. Это даёт возможность найти критические числа Марангони путём аналитических вычислений в системе Maple, показывающих, что числа Марангони зависят от геометрии контейнера и физических параметров жидкости. Данная зависимость, а также асимптотические формулы имеют слишком громоздкий вид и поэтому также не приводятся.
При получении асимптотических оценок были использованы следующие выражения: — при m ^ 0 (физически это означает, что рассматривается плоский слой, т. е. a ^ то)
Ai ~ Pri/6G1/6m1/3 + 1
m
5/3
2 Pri/6Gi/6'
+
f
+
4
Л2 ~ -1 Pr1/6G1/6m1/3 - 1 m65/31/6, 2 2 4 Pr1/6G1/6
А3 ~ — ^Pr1/6G1/6m1/3 + — m5/3 А3 2 G m + 4 Pr1/6G1/6 '
2 3 2
^C ^C ^C
sin x ~ x, cos x~ 1--, sinh x~x +--, cosh x~ 1+--;
2' 6 ' 2 '
— при m ^ œ (в данном случае либо n ^ œ, либо радиус цилиндра слишком мал). Здесь sinh x ~ ex/2, cosh x ~ ex/2.
Теперь рассмотрим конкретную заполняющую сосуд жидкость — трансформаторное масло. Физические параметры таковы: р0 = 0.86 • 103 кг/м3, v = 18.49 • 10-6 м2/с, х = 1.21 • 10-5 м2/с, k = 0.63519 • 10-4 кг • м/(с3 • К), в = 0.7 • 10-3 К-1, œ = 0.0022 Н/(м • К), а = 3.81 • 10-2 Н/м.
На рис. 2 приведены графики зависимости числа Марангони от m = 5nh/a при 1 < m < 10, We = 104, Bi = 2. Физически это означает следующее. Например, в точке (1.92674905565, 3.203963483 • 105) график функции достигает своего минимума, и если принять n = 1, т. е. = 2.4048255577, то исходя из формул для m и M получится, что при таком отношении высоты слоя жидкости к радиусу цилиндра, а именно, h/a = 0.8012011722, критическая разность температур во1 — ©о2 = 27.23325641 К. При n = 5 (S5 = 14.9309177086) данная разность температур будет достигнута уже при h/a = 0.1290442485.
Рис. 2. Графики зависимости числа Марангони от m при We = 104 и Bi = 2
Случай невесомости. Предположим, что д = 0. Тогда решение уравнения (20) будет следующим:
в = (В^'2 + (Н— 4Юг +Н) С°8Ь
+ (+ 2т (Н — ш)г +Яб)(22)
где Н.1, г = 1, 2,..., 6, — некоторые неизвестные постоянные. Теперь функция Р из (21) примет вид
F = -1
Pr
—2z + H3 ) cosh mz + ( —+ H4 ) sinhmz 2m J \2m
(23)
Используя условия на свободной границе, с учётом (22), (23) получим:
— динамическое условие 3т2 Р — Р^ = т4ШеЖ0 —
Н4 + тШеРгЖо = 0; (24)
2 2 М
— условие касательных напряжений т2Р + Р^ = т2—- (N0 + Д) —
Рг
Н + 2т2 Н3 — т2 МН5 — т2МЖ0 = 0; (25)
— условие теплового контакта + Б1 (Д + N0) = 0 —
1 ГН2 + -1- Н4 + Б1Н + тН6 + Б1Ж0 = 0; (26)
8m3 2m — кинематическое условие F = 0 —
H3 = 0. (27)
Граничные условия на нижнем основании цилиндра для скоростей и температуры
F (-1) = Fz (-1) = 0, D (-1) = 0 примут вид HH
— sinh m--cosh m + H3 cosh m - H4 sinh m = 0, (28)
2m 2m
1 / sinh m\ 1 / cosh m\ - cosh m +--Hi +— sinh m +--H 2-
2 у m J 2 \ m J
-mH3 sin m + mHi4 cos m = 0, (29)
1 / , sinh m \ 1 ( 1 . п \ H3 . .
cosh m--H i + -—- — cosh m - sinh m H2 + -— sinh m-
8m2 \ m J 8m2 \ m у 2m
H
--cosh m + H5 cosh m - H6 sinh m = 0. (30)
2m
Система, полученная из условий (24)-(30), будет являться алгебраической относительно постоянных. Нетривиальное решение системы уравнений существует тогда и только тогда, когда её определитель, составленный из коэффициентов при неизвестных, равен нулю. Это даёт возможность найти критические числа Марангони. Аналитические вычисления в системе Maple показывают, что
8m (Bi sinh m + m cosh m) (m — sinh m cosh m) mm — Q ZT1 • (31)
m3 cosh m — sinh m — 8m3 cosh m (PrWe)
Замечание. Если a ^ ж, n ^ ж таким образом, что m — h5n/a ^ m0 — const, то выражение (31) в точности совпадает с числом Марангони для бесконечного слоя [3].
Для (31) найдены следующие асимптотические формулы:
— при m0 ^ ж (sinhm ~ em0/2, coshm ~ em0/2)
M - 8m0(m0 + Bi); (32)
— при m0 ^ 0 (sinhm ~ m0 + m0/6, coshm ~ 1 + m0/2)
2
M --PrWe (Bi + 1) m2; (33)
3
— в случае недеформируемой свободной поверхности
M - 64(Bi2+1). (34)
m0
Далее, как и для g — 0, рассмотрим конкретную жидкость, заполняющую цилиндрический контейнер, — трансформаторное масло, физические параметры которого приведены выше. Графики зависимости критического числа Марангони от числа m0 для этого случая представлены на рис. 3. Анализируя полученные графики и используя
Рис. 3. Графики зависимости критического числа Марангони от m0: 1 -We = 104, Bi = 0, 2 - We = 104, Bi = 2, 3 - We = 107, Bi = 2, 4 - We = ж, Bi = 2
формулы для числа Марангони при т0 для случая п = 1, получим следующую зависимость критической разности температур от чисел Вебера и Био:
Кривая 1 Кривая 2 Кривая 3 Кривая 4
Точка
минимума (1.99,79.534) (2.4,150.614) (2.4,150.6844121) (2.4,150.6844822)
К/а 0.8275 0.998 0.998 0.998
©01 — ©02 0.00 07 К 0.001317235403 К 0.00131785121 К 0.001317851824 К
Таким образом, полученные данные позволяют сделать вывод: зная заранее геометрию контейнера и физические параметры находящейся в нём жидкости, можно определить критическую разность температур, при которой в цилиндрическом контейнере возникнет конвекция.
Автор выражает благодарность д-ру физ.-мат. наук, профессору В.К. Андрееву за постановку задачи и ценные советы при проведении настоящего исследования.
Список литературы
[1] Андреев В.К., Захватаев В.Е., Рявицкий Е.А. Термокапиллярная неустойчивость. Новосибирск: Наука, 2000. 31 с.
[2] Аврамовица М., Стиган И. Справочник по специальным функциям с формулами, графиками и математическими таблицами. М.: Наука, 1979.
[3] Рявицкий Е.А. Колебательная термокапиллярная неустойчивость равновесия плоского слоя в присутствии поверхностно-активного вещества // Изв. РАН. МЖГ. 1993. № 1. С. 6-10.
Поступила в 'редакцию 14 августа 2013 г., с доработки — 23 сентября 2013 г.
