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

2011

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Управление, вычислительная техника и информатика

№ 3(16)

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

УДК 517.958

З.Н. Есина, М.Р. Корчуганова, В.В. Мурашкин

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ФАЗОВОГО ПЕРЕХОДА ЖИДКОСТЬ - ТВЕРДОЕ

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

Ключевые слова: равновесие жидкость - твердое, избыточная энергия Гиббса, минимизация по параметру сольватации, параметр ассоциации.

Математическое моделирование фазового равновесия жидкость - твердое для бинарных и многокомпонентных смесей имеет значение для многих областей промышленности, таких, как химическая, пищевая, фармацевтическая и др. [1-4]. Основное применение методы математического моделирования находят при решении задач выбора составов и расчете термодинамических параметров антифризов, теплоаккумуляторов, теплозащитных покрытий, моделирования состава нефти, ректификации с использованием разделяющих агентов. Изучение равновесного состояния важно для всех процессов, где существует взаимодействие между твердым телом и жидкой фазой. Для таких процессов термодинамические свойства иногда трудно измерить экспериментально, особенно при экстремальных условиях (например, высокое давление, низкая температура и т.д.), и поэтому необходимо математическое моделирование. Разработка теоретических методов, позволяющих выбрать оптимальный состав антифризов и теплохладоаккумуляторов с требуемой температурой и теплотой плавления, является актуальной задачей, представляющей научный и практический интерес.

Равновесие жидкость - твердое может быть рассчитано многими способами [5, 6]. Термодинамические модели для рассмотрения фазового равновесия растворов могут быть классифицированы на две категории: модели коэффициента активности и модели уравнения состояния. Термин «прогнозирующие термодинамические модели» характеризует типы моделей, которые могут описывать фазовое равновесие при условии, что молекулярные структуры или физические свойства чистых компонентов в смеси известны. Прогнозирующие термодинамические модели также можно классифицировать на две категории: модели с использованием экспериментальных данных, например модель «ЦЖРАС» [5], в которой параметры взаимодействия найдены из экспериментальных данных и прогноз термодинамических свойств делается на основе существующих параметров для атомных групповых вкладов молекул, присутствующих в смеси; и модели без ис-

пользования экспериментальных данных (или априори прогнозирующие модели), например модель «С08М0-Я8», в которой требуются только специфические атомные параметры и прогноз термодинамических свойств делается на основе квантовых химических вычислений, что дает необходимую информацию для оценки молекулярных взаимодействий в жидкостях. Кривые фазового равновесия могут дать информацию о характеристиках чистых компонентов и смеси. Достоверное определение фазового равновесия необходимо для корректного прогноза образования фаз и их композиций.

1. Минимизация избыточной свободной энергии Гиббса

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

Избыточная свободная энергия Гиббса ОЕ является разностью между свободной энергией Гиббса реальной и идеальной систем. Для идеальных систем, где отсутствует взаимодействие между молекулами, избыточная свободная энергия Гиббса равна нулю. Таким образом, избыточная свободная энергия Гиббса является мерой взаимодействия в системе, причем она может быть положительной или отрицательной, что связано с видом взаимодействия.

Если за параметры состояния принимаются мольные доли компонентов в смеси х,, температура Т и давление Р, то уравнение состояния фазы многофазной, многокомпонентной системы имеет вид [7]

ёОЕ = -8 ёТ + V ёР + Ец ёх, , где 8 - энтропия, V - объем, ц - химический потенциал компонента.

Условием равновесия системы является минимум внутренней энергии. Из этого условия следует, что ёи = 0. Если температура и давление постоянны, то для равновесия системы, представляющей собой одну фазу, должны выполняться равенства

ёТ = 0, ёР = 0, йО = 0. (1)

В случае двухфазной системы для равновесия необходимо, чтобы были равны химические потенциалы фаз ц' = ц" и выполнялись условия (1). Фазовая диаграмма жидкость - твердое бинарной системы может быть получена минимизацией свободной энергии Гиббса. Избыточная свободная энергия Гиббса бинарной системы дается выражением [7]

АОе = КТ Ех, 1п у, ,

где у,- - коэффициент активности компонента, , = 1,2; х, - мольная доля г-го компонента.

Разность уравнений состояния бинарной системы, для реальной и идеальной равновесных фаз, можно представить в виде [7]

