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

УДК 532, 516.5

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЕСТЕСТВЕННОЙ КОНВЕКЦИИ В ВЕРТИКАЛЬНОМ КАНАЛЕ С СИСТЕМОЙ НАГРЕВАТЕЛЕЙ

© 2010 г. Е.И. Калинин, А.Б. Мазо

Казанский (приволжский) федеральный университет, Kazan (Volga Region) Federal University,

ул. Кремлевская, 18, г. Казань, Р. Татарстан, 420008, Kremlyovskaya St., 18, Kazan, Republic ofTatarstan, 420008, public.mail@ksu. ru public. mail@ksu.ru

Предложен численный метод расчета нестационарной естественной термоконвекции жидкости в вертикальном канале, основанный на решении системы уравнений Навье-Стокса в переменных функция тока - завихренность. Метод использован для исследования развития конвективного течения около системы из трех профилированных нагревателей, размещенных в канале.

Ключевые слова: естественная конвекция, численное моделирование, уравнения Навье-Стокса, переменные функция тока - завихренность, вертикальный канал, система нагревателей.

A numerical method for simulating non-stationary natural convection in vertical channel, based on numerical solving of Navier—Stokes equations in stream function — vorticity formulation is developed. The method is used to investigate the evolution of connvective flow around a system of three profiled heaters, located in a channel.

Keywords: natural convection, numerical simulation, Navier-Stokes equations, stream function - vorticity variables, vertical channel, system of heaters.

Термоконвективные течения жидкости в канале широко распространены в теплотехнике и оказывают значительное влияние на интенсивность теплообмена. В большинстве теоретических и численных работ по естественной конвекции рассматриваются горизонтальные каналы с различными режимами нагрева и охлаждения стенок [1—3]. Немногочисленные работы по течениям в вертикальном канале посвящены исследованию стационарной конвекции около единственного нагревателя в форме круга [4-6], квадрата [7, 8] или пластины [9]. Вместе с тем нестационарные конвективные течения около нескольких нагревателей остаются изученными недостаточно. Отличительной особенностью свободной конвекции в вертикальном канале является то, что расход жидкости заранее не известен и меняется со временем. В ряде случаев это приводит к смене режимов обтекания нагревателей, отрыву потока, образованию вихревых дорожек за телами и формированию сложных нестационарных течений в дальнем следе. Математическое моделирование таких течений предполагает разработку адекватных методов расчета термогидродинамических полей.

В настоящей статье предлагается метод численного решения задачи Навье-Стокса в переменных ш

(функция тока), с (завихренность), T (температура). Для определения граничных значений функции тока на обтекаемых поверхностях применяются интегральные условия Пирсона [10, 11]. Построенная численная схема используется для моделирования естественной термоконвекции в плоском вертикальном канале около системы 3 нагревателей, имеющих форму крылового профиля NACA 0040.

Постановка задачи

Рассматривается конвективное течение жидкости в канале D = [0, L] х [0, H] с помещенными в него N

нагревателями с границами у, i = 1, N (рис. 1). Во входном сечении yin безразмерная температура жидкости равна нулю, а поверхности нагревателей поддерживаются при постоянной температуре T = 1; верхняя у0 и нижняя yN+l стенки канала теплоизолированы. Термоконвективное течение развивается из состояния покоя, в котором начальная температура жидкости равна нулю.

жении Буссинеска в переменных у,ю,Т имеет вид [6]

Рис. 1. Схема течения в вертикальном канале с нагревателями

В декартовых координатах х, у (ось х направлена против вектора ускорения гравитации g) определяющая система безразмерных уравнений Навье-Стокса и конвективной теплопроводности в прибли-

V7 1 Л

--ьv V® =—Аю + F, F =

3t Re

3Т 3y 9

-Ау = ю, (1)

3v

ю =--

u =

v =

дТ , Y7T 1 AT

— + v -VT = — AT 3t Pe

du

dx ду J ду dx

где t - время; v = (u, v) - вектор скорости жидкости; Re и Pe - числа Рейнольдса и Пекле.

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

формулами Gr = Re2, Pr = Pe/Re.

Граничные условия во входном сечении канала Л,y еуш : T = 0, ш = Pw(y)Q, с = Pa(y)Q. (2)

Здесь функции Pw, Pc выражаются через заданный нормированный профиль входной скорости

~(у), Pw (y) = í~(y)dy, Pa=- , Q(t) (расход

0 dy

жидкости в канале) подлежит определению.

Условия прилипания на твердых стенках имеют вид

х,yеу : T = 1, ^f = -V,, ш = C(t), i = 1N (3)

дп

dT дш

х, y еуу+1: — = 0, ^ = 0, ш = C (t), i = 0, N +1,

