CC BY
178
УДК 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 (Z1 )= Z H0zk [Z zk/T0k - R Е zk ln z* I , Z2 =1 - zn (1)
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(1 - T1/ T)/RT1 - ln z1, ln Y^iuid = H02(1 - T02/T)/RT02 - lnz2. (2)
Находя логарифмы коэффициентов активности в точке эвтектики ln yle, ln у2е по формулам (2) и применяя термодинамически согласованную модель Редлиха-Кистера, рассчитаем коэффициенты ряда Редлиха-Кистера по формулам:
в = ln к / z2 -1/ 2(ln Ye / zi - ln / z2, К, - z2,). c = 1/ 2(ln yu / zL - ln f2, / z). (3)
Найдем состав раствора на ветвях ликвидуса в зависимости от температуры по формулам: 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)].
Из условия термодинамической согласованности построенной модели определяется коэффициент ассоциации к в жидкой фазе, имеющий смысл числа молекул, объединившихся в кластер. Если пренебречь влиянием на равновесие разности теплоемкостей жидкой и твердой фаз, то уравнение термодинамического равновесия между твердой и жидкой фазами можно записать в виде:
= (у 1щЫ / ум) гкрм ехр[- нк (1/Тк -1/т)/(4)
где ук911“1, Ук°ш - коэффициенты активности в жидкой и твердой фазах; к = 1, 2 - номер компоненты. В точке эвтектики кривые ликвидуса и солидуса имеют минимум. В случае образования эвтектики следует использовать условие одинакового состава жидкой и твердой фаз в точке эвтектики: гье9и11 = г1°ш. Из уравнения (4) следует:
У1°ш = У19и1 ехр[- нк (1/Т0к -1/Т)/я} (5)
По значениям 1п у^Лп У2°е11й, рассчитанным по формулам (5) в точке эвтектики, мож-
но найти коэффициенты ряда Редлиха-Кистера и определить кривые солидуса. Проведение процедуры термодинамического согласования коэффициентов активности компонентов позволяет вычислить коэффициенты ассоциации в жидкой фазе. При известной модели ассоциации в твердой фазе, сведения о которой можно получить из спектрального анализа, возможно более точное моделирование фазового равновесия жидкость-твердое тело. Энтальпия плавления в эвтектической точке рассчитывается по формуле:
не = [(к1 + Г1 )х1Н 1 + (к2 + Г2 )х2н02 ] , (6)
где к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]. На основе данных о взаимной растворимости компонентов можно определить
коэффициенты ассоциации в твердой фазе.
Расчеты показали, что модель РСБЛБ позволяет прогнозировать основные термодинамические характеристики фазопереходных материалов.
Вычислим эффективные значения объемной теплоемкости (рСр) , плотности рг , температуропроводности аг для г-й фазы, удельной теплоты плавления Не в точке эвтектики и теплопроводности Д:
(рСр) = Е ркс№ = I С^М,/Е х^М, / рк, г = 1,2, (7)
2 , , 2 , 2 хкМ 2 2
Л = Е рк?к = Ехкмк/Е= Е ХкеМк/Е хквМк / рк, (8)
к =1 к =1 к=1 р; к=1 к =1
а, = V (рС_)г, (9)
2 2 2 2
Не = ЕНок-к;/Е-кМк = ЕНоЧе/ Е-кМк • (10)
к=1 к =1 к =1 к =1
л=Е =Е 4к-Мк / ркі Е -кЄМк / р\, (11)
к=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], (12)
дп (1 -|) дп 1 -1
дв2 _ а2 д 2в2 п _ 1 - У , пє [0,1], (13)
дп І дп І
двз + аз_двз _ аз д2вз п _ У1 - У (14)
дт ( У1 -1)[У1 + п( У -1)] дп (У1 -1)2 дп2’ У1 -1 ’
тт dl к1 дв1 / \ к2 дв2 / \ & dl
apH^ _^^ (т,0) + (т,0), І _^~. (15)
dT 1 -| дп І дп dт
Г раничные и начальные условия:
вх (т,0)_ в2 (т,0)_ 1, ^ (т,1)_ 0, (16)
дп
в, (т,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
