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

2011

ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА

Математика. Механика. Информатика Вып.1(5)

МЕХАНИКА. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

УДК 532.5:536.2:519.6

Численное исследование конвективного вихря в неоднородно нагретой вращающейся области

Д. А. Винокуров1, Е. Л. Тарунин2

Пермский государственный университет, Россия, 614990, Пермь, ул. Букирева, 15 1df6.626@gmail.com; 89058621711 2taranm@psu.m; (342) 2-237-10-31

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

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

Введение

Ярким примером природного явления с вихревым течением газа является тропический циклон. С целью изучения тропического циклона экспериментаторы в Перми выполнили серию лабораторных исследований конвекции во вращающейся цилиндрической кювете диаметром 30 см [1—5]. Эти эксперименты были начаты Богатырёвым Г.П. и продолжены его коллегами. В качестве рабочей жидкости использовались трансформаторное масло и его смесь с авиационным керосином. Конвективное течение возникало за счет нагрева центральной части дна кюветы. Наибольшее значение

разности температур достигало 35°. В некоторых опытах азимутальная скорость вращения

© Д. А. Винокуров, Е. Л. Тарунин, 2011

жидкости превосходила скорость твердотельного вращения примерно в 10 раз.

Численная трехмерная модель, схожая с лабораторной по объекту и предмету исследования, описывается в работе [6]. Рассматривалась задача о зарождении вихря во вращающейся кювете при неравномерном подогреве ее снизу. Представлены результаты расчета режима формирования интенсивного вихря из крупномасштабной конвективной ячейки. Полученный вихрь совершал циклоническое вращение, на порядок превосходившее скорость вращения кюветы, и обладал характерными чертами тайфунов в земной атмосфере. Дополнительно были исследованы стадии перестройки течений в вихре после отключения нагрева.

В работе [7] с помощью пакета ANSYS CFX исследовалось похожее явление при

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

В данной работе методами математического моделирования также исследуется конвективное течение во вращающейся кювете. Постановка задачи выбиралась близкой (по геометрии и значениям безразмерных параметров) к лабораторным экспериментам [4].

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

циклона называют в метеорологии "глазом

тайфуна", а при отсутствии вращения течение называют "костром Каина" [8].

1. Модель

Осесимметричное конвективное движение вязкой жидкости в приближении Бусси-неска [9] описывалось следующей системой уравнений в безразмерных переменных:

дф _ 1 ду дф ду дф 1 ду

д г д% дг дг д2 г д2 Ф

(1)

Л 1 ^ дТ

н------------) + Аф —- -ф + Gr-----,

дz г дг

л 2 ду г -ф+Ау-----------------= 0,

г дг

дУа _ К дУ дУа дУ дУа +

------------(--------------------------н

дt г ді дг дг ді

Уа ду Уа

+ ) + Ауа

г ді г

дТ _ 1,ду дТ ду дТ 1

-----------(-----•-----------• —) н---------А! ,

дї г дz дг дг дz Рг

где V = (уг , Уа, ) - вектор скорости, Т -

(2)

(3)

(4)

д2

1 д д2

температура, А = —- н-----------------+-.

дг г дг ді'

скорости ф определен формулой

вихрь

ф

_с^2 _ 1 д2у 1 д2у 1 ду

+-ІГ’ (5) г дг

а функция тока у введена следующим образом:

дг ді

г ді г дг

1ду 1ду

V =-----—, V =------—.

г г ді 1 г дг

(6)

Г еометрия рассматриваемой области

изображена на рис. 1. Дно и боковые границы области полагались твердыми и непроницаемыми. Верхняя граница считалась свободной, но не деформируемой. Повышенная температура задавалась на дне полости при г < а - R. Отсчет температуры осуществлялся от заданной температуры на верхней границе. Остальные части полости были теплоизолированными. Задача содержит три безразмерных критерия подобия - числа Грасгофа, Прандтля и Рейнольдса:

Ог =

ЕРЛТН3

V 2 Н 2П

, Рг = —, ке =------------------.

X

V

для

В качестве единиц обезразмеривания расстояния, скорости, времени, температуры и давления были выбраны соответственно: Н,

у/Н , Н 2/у, АТ, ру2/Н2, где Н - высота слоя жидкости, АТ - максимальный перепад температур.

О

тттЬ

л

а-П

Я1

Ь

Рис. 1. Геометрия исследуемой области

2

V

а

Вихрь скорости на твердых границах определялся по формуле Тома, на свободной верхней границе и оси он полагался нулевым:

2_

tír

P,0 = Vi,i, i = 1, Nr - 2,

^-1,j = - h> WNr-2,j , j = 1, Nz - 2, P\r=0= P\z=1=0

(7)

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

«\г =0. (8)

Граничные условия для температуры:

т\ =-г Т | = СГ | =о