дп дп

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

На выходной границе канала yout ставятся «мягкие» условия [10], обеспечивающие свободный перенос температуры и завихренности за пределы расчетной области.

Для определения функций Ci (t) постановка задачи дополняется нелокальными граничными условиями [11]. При постоянных по времени и пространству скоростях V они будут иметь вид

ííf- Re

7, V 3П

ds = Re[p], i = 0, N +1

(4)

где 5 - касательная к границе; [р] - скачок давления вдоль контура у. На замкнутых контурах [р] = 0, на стенках канала - перепаду давления во входном и выходном сечениях.

Не нарушая общности, можно положить значение функции тока на одной стенке равным нулю, Сн+1 = 0 . Тогда величина С0 будет равна расходу 2 жидкости в канале.

Требуется найти решение системы уравнений (1) с граничными условиями (2), (3), и дополнительными нелокальными соотношениями (4).

Метод численного решения

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

о ——Ао = оа—т v-V о+rF, Re

— Ау = о.

T--AT = T — rvV-T .

Pe

(5)

(6) (7)

Здесь использовано стандартное обозначение

V

/(/) = /(/ - г). На каждом временном слое уравнения (5), (6) решаются с граничными условиями (3) на твердых стенках у1, причем условие Дирихле для у используется для уравнения (6), а с помощью условий Неймана формулируются граничные условия 1-го

рода для а вида х, у е/: а=~, ' = 0, N +1.

Для отыскания функции а решается линейная задача [12]

ду

а = —Ау, x,y е п: = —V, i = 0,N +1. дп

(8)

Функция тока у на текущем временном слое определяется набором констант С. Поэтому и завихренность а будет зависеть от этих констант. Задачи (5), (6), (8) линейны, следовательно, можно записать линейное представление искомых функций через С

W = W + T,W'Ct, о = о +2оiCi

(9)

в котором вспомогательные функции у*, у', а , а' подлежат определению, а множители С будут най-

дены из системы уравнении

ZliCj = I* + Re[p], i = 0N,

j=0

(10)

Ij =¡0ds, I* = j

n дп n

(

до дп

+ Re gs

ds,

Здесь 8/ - символ Кронекера. Граничные условия для вспомогательных функций в выходном сечении те же, что и для искомой.

Вспомогательные задачи (11), (12) не содержат

констант С ; кроме того, уравнения для функций у' однородны и могут быть решены один раз на шаге инициализации.

Аналогично, подставляя (9) в уравнение переноса

вихря (5), сформулируем задачи для вспомогательных

*

функций * , ' :

А а * = Ъ, А а' = 0, г = , (13)

А = Е--А, Ъ =1 Е-гуУ| а + гF,

Яе I )

где Е - единичный оператор.

Граничные условия для уравнений (13) получаются в результате подстановки (9) в исходные граничные условия (2), (3) для и приравнивания коэффициентов при С :

x, y е п : о * = ю*, о j = юj, i, j = 0, N, (14)

X У е Пп : о * = а о j = &0ра, j = 0, N,

ю* ~i

где о , о определятся как решения задач

ю* * ду

ю = —Ау , x,yenj : —— = —Vj дп

(15)

а' = —Ау', x,y е п : —— = 0 .

ду

дп

полученной с помощью подстановки представления (9) в интегральные условия (4).

Функции у*, у' находятся из решений вспомогательных задач, сформулированных в результате подстановки представления функции тока из (9) в линей* N ' V

ное уравнение Пуассона (6) - А у - С^ = а .

'=0

Поскольку правая часть полученного уравнения не зависит от С , его можно разбить на N + 2 уравнения следующим образом:

-Ау* = а, -А у' = 0, ' = . (11)

Граничные условия Дирихле для уравнений (11) получаются подстановкой разложения (9) в (2) и (3):

» N _

х,у е / : у + £ у= Сг-, ' = 0,N,

1=0

» N

х у е7ш : у + £ = Р (у)с0.

j=0

Приравняв члены при С , получим

х,уе/: у* = 0, у =8, ', j = , (12)

