CC BY
57
УДК 532.516.013.4:536.25
Вестник СПбГУ. Сер. 1. 2012. Вып. 4
КОНВЕКЦИЯ В ПЛОСКОМ СЛОЕ С МАКСИМУМОМ ПЛОТНОСТИ*
Д. В. Кузнецова1, И. Н. Сибгатуллин2
1. НИИ Механики МГУ им. Ломоносова, мл. науч. сотр., morven9@yandex.ru
2. НИИ Механики МГУ им. Ломоносова,
канд. физ.-мат. наук, ст. науч. сотр., sibgat@imec.msu.ru
Введение. Слой жидкости под действием перепада температур в поле силы тяжести может быть устойчивым или неустойчивым к возникающим возмущениям в зависимости от соотношения между силами вязкости и плавучести.
В природе часто встречаются ситуации, когда конвективно устойчивые и неустойчивые слои соприкасаются, вследствие чего происходит их взаимодействие. При взаимодействии устойчивой и неустойчивой областей в жидкости возмущения из неустойчивой области могут проникать в устойчивую. В таком случае это явление называется проникающей конвекцией.
Наиболее распространенный пример проникающей конвекции на Земле — это конвекция в пресной воде. При атмосферном давлении и температуре 3.98°С пресная вода имеет максимум плотности, поэтому слой воды, содержащий эту точку, может иметь устойчивую и неустойчивую области. В соответствии с экспериментальными данными зависимость плотности воды от температуры может быть аппроксимирована квадратичной функцией с максимумом при 4°С в диапазоне температур 0—10°С.
Если точка максимума плотности находится внутри слоя и движение отсутствует, слой можно условно разделить на две горизонтальные области. Граница между ними определяется вертикальной координатой, для которой частицы жидкости имеют температуру 4°С. Таким образом, для нижней области наиболее тяжелая жидкость находится сверху, а для верхней области — снизу. Соответственно, нижняя область может быть неустойчивой, а верхняя является устойчивой. В дальнейшем будем называть эти области соответственно статически неустойчивой и статически устойчивой.
Обзор некоторых задач конвекции с максимумом плотности можно найти в [1]. В указанных работах используются различные аппроксимации зависимости плотности от температуры, в частности, кусочно-линейная.
Задача о конвекции в слое воды с квадратичной зависимостью плотности от температуры была впервые математически сформулирована в статье [2], где рассматривались линеаризованная и слабо нелинейная постановки. Была также предсказана возможность существования субкритической бифуркации перехода от статического (кондуктивного) состояния к стационарному. Далее, стационарные и периодические режимы исследовались для этой или аналогичной постановки задачи о конвекции в горизонтальном слое воды [3-8]. Задача также исследовалась экспериментально, например, в [9], где с помощью интерферометрии изучалось возникновение конвекции в воде и формирование льда на одной из границ слоя, и в [10], где с помощью шли* Работа выполнена при частичной финансовой поддержке РФФИ (грант №12-01-31460 мол_а). Доклад на Международной конференции по механике «Шестые Поляховские чтения» 31 января — 3 февраля 2012 г., Санкт-Петербург, Россия.
© Д.В.Кузнецова, И. Н. Сибгатуллин, 2012
рен-метода и теневого метода были показаны как двумерные валы, так и различные трехмерные конфигурации. В [11] изучалось влияние положения точки максимума плотности в слое на критические значения числа Рэлея и амплитуды течений в прямоугольных ячейках. Для исследования решений, возникающих после потери устойчивости стационарным течением, использовался метод Ляпунова—Шмидта. В книге [12] рассмотрена проникающая конвекция в двумерном вертикальном слое, для которой найдены аналитические решения в случае модифицированной зависимости плотности от температуры. Некоторые результаты, касающиеся горизонтальных масштабов для стационарных и периодических решений, областей гистерезиса, квазипериодических решений и перемежаемости, содержатся в [13-16].
Наличие инверсии плотности в жидкости приводит к принципиально новым типам конвективного движения, влияние которых необходимо учитывать, например, при изучении термогидродинамических процессов в водных объектах. Как двумерные, так и трехмерные нестационарные течения при наличии инверсии плотности в настоящее время изучены недостаточно полно. В данной работе исследуются переходные режимы в конвекции с максимумом плотности в горизонтальном слое воды в двумерной постановке. В дальнейшем планируется изучение устойчивости двумерных движений в трехмерной постановке и описание особенностей трехмерных режимов.
Постановка задачи. Рассматривается горизонтальный слой воды высоты Н с границами г = 0 и г = Н (г — вертикальная координата). Границы являются свободными, непроницаемыми, изотермическими со значениями температур ТЬ на нижней границе и Ти на верхней границе. Предполагается, что тш(Тъ, Ти) < Т4 < тах(Ть, Ти), тогда точка максимума плотности Т4 = 4°С будет находится внутри слоя.
Уравнение состояния принимается аналогично [2]:
р = Р4 (1 - «4(Т - Т4)2) , (1)
где р — плотность, Т — температура, р4 —плотность при 4°С, Т4 = 4°С, «4 = 7.68 • 10-6 (°С)-2.
Переменные, описывающие движение, представляются в виде суммы функций в статическом состоянии и отклонений от этого состояния: / = /т + /о(г) + /', где /т — некоторое постоянное значение, /т + /о(г) —распределение в статическом состоянии и /' — отклонение.
Для системы уравнений Навье—Стокса, неразрывности и энергии, принимается приближение Буссинеска [17, 18]. Основными предположениями, позволяющими использовать приближение Буссинеска, являются следующие: 1) высота слоя мала по сравнению с характерным вертикальным масштабом функций, описывающих движение; 2) возникающие отклонения плотности и давления от статического распределения по порядку не превосходят соответствующих вариаций в статическом состоянии; 3) диссипация мала по сравнению с остальными членами в уравнении энергии.
Обезразмеривание проводится следующим образом: Т = (Ть — Ти) Т*, р = р4р*, V = (к/Н) V*, г = (Н2/к) г*, р = (р4к^/Н2)р*, (ж, у, г) = Н(ж*,у*,г*), где звездочки обозначают безразмерные переменные. Здесь V — скорость, к — коэффициент температуропроводности, г — время, р — давление, V — кинематическая вязкость, (ж, у, г) —пространственные координаты.
Поскольку функции в статическом состоянии удовлетворяют системе уравнений в приближении Буссинеска, используются безразмерные уравнения для возмущений
(штрихи и звездочки опущены) температуры и компонент скорости V = {и, V, т}:
ду й
— + (V • У)у = -аУр + <тДу + (7^(Г+ 2\-2г) е2,
дТ
— + (V • V) Т - ю = АТ,
= 0,
где используются безразмерные параметры
vк Л3 ' Ть — Т
и
Здесь а — число Прандтля, й — число Рэлея, Л определяет положение точки максимума плотности внутри слоя, и используются размерные параметры V — кинематическая вязкость, к —коэффициент температуропроводности, д — ускорение силы тяжести.
Коэффициенты вязкости и температуропроводности предполагаются постоянными и равными значениям для воды при 4° С. Тогда число Прандтля постоянно и равно а = 11.5968. Высота слоя в вычислениях фиксирована и равна 0.1 м, разность температур на границах варьируется. Точка максимума плотности в статическом состоянии находится посередине слоя г = 1/2, т. е. Л = 1/2.
Исследуется двумерная постановка при отсутствии горизонтального потока массы. Решение предполагается периодическим по пространству, что соответствует стационарным и периодическим по времени течениям при горизонтальном масштабе, намного большем вертикального.
Для выбора горизонтального масштаба ячейки периодичности проводились расчеты на больших горизонтальных длинах расчетной области при условии отсутствия касательных напряжений на вертикальных границах. Рассматривались значения длин от 2 х Ьо до 2 х 20Ьо, где Ьо = л/2—половина длины ячейки периодичности в классической конвекции по линеаризованной теории. В качестве начальных условий брались хаотически заданные по пространству возмущения. В расчетной области с большой длиной решение имело периодическую по пространству структуру. В качестве длины расчетной области для изучения перехода к стохастическим режимам был выбран пространственнный период решения, периодического по времени, равный 2Ьо. (Горизонтальный пространственный период стационарного по времени решения имеет длину в два раза меньше.)
Заметим, что данная постановка соответствует классической постановке задач конвекции с граничными условиями Рэлея. Это позволяет описать отличия от других типов конвективных течений, например, конвекции Рэлея—Бенара, двухдиффузион-ной конвекции. Подробное описание основных различий между конвекцией Рэлея— Бенара и проникающей конвекцией с максимумом плотности содержится в [15, 16].
Компоненты скорости V = {и, т} и температура Т представляются по методу Галеркина в виде части ряда Фурье с коэффициентами, зависящими от времени.
С учетом граничных условий разложение имеет вид: и = ^о итп(£)
эт(пшах)соБ(ппг), т = ^т=о Еп=о №'тп(г) сов(пшах) вт(ппг), Т = ^Е ттп(£) сов(птах) 8ш(ппг). Здесь а — волновое число, характеризующее горизонтальный масштаб ячейки периодичности, который равен 2/а.
Из-за наличия симметрий в тригонометрическом разложении численный расчет производится в области [0,1/а] х [0,1], а полная ячейка периодичности равна [—1/а, 1/а] х [0,1]. Далее результаты будут показаны только для расчетной области.
Разложение в ряд Фурье позволяет получить систему обыкновенных дифференциальных уравнений, для которой ставится задача Коши. Для численного интегрирования использовался метод Булирша—Штера (ВиИгесЬ—81юег) [19]. Нелинейные члены вычислялись с помощью псевдоспектрального метода [20], для чего использовались библиотеки быстрого преобразования Фурье [21].
Области гистерезиса. Рассмотрим эволюцию решений с возрастанием числа Рэлея. Для выбранной длины расчетной области 2Ьо было найдено несколько видов решений. Эти решения имеют различное число конвективных структур в расчетной области и поэтому отличаются своим характерным масштабом.
Решения могут быть показаны с помощью построения зависимости осредненного по времени числа Нуссельта от разности температур на границах Мит (ДТ), где ДТ = Ть — Ти и Мит = (дТ/дг}2=о. Здесь угловые скобки обозначают осреднение по ж в ячейке периодичности.
Рис. 1. Зависимость осредненного по времени числа Нуссельта от надкритичности (разности температур на границах или числа Рэлея): 1 — статический или стационарный режим; 2 — периодический режим; 3 — периодический-2 режим; 4 — квазипериодический режим; 5 — стохастический режим (Ь/Ьо = 2.0, М = 128 ... 512, N = 64 ... 256). Римскими цифрами I, II, III обозначены ветви гистерезиса.
Эта зависимость показана на рис. 1 с соответствующей шкалой для числа Рэлея вверху.
Римские цифры I, II, III обозначают ветви гистерезиса. Ветвь II является основной, так как на ней осуществляется переход к хаосу. На ветвях I и III сохраняется стационарный режим. (Для ветви III здесь приведены некоторые предварительные расчеты.)
При Д = Др = 3870.6 (ДТ = 0.59) на ветви II происходит бифуркация Пуанкаре— Андронова—Хопфа. Далее эта ветвь делится на две части, соответствующие перио-
дическому и стационарному режимам. Стационарное решение может быть получено, если взять предыдущее стационарное решение в качестве начального условия для расчета с более высокой надкритичностью. Оно является устойчивым по отношению к произвольным возмущениям ненулевых гармоник, но неустойчивым по отношению к возмущениям, являющимся произведением периодического решения для той же надкритичности и произвольных случайных чисел.
Стационарное решение ветви II становится неустойчивым к любым возмущениям после Д = Дц = 11119 (ДТ = 1.0), и после этого на ветви II существует только периодическое решение.
С возрастанием надкритичности после второй бифуркации Хопфа (бифуркации Неймарка—Сакера), которая является субкритической [16], возникает периодическое-2 решение на торе при Д = Д2р. Периодическое и периодическое-2 решения сосуществуют до Д = Д2р', когда периодическое решение исчезает. (В работах [13-15] данная бифуркация была названа бифуркацией удвоения периода. Впоследствии было обнаружено, что она является именно субкритической бифуркацией Неймарка—Сакера с синхронизацией частот.)
Для Д = Д1 (вблизи Д = Д2р и Д = Д2р') стационарные решения на ветви I становятся неустойчивыми и переходят на периодические-2 решения на ветви II.
Переход к хаосу. При Д = Дчр = 21176 (ДТ = 1.38) происходит переход от периодического-2 решения к квазипериодическому [13, 14], которое будет обозначаться д1.
На рис. 2 показаны изотермы для стационарного (а), периодического-1 (б), периодического-2 (в) и квазипериодического (г) режимов.
Для стационарного режима выделяется ячейка периодичности с длиной Ьо. Для области [0,1] х [0, Ьо] видна симметрия относительно ж = Ьо/2, аналогично для области [0,1] х [Ьо, 2Ьо] видна симметрия относительно ж = 3Ьо/2. Далее, для периодического режима длина характерного масштаба меняется и становится равной 2Ьо, при этом возникают колебания нижних частей изотерм, закрепленных в верхней части — «хвостов». Описанная выше симметрия относительно ж = Ьо/2 и ж = 3Ьо/2 исчезает, но сохраняется симметрия относительно ж = Ьо. Для периодического-2 режима и далее квазипериодического, эта симметрия также исчезает.
С увеличением числа Рэлея после режима д1 возникает перемежаемость с характерными сильными всплесками (например, теплового потока) на фоне основного квазипериодического режима.
На рис. 3 показан пример последовательности изменения изотерм для режима с перемежаемостью: а, б — для квазипериодического участка, в, г — для стохастического всплеска. Видно, что во время стохастического всплеска картина течения полностью меняется по сравнению с квазипериодическим участком.
На рис. 4 показана зависимость числа Нуссельта от времени для этого режима. Слева представлен график для 0 < г < 1000, справа — для 620 < г < 710.
Перемежаемость исчезает после Д = 25219 (ДТ = 1.506), и возникает новый квазипериодический режим д2. При возникновении этого режима среднее значение числа Нуссельта резко падает.
Режим д2 существует до Д = 25427 (ДТ = 1.5122). При дальнейшем увеличении надкритичности снова возникает режим с перемежаемостью, но на фоне режима д2.
Рис.2. Последовательность изотерм: а — для стационарного режима при Ть = 4.25, Ти = 3.75, М = 128, N = 64; б — периодического-1 при Ть = 4.375, Тц = 3.625, М = 128, N = 64; в — периодического-2 при Ть = 4.66, Ти = 3.34, М = 128, N = 64; г — квазипериодического Ть = 4.73, Тц = 3.27, М = 128, N = 64.
Заключение. Исследовались особенности проникающей конвекции с максимумом плотности в горизонтальном слое воды. Для выбранной длины ячейки периодичности при условии отсутствия потока через вертикальные границы показано существование двух основных ветвей гистерезиса, на которых решения отличаются характерным масштабом. Описан переход к хаосу через квазипериодические режимы и перемежаемость для длины расчетной области, равной горизонтальному масштабу для периодического по времени режима. Переход к хаосу происходит на одной из ветвей гистерезиса.
Эволюция движения описывается следующей последовательностью. Когда стационарный режим становится неустойчивым, возникает периодическое движение. Далее происходит удвоение периода. При дальнейшем увеличении надкритичности возникает квазипериодический режим. После этого на фоне основного квазипериодического режима возникает перемежаемость. В интервале чисел Рэлея, соответствующем режимам с перемежаемостью, существует окно квазипериодичности. Далее, перемежаемость возникает на фоне нового квазипериодического движения.
Рис. 3. Последовательность изотерм для режима с перемежаемостью: а, б — квазипериодический участок этого режима; в, г — участок, где наблюдается стохастический всплеск (ТЬ = 4.75, Ти = 3.25, М = 128, N = 64).
Рис. 4- Зависимость числа Нуссельта от времени для режима с перемежаемостью. Слева — график при 0 < I < 1000, справа —при 620 < I < 710 (Ть = 4.75, Ти = 3.25, М = 128, N = 64).
Открытыми остаются вопросы устойчивости полученных решений в трехмерном случае, а также при наличии горизонтального потока массы как в двумерном, так и в трехмерном случае.
Благодарности. Вычисления производились на суперкомпьютерных системах СКИФ МГУ имени М. В. Ломоносова и технологической платформе ИСП РАН UniHUB.
Литература
1. Гетлинг А. В. Конвекция Рэлея—Бенара. Структуры и динамика. М.: Эдитори-ал УРСС, 1999. 248 с.
2. Veronis G. Penetrative convection // Astrophys. J. 1963. Vol. 137. P. 641-663.
3. Musman S. Penetrative convection // J. Fluid Mech. 1968. Vol. 31. P. 343-360.
4. Блохин А. С., Блохина Н. С. Начало конвекции в жидкости вблизи температуры инверсии плотности // Докл. Акад. Наук СССР. 1970. Т. 193, №4. С. 805-807.
5. Блохин А. С., Блохина Н. С., Макаева О. С. Самовозбуждающиеся колебания в жидкости при развитой конвекции // Докл. Акад. Наук СССР. 1973. Т. 210, №1. С. 75-78.
6. Moore D. R., Weiss N. O. Nonlinear penetrative convection // J. Fluid Mech. 1973. Vol. 61(3). P. 553-581.
7. Blake K. R., Poulikakos D., Bejan A. Natural convection near 4 C in a horizontal water layer heated from below // Phys. Fluids. 1984. Vol. 27(11). P. 2608-2616.
8. Tong W., Koster J. N. Penetrative convection // Warme- und Stofftubertragung. 1993. Vol. 29. P. 37-49.
9. Tankin R. S., Farhadieh R. Effects of thermal convection formation of ice // Int. J. Heat Mass Transfer. 1971. Vol. 14. P. 953-960.
10. Large E. D. An experimental investigation of penetrative convection in water near 4C, PhD thesis. The Ohio State University. 2010.
11. Надолин К. А. Конвекция в горизонтальном слое жидкости при инверсии удельного объема // Изв. РАН. МЖГ. 1989. №1. С. 43-49.
12. Гершуни Г. З., Жуховицкий Е. М., Непомнящий А. А. Устойчивость конвективных течений. М.: Наука. Гл. ред. физ.-мат. лит. 1989. 320 с.
13. Кузнецова Д. В., Сибгатуллин И. Н. Переходные режимы проникающей конвекции в плоском слое воды // Докл. Акад. Наук. 2011. Т. 438, №1. С. 47-50.
14. Кузнецова Д. В., Сибгатуллин И. Н. Переходные режимы проникающей конвекции в плоском слое // Вестник Нижегородского ун-та им. Н.И.Лобачевского. 2011. №4, ч. 3. С. 907-909.
15. Кузнецова Д. В. Конвекция в плоском слое воды при наличии максимума плотности // Тепловые процессы в технике. 2012. №3. C. 105-110.
16. Kuznetsova D. V., Sibgatullin I. N. Transitional modes of penetrative convection in plane layer // Fluid Dyn. Res. 2012. Vol.44. 031410.
17. Надолин К. А. О проникающей конвекции в приближении изотермически несжимаемой жидкости // Изв. РАН. МЖГ. 1996. №2. С. 40-52.
18. Spiegel E. A., Veronis G. On the boussinesq approximation for a compressible fluid // Astrophys. J. 1960. Vol. 131. P. 442-447.
19. Press W.H., Teukolsky S.A., Vetterling W. T., Flannery B.P. Numerical Recipes: The Art of Scientific Computing / 3rd Edition. New York: Cambridge University Press, 2007. 1235 pp. + xxi
20. Флетчер К. Численные методы на основе метода Галеркина. М.: Мир, 1988. 352 c. (Fletcher C. A. J. 1984, Computational Galerkin methods. New York: Springer-Verlag.)
21. http://fftw.org (дата обращения: 12.03.2012).
Статья поступила в редакцию 26 июня 2012 г.
