Математическое моделирование нестационарных режимов теплопереноса в замкнутом двухфазном цилиндрическом термосифоне в условиях конвективного теплообмена с внешней средой

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

2011 Математика и механика № 1

УДК 536.24

Г.В. Кузнецов, М.А. Аль-Ани, М.А. Шеремет

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НЕСТАЦИОНАРНЫХ РЕЖИМОВ ТЕПЛОПЕРЕНОСА В ЗАМКНУТОМ ДВУХФАЗНОМ ЦИЛИНДРИЧЕСКОМ ТЕРМОСИФОНЕ В УСЛОВИЯХ КОНВЕКТИВНОГО ТЕПЛООБМЕНА С ВНЕШНЕЙ СРЕДОЙ

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

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

Замкнутый двухфазный термосифон (ЗДТ), который является по сути гравитационной бесфитильной тепловой трубой, использует испарение и конденсацию рабочей жидкости внутри теплообменника для переноса энергии. В отличие от обычной тепловой трубы, использующей капиллярную силу для возвращения жидкости в испаритель, ЗДТ использует гравитацию для возвращения конденсата. У него более простая конструкция, меньшее тепловое сопротивление, более высокий КПД и более низкая стоимость изготовления. Благодаря этим достоинствам ЗДТ широко используется во многих сферах, таких как промышленная регенерация тепла, охлаждение электронных компонентов и лопаток турбин, солнечные системы отопления [1-3].

Большинство проведенных исследований [4-8] направлено на выявление характерных свойств термосифонов и анализ области их применения. Так, например, результаты изучения влияния коэффициента наполнения рабочей жидкости на параметры стационарных режимов переноса тепла двухфазного замкнутого термосифона на основе многопараметрической математической модели представлены в [8]. Рассмотрены три варианта взаимодействия пленки жидкости и парового канала. Суммарная интенсивность теплопередачи жидкостного резервуара, обусловленная естественной конвекцией и пузырьковым кипением, определяется на основе комбинации их эффективных (рабочих) площадей и коэффициентов теплопередачи.

Предложены новые соотношения эффективной площади, основанные на экспериментальных результатах других работ [4-7]. Проведены натурные исследования двух различных геометрических конфигураций ЗДТ с азотом в качестве рабочей среды, которые показали достаточно хорошее согласование с теоретическими расчётами. Основной вывод работы связан с определением диапазона изменения коэффициента наполнения, в котором ЗДТ может оставаться устойчивым и эф-

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

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

Исследование характеристик умеренных режимов работы ЗДТ на основе одномерной модели течения пара с использованием корреляционных соотношений, определяющих эффекты кольцевого режима двухфазного потока, проведено в [10]. Коэффициент теплопередачи на поверхности пленки жидкости определяется как в ламинарном, так и в турбулентном приближениях. Математическая модель включает в себя дифференциальные законы сохранения массы, количества движения и энергии в паровом канале и в пленке жидкости. Одним из основных недостатков представленной модели является отсутствие влияния окружающей среды вследствие пренебрежения теплопроводностью ограждающих стенок, что может приводить к отличным от реального физического процесса результатам [11]. Необходимо отметить, что высокая интенсивность теплопереноса при наличии механизмов фазового перехода в термосифоне приводит к существенному повышению роли теплопроводности твердой оболочки даже в случае тонких стенок, выполненных из материалов с высокой проводимостью.

Математическое моделирование переходных режимов функционирования термосифона на основе плоских уравнений ламинарного потока сжимаемого идеального газа проведено в [12]. В начальный момент времени пренебрегалось наличие некоторого объема жидкости в зоне испарителя, что может значительно изменить характеристики переходных процессов в термосифоне. Вместо этого в модели рассматривался критический объём наполнения, обусловленный наличием достаточного количества рабочей жидкости для увлажнения стенок термосифона. Необходимо отметить, что наличие критического наполнения предполагает зависимость массы рабочей жидкости от тепловой нагрузки, что несколько не согласуется с условиями работы термосифонов.

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

