Моделирование двумерной конвекции Рэлея - Хэдли Текст научной статьи по специальности «Математика»

УДК 551.513

МОДЕЛИРОВАНИЕ ДВУМЕРНОЙ КОНВЕКЦИИ РЭЛЕЯ - ХЭДЛИ

© 2012 г С.А. Сухов

Ставропольский государственный университет, Stavropol State University,

ул Пушкина, 1, г. Ставрополь, 355009, Pushkin St., 1, Stavropol, 355009,

[email protected] [email protected]

Приводится система нелинейных дифференциальных уравнений в приближении Обербека—Буссинеска, описывающая двумерную конвекцию в замкнутой полости в присутствии горизонтального и вертикального градиентов температуры. Применен конечно-разностный численный метод для решения системы. Получены картины предельных стационарных режимов при различных числах Рэлея. Исследованы интегральные характеристики (число Нуссельта) в зависимости от числа Рэлея.

Ключевые слова: свободная конвекция, циркуляция Хэдли, число Рэлея, число Хэдли, число Нуссельта, вихрь, возмущение температуры.

The system of nonlinear differential equations for the Oberbec—Boussinesq approximation that describes the two-dimensional convection in a closed domain in the presence of horizontal and vertical temperature gradients is represented in the paper. The finite-difference numerical method is used to solve the system. The marginal steady-state conditions patterns are obtained for different Rayleigh numbers. Integral characteristics (Nusselt number) had been investigated subject to Rayleigh number.

Keywords: freee convection, Hadley circulation, Rayleigh number, Hadley number, Nusselt number, curl, temperature disturbance.

В гидродинамике в основном исследуются 2 вида конвекции. Каждый представляет собой атмосферный феномен, который обычно наблюдается в широких пространственных и временных масштабах. Первый возникает, когда вязкая среда нагревается снизу; как только амплитуда нагрева превышает некоторое критическое значение, среда самоорганизуется в систему отдельно вращающихся вихрей (конвекционных ячеек) [1, 2]. Это естественная конвекция, или конвекция Рэлея-Бенара. Ее примером может служить образование облаков, находящихся на вершине пограничного слоя атмосферы (1-1,5 км), в виде «облачных улиц» или пчелиных сот.

Второй тип конвекции возникает благодаря нагреву по горизонтали. В лабораторном эксперименте, в котором одна стенка вращающегося кольцевого канала нагревается, в то время как другая стена охлаждается, развивается двумерная ячейка в радиально-вертикальной плоскости [3, 4]. Прототипом такого течения является циркуляция Хэдли [5], наблюдаемая в северном и южном полушариях в поясе от экватора до 20-30 ° широты. Таким образом, можно назвать движения, поддерживаемые неоднородным нагревом по горизонтали, конвекцией Хэдли. Циркуляция Фе-релла, вообще говоря, не являются термической конвекцией, а возникает вследствие саморегуляции движений атмосферы.

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

мере, малые температурные градиенты всегда присутствуют в обоих направлениях.

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

Рассмотрим задачу о влиянии горизонтального градиента температуры на устойчивость и форму конвективных ячеек в жидкости.

Постановка задачи

Рассмотрим ограниченную область

О = {(х, 2),0 < х < Ь,0 < 2 < Н} с аспектным отношением А = Н/Ь, где Н и Ь - высота и ширина прямоугольной полости О .

Представим поля величин в виде суммы основного состояния и возмущения р(х, г, /) = р(г) + ~(х, г, /),

Т (х, г, /) = Т (г) + ~(х, г, /), р(х, г, /) = р + р(х, г, /) = р - ррб(х, г, /), где р(г) = ро - - гидростатическое давление; р -среднее значение плотности; р - коэффициент объемного расширения. Потребуем, чтобы |~>| << |р|,

|р| << |р|, |б| << |Т|. В невозмущенном состоянии отсутствует всякое движение жидкости, т.е. V = 0.

Отклонение температуры от линейного профиля

Т (х, г) = То +ДгТ - ^гТ 2 (То - начальная темпера-

тура на верхней границе;

A T H

- заданный верти-

кальный градиент температуры) положим равным

