CC BY
60
МЕХАНИКА
УДК 532.5:536.25
О ЛОКАЛЬНЫХ ЭФФЕКТАХ СЛАБЫХ ТЕРМОГРАВИТАЦИОННЫХ КОНВЕКТИВНЫХ ТЕЧЕНИЙ
И. А. Ермолаев, С. В. Отпущенников
Саратовский государственный университет Email: pers@info.sgu.ru
Методами численного моделирования исследуются особенности естественной термогравитационной конвекции малой интенсивности, возникающей в условиях микроускорений. Изучается влияние тепловых граничных условий на локальные особенности температурных полей. Показано, что величина максимума температурного расслоения не монотонно зависит от интенсивности теплоотдачи на границе области. Предложен алгоритм коррекции граничных условий для вихря скорости на твердых непроницаемых стенках.
Ключевые слова: численное моделирование, метод конечных элементов, естественная термогравитационная конвекция.
Local Effects of the Weak Thermogravitational Convective Flows I. A. Ermolaev, S. V. Otpoushchennikov
The features of the natural low intensity thermo-gravitational convection occurring in micro-acceleration condition have been investigated numerically. The effect of thermal boundary conditions on the local characteristics of temperature field has also been studied. It was shown that the value of maximum temperature stratification depends monotonically on the intensity of heat transfer at the system boundaries. The correction algorithm has been proposed for the vorticity boundary conditions on the solid impermeable walls.
Keywords: numerical simulation, finite element method, natural thermogravitational convection. ВВЕДЕНИЕ
При увеличении интенсивности естественная конвекция в неравномерно нагретой области претерпевает ряд качественных переходов между различными режимами теплообмена. Наиболее известен переход к турбулентной конвекции. Между тем, переход от режима теплопроводности к развитой стационарной конвекции также обладает рядом важных особенностей, связанных с локальными характеристиками теплообмена и температурных полей. Именно в этой области параметров возникает эффект максимума температурного расслоения [1-3], заключающийся в образовании локальных зон перегрева и переохлаждения, обусловленных неэффективностью конвективного перемешивания.
Эффекты слабой конвекции актуальны, прежде всего, в условиях микроускорений, характерных для орбитальных космических аппаратов [4,5]. Исследования этих эффектов важны при решении проблемы получения качественных полупроводниковых материалов в космических условиях [6], анализе гравитационно-чувствительных систем и процессов [7-9], разработке систем термостабилизации радиоэлектронной и иной термочувствительной аппаратуры, работающей в условиях ослабленной силы тяжести и др.
Работы, посвященные численному исследованию локальных особенностей слабой термоконвекции, представлены в обзоре [5]. Кроме того, в [10] выявлены зависимости величины локального перегрева (переохлаждения) от интенсивности конвективного течения и отношения длин сторон прямоугольной полости для одно-, двух- и трехвихревых стационарных движений, даны оценки границ режимов конвекции. В [11] получены зависимости локальных эффектов от свойств жидкости, даны оценки влияния числа Прандтля на границы режимов конвекции.
На особенности протекания конвекции существенно влияют тепловые граничные условия. При моделировании эти условия постулируют тепловой поток и поле температуры на границах области, идеализируя тем самым теплообмен в стенках и во внешней среде. В настоящей работе изучается влияние граничных условий для температуры на локальные эффекты слабой конвекции.
1. ПОСТАНОВКА ЗАДАЧИ, МОДЕЛЬ И МЕТОД РЕШЕНИЯ
Рассматривалось плоское конвективное течение вязкой, термически сжимаемой жидкости, для которой справедливо приближение Буссинеска. Жидкость заполняла двухмерную замкнутую прямоугольную область шириной L, высотой H с твердыми и непроницаемыми стенками. Для описания течения и теплообмена использовалась декартова система координат, начало которой совпадало с левым нижним углом области. Ось x была направлена горизонтально, ось y — вертикально, сила тяжести направлена вертикально вниз. В начальный момент времени поле температур считалось однородным, жидкость находилась в гидростатическом равновесии.
При решении задачи использовались нестационарные двухмерные уравнения Буссинеска [12]. В качестве масштабов расстояния, времени, скорости и температуры были выбраны H, H2/v, v/H, q0H/A. Безразмерные переменные были равны соответственно X = x/H, Y = y/H, т = vt/H2, U = uH/v, V = vH/v, 9 = Aa/(q0H). Здесь x, y — координаты, т — время, v — коэффициент кинематической вязкости, u, и — составляющие скорости в проекции на оси x, y соответственно, и = T — T0, T0 — начальная температура жидкости, A — коэффициент теплопроводности, q0 — плотность теплового потока. Безразмерные уравнения Буссинеска в переменных «вихрь скорости -функция тока - температура» записывались следующим образом:
5ш 5ф 5ш 5ф 5ш . 59 ,Л.
57 + 5Y5X — 5X5Y = Аш — Gry 5X' (1)
Аф = ш, (2)
59 + 5±5L — 55±5L = ±А9 (3)
5т + 5Y 5X 5X 5Y Pr ' ( )
где ш, ф — вихрь скорости, функция тока соответственно, Gry = gyßq00— — число Грасгофа
Av2
(Grx = 0), Pr = v/x — число Прандтля, gy — составляющая ускорения силы тяжести в проекции на ось y (gx =0), ß — температурный коэффициент объемного расширения, x — коэффициент температуропроводности.
Граничные условия для системы (1)-(3) имели следующий вид. На боковых границах во всех расчетах были заданы условия «прилипания» для функции тока и условия «адиабатической изоляции» для температуры:
X =0: ф(0, Y, т) = Щ^П = , X = L : *(L,Y,r) = 5Ф<^ = ,
на горизонтальных границах — условия «прилипания» для функции тока:
Y = X : *(X, 0,т) = .
Граничные условия для температуры варьировались
Y = H : ф^.т) = .
В начальный момент времени ш(Х,У, 0) = ф(Х,У, 0) = 6(Х,У, 0) = 0. Задача решалась методом конечных элементов (МКЭ) Галеркина [13]. Температура, вихрь скорости и функция тока аппроксимировались линейной комбинацией не зависящих от времени функций формы на линейных треугольных конечных элементах. Для временной аппроксимации использовалась неявная двухслойная схема.
Применение процедуры МКЭ Галеркина для решения системы (1)-(3) (ортогонализация невязки относительно базисных функций, определенных локально — на конечных элементах) приводит к системе из трех матричных уравнений вида
[С ](Ф) + [К ]{Ф} + } = 0,
где [С] — матрица демпфирования, [К] — матрица жесткости, {Ф} — вектор узловых значений, (^} — вектор нагрузки. Элементы матриц [С], [К] и вектора (^}: для вихря скорости:
с, = ¡8 NN ¿Б,
для функции тока:
=
ду ) е г дх \дх ) е г ду дх дх
дЫ дЫ,
¡г = -
ч а е - ^ /,
Ыг ^е
дУ
П-1 ¿Б;
ду ду
¿Б,
е„ = 0,
к.. = [ (дЫ±дЫ± + дЫгдЫ,
г] /де \ дх дх ду ду
¿Б,
для температуры:
егу = ^ ЫгЫ ¿Б,
АТ
к • =
кг] =
N. Ы, + | Ш Ыг ^
у е х
¡г = / Ыг^е ¿Б;
£) Ы, Б
х е у
¡г =
1
аТ
ЫгОе
п — 1
¿Б.
Здесь ег, — элемент матрицы демпфирования, Ат — шаг по времени, Ыг, Ы, — функции формы, Бе — площадь конечного элемента, е — индекс конечного элемента, к, — элемент матрицы жесткости, ¡г — элемент вектора нагрузки, п - индекс шага по времени. Величины с нижним индексом е усредняются по каждому конечному элементу, величины с верхним индексом п — 1 относятся к предыдущему временному шагу. Матрицы [С] и [К], а также вектор (^} формируются суммированием по всем конечным элементам. При решении системы алгебраических уравнений применялся метод Гаусса для ленточных матриц.
Расчеты проводились по конечно-элементной программе, реализующей данный алгоритм [14]. Стационарные решения были получены методом установления. Критерием установления являлось нера-
венство
где 9Г
Фг
—вкт I + \-г+1—|+\Фкг+1—Фкт | < ерв,
экстремальные значения температуры, вихря скорости и функции тока. Индекс
к — номер шага по времени, величина ерв изменялась в интервале 10 5-10 6. Шаг по времени 10 Расчеты проводились на равномерных сетках.
1-3
2. ГРАНИЧНЫЕ УСЛОВИЯ ДЛЯ ВИХРЯ СКОРОСТИ
При решении уравнений гидрогазодинамики в переменных «вихрь скорости - функция тока» граничные условия для вихря скорости на твердых непроницаемых границах не могут быть заданы в явном виде и требуют приближенного вычисления тем или иным способом. Чаще всего используются аппроксимационные формулы разного порядка точности, полученные путем разложения в ряд функции тока вблизи границы [15].
Однако использование приближенных граничных условий данного вида приводит к снижению устойчивости, ограничениям на временной шаг и замедлением сходимости. Кроме того, эффективность той или иной аппроксимационной формулы может меняться в зависимости от особенностей
— V
в
в
в
в
конкретной задачи. Одним из способов повышения устойчивости является метод нижней релаксации [16]. Однако введение релаксационного параметра, строго говоря, возможно лишь для стационарных задач. Для нестационарного режима использование релаксации приводит к дополнительной невязке в граничных условиях «прилипания». Для ее устранения обычно вводится внутренний итерационный цикл на каждом временном шаге.
Для получения приближенных граничных значений вихря скорости может производиться также экстраполяция решения за границу области [17]. Другим распространенным подходом является использование различных итерационных процедур [18], сведение решения «линеаризованной» разностной схемы для основной задачи к решению задачи Дирихле с некоторыми интегральными граничными условиями [19] и др.
В [20] показано, что локального эквивалента граничного условия для ш, вообще говоря, не существует, а должно выполняться некоторое интегральное условие вида
J шпйБ = У (^п - ад^ ¿Ь, (4)
где п — функция, определенная на всей области Б, такая, что Дп = 0, фь = а, (дф/дп)ь = Ь. В задачах конвекции для замкнутых областей, где фь = (дф/дп)ь = 0, выражение (4) упрощается
[ шпЛБ = 0. (5)
.'б
Отметим, что аналогичное условие можно получить, воспользовавшись теоремой Остроградского-Гаусса для уравнения (2)
[ ДфйБ = [ шйБ = [ ^¿Ь. (6)
'5 ¿Б ¿Ь
дп
Выполнение условия «прилипания» приводит к условию, налагаемому на поле ш вида
дф
JS шЛБ = /ь — ¿Ь = 0, что эквивалентно условию (5) при п = 1.
Условие (6) средствами МКЭ может быть использовано для корректировки поля вихря скорости на каждом временном шаге или каждой итерации. Пусть на некотором временном шаге или итерации условие (6) не выполняется:
[ шЛБ = ( ^¿Ь = е, (7)
где е — «невязка поля ш». Скорректируем поле ш так, чтобы (ш + ш') ¿Б = 0, где ш' — поправка. Нетрудно заметить, что ш' ¿Б = -е. Допуская, что величина е равномерно распределена по области Б, можно записать ш' = -е/Б, где е определяется выражением (7). При использовании МКЭ
Е
этот интеграл легко вычисляется как е = ш ¿Б = ^ шеБе, где Е — число конечных элементов
е=1
в области Б, е — индекс элемента, Б — площадь элемента. Таким образом, поправка к полю вихря
скорости на каждом временном шаге или итерации имеет вид
1 Е
ш' = - Б^2шеБе, (8)
Б е=1
Уравнения (1)-(3) решались последовательно, каждый временной шаг начинался с вычисления поля температуры (3), затем определялись граничные условия для вихря скорости по той или иной аппроксимационной формуле, например формуле Вудса [15], и решалось уравнение (1). Далее поле вихря скорости корректировалось согласно (8) и определялось поле функции тока (2).
Алгоритм тестировался на классической задаче о тепловой конвекции в квадратной области, подогреваемой сбоку [21]. Сравнение результатов с МКР и МКЭ по интегральному потоку тепла показало наибольшие отличия в 2% от МКЭ и 5.6% от МКР для От = 103 и 4% от МКЭ и 4.5% от МКР для От = 104.
3. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ
Численные расчеты проводились для жидкости, характеризуемой числом Прандтля Рг = 1, что ближе всего соответствует свойствам воздуха. Исследовались одновихревые и двухвихревые течения в двумерных прямоугольных полостях с различным относительным размером Ь/И.
Теплоперенос при числах Грасгофа менее 103 можно назвать теплопроводностным [4]. Конвективное течение в этом случае почти не влияет на теплообмен в полости, температурные поля практически соответствуют режиму теплопроводности. Изотермы параллельны горизонтальным стенкам, небольшие искажения заметны лишь вблизи нагреваемых, либо охлаждаемых границ. С ростом интенсивности выталкивающих сил усиливается температурное расслоение, растет горизонтальная составляющая теплового потока, что приводит к повышению либо к понижению температуры в некоторой локальной области, в сравнении с режимом теплопроводности. Образуются зоны «перегрева» и «переохлаждения» при росте среднего теплового потока через слой.
Развитие теплообмена с увеличением интенсивности выталкивающих сил для одновихревого течения при Ь/И = 1.5 показано на рис. 1. Относительный размер Ь/И = 1.5 соответствует максимальному температурному расслоению для прямоугольной области [10]. Здесь на нижней границе была задана теплового потока дО(Х, 0, т)/дУ = 1, на верхней границе — постоянная температура 0(Х, И, т) = 0. При таких граничных условиях существует ненулевая составляющая градиента температуры вдоль нижней границы (рис. 1, б), что дополнительно обеспечивает движение.
До Сг = 1.25 • 103 (рис. 1, кривая 2) распределение температуры по нагреваемой нижней границе линейно, в нижнем углу возникает зона локального перегрева, составляющая примерно 1/3 ширины полости.
С увеличением интенсивности конвекции усиливается горизонтальное температурное расслоение, длина зоны перегрева уменьшается, максимум температуры растет, достигая наибольшего значения при Сг « 2 • 103 (рис. 1, кривая 4). С дальнейшим ростом числа Грасго-фа уменьшаются как размер зоны перегрева, так и величина максимума температуры (рис. 1, кривые 4-6). Можно отметить, что при Сг ~ 2.5 • 103 и более температура всех точек границы изменяется с ростом критерия Сг по одному закону (линейно). При Сг ~ 3• 103 величина перегрева становится равной нулю, кривая 6 соответствует режиму развитой стационарной конвекции.
Эффекту локального перегрева в полости, подогреваемой снизу, соответствует эффект локального переохлаждения в полости при оттоке тепла сверху. На рис. 2 показаны распределения температуры при двухвихревом конвективном течение в полости с Ь/И = 2.5 при оттоке тепла с верхней границе 39(Х,И,т)/ЗУ = —1 и изотермической нижней стенкой 9(Х, 0,т) = 0.
Здесь также существует ненулевой температурный градиент вдоль верхней границы. Аналогично результатам на рис. 1 с увеличением числа Сг следует рост величины и уменьшение размеров участков переохлаждения (рис. 2, кривые 1-3) в углах верхней границы. Затем следует снижение как величины, так и длин участков переохлаждения (рис. 2, кривые 3-5). При Сг ~ 2.25 • 103 температура всех точек верхней границы изменяется с ростом Сг линейно, кривая 6 соответствует режиму развитой конвекции.
Рис. 1. Развитие теплообмена: а — распределения температуры по нагреваемой стенке: 1 — режим теплопроводности, 2-6 — Сг = 1.25 • 103, 1.5 • 103, 2.0 • 103, 2.5 • 103, 3.0 • 103; б — изотермы поля температур и в линии тока для Сг = 2.5 • 103
Наиболее интересными представляются результаты исследования влияния более общих тепловых граничных условий (условий Ньютона). Эти условия соответствуют некоторому произвольному термосопротивлению границ. Так, на верхней стенке было задано общее условие теплоотдачи дО(Х,Н,т)/дУ = ЬО, где Ь = 0.1, 1, 10, 100 — безразмерный коэффициент, на нижней границе — плотность теплового потока дО(X, 0, т)/дУ = 1, величина Ь/Н = 1. При Ь = 100 верхнюю стенку можно считать изотермической, что соответствует бесконечной теплопроводности границ. Случай Ь = 0.1 соответствует практически теплоизолированной границе. Течение одновихревое.
Здесь также образуется зона перегрева (на нагреваемой стенке) и переохлаждения (на охлаждаемой стенке), температуры которых соответственно выше и ниже экстремальных температур режима теплопроводности. Однако в этой задаче зона переохлаждения сохраняется во всем исследуемом диапазоне чисел Грасгофа (рис. 3, а). Следует отметить также, что величина переохлаждения, определяемая как 0т/0т,о, где 0т — минимум температуры при текущем значении критерия От, 0т,о — минимум температуры при От = 0, заметно больше величин перегрева на рис. 1 и рис. 2.
Изменения параметра Ь, имеющего физический смысл величины, обратной безразмерному термосопротивлению границы, существенно влияют на градиент температуры вблизи верхней стенки. Градиент температуры на вертикальных границах с ростом числа Ь также увеличивается, температурное расслоение по горизонтали более выражено при низких значениях этого коэффициента.
Зависимость 0т/0т,о от параметра Ь немонотонна и имеет экстремум при Ь =10 (рис. 3, б). Это значение обеспечивает относительный минимум температуры 0т/0т,о = 0.585. Данное сочетание тепловых граничных условий моделирует, фактически, многослойную систему, связанную тепловым взаимодействием, где Ь — параметр взаимодействия. Очевидно, это взаимодействие нелинейно, чем и объясняется существование экстремума на рис. 3, б.
Рис. 2. Развитие теплообмена: а — распределения температуры по охлаждаемой верхней: 1 — режим теплопроводности, 2-6 — От = 1.5 • 103, 1.75 • 103, 2.0 • 103, 2.25 • 103, 2.5 • 103, Ь/И = 2.5; б — изотермы поля температур, в — линии тока для От = 2.5 • 103
Рис. 3. Развитие теплообмена: а — изменения относительного минимума температуры с ростом числа От при Ь = 0.1,1,10,100; б — зависимость величины переохлаждения от коэффициента теплоотдачи при От = 104
ЗАКЛЮЧЕНИЕ
Изменения величины относительного температурного расслоения с ростом интенсивности конвективного течения качественно различны для границ с заданной температурой, тепловым потоком (оттоком) и общим условием теплоотдачи. На границе с заданным тепловым потоком (оттоком) температурное расслоение проходит через максимум, затем снижается до нуля, далее изменяясь линейно. На стенке с общим условием теплоотдачи зона переохлаждения (перегрева) сохраняется во всем исследуемом диапазоне чисел Грасгофа. Зависимость 0т/0т,о от параметра Ь не монотонна и имеет экстремум при Ь = 10.
Библиографический список
1. Полежаев В. И. Нестационарная ламинарная тепловая конвекция в замкнутой области при заданном потоке тепла // Изв. АН СССР. Механика жидкости и газа. 1970. № 4. С. 109-117.
2. Авдуевский В. С., Полежаев В. И. Некоторые особенности естественной конвекции жидкостей и газов // Избранные проблемы прикладной механики. М. : ВИНИТИ, 1974. С. 11-20.
3. Полежаев В. И. Эффект максимума температурного расслоения и его приложения // Докл. АН СССР. 1974. Т. 218, № 4. С. 783-786.
4. Полежаев В. И., Белло М. С., Верезуб Н. А., Дубовик К. Г., Лебедев А. П., Никитин С. А., Павловский Д. С., Федюшкин А. И. Конвективные процессы в невесомости. М. : Наука, 1991. 240 с.
5. Полежаев В. И. Конвекция и процессы тепло- и мас-сообмена в условиях космического полета // Изв. РАН. Механика жидкости и газа. 2006. № 5. С. 67-88.
6. Zemskov V. S., Raukhman M. R., Shalimov V. P., Goncharov V. A. Peculiarities of inhomogeneities and heat/mass transfer during directional crystallization under low and normal gravity conditions // Single crystal growth and heat & mass transfer : proceedings of First Intern. Conf. Obninsk, 2003. Vol. 2. P. 717-726.
7. Полежаев В. И. Режимы микроускорений, гравитационная чувствительность и методы анализа технологических экспериментов в условиях невесомости // Изв. РАН. Механика жидкости и газа. 1994. № 5. С. 22-36.
8. Асеведо Х., Ермаков М. К., Зыков С. Г., Комаров М. М., Либерман Е., Никитин С. А., Полежаев В. И., Рябуха С. Б., Сазонов В. В., Стажков В. М. Микроускорения на орбитальной станции «Мир» и оперативный анализ гравитационной чувствительности конвективных процессов тепло-массопереноса // Космические исследования. 1999. Т. 37, № 1. С. 86-101.
9. Сазонов В. В., Юферев В. С. Тепловая конвекция, вызванная квазистатической компонентой поля микроускорений орбитальной станции «Мир» // Изв. РАН. Механика жидкости и газа. 2000. № 3. С. 39-45.
10. Ермолаев И. А., Жбанов А. И., Отпущенников С. В. Исследование режимов малоинтенсивной конвекции в прямоугольной полости с тепловым потоком
на границе // Изв. РАН. Механика жидкости и газа. 2008. № 3. С. 3-11.
11. Ермолаев И. А., Жбанов А. И., Кошелев В. С., Отпущенников С. В. Исследование влияния числа Прандтля на локальные свойства малоинтенсивной конвекции в подогреваемой снизу прямоугольной области // Теплофизика высоких температур. 2011. Т. 49, № 4. С. 589-593.
12. Гершуни Г. З, Жуховицкий Е. М. Конвективная устойчивость несжимаемой жидкости. М. : Наука, 1972. 392 с.
13. Зенкевич О., Морган К. Конечные элементы и аппроксимация. М. : Мир, 1986. 288 с.
14. Ермолаев И. А., Жбанов А. И., Кошелев В. С. Решение двумерной нестационарной задачи тепло- и массо-переноса методом конечных элементов // Вопросы прикладной физики. Саратов : Изд-во Сарат. ун-та, 2002. Вып. 8. С. 60.
15. Тарунин Е. Л. Вычислительный эксперимент в задачах свободной конвекции. Иркутск : Изд-во Иркутск. ун-та, 1990. 225 с.
16. Тарунин Е. Л. О выборе аппраксимационной формулы для вихря скорости на твердой границе при решении задач динамики вязкой жидкости // Численные методы механики сплошной среды. 1978. Т. 9, № 7. С. 97-111.
17. Полежаев В. И., Грязнов В. Л. Метод расчета граничных условий для уравнений Навье-Стокса в переменных «вихрь, функция тока» // Докл. АН СССР. 1974. Т. 219, № 2. С. 301-304.
18. Вабищевич П. Н. Разностные схемы для задач гидродинамики в переменных «функция тока - вихрь скорости» // Докл. АН. 1996. Т. 346, № 4. С. 442-444.
19. Бабенко К. И., Введенский Н. Д. О численном решении краевой задачи для уравнений Навье-Стокса // Журн. вычисл. мат. и мат. физ. 1972. Т. 12, № 5. С. 1343-1349.
20. Флетчер К. Вычислительные методы в динамике жидкости : в 2 т. Т. 1 / пер. с англ. М. : Мир, 1991. 504 с.
21. Полежаев В. И., Бунэ А. В., Верезуб Н. А., Глуш-ко Г. С., Грязное B. Л., Дубовик К. Г., Никитин C. А., Простомолотов А. И., Федосеев А. И., Черкасов С. Г. Математическое моделирование конвективного тепло-и массообмена на основе уравнений Навье-Стокса. М. : Наука, 1987. 271 с.
