Термодинамика фазовых переходов и перенос тепла в аккумуляторе

УДК 517.958

З.Н. Есина, Е.Е. Зеленский, М.Р. Корчуганова

ТЕРМОДИНАМИКА ФАЗОВЫХ ПЕРЕХОДОВ И ПЕРЕНОС ТЕПЛА

В АККУМУЛЯТОРЕ

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

Тепловые аккумуляторы, фазовые переходы жидкость-твердое, фронт фазового перехода, время затвердевания

Z.N. Esina, Е^. Zelensky, M.R. Korchuganova

THERMODYNAMICS ОF PHASE CHANGES AND CONDUCTION OF HEAT IN ACCUMULATOR

In operation the comprehensive approach to calculation of a heat accumulator is implemented on the basis of binary mixtures permitting on temperature and melting heat of pure builders to calculate the performances in a point of an eutectic of an intermixture, and also time-space fields of temperatures current standing of front of phase change and setting time (melting) of a material in a cylindrical heat accumulator.

Heat accumulators, phase changes fluid-solid, front of phase change, setting times

Важным направлением альтернативной энергетики является использование солнечной энергии, ветра и источников тепла природного и техногенного происхождения. Рациональное использование невозобновляемых источников энергии, таких как нефть, уголь, природный газ, является одной из главных задач теплоэнергетики. Немаловажное значение имеет также накопление и преобразование энергии с помощью теплоаккумуляторов и хладоакку-муляторов, основанных на принципе фазовых переходов. Тепловые аккумуляторы (ТА) на основе фазовых переходов твердое тело-жидкость перспективны тем, что они имеют высокую плотность хранения энергии при постоянной температуре Tn выделения (накопления) тепла, находящийся в широком диапазоне изменения. Использование фазопереходных материалов (ФПМ) с низкой (жирные кислоты, парафины) и высокой температурой плавления (соли цветных металлов) позволяет создать ТА для аккумулирования тепловой энергии, выделяющейся при работе различных технических устройств и солнечной энергии [1]. В работе [1] приведена основная группа солей, которые служат как первичные материалы для разработки перспективной модели, предназначенной для применения в области высоких температур и бинарные, тройные и многокомпонентные системы на основе этих солей. В [2] рас-

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

В данной работе рассматриваем результаты моделирования фазового равновесия жид-кость-твердое в этих системах, выполненного методом минимизации избыточной энергии Г иббса по параметру ассоциации и применения разработанной нами программы PCEAS (фазовые переходы в эвтектических и азеотропных системах). В работе реализован комплексный подход к расчёту ТА на основе бинарных ФПМ, позволяющий по температуре и теплоте плавления чистых компонентов рассчитать характеристики в точке эвтектики смеси, а также пространственно-временные поля температур, текущее положение фронта фазового перехода Z(t) и время затвердевания (плавления) ФПМ в цилиндрическом ТА.

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

В [3] рассматривалась возможность расчета параметров эвтектики на основе минимизации избыточной энергии Г иббса системы по параметру сольватации ц = /и1 / ц2, характеризующему соотношение числа молекул в молекулярных соединениях разнородных чистых компонентов. Минимизация избыточного потенциала Гиббса по параметру сольватации Л = (x1 z2/(x2z1)) реального раствора приводит к уравнению Бернулли относительно T (z1), решение которого имеет вид:

2(2 2 Л-1

T

k =1 V* =1 k =1 J

где R - универсальная газовая постоянная; z1 = x1 /(x1 + ¡1 x2), z2 = x12 /(x1 / ¡1 + x2) - эффективные мольные доли компонентов бинарной смеси; xi, x2 - мольные доли компонентов рас-

твора; Hk - молярная теплота плавления k-го компонента, k = 1,2 .

Из условия наличия минимума кривой T(z1) в точке эвтектики получим трансцендентное уравнение для определения эффективного эвтектического состава:

н 1H2 (1/г;- -1/т0 )=r[h 0 ln(1 - ^)- H 0 ln zu ].

Сравнение вычисленных величин zke (k = 1, 2) в точке эвтектики с экспериментальными значениями xke позволяет найти л, т.е. установить степень неидеальности раствора и получить термодинамически согласованную кривую ликвидуса. Подставив z1e в (1), получим температуру плавления смеси Te в точке эвтектики.

Уравнение (1) позволяет моделировать фазовый переход жидкость-твердое при P = const для неэвтектических и эвтектических систем. Для эвтектических систем состав раствора (z1e, z2e) и температура плавления Te в точке эвтектики определяются из условия dT / dz = 0 . Логарифмы коэффициентов активности компонентов в жидкой фазе запишем в виде:

ln^1i?“id = Н1

Находя логарифмы коэффициентов активности в точке эвтектики ln yle, ln у2е по формулам (2) и применяя термодинамически согласованную модель Редлиха-Кистера, рассчитаем коэффициенты ряда Редлиха-Кистера по формулам:

в = ln к / z2 -1/ 2

Найдем состав раствора на ветвях ликвидуса в зависимости от температуры по формулам: ln v1 = H0(1 - T1 / T)/ RT1 - ln Y1, ln w2 = H02 (1 - T02 / T)/ RT02 - ln у 2,

где V! - мольная доля 1-го компонента раствора на левой ветви кривой ликвидуса; w2 - мольная доля 2-го компонента раствора на правой ветви кривой ликвидуса; 1п у1, 1п у 2 - логарифмы коэффициентов активности по модели Редлиха-Кистера, рассчитанные по формулам:

1п 71 = г2 [в + С(3*1 - г2)], 1п у2 = г\ [В + С^ - Зг2)].

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