(-АНе/КТ 2)ёТ + (А^/К^Р = £х, ё1п у, , где АНе - энтальпия смешения; АРе - избыточный объем; Р - давление раствора; Т - абсолютная температура; К - универсальная газовая постоянная.

Если происходит образование соединения чистых компонентов в растворе, то суммарная молярная масса ,-го компонента с учетом числа молекул данного вида X,, входящих в соединение, может быть рассчитана по формуле: ц/ = X,- ц,-, где ц - молярная масса чистого компонента. Среднее соотношение числа молекул в сольватах чистых компонентов X = Х1/Х2 характеризует устойчивую структуру раствора. С учетом изменения молярной массы эффективные мольные доли компонентов бинарной смеси вычисляются по формуле [7] = х1/(х1+ >х2),

22 = х2/(х^ + х2). В том случае, если из эксперимента известны данные о составе раствора х1э, х2э в точке экстремума температуры, можно найти отношение коэффициентов X = X1/X2 = (г2э/г1э)/(х2э/х1э). Одним из методов определения среднего параметра сольватации компонентов в растворе является вычисление Xср по формуле средней арифметической взвешенной, где в роли весов выступают концентрации одного из компонентов в точках возможных эвтектик:

Xср = Е X ] хэ ].

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

ё(1п х,- у)/ёТ = АН, /КТ 2, (2)

где АН, - парциальная молярная теплота растворения ,-го компонента в насыщенном им растворе.

Теплота перехода компонента, образующего твердую фазу, в жидкое состояние АН, складывается из теплоты плавления АН™ и дифференциальной теплоты разбавления компонента АН,раз6 до концентрации, соответствующей насыщенному раствору [7]. Теплота разбавления зависит от характера взаимодействия компонент в растворе и может быть как положительной, так и отрицательной. Поскольку теплота разбавления для большинства систем неизвестна, то расчет по уравнению (2) возможен лишь для идеальных систем, для которых у,- = 1 и АН, = АН,пл:

ё( 1п х{ у,)/ёТ = АНГ /КТ 2

Теплота плавления АН™ зависит от температуры:

АН,пл = АН,0 + (с,ж - с,т)(Т - Т,0),

где АН,0 - теплота плавления чистого ,-го компонента при температуре Т,°; с,ж - с^ - молярные теплоемкости ,-го компонента при постоянном давлении в жидком и твердом состояниях.

Отличие коэффициента X от единицы свидетельствует о наличии отклонения от идеальности в бинарной системе и необходимости перехода к эффективным мольным долям. Для реальных систем, в которых возможно образование молекулярных соединений, называемых сольватами, необходимо в уравнении (2) сделать переход к эффективным мольным долям:

ё (1п г, у,)/ёТ = АН, /КТ2,

где АН, = ДНш + АН,раз6 = АН,0 +(с,ж - с,т)(Т - Т°) + АН?336; АН,0 - энтальпия плавления и Т,0 - температура плавления компонента, образующего однокомпонентную фазу.

В том случае, если теплоемкости в жидком и твердом фазовом состояниях различаются незначительно и диапазон температур плавления невелик, пренебрегая теплотой разбавления, парциальные коэффициенты активности можно представить в виде 1п у,- = ДHІ 0(1 - Т,0 /Т)/(К Т°) - 1п г, + 1п Е, (г,), где Е, (г,) - произвольные функции от состава. При постоянном давлении и малом интервале температур плавления избыточную энергию Гиббса как функцию от эффективной мольной доли компонента находим из уравнения

АОе = ^[ДНЛт/ Т!0 - 1)] + г2[ДН2°(Т/Т1° - 1)] - КТ(г1 1пг1 + г2 1пг 2) + Е^). (3) Здесь Е^) = г1Е1(г1) + (1-г1)Е2(г1) - функция, которая выбирается из условия термодинамической согласованности модели по методу Херингтона и Редлиха -Кистера [7]:

1

I Щ / 72= °.

0

Согласно этому уравнению, должны быть равны площади, ограниченные кривой 1g у1/у2 и осями координат.

Минимизация избыточной энергии (3) по внутреннему параметру X приводит к уравнению Бернулли:

^Т/ё^ + /1(Ц)Т = /2(21)Т2, (4)

