МЕХАНИКА
УДК 532.5.013.4.536.25+519.6
Д.А. Фоминский, А.Н. Шарифулин
ЧИСЛЕННОЕ ОПРЕДЕЛЕНИЕ ГРАНИЦ СУЩЕСТВОВАНИЯ АНОМАЛЬНОГО КОНВЕКТИВНОГО ТЕЧЕНИЯ В НАКЛОНЯЕМОМ ЦИЛИНДРЕ
D.A. Fominskiy, A.N. Sharifulin
Perm National Research Polytechnic University, 29 Komsomolsky Av., Perm, 614990, Russia.
NUMERiCAL DETERMiNATiON OF THE BORDERS FOR EXiSTENCE OF ANOMALOUS CONVECTiVE FLOW iN A CYLiNDER TiLTED
Статья посвящена изучению бифуркаций стационарных режимов конвекции в замкнутой наклоняемой цилиндрической полости при гармоническом распределении температуры на стенках. Методом сеток получены бифуркационные диаграммы и кривые для трех значений числа Прандтля, соответствующих воздуху, воде и трансформаторному маслу.
КОНВЕКЦИЯ. БИФУРКАЦИЯ. АНОМАЛЬНОЕ ТЕЧЕНИЕ. ЦИЛИНДРИЧЕСКАЯ ПОЛОСТЬ.
The article is devoted to study of bifurcations of stationary modes of convection in a closed, tilted cylindrical cavity at harmonic temperature distribution on boundary. Using finite difference method the bifurcation diagrams and curves for three values of Prandtl number corresponded to air, water and transformer oil were obtained.
CONVECTION. BIFURCATION. ANOMALOUS FLOW. CYLINDRICAL CAVITY.
Интерес к бифуркациям стационарных режимов конвекции в замкнутой наклоняемой полости связан с надеждой, что через понимание закономерностей их бифуркаций в лабораторных моделях будет достигнуто понимание механизмов зарождения сложных конвективных явлений в атмосфере и мантии Земли. Наклон полости используется для лабораторного моделирования, например внезапного зарождения тропических циклонов [1] и парадоксального переноса загрязнений воздуха из расположенного на уровне моря города в горные районы с дорогой недвижимостью [2].
Стационарные надкритические режимы конвекции в замкнутой полости часто имеют форму валов, т. е. вихрей с горизонтальной осью. В них частицы жидкости движутся вдоль плоских линий тока вокруг точки пересечения этой плоскости с осью вихря. Поэтому такие режимы можно назвать квазидвумерными. Это обстоятельство позволяет авторам работы надеяться, что исследование плоских двумерных течений, т. е. бесконечно вытянутых горизонтальных вихрей в абстрактном бесконечном цилиндре при гармоническом распределении температуры на стенке, поможет в понимании
закономерностей бифуркаций стационарных режимов конвекции в полостях лабораторных экспериментов.
Поясним, что наклон при малых докри-тических значениях числа Рэлея приводит к формированию горизонтального вихря, причем с таким направлением циркуляции, что воздух около более нагретой вертикальной стенки полости движется вверх. Это вихрь с нормальным вращением. Однако при значениях числа Рэлея, превышающих критическое значение, наряду с таким нормальным вихрем возможен и вихрь с обратным направлением циркуляции, т. е. аномальный. В таком вихре жидкость вдоль нагретой стенки движется вниз. Существование аномального горизонтального вихря в цилиндрической полости при малых отклонениях направления подогрева от вертикального было обнаружено в работе [3]. Произвольные отклонения были рассмотрены в статье [4] для значения числа Прандтля Рг = 1.
Исследованию конвекции в горизонтальном цилиндре круглого сечения при гармоническом распределении температуры на стенке посвящена работа [5]. Такая конвекция представляет интерес в связи с возможностью исключить влияние геометрии полости и проследить за структурой конвективного движения при изменении ориентации нагрева по отношению к направлению силы тяжести.
Цель настоящей работы - определить теоретическую зависимость от числа Рэ-лея критического угла наклона, до достижения которого возможно существование аномального вихря, для значений числа Прандтля Рг = 7 и 14.
В связи с целью исследования было использовано численное решение полных уравнений тепловой конвекции в приближении Буссинеска, для указанных значений числа Прандтля.
Постановка задачи
Рассмотрим полость, имеющую форму бесконечного кругового цилиндра (рис. 1). Полость окружена массивом, температуропроводность которого много больше температуропроводности х для жидкости. Коэф-
фициент линейного расширения жидкости Р и ее кинематическая вязкость V постоянны. Здесь введены обозначения: у — единичный вектор, связанный с постоянным градиентом температуры V Т0, который устанавливается в полости (в отсутствие движения жидкости) соотношением
От равновесного направления подогрева у отсчитываются все угловые координаты, п — единичный вектор, указывающий направление вверх и связанный с ускорением свободного падения соотношением
g = - ^п.
Вектор g составляет угол а с равновесным направлением подогрева у.
Введены декартова (х, у, £) и цилиндрическая (г, 9, £) системы координат, связанные между собой соотношениями: х = гео89, у = гэт9, £ =
Введем масштабные множители: 0 — характерная разность температур; Я — радиус цилиндра; V — коэффициент кинематической вязкости.
Будем искать плоские решения задачи. В этом случае векторные поля завихренности и функции тока будут иметь отличными от нуля только компоненты:
Ф = (0,0,ср),ч/ = (0,0,¥). (1)
Рис. 1. Геометрия задачи о свободной тепловой конвекции в цилиндрической полости (см. обозначения в тексте)
Тогда уравнения тепловой конвекции в безразмерной форме запишутся в виде [6]:
дф
— - г(Лг (у Хф) = = -гсЛДУ х Ф) + Ог(УГ х п)г;
го1г(Ух\|/) -ф = 0;
— + а1у(уг) = — дт. дt Рг
(2)
(3)
(4)
Безразмерные критерии подобия — число Грасгофа Ог и число Прандтля Рг имеют вид:
Ог =
_ gв®R3
,Рг = -.
(5)
V х
Скорость течения связана с полем функции тока следующим образом:
у = |1 ^,,0 |.
'г дЭ дг 1
(6)
Компоненты единичного вектора п, задающего направление вверх, в цилиндрической системе координат имеют вид (см. рис. 1):
п = (ео8(а - Э), 81п(а - Э), 0). (7)
Из условия непротекания через границу полости (уг = 0) и прилипания на ней (ув = 0) получаем из соотношения (6) граничные условия для функций тока:
1 дд при г = 1 д = —— = 0.
дг
(8)
Будем полагать, что массив вокруг полости имеет коэффициент теплопроводности много больший, чем у жидкости. Тогда граничное условие для температуры можно записать в виде:
при г= 1 Т = со8(Э - а).
(9)
Метод решения
Решение задачи получалось методом сеток. Использовался итерационный метод решения стационарных уравнений, полученных из уравнений (2) — (4):
1 (дф дд дф ддЛ г ^дг дЭ дЭ дг) (10)
= Дф + Ог ^ з1п(а-Э) С0^-Э) I дг дЭ г
Дд + ф = 0; (11)
1 (дТ Эд-^12) г I дг дЭ дЭ дг ] Рг ,
оператор Ла-
л д2 1 д 1 д2
где Д = —2 +--+ —2 —:
дг г дг г дЭ' пласа.
Введем сетку с постоянными шагами Дг и Д9* по радиальной и угловой координатам, соответственно:
г= (I -1)Дг, Э, = (к -1)ДЭ*, ь = 1, 2, ..., I, к = 1,2, ..., К;
Дг = 1 / (I - 1), ДЭ* = 2п / (К - 1), ДЭ = 2 81п(ДЭ* / 2).
Конечно-разностный аналог системы (10) — (12) для внутренних узлов сетки, за исключением центрального, запишется в виде:
1
(13)
— (ф.д. -фд) = Бф +
(14)
+ Ог
Т ап(а - Эк) - Т.
со8(а-Эк)
Од + ф = 0;
— (Т д. - Т. д.) = БТ. г
(15)
Э Э г (16)
Здесь и далее использованы обозначения А.А. Самарского [7]:
Лк = /(г, Эк), /0 = (/+1 - /-1) / 2Дг,
г
/0 = (/к+1 - /к-1) / 2ДЭ;
Э
/ээ = (/+, -2/ + Л-,)/ДЭ2, (17) и = (/+1 - 2/ + / _,)/ Дг2.
Б — конечно-разностный аналог оператора Лапласа:
Б/ = /. +1 /. + /-.
гг ^ г г . ЭЭ
Выделяя в каждом из уравнений (14) — (16) значения в центре шаблона, получаем выражение для вычисления сеточных функций в центральном узле шаблона:
/; ,к
а/-1, к + Ь/+1, к + / к-1 + ^ л
+ е
,(18)
2 / Дг2 + 2 / {г. ДЭ)2 где / — одна из сеточных функций: ф, Т или д.
Коэффициенты a, b, c, d для равенств (14) определяются соотношениями:
a =
J_ Дг2
1 1
--+--v 0;
2г. Дг 2г. Дг s
i 1 1
b = -т + -
1
Дг2 2г Дг 2г Дг V
c=
1
1
V о
(г- Д$)2 2г Д$ о
d=
1
1
■ V о
(Г- ДЗ)2 2г Д^ т°"
Последние слагаемые в аналогичных выражениях для уравнения (16) домно-жены на число Прандтля, а в соотношениях, аппроксимирующих уравнение (15), они отсутствуют. Коэффициент е для (16) равен нулю, для (15) он равен значению завихренности, а для уравнения (14) определяется соотношениями:
еф = Ог[Тг 81п(а - З,) + Г8 со8(а - З,) / г-]. (20)
При получении расчетных формул для центральной точки, уравнения (14) — (16) записывались в инвариантной форме, и к ним применялся метод контрольного объема.
Граничные условия (8) аппроксимировались с использованием формулы Тома для завихренности, и в разностной форме записывались в виде:
, к
VlJC = 00, Ф1 k = --
Дг2
(22)
Граничное условие для температуры на границе (9) записывалось в виде:
Т1, = сой(З, -а). (23)
При решении уравнения Пуассона (15) для функции тока у использовался метод верхней релаксации с параметром релаксации ш = 1,5, а при решении уравнений (14) и (16) полагали параметр релаксации равным единице.
Разностная схема аппроксимирует дифференциальные уравнения со вторым порядком
0(Дг2 + (г ДЗ)2).
К ее особенностям следует отнести использование нетрадиционного выражения для шага по углу (13) в разност-
ных соотношениях, которые аппроксимируют производные по углу; такое выражение позволяет учесть физические особенности задачи в цилиндрической системе координат [9].
Опишем процедуру решения.
Шаг 1. Задаем начальные условия для цилиндрической полости:
Ti, = r¡ cos | k
K
=- б4(r> 2 -1)2'°гsin a;
Ф; ,k = 8(2 ri2 - 1)^г sin a.
Шаг 2. Вычисляем в центре ф5, T5,
Шаг 3. Считая ф5, T5, и F5 известными во всех узлах сетки, осуществляем одну итерацию для нахождения ф5+1, T^1 во внутренних узлах. Итерационным путем находятся значения ys+1. Новые граничные значения для ф определяются по формуле Тома (22).
Шаги 2, 3 повторяются до получения установившихся значений.
Использовалась экстраполяция по параметрам Gr и a. При наличии двух состояний с параметрами (Gr1, a1) и (Gr2, a 2), для получения решения с Gr и a, расположенными достаточно близко к этим двум, в качестве начального бралось состояние, полученное линейной экстраполяцией во всех узлах сетки из этих двух имеющихся состояний.
Результаты расчетов
Результаты, приведенные на рис. 2, дают представление об интенсивности конвекции и при произвольных направлениях внешнего градиента температуры. Представлены бифуркационные диаграммы, отражающие зависимости экстремального значения функции тока vm от угла a между направлениями единичного вектора n (вверх) и подогрева у для нескольких значений числа Грасгофа. Кривая 1 соответствует докрити-ческому (при a = 0°) значению числа Грасгофа, поэтому при a ^ 0° имеем ^ 0. На рис. 2 (на вставках) также приведены линии тока и изотермы для Gr = 30 и углов, соответствующих подогреву строго сбоку. Вид-
Рис. 2. Бифуркационные диаграммы для Рг = 7 и различных значений числа Грасгофа:
Ог = 30 (1), 65 (2), 100 (3), 200 (4). На вставках представлены также изотермы (I, III) и функции тока (II, IV) для Ог = 30 при двух значениях угла а: — 90° (I, II) и +90° (III, IV)
но, что линии тока практически круговые и у наиболее нагретой части стенки жидкость движется вверх. Такое течение естественно назвать нормальным, для него движение частиц жидкости при значениях —180° < а < 0 осуществляется против часовой стрелки, а при 0 < а < 180° — по часовой стрелке. Течение 1, представленное на рис. 2, является нормальным, так как направление циркуляции осуществляется по часовой стрелке. Неустойчивые участки кривых 2 — 4 не отмечены на рисунке ввиду невозможности расчета этих участков данным методом.
Максимум интенсивности тепловой конвекции при Ог = 10 достигается при а = —90° и а = 90° (этот случай не показан на рис. 2). Как видно из рисунка, увеличение числа Грасгофа приводит к их смещению в сторону угла а = 0. Функция ут (а) является двоякопериодической.
Напомним, что при подогреве снизу критическое число Рэлея равно Кас = 406,8 [5]. Поэтому критическое значение числа Грасгофа Ог, = Кас/Рг, соответствующее расчету, равно 57,1. Значения же Ог = 65, 100 и 200 — надкритические; при а ^ 0 им соответствуют конечные значения ут. Кроме того, при этих надкритических значениях числа Ог вблизи угла а = 0 имеется неоднозначность. Как видно из рис. 2, для
углов наклона от —ас до +ас существует два устойчивых состояния движения, близких по форме, но отличающихся друг от друга как интенсивностью, так и направлением циркуляции. Зависимость критического угла ас от числа Рэлея принято называть бифуркационной кривой.
На рис. 3 представлены бифуркационные кривые для трех значений числа Прандтля. Все кривые имеют максимум при Яа = (2 — 3)Кас. С увеличением числа Рэлея присутствует тенденция к насыщению. Максимальное значение критического угла слабо зависит от числа Прандтля, а его асимптотическое значение уменьшается с его увеличением.
Рис. 3. Бифуркационные кривые для значений числа Прандтля Рг = 0,7 (1); 7 (2); 14 (3) (см. пояснения в тексте)
Таким образом, в результате численного исследования получено семейство бифуркационных кривых, отражающее влияние теплофизических свойств на область существования аномального режима конвекции в цилиндрическом горизонтальном канале.
1. Показано, что аномальное конвективное течение может существовать при значительных (до 10е) отклонениях от условий подогрева снизу.
2. При плавном изменении угла наклона переход от нормального течения к аномальному происходит плавно, а от аномально-
го к нормальному скачкообразно. Глубина гистерезиса, наблюдающегося при таких переходах, отражающаяся бифуркационной кривой, максимальна в области двух-трех надкритичностей по числу Рэлея и уменьшается с его увеличением.
3. Проанализировано влияние числа Прандтля на ширину интервала гистерезиса. Показано, что максимальное значение интервала гистерезиса слабо зависит от числа Прандтля. С увеличением числа Рэлея, после достижения своего максимального значения, глубина гистерезиса уменьшается тем быстрее, чем больше число Прандтля.
СПИСОК ЛИТЕРАТУРЫ
1. Шарифулин, А.Н. Лабораторное моделирование нелокального возникновения тропического циклона [Текст] /А.Н. Шарифулин, А.Н. Полудницин, А.С. Кравчук // Журнал экспериментальной и теоретической физики. — 2008.
- Т. 134. - № 6. - С. 1269-1273.
2. Princevac, M.H. A criterion for the generation of turbulent in abatic flows [Text] / M.H. Princevac, J.S. Fernando // Phys. Fluids. - 2007.
- Vol. 19. - P. 105102.
3. Шарифулин, А.Н. Вибрационная конвекция в цилиндрической полости в невесомости при произвольных ориентациях направления подогрева [Текст] / А.Н. Шарифулин // Конвективные течения. - Пермь: Пермск. пед. ин-т, 1981. - С. 22-29.
4. Никитин, А.И. О бифуркациях стационарных режимов тепловой конвекции в замкнутой полости, порождаемых особенностью типа
сборки Уитни [Текст] / А.И. Никитин, А.Н. Шарифулин // Процессы тепло- и массопере-носа вязкой жидкости. - Свердловск: УНЦ АН СССР, 1986.- С. 32-39.
5. McHugh, J.P. The onset of convection in horizontal cylinders [Text] / J.P. McHugh // Quarter^ of applied mathematics. - 2000. - Vol. LVIII.
- № 3. -P. 425-436.
6. Чернатынский, В.И. Численное исследование конвекции в горизонтальном цилиндре кругового сечения [Текст] / В.И. Чернатынский // Гидродинамика. - 1974. - № 7. - С. 65-82.
7. Самарский, А.А. Введение в теорию разностных схем [Текст] / А.А. Самарский.
- М.: Наука, 1971. - 552 с.
8. Тарунин, Е.Л. Вычислительный эксперимент в задачах свободной конвекции [Текст] / Е.Л. Тарунин. - Иркутск: Изд-во Иркут. ун-та, 1990. -228 с.
ШАРИФУЛИН Альберт Нургалиевич — кандидат физико-математических наук, доцент кафедры прикладной физики Пермского национального исследовательского политехнического университета.
614990, г. Пермь, Комсомольский пр., 29 (342) 239-12-83 а1Ъе11.8Ьаг1М1п@§та11.ги
ФОМИНСКИЙ Дмитрий Александрович — аспирант кафедры прикладной физики Пермского национального исследовательского политехнического университета. 614990, г. Пермь, Комсомольский пр., 29 (342) 239-12-83 га91Ш;@уапдех.ги
© Санкт-Петербургский государственный политехнический университет, 2013