е(х, t) = е(х, t) + АхТ х, где х - горизонталь-Ь Ь

ный градиент температуры; е(х, г, ^ - возмущение

температуры среды на сообщенный ей нагрев.

Для условий сухой атмосферы выражение

A„ T

AT

H ,

имеет смысл вертикального градиента по-

50,

■ = ay = (y a-у),

тенциальной температуры © р ,

р дг

где у а - сухоадиабатический градиент температуры. Тогда поле температуры можно переписать в виде А, Т

T (x, z, t ) =

To + A, T--

A xT

-z +—x— x

+ e(x, z, t).

н ь

Из уравнения неразрывности V • V = 0 следует, что

, дш дш

можно ввести функцию тока ш , и =--—, ^ = —- .

дг дх

Тогда уравнения мелкой конвекции в форме Обер-бека - Буссинеска запишутся в виде

д„ 2 д(ш, V 2у) 0 д9 _ А хТ л —+ —г-^-Ря--—--vV4ш = 0,

5t

d(x, z)

dx

L

5e+5iv,el-CM AT-CMAZ-Kv2e = o, (i)

5t d(x, z) dx H dz L

где

d(a, b)_f da db db da

- якобиан; к и v - ко-

д(х, г) ^дх дг дх дг _ эффициенты температуропроводности и кинематической вязкости соответственно.

Зададим жесткие граничные условия для Б . Вертикальная и горизонтальная компоненты скорости обращаются в нуль, т.е.

М x=0,l = М z=0,1 =

Cy Hz

z=0,1

dy dx

= 0.

x=0,1

x = Lx , z = Hz * ,

t = (LH / к)/*,

*

У = Kf ,

0 = (kv/gpH3 )э*. Применяя (2) к (1), получим

(2)

J * V, * 2 * C М , v М

d * 2 * I

TTV М +-J * * \

dt dlx , z I

* / * * \ * *

- Ra ^Мг-Ha ^Мг-1 v* 2e * = 0. (3)

dt dlx , z )

„de Pr *4 *

* *! - Pr —r- Ha - — V М = 0,

dlx , z ) dx A

dx Pr cz A

Здесь Ra =

gpA zTH 3

Ha =

gPAxTH 3

Pr = V -к

kv k2

числа Рэлея, Хэдли, Прандтля. Безразмерные параметры Ra и Ha характеризуют тепловое воздействие по вертикали и горизонтали соответственно.

*

В безразмерных переменных область D стано-* If * * \ * * I

вится квадратной D = Ix , z 10 < x < 1,0 < z < 1i.

При Ha Ф 0 система (3) неоднородна и не имеет тривиального решения. Поэтому при появлении незначительного нагрева по горизонтали невозможно сохранение покоящегося состояния среды. Из систе-

Hal

мы (3) видно, что при Ra >> J—- преимущественно

Pr

наблюдается классическая рэлеевская конвекция, при

Ha |Ha|

Ra << J—- - циркуляция Хэдли. Если Ra Ф-—- Ф 0 , Pr Pr

получаем смешанную конвекцию Рэлея - Хэдли. Таким образом, модель (3) более полно описывает наблюдаемые конвективные течения в жидкости.

Численный метод

Рассмотрим плоское конвективное движение в

*

квадратной области D с аспектным отношением |Ha|

A = 1 при Ra =J—1. Для исследования надкритиче-Pr

ских движений обратимся к нелинейным уравнениям нестационарной конвекции. В отличие от (3) эти уравнения удобно записать, не выделяя из температуры и давления равновесные части.

В переменных функция тока - температура система уравнений нестационарной конвекции запишется в виде (далее звездочки опущены)

dV2M 1 - + ■

dt

Pr

dy dV2y dy dV2M^ dz dx dx dz

= V4m- Ra CT,

dx

Возмущения температуры на границах исчезают

е(0, г, Л = е(ь, г, t )= 0; 0 < г < Н ] е(х,0, t ) = е(х, Н, /) = 0; 0 < х < ь\

Введем безразмерные переменные (отмечены звездочкой):

PrdT+[ CMCT_ CMCT | =у2т

dt I dz dx dx dz I

(4)

Для численного решения предпочтительно ввести новую переменную О = -Аш (О имеет смысл проекции вихря скорости на направление, перпендикулярное плоскости движения) и добавить возмущение

температуры в безразмерном виде е = (А гТ )е* . Тогда система (4) запишется в виде

5Q 1 [ 5м SQ 5м 5Q1 „ dT -=--1 —-----l + V2Q + Ra—

dt Pr I 5z dx dx dz I dx

5T _ 1 [ 2T _ I 5М CT 5М CT

dt Pr

5z dx dx dz

(5)

Q = -Ay, e = t - rtop -1+z -

Ha PrRa

где Top =

A zT

- некоторое значение безразмерной

температуры на верхней горизонтальной границе при числе На = 0. В дальнейших расчетах Т0р =1.

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

I ду

Ч z=о = dy

z=0

= 0, T о= 2"*, Ч ,= —

h=о =i dz

= 0.

T\z=! = 1 - * , 4х=0,1 =

dy дх

= 0, T

х=0,1

х=0

z=1

= 2 - z ,

T

х=1

= 1 - z .

(6)

Задав некоторое начальное распределение вихря скорости и температуры, при помощи уравнений (5) и краевых условий (6) можно проследить за эволюцией этого начального распределения и, в частности, получить предельный стационарный режим.

Для численного решения задачи применим метод конечных разностей. Введем пространственно-временную сетку

Г х, = ¿И, 1

\ г, к = 0,1,...М, И = —,

\гк = кИ М

tn = т,

n = 0,1,...N , т = —

N

и обозначение /(хг, гк ¿п ) = /"к .

Заменяя производные по времени односторонними разностями, а производные по координатам - цен-

df

тральными,

dt

s-n+1 _ у-П

Ji,k Ji,k

+ 0(т),

i,k

df

дх

fi+1,k fi-1,k , 2)

2h

Qk=^-k+■№+Ra T+1,k - T-u

+ O\h ), получим

Ra fj-r n _ rpn )_

2h

—L2 [Ч k+1 -<k-1 -Q-u )-

4Prh : : А : : /

-4+u -4?-1,k htk+1 -ou-1)

rpn + 1 _ rpn , T i,k = Ti, k +

(7)

(8)

W - А k++1 -Чп,+-1 )T+1,k - T--1 k)

Pr [ 4h

-(ч?+и I^Tk+1- t"-1

en!k = T^ -Top -1+khih.

^ PrRa

(9)

Используя формулу д /

дх2

fi+1,k +fi-1,k - 2fink

ho(h2),

i,k

аппроксимируем оператор

¥ fхх + fzz

следующим

fk =

fi

+1,k + fi-1,k + fi "k+1 + fl\k-1 - 4fi

Приведенная схема (7) - (9) имеет погрешность аппроксимации порядка о(т + И2).

Граничные условия для функции тока и температуры в конечно-разностной форме имеют вид

V 0, к = Ум ,к = Уг,0 = Уг,М = 0 Т0, к = 2 - кИ ,ТМ,к = 1 - кИ

Т,0 = 2 - ¿И,Т,м = 1 - ¿И

Значение вихря на стенке получается из условия Вудса 2-го порядка точности 3(^+1 1.

h2

'4+1

+ o(h2).

Лапласа образом:

Пк).

Уравнение Пуассона (8) решалось итерационным методом Либмана (Гаусса - Зейделя).

ш(«+1) = I (>) +ш(«) +ш(«) +ш(«) + И2^,)

Уг,к = +1,к + -1,к + Уг,к+1 + Уг,к-1 + И Ь%к/.

Шаг по времени т выбирался таким образом, чтобы при дальнейшем его уменьшении поля величин оставались неизменными. Основные расчеты проводились на сетке 20 х 20. Вычисления на более грубой сетке 16 х16 и мелкой - 40 х 40 показали достаточную точность полученного численного решения. При На = 0 полученные поля функции тока и температуры совпали с имеющимися расчетными данными в [6] при Яа = 5300 , Яа = 8000, Яа = 60000 . Расчетные графики выполнены на сетке 40 х 40.

Для поиска стационарных режимов задавалось начальное распределение вихря скорости и затем отслеживалась его эволюция на безразмерном временном отрезке [0; 1] .Такого промежутка было достаточно для установки предельного течения. Число Прандтля во всех расчетах принималось равным 4,0 (для воды Рг = 4,8). Начальное распределение температуры было выбрано в виде Т = 2 - х - г.

Обсуждение результатов

В ходе численных расчетов показана однозначная зависимость направления циркуляции от знака числа На . Задаваемое начальное распределение вихря скорости в виде локального вихря или его отсутствия ^(х, г,0) = 0 оказалось несущественным для достижения предельного стационарного режима (рис. 1). Развиваясь, начальное возмущение с циркуляцией, обратной предельному стационарному режиму, меняло первоначальное направление вращения на противоположное, и после переходного процесса, продолжительность которого примерно составляла 0,2 безразмерного времени, устанавливалось предельное состояние.

С ростом числа Рэлея увеличивалась область изотермического ядра с температурой Т = 1, что объясняется конвективным перемешиванием. Об интенсивности течения можно судить по максимальным значениям функции тока (рис. 1а).

n

х

n

k

n

+

2

h

в г

Рис. 1. Предельный стационарный режим в квадратной области при Яа = 60000, На < 0. Представлены изолинии: а - функции тока; б - вихря скорости; в - возмущения температуры; г - полной температуры

Проводились исследования и в условиях колеба- общего центра (рис. 2а). Теплые потоки жидкости, тельного режима конвекции при больших числах Рэ- вытесненные холодным потоком, поднимались на-лея. При Яа = 600000 была обнаружена система из верх и перемещались в область пониженной темпера-

двух интенсивно вращающихся вихрей относительно туры, где затем диссипировали.

Рис. 2. Мгновенные поля: а - функции тока; б - полной температуры при Яа = 600000 , На < 0 , t = 0,26

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

Основной интегральной характеристикой конвективного теплообмена является среднее число Нус-сельта, служащее мерой интенсивности конвективного движения.

Для задачи (5), (6) число Нуссельта в момент /" рас-

считыгалось по формуле = -Т(х,0,^- Т(0,и")к .

0 0

На рис. 3 сравниваются результаты расчетов по установлению во времени среднего числа Нуссельта для различных чисел Яа . Интенсивность ного переноса тепла увеличивается с ростом числа Яа . При достижении Яа « 2 •Ю5 наблюдается колебательный режим; при дальнейшем увеличении Яа растет частота колебаний.

Рис. 3. Влияние числа Рэлея на установление среднего числа Нуссельта во времени

Результаты численных расчетов [7] показывают, что при достаточно больших числах Рэлея изменение интегральных характеристик, прежде всего Nu, слабо зависит от структуры движения. Зависимость Nu (Ra)

исследуется практически во всех работах как основная характеристика конвекции. Типичной для большинства работ является аппроксимация этой зависимости в виде Nu = aRa^.

Для задачи (5), (6) зависимость числа Нуссельта от числа Рэлея с графической точностью представляется формулой

Nu = 0,414 • Ra0,249 при 5 • 102 < Ra < 105 . (10) На рис. 4 изображены полученная зависимость числа Нуссельта от числа Рэлея (4-1) и эмпирические точки (Ra, ,Nu ,) (4-2).

Эмпирический закон O Tool и Silveston [8] при Pr = 10 и жестких границах имеет вид

Nu = 0,126• Ra0,305. Эксперименты по конвекции в

0 3

воде и воздухе [9, 10] показывают закон Nu ~ Ra , , а расчеты по двумерной модели бесконечного числа Прандтля со свободными граничными условиями -

Nu~Ra0,301 [11]. В других экспериментальных работах по конвекции в воде [12, 13] получены законы

Nu~Ra0,293 и Nu ~ Ra0,278 соответственно. В [14]

13

для турбулентного режима получено Nu = 0.104 • Ra' .

о I-1-1-1-1-

0 2xl04 4xl04 6xl04 SxlO4 lxlOJ

Рис. 4. Зависимость числа Нуссельта от числа Рэлея. 1 - по соотношению (10), 2 - эмпирические точки

В [15] по конвекции в горизонтальном слое со свободными граничными условиями в зависимости числа Нуссельта от надкритичности выделено 3 участка:

№ = 0,454 • Яа0,261 при 3,3 •Ю3 < Яа < 1,32 -105 ;

1,32 • 105 < Яа < 3,95 • 105 - переходной участок без вы-

0 302

раженного степенного закона, и № = 0,172 • Яа , при

3,95 -105 < Ra < 2,24 -107. Можно заключить, что в условиях ламинарной конвекции при увеличении промежутка аппроксимации величина b в степенном законе

Nu = aRab стремится к 0,3 .

Как указано в [15], данные расчетов с жесткими граничными условиями лучше соответствуют экспериментальным данным. Добавление горизонтальной неоднородности температуры вместе с жесткими граничными условиями усилило устойчивость конвективных течений, что является важным фактором их реализуемости в реальной жидкости (или атмосфере).

Вывод

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

Работа выполнена под научным руководством доктора физико-математических наук Роберта Гурге-новича Закиняна.

Работа выполнена в рамках реализации ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг. (№ 02.740.11.0739).

Литература

1. Гетлинг А.В. Конвекция Рэлея - Бенара. Структура и

динамика. М., 1999. 247 с.

2. Chandrasekhar S. Hydrodynamic and hydromagnetic sta-

bility. Oxford, 1961. 652 p.

Поступила в редакцию_

3. Riehl H., Fultz D. Jet stream and long waves in a steady

rotating-dishpan experiment: Structure of the circulation // Quart J. R. Met. Soc. 1957. № 356. P. 215-231.

4. Riehl H., Fultz D. The general circulation in a steady rotat-

ing-dishpan experiment // Quart J. R. Met. Soc. 1958. № 362. P. 389-417.

5. Лоренц Э.Н. Природа и теория общей циркуляции атмо-

сферы. Л., 1970. 259 с.

6. Гершуни Г.З., Жуховицкий Е.М. Конвективная устойчи-

вость несжимаемой жидкости. М., 1972. 392 с.

7. Математическое моделирование конвективного тепло-

массообмена на основе уравнений Навье-Стокса / под общ. ред. В.С. Авдуевского. М., 1987. 272 с.

8. Denton R.A., Wood I.R. Turbulent convection between two

horizontal plates // Intern. J. Heat and Mass Transfer. 1979. Vol. 22. P. 1339-1346.

9. Rossby H.T. A study of Benard convection with and without

rotation // J. Fluid Mech. 1969. Vol. 36. P. 309-335.

10. Fitzjarrald D.E. An experimental study of turbulent convec-

tion in air // J. Fluid Mech. 1976. Vol. 73. P. 693-719.

11. Malevsky A.V., Yuen D.A. Characteristics-based methods

applied to infinite Prandtl number thermal convection in the hard turbulent regime // Phys. Fluids. A. 1991. Vol. 3, № 9. P. 2105-2115.

12. Garon A.M., Goldstein R.J. Velocity and heat transfer mea-

surements in thermal convection // Phys. Fluids. 1973. Vol. 16, № 11. P. 1818-1825.

13. Chu T.Y., Goldstein R.J. Turbulent convection in a horizon-

tal layer of water // J. Fluid Mech. 1973. Vol. 60, pt 1. P. 141-159.

14. Турбулентная естественная конвекция в вертикальном

слое / С.С. Кутателадзе [и др.] // Теплофизика высоких температур. 1977. Т. 15. № 3. С. 545-553.

15. Палымский И.Б. Численное моделирование двумерной

конвекции при высокой надкритичности // Успехи механики. 2006. № 4. С. 3-28.

14 ноября 2011 г.