Рассматривается краевая задача естественной конвекции в замкнутой цилиндрической области с теплопроводными стенками конечной толщины (рис. 1) при наличии локальных участков испарения и конденсации, отражающих зоны фазового перехода. В отличие от [13] в настоящей работе учитывается конвективный теплообмен с окружающей средой на границе г = Ь2, отражающий реальные условия передачи энергии от нижних слоев теплообменника к верхним. Учет конвективного теплообмена с окружающей средой на верхней границе области решения позволяет оценить эффективность и производительность термосифона в реальных условиях работы.

При математическом моделировании предполагается, что внешние поверхности вертикальных стенок являются адиабатическими, а на нижней границе г = 0 поддерживается постоянная температура [13]. Рассматриваемая геометрия задачи и граничные условия позволяют исключить влияние азимутальной координаты и проанализировать процесс переноса массы, импульса и энергии в осесимметрич-

ной постановке. При проведении вычислительных экспериментов предполагалось, что теплофизические свойства материала стенок, пара и жидкости не зависят от температуры, а режим течения является ламинарным. Пар считался вязкой, теплопроводной, ньютоновской жидкостью, удовлетворяющей приближению Буссинеска. Для исследования гидродинамических режимов в пленке жидкости использовался аналитический подход, изложенный в [14].

Математическая модель была сформулирована с учетом следующих допущений:

- конвективный теплообмен с окружающей средой на внешних границах вертикальных стенок считался пренебрежимо малым по сравнению с интенсивной продольной теплопроводностью;

- течения в паровом канале и в пленке жидкости рассматривались в ламинарном приближении;

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

Принимая во внимание сделанные допущения основные уравнения переноса массы, импульса и энергии в безразмерных переменных «функция тока - завихренность - температура» в цилиндрических координатах в системе «паровой канал-пленка жидкости» примут вид [13]

0

Рис. 1. Область решения задачи: 1 - пар; 2 - пленка жидкости; 3 - металлический корпус; 4 - поверхность испарения; 5 - поверхность конденсации

У2Т1 - — ^ = - ЯЦ;

Я ЗЯ

5©1 + 5Ц0) + 3(^0!) = 1

5т 5Я д1 ^Яа, • Рг1

2 и 0,

У20, —^ •

Я

2 РЯг -Ф „2

(г1 + г2 )/2

4^2

Я2 + С, 1п Я + С2, Т 2 = | У2 ёЯ ;

Г1І 22

д02

1

гУ202.

дт л/^ • РГ2

Уравнение теплопроводности в стенках термосифона:

503

3 =У203.

5Бо

Здесь Яа, = gz рДТг^/va1, Яа2 = gz РДГ23/va2 - числа Рэлея в паровом канале и в пленке жидкости; V - кинематический коэффициент вязкости; а - коэффициент температуропроводности; Рг, = у,/а,, Рг2 = у2/а2 - числа Прандтля пара и жид-

„2 1 2 а3/0

кости; У =—— | Я _— |+-------г - безразмерный оператор Лапласа; го = ^р-

я 5Я ^ 5Я у1 5г2 2

число Фурье, ф - градиент давления. Постоянные С1, С2 определяются из граничных условий.

Число Нуссельта, отражающее интенсивность теплообмена на границе раздела

50

пар - твердая стенка 2 = 2, /г2, вычисляется по формуле: Ми = Г ёЯ .

„ д2

Безразмерные граничные условия для уравнений (1) - (6) получим в виде

Я = 0, 0 < г< Ь,1х2, -------= 0, ¥ = 0;

г' 2 дЯ

Я = ^-, 0<г< Ь7/22, — = 0; г’ г' 2 дЯ

^2 2

501 = х 502

дЯ 2,1 дЯ

0, = 02,

дТ

дЯ

5Т

дЯ

2

5Т

302 503

--2 = Х32—3, 0? =03, -------

дЯ ’2 дЯ дЯ

= Т = 0;

г = о, о <я < Ь-, 0 = 0Й;

22

г = 0 <Я<^-, — = Бі(0е-0);

г2 г2 д2

2

г, г 50

г =^, о <я <^-

50,

2

2

= Я.,,3 ^-1+д „

дЯ дЯ

исп> 03 0,;

2, + 22 Г, 50, 503 ^

7 =----’ 0<Я<~ЙЯ~ = хзл— -<2кон^

22 22 дЯ дЯ

кон, 0 =03,