\г<аЯ, ¿=0 ’ 1г>=ой, ¿=0 1г=R ’

СП СП

T\z=i = 0, 8T

— \r=0=WrrTo +— (T + Tzz (1 +

(9)

Sí

Pr

1

+ -• Pr-hz- \Wtf \)).

Граничные условия для азимутальной скорости:

Va \z=0= Va \r=R = 05^ Re • r,

dv„

\ z=1 = Va \ r=0 = 0.

(10)

Sn

Начальные условия соответствовали непрогретой покоящейся жидкости:

р(0) =w(0) = v (0) = T(0) =0.

11)

Эволюционные уравнения системы (1,3-4) решались по явным двухслойным схемам с использованием направленных разностей в конвективных слагаемых (центральные разности приводили к неустойчивости схемы). Порядок погрешности аппроксимации схемы O(z + h). Уравнение Пуассона для функции тока (2) на каждом шаге по времени решалось итерационным методом последовательной верхней релаксации [10]. Шаг по времени выбирался из соображений устойчивости используемой разностной схемы, но с некоторым запасом.

Пространственная сетка была равномерной. Число используемых узлов сетки по горизонтали и вертикали изменялось от 100 х 20 до 300 х 60. Оценка погрешности по кинетической энергии с использованием идеи Рунге-Ромберга показывает, что на сетке 100 х 20 относительная погрешность составляет

7%, а при 300 х 60 - 2%. В основном для расчетов использовалась сетка 100 х 20 .

В ходе вычислений определялись различные интегральные характеристики течения: кинетическая энергия (полная, меридиональная, вращательная, циклоническая и антицикло-ническая), максимальная скорость, максимальная функция тока, максимальное отклонение азимутальной скорости от твердотельного вращения. Расчет останавливался при установлении характеристик течения. Безразмерное время установления решения в зависимости от параметров менялось от 0.5 до 10.

2. Результаты

2.1. Формирование интенсивного циклонического вихря

Основными изменяемыми параметрами были числа Грасгофа (Gr <105) и Рейнольдса (Re < 100). Число Прандтля изменялось от 60 до 180. Указанные значения параметров близки к тем, что использовались в лабораторных экспериментах. Отметим, что полноценное сравнение результатов расчета с результатами лабораторных экспериментов едва ли возможно.

Во-первых, в расчетах не учитывается зависимость параметров жидкости от температуры. А при максимальной разности температур в 35° следует ожидать изменения числа Прандтля более чем в два раза.

Во-вторых, в эксперименте на верхней границе происходил теплообмен с воздухом, а в расчетах температура на верхней границе была фиксированной.

Были получены два типа характерных картин установившихся течений. В первом случае циклоническая и антициклоническая зоны примерно одинаковы по максимальным значениям азимутальной скорости, во втором - циклон вытесняет антициклон из центральной зоны и значительно превосходит его по максимальной скорости. Второй случай характерен для больших чисел Грасгофа. Наглядной характеристикой интенсивности циклона и антициклона является относительное отклонение азимутальной скорости от скорости твердотельного вращения: V — V

и а = —-------, где V = г^е/2 - азиму-

s

тальная скорость твердотельного вращения.

-0.5 -0.3 -0.2 0.0 0.2 0.3

1.0

0.8

0.6

0.4

0.2

0.0

Рис. 2. Изолинии у (а) и и а (Ь) для случая Ог = 5000, Рг =100, Яе = 30

4 5

-20.0 -16.0 -12.0 -8.0 -4.0 0.0

1.0

Рис. 3. Изолинии у (а) и и а (Ь) для случая Ог = 105, Рг =100, Яе = 30

Вид линий тока и и а для первого случая представлен на рис. 2, для второго - на рис. 3. Отметка на дне около г = 1.66 показывает границу нагреваемой зоны. Как видно из изолиний иа, в нижней части вблизи оси образуется зона циклона - жидкость движется быстрее, чем при твердотельном вращении (и а >0), а в остальной части - зона антициклона, жидкость движется медленнее (иа < 0). Черной сплошной линией на рисунках обозначалась граница между циклоном и антициклоном (изолиния и а = 0).

В лабораторных экспериментах [4] граница возникновения интенсивного циклонического вихря определялась следующим образом. Использовался поплавковый датчик, состоящий из трех плексигласовых стаканчиков. Ось вращения датчика совпадала с осью

С

-1.5 -1.2 -0.9 -0.6 -0.3

Ь

а

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