=

где ук911“1, Ук°ш - коэффициенты активности в жидкой и твердой фазах; к = 1, 2 - номер компоненты. В точке эвтектики кривые ликвидуса и солидуса имеют минимум. В случае образования эвтектики следует использовать условие одинакового состава жидкой и твердой фаз в точке эвтектики: гье9и11 = г1°ш. Из уравнения (4) следует:

У1°ш = У19и1 ехр[- нк

По значениям 1п у^Лп У2°е11й, рассчитанным по формулам (5) в точке эвтектики, мож-

но найти коэффициенты ряда Редлиха-Кистера и определить кривые солидуса. Проведение процедуры термодинамического согласования коэффициентов активности компонентов позволяет вычислить коэффициенты ассоциации в жидкой фазе. При известной модели ассоциации в твердой фазе, сведения о которой можно получить из спектрального анализа, возможно более точное моделирование фазового равновесия жидкость-твердое тело. Энтальпия плавления в эвтектической точке рассчитывается по формуле:

не = [

где к1, к2 - коэффициенты ассоциации компонентов в жидкой фазе, которые определяются

методом термодинамического согласования коэффициентов активности; г1, г2 - параметры, учитывающие комплексообразование в твердой фазе. Если при растворении компонента в жидкой смеси энтальпия растворения положительна, то г > 0, в противном случае г < 0. Для расчета коэффициента ассоциации компонента второго компонента в твердой фазе s2 = к2 + г2, при известных энтальпии плавления чистых компонент и раствора, а также коэффициенте ассоциации первого компонента в твердой фазе s2 = к2 + г2, можно использовать выражение, следующее из (6):

к2 + Г2 = Не /(Х2Но2) - (к + Г1)Х,Н 1 /(Х2Но2).

Соли цветных металлов и их композиции успешно применяют для накопления энергии при высоких температурах. В интервале температур от 250 до 1680° С, большой интерес представляют неорганические соли с теплотой плавления от 68 до 1041 Дж/г.

В табл. 1 представлены экспериментальные и расчетные данные для состава, температуры и теплоты плавления в бинарных системах на основе солей металлов.

В табл. 2 приведены экспериментальные и расчетные эвтектические параметры для различных композиций насыщенных жирных кислот.

Как можно видеть из табл. 2, экспериментальные и расчетные величины состава и температуры плавления хорошо согласуются. Измеренные и расчетные значения теплоты плавления различаются в большей степени. Последнее возможно связано со значительным разбросом опытных значений Н0 [1], либо с ассоциацией компонентов в твердой фазе.

Таблица 1

Термофизические свойства бинарных систем на основе солей металлов

Раствор Эксперимент [1] Расчет по модели минимизации избыточной энергии Гиббса

Состав эвт., Мол Температура пл., С° .л п Р И а 5 п е Т .т в І ° Ба о С Температура пл., С° .л п а 5 п е Т .оц ,е ° га -і* ё ^ Ж 5 11 >§ т .ц ,е оз ОГО ^ ГО ~ _о Ч СП & ш 11 го “ « о т ш

NaF-NaBr 0,270 0,730 642 360 0,292 0,708 640,56 369,35 4/5 8/5

NaF-LiF 0,400 0,600 652 816 0,531 0,469 642,39 806,97 10/11 20/9

NaF-NaCl 0,335 0,665 675 572 0,342 0,658 675,03 568,87 7/8 14/7

NaF-CaF2 0,680 0,320 810 600 0,658 0,342 794,1 601,88 16/15 32/7

NaF-MgF2 0,750 0,250 832 627 0,762 0,238 896,25 741,72 9/4 18/1

^-иа 0,500 0,500 487 344 0,573 0,427 464,4 346,90 5/7 10/2

KF-KBr 0,400 0,600 576 315 0,397 0,603 566 316,47 23/24 46/22

KF-KCl 0,450 0,550 605 407 0,430 0,570 607,78 407,77 23/24 46/20

KF-CaF2 0,850 0,150 780 440 0,751 0,249 747,09 447 7/6 14/2

LiF-LiСl 0,736 0,264 485 403 0,714 0,286 501,52 419,86 19/10 38/1

LiF-MgF2 0,670 0,330 746 947 0,630 0,370 761,98 1004,99 9/5 18/5

LiF-CaF2 0,815 0,195 769 820 0,753 0,247 745,41 822,15 14/11 28/8

Таблица 2

Экспериментальные и расчетные эвтектические параметры для различных смесей жирных кислот

Тепло- аккумулирующая композиция Массовые концентрации Т, к Н е, кДж / кг

Экспе- римент Расчет є=Дї/-, % Экспе- римент Расчет * = со Экспе- римент Расчет * = со

Лауриновая кислота 0,6250 0,6245 0,08 305,75 305,03 0,23 156 182,32 -16,87

Миристиновая кислота 0,3750 0,3755 -0,4

Лауриновая кислота 0,7550 0,7720 -2,25 310,45 310,18 0,08 171 180,79 5,73

Стеариновая кислота 0,2450 0,2280 6,94

Капроновая кислота 0,7520 0,7524 -0,05 295,26 297,74 -0,84 153 161,19 5,35

Пальмитиновая кислота 0,2480 0,2470 0,16

Как показано на примере систем на основе солей цветных металлов, учет ассоциации в твердой фазе позволяет уточнить расчетные величины энтальпии плавления в точке эвтектики Не [4]. На основе данных о взаимной растворимости компонентов можно определить

коэффициенты ассоциации в твердой фазе.

Расчеты показали, что модель РСБЛБ позволяет прогнозировать основные термодинамические характеристики фазопереходных материалов.

Вычислим эффективные значения объемной теплоемкости (рСр) , плотности рг , температуропроводности аг для г-й фазы, удельной теплоты плавления Не в точке эвтектики и теплопроводности Д:

2 , , 2 , 2 хкМ 2 2

Л = Е рк?к = Ехкмк/Е= Е ХкеМк/Е хквМк / рк,

к =1 к =1 к=1 р; к=1 к =1

а, = V

2 2 2 2

Не = ЕНок-к;/Е-кМк = ЕНоЧе/ Е-кМк •

к=1 к =1 к =1 к =1

л=Е =Е 4к-Мк / ркі Е -кЄМк / р\,

к=1 к=1 к=1

где - изобарная удельная теплоемкость; М - молярная масса; ф - объемная доля; нижний индекс - номер фазы (і = 1 - твердая фаза), верхний индекс - номер компонента в г-й фазе.

Кроме того учтено, что в точке эвтектики —ке = -2е = -ке. Аддитивность формул (7)-(10) для

(рСр) , рг, температуропроводности аг для г-й фазы и удельной теплоты плавления Не следует из аддитивности теплоемкости, плотности и теплоты плавления. Коэффициент теплопроводности Д не является аддитивной величиной, но при параллельном включении элементов обобщенная теория проводимости даёт аддитивную формулу (11), которая вследствие простоты применяется даже для изотропных сред - жидких растворов при небольших различиях между А и А [5].

Задача Стефана для цилиндрического случая после введения безразмерных переменных и применения метода выпрямления фронтов (метода функциональных преобразований) сводится к решению системы уравнений:

дв1 —1 + дт

П -1 & а1

W -(1 - І)[| (1 - п)+п]

дв2

—2 + дт

І(1 - ")+fb

М=-Ja^. n=f!, n.fei],

дп (1 -|) дп 1 -1

дв2 _ а2 д 2в2 п _ 1 - У , пє [0,1],

дп І дп І

двз + аз_двз _ аз д2вз п _ У1 - У

дт ( У1 -1)[У1 + п( У -1)] дп (У1 -1)2 дп2’ У1 -1 ’

тт dl к1 дв1 / \ к2 дв2 / \ & dl

apH^ _^^

dT 1 -| дп І дп dт

Г раничные и начальные условия:

вх

дп

в, (т,1) _ вз (т,1), -А- дв*- (т,1) _ -Л- (т, 1),

1 -1 дп У1 -1 дп

в з (т,0)_ в , в.(т, п)_ в0, 1 (0)_ 1.

Безразмерные переменные вводились по формулам

т _ аЧ2, У _ г/г„ | _ |/Г1, в. _ Т{/Т0, рк _ рк/р\, р _ р/р1,

к, _ АД1, Cp, _ Cr,/Ch, He _ (H,M 1 )/H0, a _ HJ/(m.C'7).

При записи формул (12)-(16) черта над безразмерными величинами для удобства была опущена.

Данные для цилиндрического ТА на основе бинарной смеси жирных кислот малодоступны, поэтому расчеты проводились для ТА с однокомпонентным ФПМ (парафин RT 60, Rubitherm GmbH), теплофизические характеристики и физическое моделирование которого описаны в [1]. Цилиндрическая вертикальная алюминиевая труба различного внутреннего радиуса: r1 _ 0,002 м, r1 _ 0,0з м, r1 _ 0,04 м, длиной l _ 0,6 м и толщиной стенок

h _ 0,00з м с расплавом парафина помещалась в водяную баню с температурой 7* > 7^. Затем температура 7* быстро понижалась до постоянной величины 7* < 7ц. Дно и крышка трубы теплоизолированы, что наряду с r1 << l , h << r1 обеспечивало одномерность процесса теплопроводности. Значения 701 , H 01 и теплофизические свойства RT 60 в твердой и жидкой фазах приведены в [1]. В этом случае в (12)-(16) необходимо положить 01 _ А1 _ р1 _ 1, к2 _ А1 /А _ 1, а2 _ а2/а1 _ 0,87, 0з _ а1 /0^ _ 758, Аз _ А/А _ 1020, y1 _ 1 + h/r1 _ 1,1, а _ 1,з8.

Краевая задача (12)-(16) решалась на основе солвера PDEPE среды программирования Matlab, сводящего уравнения в частных производных к системе ОДУ. Для решения последней применялась функция ode 15 s, реализующая метод Гира.

Расчеты показали, что для алюминиевой оболочки можно пренебречь её тепловым сопротивлением и значительно упростить постановку задачи, заменив граничные условия сопряжения в (16) условием в1 (т, 1) = , что практически не сказывается на результатах. Данный вы-

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

В табл. 3 приведено сравнение экспериментальных гэксп и численных гчисп значений интегральной характеристики ТА - времени затвердевания ФПМ для разных внутренних радиусов трубы.

Таблица 3

Сравнение опытных и расчетных времен затвердевания

r1, м 0,02 0,03 0,04

^эксп, час 5,75 12,92 23,05

tчислj час 5,83 13,125 23,33

Uисл5 час 12,44 28,0 49,78

Относительная ошибка не превышает 2%. В последней строке таблицы приведены расчетные времена затвердевания эвтектической смеси лауриновой и стеариновой кислот.

ЛИТЕРАТУРА

1. Kenisarin M. High-temperature phase change materials for thermal energy storage / Renewable and Sustainable Energy Reviews. 2010. №14. P.955-970.

2. Kenisarin M. Solar energy storage using phase change materials / M. Kenisarin, K. Mahkamov // Renewable and Sustainable Energy Reviews. 2007. № 11. P.1913-1965.

3. Есина З.Н., Корчуганова М.Р., Мурашкин В.В. Математическое моделирование фазового перехода жидкость - твердое // Вестник Томского государственного университета. Управление, вычислительная техника и информатика. 2011. №3 (16). С.13-23.

4. Есина З.Н., Корчуганова М.Р., Мурашкин В.В. Эвтектические смеси для теплоаккуму-ляторов // Ресурсосбережение в химической технологии. Сборник трудов международной конференции. СПб : СПбГТИ (ТУ), 2012. С.118-121.

5. Дульнев Г.Н., Заричняк Ю.П. Теплопроводность смесей и композиционных материалов. Л.: Энергия, 1974. 263 с.

Есина Зоя Николаевна —

кандидат технических наук, доцент кафедры «Вычислительная математика» Кемеровского государственного университета

Зеленский Евгений Евгеньевич —

кандидат физико-математических наук, доцент кафедры «Вычислительная математика» Кемеровского государственного университета

Корчуганова Маргарита Рашидовна -

старший преподаватель, соискатель кафедры «Вычислительная математика» Кемеровского государственного университета

Статья поступила в редакцию 29.02.12, принята к опубликованию 12.03.12