где X21 = X2/Х1 - относительный коэффициент теплопроводности; ц21 = ц2/ц -относительный коэффициент динамической вязкости; 2исп, Qкон и wисп, wкон -безразмерные теплоты и скорости испарения и конденсации, Ы = ыЬг/Х3 - число Био.

Сформулированная краевая задача с соответствующими начальными и граничными условиями решена методом конечных разностей [13 - 17].

Разработанный метод решения был протестирован на модельной задаче. Рассматривалось плоское ламинарное конвективное течение вязкой теплопроводной жидкости [18] в замкнутой цилиндрической полости. На вертикальных стенках и на верхней горизонтальной стенке поддерживался постоянный тепловой поток. Нижняя горизонтальная стенка являлась адиабатической.

На следующих рисунках показано сравнение линий тока и полей температуры при различных значениях числа Рэлея с численными результатами [18].

¥

©

1т^

ОЩ 0,5Ш1^;:р1Щ

' п "11114\> N Л':/-',"Iм I

0----------------------------------

-1 -0,5 0 0,5 1 -1 -0,5 0 0,5 1

Рис. 2. Линии тока ¥ и поля температуры © при Рг = 0,7, Яа = 103: результаты [18] - a, полученные результаты - Ь

¥

©

а 0,5

-0,5 0

0,5 1 -1 -0,5 0

0,5 1

I

б 0,5 |Чу.-"

I , I • 1 ;•  } I I ' '| ' I  V Ч- / ' I , ' I Г •' Г / ^ *

о п;к;ОЭ

// ^Хуги1.1

0 I Т, . ' _I_,_I__ ' . I

-0,5 0

0,5 1 -1 -0,5 0

0,5 1

Рис. 3. Линии тока ¥ и поля температуры © при Рг = 0,7, Яа = 10 результаты [18] - a, полученные результаты - Ь

Результаты, представленные на рис. 2 - 4, наглядно показывают, что созданная математическая модель и используемый численный алгоритм решения приводят к достаточно хорошему согласованию с результатами других авторов.

Анализ полученных результатов

Численные исследования проводились для термосифона в форме цилиндра со стальными стенками, в качестве рабочей жидкости рассматривалась вода. Были выбраны следующие геометрические характеристики термосифона: высота -100 мм, радиус - 25 мм, толщина стенок - 2 мм.

На рис. 5 - 8 показаны характерные линии тока, поля скорости и температуры при различных числах Рэлея. Увеличение Яа приводит к интенсификации конвективного теплопереноса в паровом канале. Наблюдается рост объемов испаряемой жидкости при условии, что скорость движения в пленке жидкости является постоянной величиной. Необходимо отметить, что в зоне конденсации происходит разделение потока - некоторое количество пара конденсирует и образует течение в пленке жидкости, а оставшаяся часть отражает формирование циркуляционной зоны в паровом канале. При Яа = 103 (рис. 5) основным механизмом переноса энергии является теплопроводность, поэтому наблюдаются незначительные скорости циркуляции. Время передачи тепла из зоны испарения в зону конденсации достаточно большое, что несоизмеримо с реальными условиями использования термосифонов.

Увеличение температурного напора в 10 раз (рис. 6) приводит к формированию рециркуляционной зоны в нижней части теплообменника. Последнее отражает образование застойной тепловой области, препятствующей интенсивному переносу энергии в верхние слои термосифона. Дальнейшее увеличение числа Рэлея (рис. 7, 8) проявляется в смещении вихря в вертикальном направлении к зоне конденсации, что также может приводить к снижению КПД анализируемого объекта.

На рис. 9 показаны профили температуры в сечении R = 0 при различных значениях числа Рэлея. Увеличение перепада температуры в паровом канале приводит к росту температуры на оси цилиндра. Наиболее существенные изменения наблюдаются при 104 < Яа < 105, что объясняется как формированием термического факела, так и устойчивой рециркуляционной зоны в паровом канале (рис. 7).

* * н н у /

t * ' ? ^ . г+м

тМ + '? ! Г

жж?м+*

Г I I м м

^Л’\\*^’\*-Л

0 0,1 0,2

Рис. 5. Линии тока (а), поля скорости (б) и температуры (в) при Яа = 103

б в

' 1 « г ? г ? ?

* + + *Г Г* + М*7 ++ 1^^ , т

^ч

0 0,1 0,2