изменения угла наклона графика (Оа (Ог) или (Оа (Яе) считалась точкой образования интенсивного вихря.

В данной работе была сделана попытка осреднять угловую скорость аналогично лабораторным опытам. Однако на большинстве графиков а>а (Ог) не было заметно изменений угла наклона, аналогичных изменениям на графиках, полученных лабораторным путем. В то же время на рисунках и был отчетливо виден интенсивный вихрь вблизи оси.

О 20 40 60 8В 100 120

С1

0 20 40 50 80 100 120

0 20 40 50 80 100 120

Сг х 10'3

ъ

Рис. 4. Рг = 60, Яе = 20. Ект (Ог)

(а), Екг (Ог) (Ь, 1) и Ек (Ог) (Ь, 2)

Поэтому были испробованы другие способы определения границы режимов.

Во-первых, для определения точки возникновения вихря использовались графики за-

висимости кинетических энергий от числа Грасгофа - энергии меридионального движения Ект (рис. 4а), вращательного движения Екг

(рис. 4Ь, 1) и полной кинетической энергия Ек

(рис. 4Ь, 2). Ект монотонно увеличивается с

ростом числа Грасгофа. Екг и Ек отражают

более сложное поведение. Сначала полная энергия падает с ростом числа Грасгофа. Это может быть объяснено следующим образом.

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

Ре

Рис. 5. Граница существования интенсивного циклонического вихря (1, 2 — вычислительный эксперимент, 3 — лабораторный эксперимент [4])

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

после некоторого критического числа Ог начинает линейно возрастать.

Для нескольких чисел Рейнольдса каждым из двух способов были определены

Ог* и в результате были установлены границы существования интенсивного циклонического вихря. Граница, определенная по кинетическим энергиям, изображена на рис. 5, линия 1. Граница, определенная по площади циклонических

течений на поверхности, изображена на рис. 5, линия 2. Для сравнения линией 3 показан результат, полученный в лабораторном эксперименте Богатырева Г. П. [4]. Зона существования интенсивного циклонического вихря расположена выше этих границ. Как видно, правые ветви обеих линий, полученных в вычислительном эксперименте, качественно схожи с результатом лабораторных экспериментов.

Однако получить значения на левой ветви не получилось, хотя в лабораторном эксперименте они были.

2.2. Перестройка течения на оси цилиндра

Обсудим результаты исследования перестройки течения на оси цилиндра. При определенных значениях Ог и Яе движение жидкости меняется с восходящего на нисходящее. Зависимость вертикальной компоненты скорости от радиуса на уровне г = 0.75 изображена на рис. 6а для значений Ог = 7000 , Рг = 60, Яе = 80 . На рисунке отчетливо видны большие значения скорости, направленной вниз. Для сравнения на рис. 6Ь показан случай восходящего течения при других значениях чисел Ог и Яе .

Подобные течения существуют только при определенных значениях параметров. Граница существования этих течений показана на рис. 7, она определялась по знаку скорости Уг1 = Уг (г = 0, г = 0.75). Выше границы раздела реализуется течение с восходящей струей на оси цилиндрической области. Как видно по положению границы, нисходящее течение соответствует большим числам Рейнольдса. Вихрь (обсуждался в предыдущем разделе) при этих значениях параметра значительно ослабевает из-за большой скорости вращения кюветы, интенсивного циклонического вихря не образуется.

Для определения границы рассматриваемых режимов течения требуется большое число вычислительных экспериментов с малым шагом изменения параметров задачи. Затруднения связаны также с отсутствием стационарного течения при больших значениях чисел Грасгофа и Рейнольдса. Нестационарность часто проявляется в периодических колебаниях характеристик одного типа течения (восходя-

Не

Рис. 7. Граница перестройки течения у оси кюветы

щего или нисходящего).

Однако при некоторых значениях параметров со временем периодически меняется и тип течения. Пример колебательного режима с перестройкой течения на оси изображен на рис.8. Как видно, вертикальная компонента скорости Ул меняет знак на каждом периоде коле-

г>2

а

г;2

ь

Рис. 6. Уг(г) на уровне г = 0.75 для случая Ог = 7000 , Рг = 60,

Яе = 80 (а) и для случая

Ог = 28000 , Рг = 60, Яе = 40 (Ь) баний ~ 2.3 единиц безразмерного времени. На слайдах изотерм и функции тока в этом случае видна серия неустойчивых поперечных валов, возникающих в пограничном слое вблизи границы подогрева. Эти валы "проталкивают" нисходящее течение к оси цилиндра.

Рис. 8. ул() при Ог = 28000 , Яе = 160, Рг = 60

1.0

0.8

0.6

0.4

0.2

0.0,

т

-2.0 -1.0 0.0 1.0 2.0 3.0

а

-0.6 0.0 0.6 1.2 1.8 2.4

Ь

Рис. 10. Ог = 28000 , Яе = 160,

Рг =60, г = 4.09; а - у, Ь - и а

Рассмотрим структуры азимутальных и меридиональных течений для двух моментов времени, соответствующих положительной и отрицательной скоростям Уг1 на графике колебаний на рис. 8. Для восходящего течения и г = 2.38 на рис. 9 показаны линии тока (рис.9а) и поле и а (рис. 9Ь). Для нисходящего течения

и г = 4.09 аналогичные графики приведены соответственно на рис. 10а и рис.10Ь.

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

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

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

гт

-3.0 -2.0 -1.0 0.0 1.0 2.0

а

-0.8 0.0 0.8 1.6 2.4 3.2

Ь

Рис. 9. Ог = 28000 , Яе = 160 Рг = 60, г = 2.38 . а - у, Ь - и а

Вычислительные эксперименты показали, что граница режимов понижается при увеличении размеров теплового пятна на дне кюветы. Этот эффект понятен с точки зрения конвективной устойчивости равновесия подогреваемой снизу жидкости даже без вращения [9].

10

15

20

25

Рис. 11. Уг1(Ог) для Яе = 70, Рг =60

Перейдем к описанию вычислительных экспериментов, позволивших выяснить характеристики гистерезисных эффектов смены режима течения. На рис. 11 показаны зависимости Уг1 от числа Грасгофа. Линия 1 соответствует изменению числа Грасгофа от 2000 до 22000 с шагом 2000 с использованием метода продолжения по параметру. Линия 2 соответствует изменению в обратном направлении. Жирные линии соответствуют средним значениям Уг1 по времени, штриховые линии показывают минимальные и максимальные значения в процессе колебаний.

Как видно, амплитуда колебаний возрастает с ростом числа Грасгофа. При значениях Ог < 5000 колебания практически не заметны.

Заключение

В ходе вычислительных экспериментов получены следующие результаты:

1. Частично подтверждена граница существования интенсивного циклонического вихря, выявленная в лабораторных экспериментах.

2. Определена граница перестройки течения жидкости на оси цилиндра со сменой направления движения.

3. Выявлен колебательный характер течения жидкости при больших значениях чисел Грасгофа и Рейнольдса.

4. Выявлен эффект гистерезиса для перестройки течения.

Список литературы

1. Богатырев Г.П. Возбуждение циклонического вихря, или лабораторная модель тро-

пического циклона // Письма в Журнал экспериментальной и теоретической физики. 1990. Т.51, № 11. С.557-559.

2. Богатырев Г.П., Попова Э.В. Исследование поля скорости в лабораторной модели тропического циклона // Вестник Перм. ун-та. Сер. Физика. 1994. №2. С.141-149.

3. Богатырев Г.П., Смородин Б.Л. Физическая модель вращения тропического циклона // Письма в Журнал экспериментальной и теоретической физики. 1996. Т.63, №1. С.25-28.

4. Богатырев Г.П., Колесниченко И.В., Левина Г.В., Сухановский А.Н. Лабораторная модель процесса образования крупномасштабного спирального вихря в конвективно-неустойчивой вращающейся жидкости // Изв. РАН. Физика атмосферы и океана. 2006. Т.42, №4. С.60-466.

5. Batalov V., Sukhanovsky A., Frick P. Laboratory study of differential rotation in convective rotating layer // Geophysical and Astrophysical Fluid Dynamics. 2010. Vol.104, №4. P.349-368.

6. Иванов М.Ф., Поварницын М.Е. Численное моделирование эволюции интенсивных конвективных вихрей тайфунного типа во вращающейся жидкости // Изв. РАН. МЖГ. 2003. № 3. С. 69-77.

7. Сухановский А.Н. Формирование дифференциального вращения в цилиндрическом слое жидкости // Вычислительная механика сплошных сред. 2010. Т.3, № 2. С. 103-115.

8. Тарунин Е.Л., Шарапова А.М. Ветвление осесимметричного конвективного течения // Изв. РАН. МЖГ. 2010. № 2. С.23-33.

9. Гершуни Г.З., Жуховицкий Е.М. Конвективная устойчивость несжимаемой жидкости. М.: Наука, 1972. С.392.

10. Тарунин Е.Л. Вычислительный эксперимент в задачах свободной конвекции. Иркутск: Изд-во Иркут. ун-та, 1990. С.228.

Numerical investigation of convective vortex in locally heated rotating vessel

D. A. Vinokurov1, E. L. Tarunin2

Perm State University, Russia, 614990, Perm, Bukireva st., 15 'df6.626@gmail.com; 89058621711

2tarunin@psu.ru; (342) 2-237-10-31

The convection in rotating cylindrical vessel was studied with using methods of mathematical modeling. The original formulation of the problem corresponded to the laboratory setup made in Perm for a typhoon modeling. Conditions of the intensive cyclonic vortex formation near a cylinder axis and conditions of flow restructuring were studied. Borders of concerned flows existence were determined in a plane of Grashof and Reynolds dimensionless numbers.

Key words: convection; Boussinesq approximation; differential rotation; tropical cyclone; flow restructuring, two-field method.