где/1(г1) = а/(АН10/К - а 21), /2^) = -[в + 1п (1-21)/г1]/(АН10/К - а 21), а = (АН10 - АН20)/К, в = АН20/КТ20 - АН10/К Т10.

Решение уравнения (4) имеет вид

Т(21) = [АН1021 + АН20(1-21)]/

/{ДН^Т0 + ДН20(1-21)/Т20 - К^1 1п 21 + (1-21) 1п(1-21)]}. (5)

Таким образом, минимизация избыточной энергии Гиббса АОе по параметру сольватации X позволяет найти уравнение, моделирующее Т2-диаграмму фазового равновесия. Эта диаграмма представляет собой сечение при Р = сош! трехмерной РТ2-диаграммы состояния бинарной системы.

Из условия экстремума функции (5) получено алгебраическое уравнение

1п^(1 -^ )ЛH1°/К .( уш2/К^ = ДН^дн0 (1/( -1/Т°)к2. (6)

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

с л н О Л (

3 + АИ

аи

2 У

з - АИі

О

аи.

2 У

3 +-

аи,

аи 0

О А2 ( /

3 —

аи,

аи о

2АИ,0

л

1 1

О

V т2

710

11

3

3 —

аи1

аи2о

о Л

(7)

Мольная доля компонента в растворе не превышает единицы, поэтому выбираем решение 2е (0,1). Подставляя 21 = 2\ в функцию (5), можно найти температуру плавления в эвтектической точке.

В общем случае, если известна парциальная молярная теплота растворения /-го компонента в насыщенном им растворе АИ,- , то в уравнениях (3) - (7) следует заменить энтальпию плавления чистого компонента АИ,0 на АИ,-.

2. Численные результаты

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

Т аблица 1

Расчет состава, температуры плавления эвтектики и параметров сольватации и ассоциации в растворах гликолей

Компоненты системы Состав эвтектики хД х2э масс Температура эвтектики Т, °С Энтальпия плавления в точке эвтектики ДЯплэ, кДж/моль Коэффициент сольватации X Коэффициент ассоциации к

Этиленгликоль 0,6977 -62,7 8,306 0,97 0,98

Вода 0,3023

Диэтиленгликоль 0,7306 -57,22 8,909 0,96 0,875

Вода 0,2694

Триэтиленгликоль 0,6871 -48,81 9,897 0,94 0,75

Вода 0,3129

Этиленгликоль 0,4364 -38,64 12,455 1 1,037

Диэтиленгликоль 0,5636

Триэтиленгликоль 0,4688 -29,87 15,475 1 0,917

Диэтиленгликоль 0,5312

На рис. 1 приведена кривая, моделирующая среднюю температуру плавления системы (этиленгликоль - вода); точками, соединенными ломаными линиями, показаны возможные в данной системе эвтектики. На рис. 2 приведены кривые ^ (уі/у2) для системы (этиленгликоль - вода), полученные в результате предварительного моделирования и проведения процедуры термодинамического согласования. На рис. 3 приведены кривые 1п у1, 1п у2 для системы (этиленгликоль - вода), полученные в результате проведения процедуры термодинамического согласования.

Из условия термодинамической согласованности построенной модели определяется коэффициент ассоциации к = к1 / к2, имеющий смысл отношения числа молекул компонента А, объединившихся в кластер до образования раствора или непосредственно в данном растворе, к числу молекул в кластере компонента В. Это позволяет построить кривые ликвидуса, близкие к экспериментально найденным (рис. 3).

Мольная доля 1-го компонента

Рис. 1. Средняя температура плавления раствора (этиленгликоль - вода); точками показаны предполагаемые эвтектики

Мольная доля 1-го компонента

Рис. 2. Зависимость lg у! Iy2 от мольной доли первого компонента в бинарной смеси (этиленгликоль - вода); кр. 1 - расчет по модели; кр. 2 - результат термодинамического согласования

0 0,2 0,4 0,6 0,8 1

Мольная доля 1-го компонента

Рис. 3. Зависимости ln уь ln у2 от мольной доли первого компонента в бинарной смеси (этиленгликоль - вода)

В том случае, если кривые 1п уь 1п у2 удовлетворяют условию термодинамической согласованности, представим их в виде ряда Редлиха - Кистера

1п У1 = 222[Б + С(32] - 22)], 1п У2 = 2]2[Б + С(2] - 3*2)] (8)

или аппроксимируем их уравнениями Ван-Лаара

1п у! = А [Б22 /(А2: + Б*2)]2, 1п у2 = Б[А21/(А21+Б22)]2.

По эвтектическим данным: 21э, 22э, 1п у:э , 1п у2э находим коэффициенты моделей Редлиха - Кистера и Ван-Лаара.

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

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

1п V] = ДЯДЯ Т:0 (1 - Г:0 /Г) - 1пу: , 1п ^2 = ДЯ20/Я Т20(1 - Т20 /Г) - 1п у2, где v1 - мольная доля 1-го компонента раствора на левой ветви кривой ликвидуса; ^2 - мольная доля 2-го компонента раствора на правой ветви кривой ликвидуса; 1п у:, 1п у2 - логарифмы коэффициентов активности по модели Редлиха - Кистера, рассчитанные по формулам (7).