Рис. 6. Линии тока

при Яа = 104

в

Z . 1 -I

0,9

0,8 -

0,7 -

0,6 -

0,5 -

0,4 -

0,3 -

0,2 -

0,1 -

(ЩЙ

11'"II' _  '.|М 1

11 ..;;У

11 \ 1.1111 11!

щ

'Л,

I

I

11.1 '.'Л 0,-\Ч 11\|

1 --чп:;чи*

т* "Щ* 1*/^!

к:'

111  • V ,

1 ..х,

щр

’И

* *

л У?

., I

I I

I I

I I \\lil ''11

\\\‘| 11 ч Ч I 1 I I , \\ ,':'1

\\\\\\;' I ’\\VAV'

:%;л '.

ЭДЩ

---. == 1гЛ 5&*4Т-"Д

0 0,1 0,2 0 0,1 0,2 0 0,1 0,2 R

Рис. 7. Линии тока

при Яа = 105

Z

1 н

0,9 0,8 -0,7 -0,6 -0,5 0,4 -0,3 -0,2 -0,1 -

ЩШт

'I III

А № |-\ /'>

Г, Ч', I м,\ :'

Г,1'.

1,1'!'

...... //:

'У;

1*1 " +

‘Ть ^

*** 11 + 4-4*

I ч:ч М1 **.

I 1 1^4

\!

1111 о I в! 1* 1 ' .

Шк\л\

мштШё

0 0,1 0,2 0 0,1 0,2 0 0,1 0,2 R

Рис. 8. Линии тока (а), поля скорости (б) и температуры (в) при Яа = 106

в

Рис. 9. Профили температуры в сечении R = 0

На рис. 10, а показана зависимость среднего числа Нуссельта от числа Рэлея на границе Z = z1lz2. Следует отметить слабоинтенсивные режимы переноса энергии при Яа < 5 104, отражающие доминирование кондуктивного механизма теплообмена. Дальнейшее увеличение температурного напора приводит к существенному повышению числа Нуссельта Ми.

Рис. 10. Зависимость среднего числа Нуссельта от числа Рэлея

Рис. 10, б демонстрирует нестационарность анализируемого процесса. Изменение среднего числа Нуссельта во времени наиболее существенно на начальном этапе 0 < т < 5 . При т> 5 наблюдается достижение термически стабилизированного режима.

Профили вертикальной компоненты скорости в различных сечениях представлены на рис. 11. Распределения на рис. 11, а еще раз подтверждают формирование

слабоинтенсивного конвективного течения при Яа = 104. По мере удаления от границы испарения наблюдается увеличение скорости течения пара, что характеризует наличие температурного напора вниз по потоку.

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

Сравнение полученных распределений изотерм и линий тока со случаем постоянной температуры, поддерживаемой на границе Z = Lz/z2 [13], отражает существенное влияние вида граничных условий на начальном временном этапе. С течением времени различия становятся незначительными, что объясняется влиянием оболочки термосифона, а также достижением установившегося режима переноса энергии. На рис. 12 представлено сравнение профилей температуры в характерных сечениях термосифона. Использование граничных условий, определяющих

Рис. ll. Профиль вертикальной компоненты скорости в различных сечениях Z = const при Ra = lO4 -

Рис. 12. Профили температуры при Яа = 105 в сечении R = 0 при т = 2 -

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

Заключение

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

анализируемых режимов течения и теплопереноса при Ra = 103 -106, Pr = 0,7, 0 < т < 100. Установлены особенности формирование термогидродинамических режимов при изменении температурного напора в анализируемом объекте. Показано, что при 0 < т < 5 наблюдаются значительные изменения среднего числа Нуссельта на поверхности испарения. Установлены существенные отличия в локальных распределениях линий тока и изотерм со случаем постоянной температуры на внешней поверхности.

ЛИТЕРАТУРА

1. Esen M., Esen H. Experimental investigation of a two-phase closed thermosyphon solar water heater // Solar Energy. 2005. V. 79. P. 459-468.

2. Desrayaud G., Fichera A., Marcoux M. Numerical investigation of natural circulation in a 2D-annular closed-loop thermosyphon // Int. J. Heat Fluid Flow. 2006. V. 27. P. 154-166.

3. Бакиев Т.А., Юсупов С.Т. Перспективы применения термосифонов в газовой промышленности // Материалы научно-технической конференции. М.:ООО «ИРЦ Газпром», 2005. С. 16-22.

4. Imura H., Ssasaguchi K., Kozai H. Critical heat flux in a closed two phase thermosyphon // Int. J. Heat Mass Transfer. 1983. V. 26. P. 1181-1188.

5. Shiraishi M., Kikuchi K., Yamarcishi T. Investigation of heat transfer characteristics of a two phase closed thermosyphon // Proc. Fourth Intern. Heat Pipe Conf. 1981. P. 95-104.

6. Ueda T., Miyashita T., Chu P.H. Heat transport characteristics of a closed two-phase thermosyphon // Trans. JSME. 1988. Part B 54. V. 506. P. 2848-2855.

7. El-GenkM.S., Saber H.H. Heat transfer correlations for small, uniformly heated liquid pools // Intern. J. Heat Mass Transfer. 1998. V. 41. P. 261-274.

8. Jiao B., Qiu L.M., Zhang X.B., Zhang Y. Investigation on the effect of filling ratio on the steady-state heat transfer performance of a vertical two-phase closed thermosyphon // Appl. Therm. Eng. 2008. V. 28. P. 1417-1426.

9. Seok-Ho Rhi An Experimental and Analytical Study on Two-Phase Loop Thermosyphons. Very Small to Very Large Systems // The Ottawa-Carleton Institute for Mechanical and Aeronautical Engineering. D.S.Thesis. Canada, 2000.

10. Reed J.G., Tien C.L. Modeling of the Two-Phase Closed Thermosyphon // ASME. J. Heat Transfer. 1987. V. 109. P. 722-730.

11. Sheremet M.A. The influence of cross effects on the characteristics of heat and mass transfer in the conditions of conjugate natural convection // J. Engineering Thermophysics. 2010. V. 19, № 3. P. 119-127.

12. Harley C., Faghri A. A complete transient two-dimensional analysis of two-phase closed thermosyphons including the falling condensate film // ASME. J. Heat Transfer. 1994. V. 116. P. 418-426.

13. Кузнецов Г.В., Аль-Ани М.А., Шеремет М.А. Численный анализ влияния температурного перепада на режимы переноса энергии в замкнутом двухфазном цилиндрическом термосифоне // Изв. Томского политехнического университета. 2010. Т. 317. № 4. С. 13-19.

14. Семенов П. Течение жидкости в тонких слоях // ЖТФ. 1944. Т. 14. № 7-8. С. 427-437.

15. Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. б1б с.

16. Самарский А.А. Теория разностных схем. М.: Наука, 1977. б5б с.

17. Пасконов В.М., Полежаев В.И., Чудов Л.А. Численное моделирование процессов тепло-и массообмена. М.: Наука, 1984. 288 с.

18. Lemembre A., Petit J.P. Laminar natural convection in a laterally heated and upper cooled vertical cylindrical enclosure // Int. J. Heat Mass Transfer. 1998. V. 41. P. 2437-2454.

Статья принята в печать 17.12.2010г.

Ruznetsov G.V., Al-Ani Maathe, Sheremet M.A. MATHEMATICAL SIMULATION OF TRANSIENT HEAT TRANSFER IN A TWO-PHASE CLOSED CYLINDRICAL THERMOSiPHON IN CONDITIONS OF CONVECTIVE HEAT EXCHANGE WITH AN ENVIRONMENT. Numerical simulation of mass transfer, momentum, and energy transport regimes in a two-phase closed cylindrical thermosiphon in conditions of convective heat exchange with an environment has been carried out. The mathematical model has been formulated in dimensionless variables such as stream function, vorticity, and temperature in cylindrical coordinates. Distributions of local thermo-hydrodynamic parameters reflecting influence of an environment have been obtained. Stages of energy transport from the evaporation zone to the condensation zone of the thermosyphon have been determined.

Keywords: two-phase closed thermosyphon, natural convection, solid walls, Newton - Richmann law, cylindrical coordinates.

KUZNETSOV Genii Vladimirovich

AL-ANI Maathe

E-mail: maathe_a@yahoo.com

SHEREMET Mikhail Aleksandrovich

E-mail: sheremet@math.tsu.ru