Содержание к диссертации
Введение
1. Краткий обзор литературы 9
1.1 Прямоугольные полости 11
1.2 Цилиндрические полости 19
2. Конвекция жидкости в квадратной ячейке с учетом процессов плавления-замерзания 25
2.1 Постановка задачи 25
2.2 Математическая модель 27
2.3 Численная реализация метода энтальпии 32
2.4 Анализ результатов в задаче с линейным распределением температуры на боковых стенках 35
2.5 Анализ результатов в задаче с адиабатическими боковыми стенками 49
2.6 Обоснование существования несимметричного решения в задачах с симметричными граничными условиями 59
Выводы к главе 2 66
3. Конвекция жидкости в круглой полости 67
3.1 Постановка задачи 68
3.2 Математическая модель 70
3.3 Метод циклической прогонки 75
3.4 Граничное условие в центре 78
3.5 Аппроксимация источниковых членов 80
3.6 Тестовые расчеты 86
3.7 Различные решения задачи свободной конвекции в круглой полости 93
Выводы к главе 3 101
4. Конвекция жидкости в цилиндрической полости 102
4.1 Постановка задачи 102
4.2 Математическая модель 104
4.3 Тестирование программы 108
4.4 Анализ полученных решений 116
Выводы к главе 4 133
Заключение 134
Литература 136
- Цилиндрические полости
- Анализ результатов в задаче с линейным распределением температуры на боковых стенках
- Граничное условие в центре
- Тестирование программы
Введение к работе
Актуальность темы. Конвективные течения играют важную роль во многих природных явлениях и технологических процессах, в частности, при протаивании и промерзании грунта при строительстве в зоне вечной мерзлоты, почти во всех способах производства энергии; являются определяющими при обогреве и вентилировании зданий. В металлургических топках, химических реакторах, теплообменниках, конденсаторах происходит тепломассообмен при налігчии естественной конвекции в жидкостях и газах. Процессы фазовых переходов при наличии естественной конвекции происходят в аккумуляторах солнечной энергии и установках по выращиванию кристаллов.
Ввиду важности этгос процессов, управление ими является актуальной задачей. Для этого, в свою очередь, необходимо понимание существа этих явлений и методов их моделирования.
Целью работы являлось численное исследование естественной конвекции и фазовых переходов в прямоугольных и цилиндрических полостях и создание пакета программ для моделирования данных процессов.
Научная новизна данной работы заключается в следующем:
Показано, что в цилиндрической полости при определенных условиях существует, по крайней мере, пять типов решений, причем три из них не имеют двумерных аналогов. Получены значения суммарного теплового потока через цилиндрическую полость, а также максимальные и минимальные значения локальных тепловых потоков.
Получены диапазоны чисел Грасгофа, при которых существуют различные типы решений в задаче плавления-замерзания жидкости в квадратной ячейке при наличии естественной конвекции.
Обнаружено, что при малых числах Грасгофа существует тип решения, при котором конвекция приводит не к уменьшению, а к увеличению доли твердой фазы по сравнению с решением без конвекции.
Показано, что в круглой полости возможно существование трех различных типов нетривиальных стационарных решений, получены диапазоны их устойчивости.
На основе метода контрольного объема создан пакет программ, позволяющий моделировать двумерные и трехмерные процессы
в круглых и цилиндрических полостях. Предложен эффективный
способ аппроксимации источниковых членов.
Практическое значение работы. Результаты данной работы позволяют лучше поїить процессы, происходящие при наличии фазовых переходов. Диапазоны чисел Грасгофа, при которых существуют решения различных типов, дают возможность оценить вероятность появления того или иного типа течения, а также связанных с ним термических режимов. Показанный способ аппроксимации источниковых членов позволяет повысить точность расчетов.
Достоверность результатов, полученных в работе, обусловлена
использованием фундаментальных физических законов,
корректностью математических постановок задач, использованием современных численных методов, сопоставлением численных расчетов с экспериментальными и численными исследованиями других авторов.
На защиту выносятся:
Результаты моделирования естественной стационарной конвекции жидкости в квадратной ячейке при наличии фазового перехода (эффект увеличения доли твердой фазы с развитием конвекции, диапазоны существования решений).
Результаты исследования естественной стационарной конвекции жидкости в круглой ячейке (существование несимметричного решения, диапазоны существования полученных решений).
Результаты моделирования естественной стационарной конвекции жидкости в цилиндрической области, в частности, решения не имеющие двумерных аналогов.
Апробация работы и публикации. Основные результаты работы докладывались на международных конференциях «Всесибирские чтения по математике и механике» (Томск, 1997), «International Symposium on Multiphase Fluid, Non-Newtonian Fluid and Physico-Chemical Fluid Flow» (Beijing, China, 1997). Обсуждение результатов проводилось на семинарах Тюменского филиала Института теоретической и прикладной механики СО РАН (Тюмень), на семинаре Института проблем механики РАН (Москва), на семинарах лаборатории математического моделирования Сургутского государственного университета (Сургут).
По результатам диссертации опубликовано шесть печатных работ, список которых приведен в конце автореферата.
Структура и объем диссертации. Диссертация состоит из введения, четырех глав и заключения, содержит 142 страницы, 82 рисунка, библиографию из 62 наименований.
Цилиндрические полости
В данном разделе рассматриваются работы, исследующие течения в замкнутых полостях конечного размера. Двумерные течения в вертикальных, горизонтальных и наклонных прямоугольных полостях широко исследованы ввиду гораздо большей простоты этого случая по сравнению с моделированием трехмерных течений. Хотя внутренние естественноконвективные течения, представляющие интерес с точки зрения практических приложений (например, в условиях роста кристаллов, при расчете солнечных коллекторов или пожаров в помещениях), обычно являются трехмерными, соответствующие двумерные расчеты зачастую дают вполне удовлетворительные результаты. Кроме того, переход от двумерного анализа к трехмерному часто требует существенных затрат времени и средств. В связи с этим большая часть информации, касающейся течений в замкнутых полостях, была получена с помощью именно двумерных расчетов, а также из сравнительных экспериментальных исследований. В последние годы, однако, число работ, в которых проводятся соответствующие трехмерные расчеты, заметно выросло. Поэтому в обзоре литературы, касающейся течений в прямоугольных полостях сначала рассмотрены работы в двумерной постановке (с учетом или без фазовых переходов), а после - работы в трехмерной постановке.
В работе [27] (S. Devahastin, Z.X. Gong, A.S. Mujumdar, N. Arai (1998)) проводилось численное исследование процесса плавления парафина в квадратной ячейке. Такие процессы являются основными в устройствах хранения тепловой энергии. Для моделирования процесса фазового перехода использовался метод энтальпии. Две соседние стенки теплоизолированы, а на двух других поддерживается температура выше температуры фазового перехода. В работе показано, что вдоль нижней подогреваемой стенки образуется система вихрей, тогда как вдоль вертикальной стенки существует всего один вихрь. Также в работе исследовалось влияние поворота ячейки на 45. В частности, показано, что при таком повороте ячейки системы из множества вихрей не образуется, а появляется только два вихря, при этом процесс плавления ускоряется. Получено, что тепловой поток на подогреваемых стенках за счет влияния конвекции увеличивается в два раза при числе Релея, равном 2.84 -106.
Экспериментальное и численное исследование замерзания чистой и соленой воды в квадратной ячейке провели Tomasz А. Kowalewski и Jerzy Banaszek, Marek Rebow, Piotr Furmanski and Tomasz S. Wisniewski (1998) [28]. Одна боковая стенка подогревалась, а другая охлаждалась, причем температура бралась ниже температуры замерзания. Точка инверсии плотности воды попадала в рассматриваемый диапазон температур. В работе получено хорошее совпадение численных и экспериментальных результатов. Показано, что из-за наличия точки инверсии плотности образуются вторичные течения.
В работе А.В. Гудзовского (1999) [29] численно изучалась естественная конвекция воды и воздуха в вытянутой прямоугольной ячейке (отношение ширины ячейки к высоте равно 0.5). На небольшом участке снизу производился подогрев, а на аналогичном участке сверху - охлаждение. Обнаружено, что в широком диапазоне чисел Gr для данных граничных условий существует более одного колебательного режима, пространственно временные структуры которых существенно различаются. Выявлено, что два метода получения режимов при заданном Gr - посредством изменения параметра и путем эволюции из состояния покоя -приводят к разным выводам о количестве возможных режимов течения.
В статье [36] (Ekkehard Holzbecher (1997)) численно рассматривалась конвекция холодной воды в квадратной ячейке. Исследования проводились в широком диапазоне температур, максимальная разность между температурой холодной и горячей стенок была 40С. В работе учитывалась зависимость вязкости от температуры. Найдены критические числа Релея, соответствующие началу конвекции при различных граничных условиях, получена диаграмма зависимости теплового потока через ячейку от числа Релея.
П.Т. Зубков и К.М. Федоров (1994) [37] в своей работе рассматривают термогравитационную конвекцию жидкости при наличии фазового перехода. Процесс фазового перехода моделировался с помощью метода энтальпии. В работе показано, что при одном и том же числе Релея возможно существование нескольких типов решений. Различные решения появляются из-за разных типов начальных условий, а их устойчивость объясняется формой границы между твердой и жидкой фазами. Все решения отличаются различным теплопереносом через ячейку и долей твердой фазы.
Анализ результатов в задаче с линейным распределением температуры на боковых стенках
Зависимость доли твердого вещества и стационарных тепловых потоков от числа Gr являются наиболее интересными зависимостями с точки зрения анализа полученных решений и переходов решений из одного типа в другой (рис. 2.5.1, 2.5.2). В отличие от задачи с линейным распределением температуры на боковых стенках, рассмотренной в предыдущей главе, в данной ситуации каждый тип решения имеет свое собственное критическое число Gr, отвечающее началу конвективного движения. При малых числах Грасгофа существует только несимметричное решение. Для него критическое число Gr начала конвекции равно 7.8-103, что меньше, чем в задаче с линейным распределением температуры на боковых стенках (1.3 104). Уменьшение первого критического числа Gr можно объяснить тем, что линейное распределение температуры на боковых стенках, рассмотренное в предыдущем пункте, по сравнению с их теплоизолированностью является более жестким удерживающим условием.
Далее, при малых Gr симметричные решения можно получить только искусственным путем (постановкой условия симметрии на вертикальной линии, проходящей через центр квадрата). Эти "нестабильные" решения (на рис. 2.5.1 и 2.5.2 показаны пунктирными линиями) имеют следующие критические числа Gr: для решения с течением жидкости у стенок вверх (решение первого типа) Gr = l-104, для решения с течением жидкости у стенок вниз (решение второго типа) Gr = 1.15 104. Начиная с Gr = 2.1 104, вплоть до Gr = l-105, существует "стабильное" решение первого типа. Область существования "нестабильного" решения первого типа - от Gr = 1 104 до Gr = 2.1 104. Стрелками показаны переходы этого решения на несимметричное решение. Заметим, что при
увеличении Gr с 9-Ю4 до 1-Ю5 наблюдается не уменьшение, а некоторое увеличение доли твердого вещества в ячейке. В то же время в этом диапазоне изменения числа Gr незначительно уменьшается тепловой поток через ячейку. Такое поведение можно понять, рассмотрев характер течения, представленный на рис. 2.5.3 и 2.5.5. Как видно из этих графиков, вверху области существуют вторичные течения, которые развиваются при увеличении числа Грасгофа. Это приводит к потерям и уменьшению теплового потока, и, как следствие, увеличивается доля твердой фазы.
Перейдем к описанию решения второго типа, при котором в центре квадрата жидкость течет вверх. В этом случае стабильное решение существует от Gr = 1.6-104 до Gr = 12 -105. Область существования "нестабильного" решения второго типа - от Gr = 1.15-104 до Gr = 1.6-104. При Gr = 1.2-105 и Gr = 1.6-104 это симметричное решение тоже становится неустойчивым и "сваливается" на несимметричное решение. Аналогично предыдущему случаю перед переходом на несимметричное решение наблюдается незначительное увеличение доли твердого вещества и, соответственно, уменьшение теплового потока, начиная от Gr = 9.3-104 до критического значения для этого решения Gr = 1.2-105. Здесь так же, как и в предыдущей ситуации, существуют вторичные течения, находящиеся теперь внизу области (рис. 2.5.7 и 2.5.9). В указанном выше диапазоне изменения числа Грасгофа вторичные течения увеличиваются, что приводит к потерям и, соответственно, уменьшению теплового потока и увеличению доли твердой фазы.
Сравним результаты для решений первого и второго типа, представленные на рис. 2.5.1 и 2.5.2. При данном Gr решение первого типа дает меньшее значение как доли твердого вещества, так и теплового потока, переносимого через ячейку, чем значение для решения второго типа. Такой необычный, на первый взгляд, результат объясняется тем, что для решения второго типа поверхность раздела твердое-жидкое значительно более выпукла вверх, чем в случае решения первого типа, где граница раздела хотя и выпукла вниз, но является почти горизонтальной. Тем самым распределение плотности теплового потока на верхней стороне квадрата в случае решения первого типа ближе к постоянному значению (рис. 2.5.4), чем для решения второго типа (рис. 2.5.6). Кроме того, длина границы раздела фаз, изображенная штриховой линией на рис. 2.5.3, для решения второго типа превышает длину границы раздела фаз решения первого типа, а через большую поверхность возможно прохождение большего потока тепла.
Наибольшему диапазону изменения Gr отвечает несимметричное решение, которое получается при начальном распределении температуры = 1 - X и существует от Gr = 7.8 10 до Gr = 5.9-105. При числах Gr больше, чем Gr = 5.9-105, это решение становится неустойчивым и получить стационарное решение не удается. Если задать начальное условие 0 = X, то получим несимметричное решение, являющееся зеркальным отображением рассматриваемого решения относительно оси симметрии, проходящей через центр квадрата. Но поскольку по переносимым тепловым потокам, структуре вихревого движения жидкости, а также по количеству твердой фазы эти решения одинаковы, то рассматривалось только одно из этих решений.
Граничное условие в центре
Как видно из уравнений (3.2.12) и (3.2.13), в полярной системе координат появляются дополнительные члены в уравнениях движения жидкости. Для корректного счета необходима аппроксимация этих членов.
Выпишем отдельно члены, которые необходимо аппроксимировать (считая, что они находятся в правых частях уравнений).
На рис. 3.5.1 представлены контрольные объемы для скоростей. Поскольку в алгоритме SIMPLER используется «шахматная» сетка для скоростей, то, как видно из рис. 3.5.1, нам неизвестно значение скорости по радиусу в точке, где вычисляется скорость по углу, и наоборот. Для аппроксимации источниковых Контрольный объем Контрольный объем для скорости Уф для скорости Vr членов необходимо знать значение скоростей в соответствующих точках. Например, для скорости V значение скорости по углу нам известно, тогда как значение скорости Vr определено в других узлах сетки. Поэтому в программе производится вычисление компонент скоростей в соответствующих точках. Один из возможных вариантов - это нахождение значения скоростей по среднему арифметическому от значений скоростей в углах контрольного объема. = ,пьХ= \ (3.5.3) где Vr - значение скорости Vr в точке, где вычисляется Уф; Уф значение скорости Уф в точке Vr; Vrnb HV - значения соответствующих компонент скорости в углах контрольных объемов. Для нахождения производных используется формула центральной разности. Однако при использовании такой аппроксимации ошибка, получаемая в решении, достаточно велика. Она возникает из-за полуторного контрольного объема для скорости Vr в центре (рис. 3.4.1). Для улучшения аппроксимации в- центре проинтегрируем выражение для источникового члена по контрольному объему и разделим на его площадь.
Здесь Rq - радиус, проведенный к центру контрольного объема, RuARAq - площадь контрольного объема, ч дУ л значение XR-Эф), производной скорости Уф по дуге, вычисленное с помощью центрально-разностной аппроксимации. Необходимо учесть, что при равномерной сетке (наиболее часто используемый вариант) Ru будет равно радиусу, проведенному к вектору скорости, во всей области, кроме центрального контрольного объема для скорости Vr, где RU=0.5-AR. Ниже представлен блок программы, в котором производится вычисление значений для источниковых членов.
В качестве теста оптимальности данной аппроксимации источниковых членов рассматривалась задача течения жидкости через полость, поскольку существует возможность сравнения численных расчетов с известным аналитическим решением. Направление протекания потока со скоростью UQ Направление протекания потока со скоростью U0 Различные варианты тестовых задач На рис. 3.5.2 показаны примеры различной постановки подобной задачи для половины и для полного круга. При этом рассматривались следующие граничные условия: Уф=и0, Vr=0, приФ = 0 a) Уф =-U0, Vr =0, при ф = 7Т (3.5.5) Уф = -U0 coscp, Vr = U0 sin (p, при R = R0 Vr=U0, Уф=0, приср = 0 b) Vr = -U0, V9 = 0, при q = к (3.5.6) Уф = -U0 sin ф, Vr = U0 cosф, при R = R0 c) Уф = -U0 зіпф, Vr = U0 соБф, при R = R0 (3.5.7) В первых двух вариантах расчет производился на разных сетках: 32x32 и 62x32 (первое число - это количество контрольных объемов по углу), в варианте с полным кругом - на сетке 62x32. Решение задачи определялось при различных источниковых членах: 1) источниковые члены без учета коррекции радиуса в центре, 2) источниковые члены с коррекцией радиуса в центре, 3) аналитическое значение источниковых членов. Полученное в результате расчетов поле скоростей сравнивалось с известным аналитическим решением для каждой узловой точки расчетной области. Результаты расчетов представлены в таблице 1. Таблица Задача Сетка Аппроксимация Ошибка вычислений АУф/Уф % AVr/Vr % половинакруга,граничныеусловия(3.5.5) 32x32 обычная 1.0091 4.5270 с коррекцией в центре 0.0387 0.9655 аналитическая 0.0408 0.9590 62x32 обычная 1.0050 4.7801 с коррекцией в центре 0.0119 0.7093 аналитическая 0.0125 0.7070 половинакруга,граничныеусловия(3.5.6) 32x32 обычная 4.0501 7.4234 с коррекцией в центре 0.3466 0.1382 аналитическая 0.3426 0.1521 62x32 обычная 4.1598 7.7781 с коррекцией в центре 0.1722 0.0343 аналитическая 0.1717 0.0646 полный круг 62x32 обычная 4.2319 4.5723 с коррекцией в центре 0.1012 0.0367 аналитическая 0.0917 0.0534 Как видно из таблицы 1, предложенная аппроксимация с коррекцией радиуса в центре дает удовлетворительный результат во всех случаях. 3.6. Тестовые расчеты
Помимо сравнения решения, получаемого численно, с аналитическим решением, проведенным в предыдущем пункте, проводилось сравнение с результатами других авторов. В качестве теста был взят один из экспериментов, выполненных в работе К.С. Cheng, М. Takeuchi "Transient natural convection of water in a horizontal pipe with constant cooling rate through 4C" (1976) [31]. В работе численно исследовалась конвекция воды в горизонтальной трубе, охлаждаемой с постоянной скоростью (см. рис. 3.6.1). Процессы, происходившие в области, считались двумерными.
В начальный момент времени вода, находящаяся в области, имеет постоянную температуру, выше температуры инверсии плотности. Далее температура на боковой поверхности уменьшалась с постоянной скоростью по следующей формуле: TW=T0-Ht (3.6.1.) где Н - скорость охлаждения стенки. Обладая максимальной плотностью при температуре 4С, вода подвергается нормальному тепловому расширению выше этой температуры и претерпевает инверсию теплового расширения при температурах ниже 4С. Зависимость плотности воды от температуры, в интервале 0-20С, можно аппроксимировать соотношением в форме Гебхарта-Моллендорфа [15]: p = Pm(l-pTiq) (3.6.2) где рт = 999.9720 кг/м3, р = 9.297173-10"6 (С)Л Tj = 4.029325(С), q= 1.894816.
Тестирование программы
В качестве теста было взято несколько экспериментов, выполненных в работе (Stefan Schneider, Johannes Straub "Laminar natural convection in a cylindrical enclosure with different end temperatures" (1992)) [33]. В работе численно и экспериментально исследовалась конвекция жидкости и газа в цилиндрической полости, охлаждаемой с одной стороны и подогреваемой с другой. В качестве одного из параметров экспериментов брался угол наклона оси цилиндра к горизонту (см. рис. 4.3.1).
Запишем в безразмерном виде систему уравнений, описывающих процессы, происходящие в жидкости, предполагая, что: 1) жидкость несжимаема и ее течение ламинарное; 2) теплофизические свойства постоянны и равны для двух фаз; 3) изменения объема в системе незначительны; 4) для жидкой фазы справедливо приближение Буссинеска. dV VV і ЯРФ+(УУ)УФ+ Гф= +v2v-5т ф R Rdq ф \ 2—T + —T R2 R2 avr + sin((p)sin(a) Gr 0 aVr+(vv)vr Уф= ap+v2vr Vr-ax r R aR r R2 2 av,R2 dq (4.3.1)
На рис. 4.3.2-4.3.5 представлены результаты тестовых расчетов в сравнении с результатами работы [33]. На рис. 4.3.2 для разных чисел Релея представлена зависимость теплового потока через цилиндр в зависимости от угла наклона как для работы [33], так и для данной работы. В качестве исследуемого вещества рассматривался воздух. На рис. 4.3.3 также показана зависимость теплового потока через цилиндр от угла наклона для гелия. Здесь производится сравнение экспериментальных и численных результатов работы [33] с тестовыми расчетами. На рис. 4.3.4 представлена зависимость теплового потока через цилиндр от угла наклона для двух веществ (воздух и глицерин) при одном и том же числе Релея, равном 5000. Поле скоростей для числа Релея, равного 5000, в случае, когда отношение высоты цилиндра к его радиусу равно единице, для работы [33] и тестового расчета представлено на рис. 4.3.5.
Таким образом, результаты тестовых расчетов показывают хорошее совпадение с численными и экспериментальными данными работы [33].
Течение жидкости в цилиндре исследовалось при различных значениях отношения высоты цилиндра к его радиусу, в частности при L/Ro=2 и L/Ro=l. В работе сравниваются стационарные решения, получаемые в двумерной постановке, с решениями, получаемыми в трехмерном случае, а также определяется количество стационарных течений в трехмерной постановке.
Отношение высоты к радиусу L/RQ=2. При данном отношении высоты цилиндра к его радиусу и Gr=5-104 в двумерной осесимметричной постановке можно получить два вида решения: с течением жидкости у стенки цилиндра вверх и с течением жидкости у стенки цилиндра вниз (см. рис. 4.4.1). Причем количество тепла, проходящее через ячейку, для этих двух решений одинаково. Однако при рассмотрении такой задачи в полной трехмерной постановке решения, соответствующего данным двумерным течениям, получить не удалось. При данном числе Грасгофа в трехмерной постановке было получено только одно решение (рис. 4.4.2, 4.4.3). Данное решение несимметрично и имеет трехмерную структуру, а также значительно превышает двумерные решения по переносимому тепловому потоку (в двумерном случае Nu=2.878, а в трехмерном Nu=3.312). При числе Грасгофа, равном 1-Ю4, в двумерном случае возможно получение только тривиального решения без течения жидкости (то есть Nu=l), а в трехмерном случае также получается решение с несимметричным течением, при этом число Нуссельта равно 1.986.
Отношение высоты к радиусу L/Pypl. В данном случае течение жидкостей исследовалось только при одном значении числа Грасгофа, равном 1-Ю4. При таком значении числа Грасгофа и отношения высоты цилиндра к его радиусу в двумерной постановке также возможно получение двух различных решений с течением жидкости. Эти решения представлены на рис. 4.4.4. Как и в предыдущем случае, тепловой поток через цилиндр для этих двух решений одинаков и равен 2.093. При решении этой задачи в полной трехмерной постановке было получено пять различных решений, данные по которым содержатся в таблице 2. Первые два из них (рис. 4.4.5-4.4.8) соответствуют полученным двумерным течениям, течение представляет собой один вихрь в виде тора. По переносимым тепловым потокам данные решения также соответствуют полученным двумерным течениям. Решение, представленное нарис. 4.4.9 и 4.4.10 является несимметричным. Для этого решения число Нуссельта равно 2.048, что является наименьшим среди всех полученных решений. Этому решению также соответствует минимальное значение локального теплового потока, оно равно 0.38 и достигается у границы цилиндра (см. рис. 4.4.10). На рис. 4.4.11 и 4.4.12 представлено симметричное решение относительно плоскости проходящей через ось цилиндра, течение имеет два основных вихря и обладает трехмерной структурой. Решение, показанное на рис. 4.4.13 и 4.4.14, также симметрично. Оно подобно предыдущему, только вращение вихрей направлено в обратную сторону. Эти два решения имеют одинаковое значение теплового потока, проходящего через ячейку и равного 2.311. Из всех полученных решений (при данном числе Грасгофа) последние два дают максимальный тепловой поток через цилиндр.