На рис. 4 можно видеть, что экспериментальные данные [8], полученные для левой и правой ветвей кривых ликвидуса, имеют хорошее согласие с расчетными кривыми.

Массовая доля этиленгликоля

Рис. 4. Зависимость температуры плавления раствора (этиленгликоль - вода) от массовой доли этиленгликоля. кр. 1, 2 - эксперимент [8], кр. 3 и 4 - расчет по модели

Для бинарных эвтектик, входящих в тройную эвтектику №2С03 - МаР - №С1, проведено математическое моделирование параметров кристаллизации. Полученные результаты (табл. 2) согласуются с данными, приведенными в [9].

Относительная ошибка между экспериментальными и расчетными значениями параметров:

ДХ/Х = [(Хэксп - Храсч)/ХэКсп]-100 %.

Т аблица 2

Расчет состава и температуры плавления эвтектики в неорганической бинарной системе

Компоненты системы Состав эвтектики, х, мол. Температура эвтектики, Т, К

Эксперимент [9] Расчет Относит. ошибка Дх/х, % Эксперимент [9] Расчет Относит. ошибка ДТ/Т, %

Ыа2С03 0,610 0,599 1,903 985,00 985,6 -0,06

ЫаР 0,390 0,401 - 2,99

Ыа2С03 0,430 0,431 -0,233 910,6 910,6 0

ЫаС1 0,570 0,569 0,175

ЫаР 0,340 0,341 -0,294 954 948,1 0,6

ЫаС1 0,660 0,659 0,152

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

Т аблица 3

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

Компоненты системы Состав эвтектики, х, мол. Температура эвтектики, Т, К

Экспери- мент Расчет Относит. ошибка Дх/х, % Экспери- мент Расчет Относит. ошибка ДТ/Т, %

Лауриновая 0,7 [3] 0,669 4,43 308,45 307,95 0,16

Миристиновая 0,3 [3] 0,331 -10,33

Лауриновая 0,849 [4] 0,846 0,35 312,31 313,6 -0,41

Стеариновая 0,151 [4] 0,154 - 2,99

Капроновая 0,848 [4] 0,891 -10,22 299,28 301,33 0,68

Пальмитиновая 0,152 [4] 0,109 28,28

Моделирование состава и температуры плавления тройной эвтектики осуществлялось на основе результатов расчета бинарных систем. Расчет состава и температуры кристаллизации тройной эвтектики проводили по алгоритму [9]. Проводится минимизация целевой функции ^(х3э), с этой целью вычисляется такое значение х3э, что Щх3э) < W(x3 < х3э);

Щх3э) = {х13 ехр[(ДЯ10/ДЯ30)1п(хэ/х13)] + х23ехр[(ДЯ20/ДЯ30)1п(хэ/х23)] + х! -1}2, где х,э - мольная доля /-го компонента в тройной эвтектике; х* - мольная доля /-го компонента в бинарной эвтектике, образованной компонентами / и к; ДИ® - молярная теплота плавления чистого /-го компонента.

Рассчитываются мольные доли первого х/ и второго х2э компонентов в тройной эвтектике:

хэ = х/3 (х3 /х33)ЛИ°/ЛИ°, / = 1,2.

По формуле

(Г3)-1 = (Т/к3)-1 + Я/ ДИ/0 1п(х/Ух*), (/, к =1, 2, 3) определяется температура кристаллизации тройной эвтектики Г. Здесь Гк - температура плавления двойной эвтектики, включающей компоненты / и к.