х, У е : у = 0, у = 8{Ру, j = 0, N.

следующих из (8).

Правые части уравнений (13) для а', ' = 0, N равны нулю, а граничные условия (14), (15) зависят только от у', которые, как отмечалось ранее, вычисляются один раз. Следовательно, и функции а' могут быть найдены на шаге инициализации. Непосредственно на каждом временном слое нужно решить зада-^ * * ^

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

Для определения I*, I/ из системы (10) необходимо найти значения нормальных производных да * / дп , 5®/ дп, j = 0, N на твердых стенках. Ниже представлен способ их вычисления в рамках метода конечных элементов (МКЭ).

Требуется определить значение производной да */ дп в узле (х1, у), лежащем на границе / . Для этого 1-е из уравнений (13) умножается на пробную функцию ^, равную единице в выделенном узле ' и нулю во всех остальных. После интегрирования по области получим

, до , *

п дп

Re

(16)

Re rl * Г ^ * __ , |

Ф' ^ —^V© -V^i — b<Pi \dD,

r d v Re j

где pP - след функции p на границе yk. Чтобы раз-

V

V

V

V

V

V

V

V

V

i=0

i=0

решить полученное равенство относительно производной до* / дп в /-м узле, интеграл по границе в левой части уравнения (16) вычисляется с помощью

квадратурной формулы |/(« XI(^ )Ь/, К - \Mids,

у 1 у

в которой 5 - узлы интегрирования, совпадающие с узлами сетки; %1 - весовые коэффициенты. В ре-

зультате

f да до* да*

дп

дп

дп

(17)

После замены левой части (16) на (17) получим

М . (18)

дп Ь /

Аналогично для определения узловых значений производной до/ / дп, у - 0, N

да ФJ

дп

-, j = 0, N,

(20)

Ф j = J^ — -Vp^dD .

К

Яе

Ь V т

Итак, алгоритм численного решения нестационарной задачи (1) состоит из шага инициализации и цикла по времени.

На шаге инициализации сначала решается N +1 задача (11), (12) для вспомогательных функций у1 , которые подставляются в соотношения (15), из которых определяются функции О , дающие значения О на границах у.. Сами функции оО находятся как решения задач Дирихле (13), (14). По формулам (20) подсчиты-ваются производные дО / дп , после чего вычисляются

коэффициенты I/ системы уравнений (10).

На каждом шаге временного цикла сначала определяется температурное поле Т как решение уравнения (7) с граничными условиями (3). Затем решается одна задача (11), (12) для вспомогательной функции

у*. Она подставляется в соотношение (15), из которого находится функция о*, определяющая граничные значения оО на жестких стенках. Сама функция оо находится как решение задачи (13), (14). По формуле (18) подсчитываются производные до / дп , после чего вычисляются коэффициенты I* правой части системы уравнений (10). Решение этой алгебраической системы дает набор констант С, который вместе с найденными вспомогательными функциями у*, *

у1, о , О позволяет найти решение у, о задачи

на текущем временном слое по формулам (9).

Для дискретизации уравнений по пространственным переменным и построения сеточных схем использовались стандартные процедуры МКЭ на неструктурированных четырехугольных сетках билинейных элементов с локальным сгущением около границ нагревателей и стенок канала. При аппроксимации конвективного слагаемого применялся ТУБ-подход [13].

Для проверки разработанного численного алгоритма проведена серия тестовых расчетов репрезентативных течений, таких как вынужденное обтекание системы тел с образованием дорожки Кармана [14], естественная конвекция в вертикальном канале около квадратного нагревателя [8], обтекание вращающихся цилиндров [15]. Во всех случаях получено вполне удовлетворительное согласование результатов расчета по описанному выше алгоритму с альтернативными аналитическими, численными и экспериментальными результатами.

Результаты расчетов

Нестационарная задача решалась с шагом по времени т — 0,002. МКЭ-сетка содержала около 40 000 элементов со средним размером в зоне сгущения к«0,01. Расчет производился в канале шириной Н — 7 и длиной Ь — 40 при Яе — 100, Ре — 70; 3 нагревателя профилированной формы располагались равномерно по ширине канала на расстоянии х — 5 от входного сечения. Периметр обтекаемых профилей равен п, что соответствует кругу единичного диаметра. Таким образом, принятые в модели числа Рей-нольдса и Пекле для тел сложной формы следует трактовать как вычисленные по диаметру эквивалентного кругового цилиндра.

На рис. 2 представлены фазы развития термоконвективного течения в канале. В начальной фазе жидкость медленно поднимается над нагревателями в виде тонких струй, образуя при этом характерные температурные «грибы», которые со временем объединяются в один (рис. 2а). Эпюра скорости в ближнем следе имеет ярко выраженные максимумы продольной скорости над телами. Температурные возмущения и вызванное ими течение на этом этапе локализованы вблизи нагревателей (рис. 3а, б, кривые 1).

Рис. 2. Поле температур при естественной конвекции от профилированных нагревателей на моменты времени г — 20 (а), 50 (б), 80 (в), 140 (г) и 360 (д)

Со временем жидкость над телами прогревается, что приводит к формированию 3 температурных факелов (рис. 2б) и резкому возрастанию расхода ( (рис. 4). При этом тяга обеспечивается, главным образом, самой нагретой жидкостью, а нагреватели, напротив, препятствуют течению. Это проявляется в виде локальных минимумов продольной скорости над телами в ближнем следе (рис. 3 а, кривая 2).

Рис. 2. Эпюры продольной скорости жидкости в сечениях: а - х = 10; б - х = 25 на моменты времени г = 20 (кривые 1), 50 (2) и 80 (3)

3-я фаза развития естественной конвекции характеризуется возрастанием местного числа Рейнольдса до величин порядка 300, потерей устойчивости симметричного течения и формированием дорожек Кармана в следе за телами. Вначале это происходит в следе за центральным профилем, а несколько позже -и за крайними (рис. 2в). Взаимодействие вихревых дорожек приводит к сложному течению в дальнем следе и квазипериодическому колебанию расхода Q (врезка при t«100 , рис. 4).

Интенсивное перемешивание вызывает выравнивание температуры поперек канала и как следствие замедление течения. Происходит затухание автоколебаний сначала в следе за крайними профилями (рис. 2в). При этом колебания расхода жидкости в канале приобретают гармонический характер, как показано на врезке г« 200 (рис. 4).

Рис. 4. Расход жидкости в канале

На заключительном этапе интенсивность вихревых структур постепенно затухает и за центральным нагревателем (врезка при г« 350 , рис. 4). Устанавливается стационарный режим конвекции (рис. 2г) с эффективным числом Рейнольдса порядка 200.

Выводы

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

сти жидкости происходит периодический срыв вихрей с обтекаемых профилей и формирование дорожки Кармана. На завершающем этапе течение замедляется и становится стационарным.

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

Работа выполнена при финансовой поддержке РФФИ, гранты 08-01-00548, 08-01-00163, 07-01-00499.

Литература

1. Piva S., Barozzi G.S., Collins M.W. Combined convection and wall conduction effects in laminar pipe flow: numerical predictions and experimental validation under uniform wall heating // Heat and Mass Transfer. 1995. Vol. 30, № 6. P. 401-409.

2. Ермолаев И.А., Жбанов А.И. Смешанная конвекция в горизонтальном канале при локальном нагреве снизу // Изв. РАН. МЖГ. 2004. № 1. С. 33-40.

3. Егоров А.Г., Мазо А.Б. Термоконвективный эффект при медленном течении вязкой жидкости в горизонтальном канале // Докл. РАН. Механика. 2007. Т. 415, № 6. С. 751-754.

4. Farouk B., Guceri S.I. Natural and mixed convection heat transfer around a horizontal cylinder within confining walls // Numer. Heat Transfer. 1982. № 5. P. 329-341.

5. Sadeghipour M.S., Razi Y.P. Natural convection from a confined horizontal cylinder: the optimum distance between the confining walls // Int. J. Heat Mass Transfer. 2001. № 44. P. 367-374.

6. Мазо А.Б. Численное моделирование свободной конвекции вязкой жидкости в канале с нагретым цилиндром // Уч. зап. Казан. ун-та. Серия физ.-мат. науки. 2005. Т. 147, кн. 3. C. 141-147.

7. Saha A.K. Free convection in a vertical channel with a built-in heated square block // Proc. of the 4th ISHMT-ASME Heat Mass Transfer Conf. and 15th National Heat Mass Transfer Conf. Pune, India, 2000. P. 533-538.

8. Khodary K., Bhattacharyya T.K. Optimum natural convection from square cylinder in vertical channel // Int. J. of Heat and Fluid Flow. 2006. № 27. P. 167-180.

9. Naylor D., Tarasuk J.O. Natural convective heat transfer in a divided vertical channel: Part I - numerical study // ASME J. Heat Transfer. 1993. № 115. P. 377-387.

10. Флетчер К. Вычислительная гидродинамика. Ч. 2. М., 1991. 552 с.

11. Glovinski R. Finite element methods for incompressible viscous flow // Handbook of numerical analysis. Vol. 9. Numerical Methods for Fluids (Part 3). Amsterdam, 2003. 1176 p.

12. Мазо А.Б., Даутов Р.З. О граничных условиях для уравнений Навье-Стокса в переменных функция тока -завихренность при моделировании обтекания системы тел // ИФЖ. 2005. Т. 78, № 2. C. 75-79.

13. Kuzmin D., Turek S. High-resolution FEM-TVD schemes based on a fully multidimensional flux limiter // J. of Comp. Physic. 2004. № 198. P. 131-158.

14. Liu J.-G., Wang C. High order finite difference methods for unsteady incompressible flows in multi-connected domains // Computers & Fluids. 2004. № 33. P. 223-255.

15. Flow characteristics of two rotating side-by-side circular cylinder / H.S. Yoon [et al.] // Computers & Fluids. 2009. № 38. P. 466-474.

Поступила в редакцию

24 ноября 2009 г.