Mathematical simulation of phase change fluid - solid..pdf Математическое моделирование фазового равновесия жидкость - твердое длябинарных и многокомпонентных смесей имеет значение для многих областейпромышленности, таких, как химическая, пищевая, фармацевтическая и др. [1- 4].Основное применение методы математического моделирования находят при ре-шении задач выбора составов и расчете термодинамических параметров антифри-зов, теплоаккумуляторов, теплозащитных покрытий, моделирования состава неф-ти, ректификации с использованием разделяющих агентов. Изучение равновесно-го состояния важно для всех процессов, где существует взаимодействие междутвердым телом и жидкой фазой. Для таких процессов термодинамические свойст-ва иногда трудно измерить экспериментально, особенно при экстремальных усло-виях (например, высокое давление, низкая температура и т.д.), и поэтому необхо-димо математическое моделирование. Разработка теоретических методов, позво-ляющих выбрать оптимальный состав антифризов и теплохладоаккумуляторов стребуемой температурой и теплотой плавления, является актуальной задачей,представляющей научный и практический интерес.Равновесие жидкость - твердое может быть рассчитано многими способами[5, 6]. Термодинамические модели для рассмотрения фазового равновесия раство-ров могут быть классифицированы на две категории: модели коэффициента ак-тивности и модели уравнения состояния. Термин «прогнозирующие термодина-мические модели» характеризует типы моделей, которые могут описывать фазо-вое равновесие при условии, что молекулярные структуры или физические свой-ства чистых компонентов в смеси известны. Прогнозирующие термодинамиче-ские модели также можно классифицировать на две категории: модели с исполь-зованием экспериментальных данных, например модель «UNIFAC» [5], в которойпараметры взаимодействия найдены из экспериментальных данных и прогнозтермодинамических свойств делается на основе существующих параметров дляатомных групповых вкладов молекул, присутствующих в смеси; и модели без ис-пользования экспериментальных данных (или априори прогнозирующие модели),например модель «COSMO-RS», в которой требуются только специфическиеатомные параметры и прогноз термодинамических свойств делается на основеквантовых химических вычислений, что дает необходимую информацию дляоценки молекулярных взаимодействий в жидкостях. Кривые фазового равновесиямогут дать информацию о характеристиках чистых компонентов и смеси. Досто-верное определение фазового равновесия необходимо для корректного прогнозаобразования фаз и их композиций.1. Минимизация избыточной свободной энергии ГиббсаПонятие термодинамического потенциала было введено Гиббсом по аналогиис потенциалом в механике. Термодинамический потенциал GEназван потенциа-лом Гиббса, или изобарно-изотермическим потенциалом. В данной работе для по-строения математической модели применяется минимизация свободной энергииГиббса по параметру сольватации ƒ, введенному как отношение числа молекулкомпонента А к числу молекул компонента B в молекулярном соединении, обра-зующемся в растворе. Модель, в которой используется минимизация по парамет-ру сольватации, может быть отнесена к прогнозирующей модели уравнения со-стояния, так как вначале предлагается модель уравнения состояния реальной сис-темы, зависящая от параметра сольватации ƒ, а также потому, что она не требуетввода данных о взаимодействии в растворе, а основана только на эксперимен-тальных данных о температуре и энтальпии плавления чистых компонентов.Представленный метод вычисления равновесия твердое тело - жидкая фаза длябинарных систем позволяет получить выражение для температуры как функциимольной доли компонента смеси. Модель также может применяться для описанияравновесия жидкость - пар и для систем твердое - твердое с образованием соеди-нений в твердой фазе. Этот метод используется при расчете фазовых равновесий втройных и многокомпонентных реальных системах, сопровождающихся сольва-тацией молекул в растворе и ассоциацией чистых компонентов. Преимуществомописываемого метода моделирования фазового равновесия является использова-ние малого числа параметров, необходимых для расчета, для идеальных систем -это температура и энтальпия плавления чистых компонентов.Избыточная свободная энергия Гиббса GE является разностью между свобод-ной энергией Гиббса реальной и идеальной систем. Для идеальных систем, где от-сутствует взаимодействие между молекулами, избыточная свободная энергияГиббса равна нулю. Таким образом, избыточная свободная энергия Гиббса явля-ется мерой взаимодействия в системе, причем она может быть положительнойили отрицательной, что связано с видом взаимодействия.Если за параметры состояния принимаются мольные доли компонентов в сме-си xi, температура T и давление P, то уравнение состояния фазы многофазной,многокомпонентной системы имеет вид [7]dGE = -S dT + V dP + ƒƒi dxi ,где S - энтропия, V - объем, ƒi - химический потенциал компонента.Условием равновесия системы является минимум внутренней энергии. Из это-го условия следует, что dU = 0. Если температура и давление постоянны, то дляравновесия системы, представляющей собой одну фазу, должны выполняться ра-венстваdT = 0, dP = 0, dG = 0. (1)В случае двухфазной системы для равновесия необходимо, чтобы были равныхимические потенциалы фаз ƒ = ƒ и выполнялись условия (1). Фазовая диаграм-ма жидкость - твердое бинарной системы может быть получена минимизациейсвободной энергии Гиббса. Избыточная свободная энергия Гиббса бинарной сис-темы дается выражением [7]ƒGE = RT ƒxi ln ƒi ,где ƒi - коэффициент активности компонентa, i = 1,2; xi - мольная доля i-го ком-понента.Разность уравнений состояния бинарной системы, для реальной и идеальнойравновесных фаз, можно представить в виде [7](-ƒHE/RT 2)dT + (ƒVE/RT)dP = Σxi dln ƒi ,где ƒHE - энтальпия смешения; ƒVE - избыточный объем; P - давление раствора;T - абсолютная температура; R - универсальная газовая постоянная.Если происходит образование соединения чистых компонентов в растворе, тосуммарная молярная масса i-го компонента с учетом числа молекул данного видаƒi, входящих в соединение, может быть рассчитана по формуле: ƒi = ƒi ƒi,где ƒi - молярная масса чистого компонента. Среднее соотношение числа молекулв сольватах чистых компонентов ƒ = ƒ1/ƒ2 характеризует устойчивую структурураствора. С учетом изменения молярной массы эффективные мольные доликомпонентов бинарной смеси вычисляются по формуле [7] z1 = x1/(x1+ ƒx2),z2 = x2/(x1/ƒ + x2). В том случае, если из эксперимента известны данные о составераствора x1э, x2э в точке экстремума температуры, можно найти отношение коэф-фициентов ƒ = ƒ1/ƒ2 = (z2э/z1э)/(x2э/x1э). Одним из методов определения среднего па-раметра сольватации компонентов в растворе является вычисление ƒср по формулесредней арифметической взвешенной, где в роли весов выступают концентрацииодного из компонентов в точках возможных эвтектик:ƒср = ƒ ƒ j xэ j.Растворимость компонента, образующего однокомпонентную фазу в много-фазной многокомпонентной смеси конденсированных сред при постоянном дав-лении описывается уравнениемd(ln xi ƒi)/dT = ƒHi /RT 2, (2)где ƒHi - парциальная молярная теплота растворения i-го компонента в насыщен-ном им растворе.Теплота перехода компонента, образующего твердую фазу, в жидкое состоя-ние ƒHi складывается из теплоты плавления ƒHiпл и дифференциальной теплотыразбавления компонента ƒHiразб до концентрации, соответствующей насыщенномураствору [7]. Теплота разбавления зависит от характера взаимодействия компо-нент в растворе и может быть как положительной, так и отрицательной. Посколь-ку теплота разбавления для большинства систем неизвестна, то расчет по уравне-нию (2) возможен лишь для идеальных систем, для которых ƒi = 1 и ƒHi = ƒHiпл:d( ln xi ƒi)/dT = ƒHiпл/RT 2.Теплота плавления ƒHiпл зависит от температуры:ƒHiпл = ƒHi0+ (сiж - сiт)(T - Ti0),где ƒHi0- теплота плавления чистого i-го компонента при температуре Ti0;сiж − сiт - молярные теплоемкости i-го компонента при постоянном давлении вжидком и твердом состояниях.Отличие коэффициента ƒ от единицы свидетельствует о наличии отклоненияот идеальности в бинарной системе и необходимости перехода к эффективныммольным долям. Для реальных систем, в которых возможно образование молеку-лярных соединений, называемых сольватами, необходимо в уравнении (2) сделатьпереход к эффективным мольным долям:d (ln zi ƒi)/dT = ƒHi /RT2,где ƒHi = ƒHiпл + ƒHiразб = ƒHi0+(сiж − сiт)(T - Ti0) + ƒHiразб; ƒHi0- энтальпия плав-ления и Ti0 - температура плавления компонента, образующего однокомпонент-ную фазу.В том случае, если теплоемкости в жидком и твердом фазовом состояниях раз-личаются незначительно и диапазон температур плавления невелик, пренебрегаятеплотой разбавления, парциальные коэффициенты активности можно предста-вить в виде ln ƒi = ƒHi0(1 - Ti0 /T)/(R Ti0) - ln zi + ln Fi (zi), где Fi (zi) - произвольныефункции от состава. При постоянном давлении и малом интервале температурплавления избыточную энергию Гиббса как функцию от эффективной мольнойдоли компонента находим из уравненияƒGE = z1[ƒH10(T/ T10- 1)] + z2[ƒH20(T/T10- 1)] - RT(z1 ln z1 + z2 ln z 2) + F(z1). (3)Здесь F(z1) = z1F1(z1) + (1−z1)F2(z1) - функция, которая выбирается из условия тер-модинамической согласованности модели по методу Херингтона и Редлиха -Кистера [7]:11 20lg  γ / γ dz1 = 0.Согласно этому уравнению, должны быть равны площади, ограниченные кри-вой lg ƒ1/ƒ2 и осями координат.Минимизация избыточной энергии (3) по внутреннему параметру ƒ приводит куравнению Бернулли:dT/dz1 + f1(z1)T = f2(z1)T2, (4)где f1(z1) = ƒ/(ƒH10/R - ƒ z1), f2(z1) = -[ƒ + ln (1-z1)/z1]/(ƒH10/R - ƒ z1),ƒ = (ƒH10 - ƒH20)/R, ƒ = ƒH20/RT20- ƒH10/R T10.Решение уравнения (4) имеет видT(z1) = [ƒH10z1 + ƒH20(1−z1)]//{ƒH10z1/T10 + ƒH20(1−z1)/T20 - R[z1 ln z1 + (1−z1) ln(1-z1)]}. (5)Таким образом, минимизация избыточной энергии Гиббса ƒGE по параметрусольватации ƒ позволяет найти уравнение, моделирующее Tz-диаграмму фазовогоравновесия. Эта диаграмма представляет собой сечение при P = const трехмернойPTz-диаграммы состояния бинарной системы.Из условия экстремума функции (5) получено алгебраическое уравнение( ) ( ) ( ) 10 20 0 0 0 0 21 1 1 2 2 1ln 1 z H/R z H/R H H 1/T 1/T R . ⎛⎜ − ƒ ⋅ −ƒ ⎞⎟= ƒ ƒ −⎝ ⎠(6)Уравнения (5), (6) позволяют найти температуру плавления эвтектики и составраствора в точках экстремума температуры в случае их существования. Оставляяв разложении в степенной ряд функции в левой части (6) члены второго порядка,получим алгебраическое уравнение второй степени. Решение уравнения имеетвидzH HH H H HH H R T T H⎛ ƒ ⎞ ⎛ ƒ ⎞=⎜ + ⎟ ⎜ − ⎟⎝ ƒ ⎠ ⎝ ƒ ⎠⎛ ƒ ⎞ ⎛ ƒ ⎞ ⎛ ƒ ⎛ ⎞ ⎞ ⎛ ƒ ⎞ ⎜ + ⎟ ⎜ − ⎟ +⎜ ⎜ − ⎟− ⎟ ⎜ − ⎟⎝ ƒ ⎠ ⎝ ƒ ⎠ ⎝ ⎝ ⎠ ⎠ ⎝ ƒ ⎠(7)Мольная доля компонента в растворе не превышает единицы, поэтому выби-раем решение z(0,1). Подставляя z1 = z1э в функцию (5), можно найти температу-ру плавления в эвтектической точке.В общем случае, если известна парциальная молярная теплота растворения i-гокомпонента в насыщенном им растворе ƒHi , то в уравнениях (3) - (7) следует за-менить энтальпию плавления чистого компонента ƒHi0 на ƒHi.2. Численные результатыВ табл. 1 приведены результаты расчета характеристик эвтектических систем,включающих гликоли и воду.Т а б л и ц а 1Расчет состава, температуры плавления эвтектикии параметров сольватации и ассоциации в растворах гликолейКомпонентысистемыСоставэвтектикиx1э, x2э массТемператураэвтектикиТ, °СЭнтальпияплавления в точкеэвтектики ƒHплэ,кДж/мольКоэффи-циентсольва-тации ƒКоэффи-циентассоциа-ции kЭтиленгликоль 0,6977Вода 0,3023 -62,7 8,306 0,97 0,иэтиленгликоль 0,7306Вода 0,2694 -57,22 8,909 0,96 0,875Триэтиленгликоль 0,6871Вода 0,3129 -48,81 9,897 0,94 0,75Этиленгликоль 0,4364Диэтиленгликоль 0,5636 -38,64 12,455 1 1,037Триэтиленгликоль 0,4688Диэтиленгликоль 0,5312 -29,87 15,475 1 0,917На рис. 1 приведена кривая, моделирующая среднюю температуру плавлениясистемы (этиленгликоль - вода); точками, соединенными ломаными линиями, пока-заны возможные в данной системе эвтектики. На рис. 2 приведены кривые lg (ƒ1/ƒ2)для системы (этиленгликоль - вода), полученные в результате предварительногомоделирования и проведения процедуры термодинамического согласования. Нарис. 3 приведены кривые ln ƒ1, ln ƒ2 для системы (этиленгликоль - вода), получен-ные в результате проведения процедуры термодинамического согласования.Из условия термодинамической согласованности построенной модели опреде-ляется коэффициент ассоциации k = k1 / k2, имеющий смысл отношения числа мо-лекул компонента А, объединившихся в кластер до образования раствора или не-посредственно в данном растворе, к числу молекул в кластере компонента В. Этопозволяет построить кривые ликвидуса, близкие к экспериментально найденным(рис. 3).Т, К12Рис. 1. Средняя температура плавления раствора (этиленгликоль - вода);точками показаны предполагаемые эвтектики0 0,2 0,4 0,6 0,8 1Мольная доля 1-го компонента210-1-2lg ƒ1/ƒ212Рис. 2. Зависимость lg ƒ1 /ƒ2 от мольной доли первого компонентав бинарной смеси (этиленгликоль - вода); кр. 1 - расчет по моде-ли; кр. 2 - результат термодинамического согласования0 0,2 0,4 0,6 0,8 1Мольная доля 1-го компонентаln ƒ ƒ 1 2 , ln 012341 2Рис. 3. Зависимости ln ƒ1, ln ƒ2 от мольной доли первого компонентав бинарной смеси (этиленгликоль - вода)В том случае, если кривые ln ƒ1, ln ƒ2 удовлетворяют условию термодинамиче-ской согласованности, представим их в виде ряда Редлиха - Кистераln ƒ1 = z22[B + C(3z1 - z2)], ln ƒ2 = z12[B + C(z1 - 3z2)] (8)или аппроксимируем их уравнениями Ван-Лаараln ƒ1 = A[Bz2 /(Az1 + Bz2)]2, ln ƒ2 = B[Az1/(Az1+Bz2)]2.По эвтектическим данным: z1э, z2э, ln ƒ1э , ln ƒ2э находим коэффициенты моде-лей Редлиха - Кистера и Ван-Лаара.Используя полученные выражения для логарифмов коэффициентов активно-сти, находим функции (рис. 4 - 7), аппроксимирующие экспериментально опреде-ленные кривые ликвидуса для систем: (этиленгликоль - вода) и насыщенныхжирных кислот.Найдем состав раствора на ветвях ликвидуса в зависимости от температуры поформулам:ln ƒ1 = ƒH10/R T10 (1 - T 10 /T) - ln ƒ1 , ln w2 = ƒH20/R T20(1 - T20 /T) - ln ƒ2,где ƒ1 - мольная доля 1-го компонента раствора на левой ветви кривой ликвидуса;w2 - мольная доля 2-го компонента раствора на правой ветви кривой ликвидуса;ln ƒ1, ln ƒ2 - логарифмы коэффициентов активности по модели Редлиха - Кистера,рассчитанные по формулам (7).На рис. 4 можно видеть, что экспериментальные данные [8], полученные длялевой и правой ветвей кривых ликвидуса, имеют хорошее согласие с расчетнымикривыми.0 0,2 0,4 0,6 0,8 1Массовая доля этиленгликоляТ 1*, К2102202302402502601234Рис. 4. Зависимость температуры плавления раствора(этиленгликоль - вода) от массовой доли этиленгликоля.кр. 1, 2 - эксперимент [8], кр. 3 и 4 - расчет по моделиДля бинарных эвтектик, входящих в тройную эвтектику Na2CO3 - NaF - NaCl,проведено математическое моделирование параметров кристаллизации. Получен-ные результаты (табл. 2) согласуются с данными, приведенными в [9].Относительная ошибка между экспериментальными и расчетными значениямипараметров:ƒX/X = [(Xэксп - Xpасч)/Xэксп]·100 %.Т а б л и ц а 2Расчет состава и температуры плавления эвтектикив неорганической бинарной системеСостав эвтектики, х, мол. Температура эвтектики, Т, ККомпонентысистемы Экспери-мент [9] РасчетОтносит.ошибкаƒx/x, %Экспери-мент [9] РасчетОтносит.ошибкаƒT/T, %Na2CO3 0,610 0,599 1,903NaF 0,390 0,401 - 2,99 985,00 985,6 -0,06Na2CO3 0,430 0,431 -0,233NaCl 0,570 0,569 0,175 910,6 910,6 0NaF 0,340 0,341 -0,294NaCl 0,660 0,659 0,152 954 948,1 0,6В табл. 3 приведены экспериментальные и расчетные параметры эвтектиче-ских точек в системах насыщенных жирных кислот.Т а б л и ц а 3Экспериментальные и расчетные значения состава и температуры плавленияэвтектики в бинарных системах насыщенных жирных кислотСостав эвтектики, х, мол. Температура эвтектики, Т, ККомпонентысистемы Экспери-мент РасчетОтносит.ошибкаƒx/x, %Экспери-мент РасчетОтносит.ошибкаƒT/T, %Лауриновая 0,7 [3] 0,669 4,43Миристиновая 0,3 [3] 0,331 -10,33 308,45 307,95 0,16Лауриновая 0,849 [4] 0,846 0,35Стеариновая 0,151 [4] 0,154 - 2,99 312,31 313,6 -0,41Капроновая 0,848 [4] 0,891 -10,22Пальмитиновая 0,152 [4] 0,109 28,28 299,28 301,33 0,68Моделирование состава и температуры плавления тройной эвтектики осуще-ствлялось на основе результатов расчета бинарных систем. Расчет состава и тем-пературы кристаллизации тройной эвтектики проводили по алгоритму [9]. Прово-дится минимизация целевой функции W(x3э), c этой целью вычисляется такое зна-чение x3э, что W(x3э) < W(x3 < x3э);W(x3э) = {x113 exp[(ƒH10/ƒH30)ln(x3 э /x313)] + x323exp[(ƒH20/ƒH30)ln(x3 э /x323)] + x3 э -1}2,где xiэ - мольная доля i-го компонента в тройной эвтектике; xiik - мольная доля i-гокомпонента в бинарной эвтектике, образованной компонентами i и k; ƒHi0 - мо-лярная теплота плавления чистого i-го компонента.Рассчитываются мольные доли первого x1э и второго x2э компонентов в трой-ной эвтектике:( ) 0 0э 3 э 3 33 3// , 1, 2. i i Hi Hxi xi x x i ƒ ƒ= =По формуле(Tэ)-1 = (Tikэ)-1 + R/ ƒHi0 ln(xiэ/xiik), (i, k =1, 2, 3)определяется температура кристаллизации тройной эвтектики Tэ. Здесь Tikэ - тем-пература плавления двойной эвтектики, включающей компоненты i и k.Т, К13 24Рис. 5. Экспериментальные данные [3], полученные для левой и правойветвей кривых ликвидуса в системе (лауриновая кислота - миристино-вая кислота) (кр. 3 и 4) и расчетные кривые (кр. 1 и 2)315320325330335Т, К0,1 0,3 0,5 0,7 0,9Мольная доля 1-го компонента1234Рис. 6. Экспериментальные данные [4], полученные для левой и правойветвей кривых ликвидуса в системе (лауриновая кислота - стеариноваякислота) (кр. 3) и расчетные кривые (кр. 1 и 2)300310320330Т, К0,1 0,3 0,5 0,7 0,9Мольная доля 1-го компонента13 24Рис. 7. Экспериментальные данные [4], полученные для левой и правойветвей кривых ликвидуса в системе капроновая кислота - пальмитино-вая кислота (кр. 3 и 4) и расчетные кривые (кр. 1 и 2)Модель, основанная на опытных данных о температуре и теплоте плавлениячистых компонентов, входящих в тройную систему, может использоваться припрогнозировании температур плавления эвтектик трехкомпонентных систем. Рас-считанные по данной методике параметры изучаемых эвтектических систем со-гласуются с результатами экспериментально определенных состава и температу-ры в точках эвтектики трехкомпонентных систем (табл. 5).Т а б л и ц а 5Расчет состава и температуры плавления эвтектикив неорганической трехкомпонентной системеСостав эвтектики, х, мол. Температура эвтектики, Т, ККомпонентысистемы Экспери-мент [9] РасчетОтносит.ошибкаƒ x/x, %Экспери-мент [9] РасчетОтносит.ошибкаƒT/T, %Na2CO3 0,370 0,330 10,8NaF 0,220 0,218 0,9NaCl 0,410 0,452 -10,2854 857,572 −0,418В табл. 6 приведены результаты расчета эвтектических параметров в тройнойсистеме насыщенных жирных кислот.Т а б л и ц а 6Экспериментальные и расчетные значения теплоты и температуры плавленияэвтектики в трехкомпонентной системе насыщенных жирных кислотТеплота кристаллизации ƒH, Дж/моль Температура эвтектики Т, КМолекулярноесоединениеЭкспе-римент[10]РасчеОтн.ошибкаƒH/H, %Абс.ошибка,ƒHЭкс-пери-ментРас-четОтн.ошибкаƒT/T,%Абс.ошибкаƒT, КМи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[Пе 1*©3 (Па Ст)]·2 154200 170606,2 10,6 16406,2 325,3 319,62 1,75 5,68ЗаключениеВ заключение отметим, что предлагаемая математическая модель может ис-пользоваться для построения кривых ликвидуса бинарных систем жидкость -твердое, при прогнозировании составов и температур плавления эвтектик бинар-ных и многокомпонентных систем, а также для уточнения формул возможныхмолекулярных соединений в изучаемой системе.

Скачать электронную версию публикации