Рис. 5. Экспериментальные данные [3], полученные для левой и правой ветвей кривых ликвидуса в системе (лауриновая кислота - миристино-вая кислота) (кр. 3 и 4) и расчетные кривые (кр. 1 и 2)

Мольная доля 1-го компонента

Рис. 6. Экспериментальные данные [4], полученные для левой и правой ветвей кривых ликвидуса в системе (лауриновая кислота - стеариновая кислота) (кр. 3) и расчетные кривые (кр. 1 и 2)

Рис. 7. Экспериментальные данные [4], полученные для левой и правой ветвей кривых ликвидуса в системе капроновая кислота - пальмитиновая кислота (кр. 3 и 4) и расчетные кривые (кр. 1 и 2)

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

Т аблица 5

Расчет состава и температуры плавления эвтектики в неорганической трехкомпонентной системе

Компоненты системы Состав эвтектики, х, мол. Температура эвтектики, Т, К

Эксперимент [9] Расчет Относит. ошибка A x/x, % Эксперимент [9] Расчет Относит. ошибка AT/T, %

Na2CO3 0,370 0,330 10,8 854 857,572 -0,418

NaF 0,220 0,21В 0,9

NaCl 0,410 0,452 -10,2

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

Т аблица 6

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

Молекулярное соединение Теплота кристаллизации AH, Дж/моль Температура эвтектики Т, К

Экспе- римепт [10] Расчет Отн. ошибка AH/H, % Абс. ошибка, AH Экс- пери- мент Рас- чет Отп. ошибка AT/T, % Абс. ошибка AT, К

Ми6(Па Ст)4 155720 152026,0 2,37 3694 320,6 320,7 0,05 0,1

[Ми6 (Пе Ст)] -2 139900 134056,5 4,18 5843,5 324,4 318,24 1,9 6,16

[Пе3 (Па Ст)] -2 154200 170606,2 10,6 16406,2 325,3 319,62 1,75 5,68

Заключение

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

ЛИТЕРАТУРА

1. Rocha S.A., Guirardello R. An approach to calculate solid-liquid phase equilibrium for binary mixtures // Fluid Phase Equilibria. 2009. V. 281. P. 12-21.

2. Tumakaka F.I., Prikhodko V., Sadowski G. Modeling of solid-liquid equilibria for systems with solid-complex phase formation // Fluid Phase Equilibria. 2007. V. 260. P. 98-104.

3. Costa M.C., KrahenbuhlM.A., Meirelles A.J.A., et al. High pressure solid-liquid equilibria of fatty acids // Fluid Phase Equilibria. 2007. V. 253. P. 118-123.

4. Costa M.C., Rolemberg M.P., Boros L.A.D., et al. Solid-Liquid Equilibrium of Binary Fatty Acid Mixtures// J. Chem. Eng. Data. 2007. V. 52. P. 30-36.

5. Domanska U., Goskowska M. Experimental solid + liquid equilibria and excess molar volumes of alkanol + hexylamine mixtures: Analysis in terms of the ERAS, DISQUAC and Mod. UNIFAC models // Fluid Phase Equilibria. 2004. V. 216. P. 135-145.

6. Khimeche K., Dahmani A. Measurement and prediction of (solid + liquid) equilibria of (al-kanediamine + biphenyl) mixtures // J. Chem. Thermodynamics. 2006. V. 38. P. 1192-1198.

7. Коган В.Б. Гетерогенные равновесия. Л.: Химия, 1968. 432 с.

8. Дымент О.Н., Казанский К.С., Мирошников А.М. Гликоли и другие производные окисей этилена и пропилена. М.: Химия, 1976. 373 с.

9. Мариничев А.Н., Турбович М.Л., Зенкевич И.Г. Физико-химические расчеты на микроЭВМ: Справ. изд. Л.: Химия, 1990. 256 с.

10. Доценко С.П. Марцинковский А.В., Данилин В.Н. // Физико-химический анализ свойств многокомпонентных систем [Электронный ресурс]. 2004. Вып. 1. URL: http//:www. kubstu.ru/th/fams/dopoln3.htm

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

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

Мурашкин Виталий Васильевич

Кемеровский государственный университет

E-mail: ezn2@rambler.ru, markarina@mail.ru, zitner@mail.ru

Поступила в редакцию 17 февраля 2011 г.