Abstract

Abstract

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

Full Text

Исследование гидродинамических явлений вязкой среды, возникающих в условиях свободной тепловой конвекции, представляет интерес для понимания различных процессов получения материалов в черной и цветной металлургии. Ранее авторами были опубликованы результаты моделирования теплового состояния реактора для получения губчатого титана в бесконвективном приближении [1]. Однако, как показано в работе [2], в экспериментах на оптической установке, учет свободной тепловой конвекции в жидкой фазе, моделируемой прямоугольной областью, оказывает значительное влияние на неравновесную тепловую обстановку. В этой же работе, а также в монографии [3] предложены более сложные тесты с намораживанием воды в прямоугольной области с адиабатическими изотермическими границами в условиях инверсии ее плотности. Возникает задача разработки математической модели теплового состояния реактора, учитывающей конвективный тепломассоперенос, и проверки адекватности этой модели и метода ее численной реализации экспериментальным данным. Постановка задачи и реализация метода решения. В работе [2] экспериментально оптическими методами определялось неравновесное тепловое состояние в канале прямоугольного сечения с непроницаемыми стенками, полностью заполненном водой. Размеры канала: ширина – 0,045 м, высота –0,073 м. Верхняя и нижняя граница принимались адиабатными, правая и левая стенки – изотермическими с температурами соответственно 16,2 и 14,1 °С. Свойства воды в заданном интервале температур описываются температурными зависимостями, представленными в табл. 1. Коэффициент кинематической вязкости задается постоянным: ν = 1,429·10–6 м2/с. Таблица 1 Соотношения для описания свойств воды Наименование параметраРазрешающее соотношениеЕдиницы измерения Коэффициент теплопроводностиl = 0,569 + 0,0016·tВт/(м·К) Теплоемкостьc = 4,212 – 0,0011·tкДж/(кг·К) Плотностьr = r0(1 + 6,3·10–3·t – 7,491·10–4·t), где r0 – плотность воды при 0 °С, r0 = = 999,841 кг/м3кг/м3 Коэффициент объемного расширенияb = (0,15·t – 0,63)·10–4К–1 В перечисленных уравнениях значения температуры t подставляются в градусах Цельсия. Расчетная схема к модели тепломассопереноса в канале прямоугольного сечения приведена на рис. 1. t(1, y) = t1 t(N+1, y) = t2 y x i = 1 i = N+1 j = 1 j = M+1 u v ψ(1, y) = 0 ψ(x,M+1) = 0 ψ(N+1, y) = 0 ψ(x, 1) = 0 Рис. 1. Расчетная схема В обозначениях работ [3, 4] формулировка задачи тепловой конвекции в прямоугольных координатах (x, y) сводится к системе трех дифференциальных уравнений: переноса тепловой энергии, завихренности и Пуассона, которые имеют вид (1) где t – температура; w – завихренность; ψ – функция тока; τ – время; а – коэффициент температуропроводности, м2/с, а = λ0/(ρс); u, ν – компоненты скорости в проекциях на оси x и y соответственно, м/с; g – ускорение свободного падения, g = 9,81 м/с2. Система уравнений (1) решалась численным методом на регулярной сетке с числом разбиений N и M в направлении координат соответственно x и y при задании следующих начальных и граничных условий: t (τ = 0) = t0, w(τ = 0) = ψ(τ = 0) = 0, t(1, y) = t1 = 16,2 °С, t(N+1, y) = t2 = 14,1 °С, ψ(1, y) = ψ(x,M+1) = ψ(N+1, y) = ψ(x, 1) = 0, где t0 – температура в начальный момент времени. Перечислим основные методы, применявшиеся при разработке математической модели и программы ее реализации. Уравнения переноса тепловой энергии и завихренности решались методом продольно-поперечной прогонки, уравнение Пуассона – методом последовательной линейной верхней релаксации. В дискретных аналогах уравнений для расчета приграничных значений завихренности использовалась добавка к подъемной силе. При аппроксимации конвективных членов в уравнениях переноса использована схема с разностями против потока, известная еще как схема с донорными ячейками или схема Роуча. На непроницаемых границах выполнялось подправление условий прилипания с применением двухконтурной методики В.Л. Грязнова и В.И. Полежаева. В ходе счета шаг по временной сетке выбирался исходя из условия Куранта: ht < hx/vmax, где hx – шаг пространственной сетки, м; vmax – максимальное значение, выбранное из массива компонент скорости. Обсуждение результатов вычислительного эксперимента. Разработана математическая модель и компьютерная программа, позволяющие анализировать неравновесные теплофизические процессы, возникающие в условиях свободной конвекции вязкой среды без изменения ее агрегатного состояния. В вычислительном эксперименте задавались граничные условия и свойства воды. Счет проводился до выхода процесса на установление с наперед заданной точностью, решение тестировалось со сгущением сетки оценкой стремления численного решения к точному. Вычислительный эксперимент проводился на сетках, перечисленных в табл. 2. Таблица 2 Расчётные сетки, используемые для тестирования математической модели N203040506070 M*3249658197114 Число узлов (N+1)´(M+1)69315502706418259788165 * Число разбиений М по вертикальной оси выбиралось исходя из условия, что ячейка сетки имеет форму квадрата. На рис. 2 представлены результаты сравнения расчетных и экспериментальных данных. В условиях свободной конвекции вода течет по замкнутому контуру вдоль стенок канала (рис. 2, а). Более нагретые струи воды выносятся в верхнюю часть камеры, откуда, постепенно охлаждаясь, опускаются вдоль правой стенки. Вблизи границ наблюдается значительная плотность линий тока. В районе левого верхнего и правого нижнего углов канала отчетливо видны характерные изогнутости течения воды. Линии тока математической модели (рис. 2, б), повторяют особенности течения реальной жидкости: несимметричность течения в угловых зонах, два вихря в средней части расчетной области, пограничные слои у стенок канала. На рис. 2, в показаны расчетные изотермы, в результате конвекции происходит значительное их искривление. Вблизи вертикальных границ образуются перегибы, повторяющие направление течения воды. В средней части расчетной области изотермы располагаются практически горизонтально со стратификацией температуры по высоте. Счет на сетке 70´114. На основании сравнения расчетных и реальных линий тока воды можно сделать вывод об адекватности математической модели, которая позволяет детально описать наблюдаемые неравновесные процессы конвективного теплообмена в расчетной области прямоугольного сечения. На рис. 3 приведен график установления процесса счета на сетке 70´114, где отображена динамика изменения максимальных значений функции тока (ψmax), скорости (νmax) и относительной погрешности по завихренности (εω). а б в Рис. 2. Тепловая конвекция жидкости в прямоугольной области: а – экспериментальные траектории движения воды; б – расчетные траектории движения воды, ψ·105, м2/с; в – расчетные изотермы, t, °C Рис. 3. График установления итерационного процесса на сетке с числом разбиений 70´114: ψmax; νmax; εω На первых временных слоях процесс установления характеризуется скачкообразным изменением параметров, величина εw принимает значения более 100 %. Количество внутренних итераций по функции тока при этом более 200. По мере того как решение выходит на установление, осцилляции по завихренности затухают и после 12-го слоя по времени εw не превышает 1 %. Число внутренних итераций после выхода на установление также снижается до 4–5 на каждом слое по времени. Максимальное число внешних итераций наблюдается в начальный период счета, на первых трех слоях по времени. Дальнейшее изменение величин ψmax и νmax стремится к минимуму, небольшие колебания наблюдаются только по εw. Шаг сетки по времени следует выбирать исходя из условия Куранта с учетом коэффициента запаса 0,5. Для сетки 70´114 величина ht находилась в диапазоне 0,25–0,31 с, при этом достигается удовлетворительная сходимость решения и устойчивость счета. На более крупных сетках коэффициент запаса можно увеличивать до 0,7–0,8. Для оценки погрешности численного решения проводились сравнительные расчеты на 6 различных сетках. Оценивалось стремление численного решения к точному при сгущении сетки. После выхода решения на установление максимальные значения, полученные на мелкой сетке, по завихренности, скорости, функции тока и температуре фиксировались, а затем сравнивались с их аналогами на более крупной сетке. При этом рассчитывалась величина погрешности: где Sn1 и Sn2 – максимальное значение функции (w, ψ, v, t) на сетках с числом узлов n1 < n2. На рис. 4 представлен график изменения величины погрешности при сгущении сетки. Наблюдается устойчивая тенденция снижения погрешности по всем искомым параметрам. Представленный график характеризует стремление полученных математическим моделированием численных решений к точному при постепенном сгущении сетки. На максимально подробной сетке с числом разбиений 70´114 (по сравнению с сеткой 60´97) достигнута величина погрешности по завихренности 3,9 %. По остальным переменным R составляет менее 2 % и в частности для максимальной скорости равна 1,6 %. [(N+1)·(M+1)] ·102 Рис. 4. Динамика уменьшения величины погрешности R со сгущением сетки: vmax; ψmax; εω, εt Таким образом, показана адекватность разработанной математической модели и метода ее реализации для описания неравновесного процесса тепломассопереноса вязкой жидкости. Параметрические расчеты для прямоугольной области показали уменьшение погрешности расчетов при сгущении регулярной сетки. Погрешность расчетов на сетке с числом разбиений N´M = = 70´114 не превышает 3,9 %. Построены графики сходимости итерационного процесса и стремления численного решения к точному. Разработанная математическая модель оказывается полезной для описания процессов тепломассопереноса в аппарате для получения губчатого титана.