Analytical approximation and numerical simulation of liquidus for non-equilibrium multi-component eutectic systems

Full Text

Введение

Целью исследования фазовых состояний химических многокомпонентных систем является не только определение их количественных характеристик, но и в идеале прогнозирование и создание новых композиций с заданными свойствами. Экспериментальное получение диаграмм плавкости двухкомпонентных и поверхностей ликвидусов трехкомпонентных систем представляет собой затратный по времени и трудовым ресурсам процесс. Результаты исследовательской работы химиков-экспериментаторов демонстрируют, что для большинства реальных систем прогнозирование диаграмм плавкости оказывается принципиально невозможным в виду бесконечного разнообразия свойств химических соединений. Тем не менее, как показывают исследования [1–9], существует целый класс химических соединений, для которых расчет, основанный на выявлении аналитических зависимостей между термическими характеристиками, позволяет вычислять характеристики эвтектик и поверхностей ликвидусов как для двухкомпонентных, так и для трехкомпонентных систем. Главным требованием к таким системам является наличие у них единственной точки эвтектики. Это требование, в частности, выполняется для трехкомпонентных систем при том условии, что входящие в их состав бинарные системы являются эвтектическими (с наличием единственной точки эвтектики). Поскольку двухкомпонентных систем, удовлетворяющих этому требованию, огромное количество, то проблема аналитического оценивания и компьютерного моделирования фазовых портретов трехкомпонентных систем остается актуальной.

Построение ликвидусов для двухкомпонентных систем с помощью численной аппроксимации

Рассмотрим вопросы построения графиков кривых ликвидуса для двухкомпонентных систем аналитическими методами. Хорошо известна система уравнений Шредера Ле Шателье:

$\left\{\begin{array}{l}\ln x=-\frac{\Delta H_1}{R}\left(\frac{1}{T_1}-\frac{1}{T}\right),\\ \ln \left(1-x\right)=-\frac{\Delta H_2}{R}\left(\frac{1}{T_2}-\frac{1}{T}\right),\end{array}\right.$ (1)

здесь $\Delta H_1$, $\Delta H_2$, $R$ – энтальпии плавления первого и второго компонентов и универсальная газовая постоянная соответственно; $x_1$, $x_2$ – мольные доли первого и второго компонентов; $T_1$, $T_2$, $T$ – температуры плавления компонентов и искомая температура термодинамического равновесия соответственно.

Система (1) позволяет определять не только мольные доли компонентов и температуру эвтектики, но и строить ветви ликвидуса для двухкомпонентных систем [1–3]:

$\begin{array}{l}T=\frac{\Delta H_1{T}_1}{\Delta H_1-\ln \left(1-x\right)R{T}_1};\\ T=\frac{\Delta H_2{T}_2}{\Delta H_2-\ln \left(x\right)R{T}_2}.\end{array}$ (2)

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

$T_{э}\approx T_1\left(1-\frac{{T_1}^2}{{\left(T_1+T_2\right)}^2}\right)$, $x_1=\frac{T_2-T_{э}}{T_1+T_2-2T_{э}}$, $x_2=1-x_1$. (3)

Здесь $x_1$, $x_2$ – мольные доли первого и второго компонентов; $T_1$, $T_2$, $T_{э}$ – температуры плавления компонентов и искомая температура термодинамического равновесия соответственно. Были произведены расчеты для 25 систем с вычислением соответствующих погрешностей. Для температуры использовалась оценка точности с помощью относительной погрешности

$\delta T=\frac{\left|T_{э}-T\right|}{T}$,

для концентраций использовалась среднеквадратическая, несмещенная оценка погрешности

$\Delta x=\frac{1}{2}\sqrt{{\left(x_1-{\tilde{x}}_1\right)}^2+{\left(x_2-{\tilde{x}}_2\right)}^2+{\left(x_3-{\tilde{x}}_3\right)}^2}$.

Здесь $T$ – вычисленные значения температуры эвтектики; $T_{э}$ – экспериментально полученное значение температуры эвтектики; $x_i$, ${\tilde{x}}_i$ – экспериментальные и вычисленные значения концентраций компонентов системы. Сравнение результатов расчетов с помощью системы (1) и по формулам (3) для 25 систем показали, что среднее арифметическое погрешностей решений для системы Шредера Ле Шателье составляют 5 % по температуре и 6 % по составу, а по формулам (3) – 6 и 7 % соответственно. В связи с тем, что погрешности численной процедуры, предложенной в работе [10], и результатов расчетов с использованием системы (1) оказались сопоставимыми, поскольку их средние значения разнятся на 1 %, была продолжена работа по исследованию зависимости температуры систем от соотношения концентраций ее компонентов.

На рис. 1 видно, что если отразить симметрично относительно ветвей ликвидуса, аппроксимированных отрезками прямых, параболу, проходящую через три точки $\left(\mathrm{0,}{T}_1\right)$, $\left(x_2,T_{э}\right)$, $\left(\mathrm{1,}{T}_2\right)$, известные из расчетов по формулам (3), то верхние ветви этих кривых будут качественно совпадать по поведению с экспериментальными линиями ликвидуса. Понятно, что погрешность такого приближения велика и необходим весовой коэффициент, который позволит приблизить значения этих ветвей к значениям ликвидуса с допустимой погрешностью, не превышающей 10–15 %. Для определения величины или порядка этого коэффициента можно использовать выражения для ветвей ликвидуса (2). Рассмотрим первое из них, считая, что $T_1$ – вещество с меньшей температурой плавления:

$T=\frac{\Delta H_1{T}_1}{\Delta H_1-\ln \left(1-x\right)R{T}_1}$.

 

Рис. 1. Диаграмма плавкости двухкомпонентной системы: 1 – исходная парабола; 2 – симметричная парабола; 3 – парабола с весовым коэффициентом; 4 – отрезок прямой $T_{э}{T}_2$

 

Воспользуемся тем, что для более тугоплавкого компонента $x<0.5$ и $\ln \left(1-x\right)\approx x$, его можно преобразовать к следующему виду:

$T=\frac{T_1}{1-\ln \left(1-x\right)R{T}_1/\Delta H_1}\approx \frac{T_1}{1-xR{T}_1/\Delta H_1}$.

Учитывая также, что для большинства химических соединений $\frac{R{T}_1}{\Delta H_1}\le 1$, можно записать дробь в правой части в виде ряда

$T\approx T_1\left[1+x\frac{R{T}_1}{\Delta H_1}+x^2{\left(\frac{R{T}_1}{\Delta H_1}\right)}^2\left[1+x\frac{R{T}_1}{\Delta H_1}+x^2{\left(\frac{R{T}_1}{\Delta H_1}\right)}^2\right]+o\left(x^5{\left(\frac{R{T}_1}{\Delta H_1}\right)}^5\right)\right]$. (4)

Аналогичную процедуру можно произвести с другой ветвью. В результате вычитания из (4) его линейной части $\tilde{T}\approx T_1\left(1+x\frac{R{T}_1}{\Delta H_1}\right)$ получаем оценку разности между линейной и квадратичной аппроксимациями этой ветви ликвидуса:

$T-\tilde{T}\approx T_1\left[x^2{\left(\frac{R{T}_1}{\Delta H_1}\right)}^2\left(1+x\frac{R{T}_1}{\Delta H_1}+x^2{\left(\frac{R{T}_1}{\Delta H_1}\right)}^2\right)\right]$.

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

$x^2{T}_1{\left(\frac{R{T}_1}{\Delta H_1}\right)}^2\le \mathrm{0,25}{\left(\frac{R{T}_1}{\Delta H_1}\right)}^2{T}_1\le \mathrm{0,25}{T}_1$.

Исходя из предположения, что $\max T_1<3500$°С и $\left(T-\tilde{T}\right)\le \left(T_1-T_{э}\right)<1000$°С, коэффициент, на который следует умножать эту разность и который обеспечит достижение максимальной относительной погрешности значений разности

$\frac{T-\tilde{T}}{T_1-T_{э}}\le \mathrm{0,06}$,

должен иметь порядок величины $x_2^2$ – квадрата мольной доли тугоплавкого компонента в системе. Запишем формулы, соответствующие параметрам полиномов первой и второй степени, проходящих через три известные точки $\left(\mathrm{0,}{T}_1\right)$, $\left(x_2,T_{э}\right)$, $\left(\mathrm{1,}{T}_2\right)$:

$k_1=\frac{T_1-T_{э}}{x_2}; b_1=T_1; T_1=b_1-k_{1}x; a=-\frac{T_{э}-\left(T_2-T_1\right)x_2-T_1}{x_1{x}_2}$;

$\tilde{T}=a{x}^2+bx+c; b=T_2-T_1-a; c=T_1; k_2=\frac{T_2-T_{э}}{x_1}$;

 $b_2=T_{э}; T_2=b_2+k_2\left(x-x_2\right)$.

Тогда ветви ликвидуса будут выражаться следующими формулами зависимости температуры от концентрации веществ:

$\begin{array}{l}T=T_1+\left(T_1-\tilde{T}\right){\left(x_2\right)}^2,0\le x\le x_2;\\ T=T_2+\left(T_2-\tilde{T}\right){\left(x_2\right)}^2,\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{x}_2<x\le 1.\end{array}$(5)

Для 8 систем проводилось сравнение полученных с помощью данной процедуры аппроксимаций ликвидусов с кривыми ликвидусов, построенными с помощью системы уравнений Шредера Ле Шателье. Результаты расчетов систем показали довольно существенные расхождения в значениях температур и соответствующих концентраций. Максимальное расхождение в значениях графиков кривых было выявлено для системы CsCl – CsF и составило 157 °С по температуре и 0,11 по мольным долям компонентов эвтектики.

Для систем Ba(ClO₄)₂–KClO₄ и BaF₂–BaSO₄ представлены диаграммы плавкости, полученные экспериментально [11, 12] и с помощью аналитических расчетов по формулам (5) и уравнений Шредера Ле Шателье (табл. 1 и 2).

 

Таблица 1. Расчет систем Ba(ClO₄)₂–KClO₄ по формулам (5) и уравнениям Шредера Ле Шателье

	$\Delta H$ (КДж/моль)	$\Delta H/R$	$T$ плавления
Ba(ClO4)2	9,87	1187,00	473
KClO4	10,77	1295,49	570
$T_{э}$ – уравнение Шредера	$x_1$ – уравнение Шредера	$T_{э}$ – формулы (5)	$x_1$ – формулы (5)
393	0,71	375	0,67

 

 Таблица 2. Расчет систем BaF₂–BaSO₄ по формулам (5) и уравнениям Шредера Ле Шателье

	$\Delta H$ (КДж/моль)	$\Delta H/R$	$T$ плавления
BaF2	17,50	2149,88	1332
BaSO4	41,00	5036,86	1580
$T_{э}$ – уравнение Шредера	$x_1$ – уравнение Шредера	$T_{э}$ – формулы (5)	$x_1$ – формулы (5)
1117	0,73	1053	0,65

 

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

 

Рис. 2. Максимальные отклонения от экспериментальных данных (кривая 1) системы Ba(ClO₄)₂–KClO₄ расчетных кривых по уравнению Шредера (кривая 2) – 28,8 °С и формулам (5) (кривая 3) – 37,5 °С

 

Рис. 3. Максимальные отклонения от экспериментальных данных (кривая 1) системы BaF₂–BaSO₄ расчетных кривых по уравнению Шредера (кривая 2) – 229 °С и формулам (5) (кривая 3) – 107 °С

 

Построение изотерм для трехкомпонентных систем с помощью линейной аппроксимации

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

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

 

Рис. 4. Линейная аппроксимация эвтектики двухкомпонентной системы KCl–KBO₂

 

Введем обозначение для концентрации эвтектики данной двухкомпонентной системы: KBO₂– $x_1$ и для KCl– $\left(1-x_1\right)$ соответственно. Температура состава KBO₂–KCl в зависимости от изменения концентрации KBO₂– $x$ будет вычисляться по формулам

$T=T_{э}+\frac{\left(T_1-T_{э}\right)x}{x_1}$(6)

или

$T=T_1-\frac{\left(T_1-T_{э}\right)x}{x_1}$.

На участке, расположенном правее точки эвтектики, концентрация компонента KBO₂ будет равна $\left(x_1+x\right)$, KCl соответственно $(1-x_1-x)$, поэтому температура состава будет вычисляться по формулам

$T=T_{э}+\frac{\left(T_2-T_{э}\right)x}{1-x_1}$(7)

или

$T=T_2-\frac{\left(T_2-T_{э}\right)\left(1-x_1-x\right)}{1-x_1}$.

Таким же образом с использованием зависимостей (6) и (7) между температурами эвтектик и концентрациями можно приближенно вычислять значения промежуточных температур на сторонах треугольника Гиббса и вдоль прямых, соединяющих вершины треугольника с точками (в том числе и эвтектическими) на противоположных сторонах. В результате оказывается возможным программно рассчитывать значения температур в любых точках внутри треугольника и строить на их основе изотермы в качестве отрезков прямых, соединяющих точки на сторонах треугольника и проходящих через секущие, связывающие вершины с точками эвтектик на противоположных сторонах. На основе предложенной методики линейной аппроксимации была разработана программа в визуальной среде программирования Delphi, позволяющая строить приближенную поверхность кристаллизации в треугольнике Гиббса.

В данной программе для расчетов используются только значения температур плавления, входящих в систему веществ. Промежуточные вычисления, однако, требуют знания значений составов концентраций и температур эвтектик бинарных систем, входящих в качестве компонентов низшей размерности в исследуемую трехкомпонентную систему, а также состава и температуры ее эвтектики. Процедура вычисления вспомогательных значений составов и температур эвтектик бинарных систем в предлагаемой программе производится по формулам из статьи авторов [10]. Так, вычисление температуры эвтектики и концентрации тугоплавкого компонента для двухкомпонентной системы (KBO₂ на рис. 4) осуществляется по формулам

$x_1=\frac{T_2}{T_1+T_2}, x_2=1-x_1$;

$T_{э}\approx T_1{x}_1+T_1{x}_1\left(1-x_1\right)=T_1{x}_1\left(1+x_2\right)=T_1\left(1-x_2^2\right)$,

уточненное значение концентрации легкоплавкого компонента

$x_1=\frac{T_2-T_{э}}{T_1+T_2-2T_{э}}$.

Здесь $T_1, T_2, T_{э}$ – соответственно температуры плавления 1-го и 2-го компонентов и температура эвтектики; $x_1, x_2$ – содержание 1-го и 2-го компонентов в мольных долях.

Для вычисления температур эвтектики и ее состава для трехкомпонентных систем использовались формулы из [18]:

$T_{э}^{1-2-3}=\frac{T_{э}^{1-2}\left(T_{э}^{1-3}+T_{э}^{2-3}+T_3\right)}{T_1+T_2+T_3}$;

$c_3=\frac{T_{э}^{1-2}-T_{э}^{1-2-3}}{T_{э}^{1-2}+T_3-2T_{э}^{1-2-3}}$.

Концентрации двух других веществ определяются по формулам

$c_2=x_2^{1-2}-x_2^{1-2}{c}_3, c_1=x_1^{1-2}-x_1^{1-2}{c}_3$,

где $T_i, T_{э}^{i-k}, (i,k=\mathrm{1,2,3})$ – температуры плавления компонентов и температуры эвтектик в двухкомпонентных системах соответственно; $T_{э}^{1-2-3}$ – температура эвтектики в трехкомпонентной системе; $x_i^{i-k}, c_i, i, k=\mathrm{1,2,3}$ – содержание i-го компонента в мольных долях в двухкомпонентных и трехкомпонентной системах.

На рис. 5 представлены поверхности ликвидуса из справочника [19], построенные с помощью вычисленных программой значений. Экспериментально полученные значения температуры эвтектики трехкомпонентной системы и концентрации входящих в нее веществ равны соответственно: 640 ºС, KCl – 56 %, KBO₂ – 30 %, K₂SO₄ – 14 %. На рис. 5, б представлены вычисленные значения мольных долей компонентов, нижняя, верхняя и усредненная оценка температуры эвтектики и изотермы, построенные с помощью разработанной программы.

 

Рис. 5. Изотермы: а – из справочника [14]; б – построенные с помощью линейной аппроксимации (погрешность составляет по температуре 1 % и по составу 1 %)

 

Для температуры использовалась оценка точности с помощью относительной погрешности

$\delta T=\frac{\left|\frac{1}{2}\left({\tilde{T}}_d+{\tilde{T}}_u\right)-T\right|}{T}$, (8)

для концентраций использовалась среднеквадратическая, несмещенная оценка погрешности

$\Delta x=\frac{1}{2}\sqrt{{(x_1-{\tilde{x}}_1)}^2+{(x_2-{\tilde{x}}_2)}^2+{(x_3-{\tilde{x}}_3)}^2}$. (9)

Здесь ${\tilde{T}}_d$ – нижняя и ${\tilde{T}}_u$ – верхняя оценки, вычисленные по программе значения температуры эвтектики; $T$ – экспериментально полученное значение температуры эвтектики; $x_i, {\tilde{x}}_i$ – экспериментальные и вычисленные по программе значения концентраций компонентов системы.

В качестве характерного примера на рис. 6 представлены экспериментально построенные изотермы и поверхность кристаллизации, рассчитанная с помощью программы для системы LiBO₂/NaBO₂/Mg(BO₂)₂.

 

Рис. 6. Изотермы: а – из справочника [14]; б – построенные с помощью линейной аппроксимации (погрешность по температуре 3 %, по составу 3 %)

 

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

В таблице 3 представлены погрешности результатов расчетов изотерм для 30 трехкомпонентных систем в сравнении с данными из [19, 20].

 

Таблица 3. Погрешности результатов расчет трехкомпонентных систем

Название веществ системы	Погрешность		Название веществ системы	Погрешность
	по температуре, %	по составу, %		по температуре, %	по составу, %
K//SO4/Cl/BO2	1	1	Ce/Na/Sn//Cl	1,2	5
Li//F/Cl/CO3	2	9	Ce/Fe/Sn//Cl	0,1	4
Na/K/Li//F	8	4	Ca/La/Na//Cl	8	7
Ba/Na/Cs//F	0,1	6	Li/Mg/Sr//F	6	14
Ca/K/Na//NO3	0,1	7	Na//CNS/SO4/Cl	5	7
Ba/Na/Ca//F	4	5	Na//CO3/C/F	1	4
Li/Na/Sr//F	3	5	Na/Cl/P2O7/SO4	4	4
Li/Mg/Na//BO2	3	3	Na//Cl/F/P2O7	5	7
K/Na/Sr//F	1	12	K//PO4/P2O7/SO4	1	9
Cs/K/Na//NO3	1	7	Ba/Ca/Na//F	7	6
Cs/K/Na//Br	4	1	H2O-Nh4NO3-KNO3	5	6
Ca/Li/Na//F	5	7	H2O-Nh4NO3-NaNO3	5	7
Ca/Li/Mg//F	2	2	H2O-NaNO3-KNO3	7	3
K/Li/Sr//F	7	2	CO(NH2)2-NaNO3-KNO3	4	6
Li4/K4/Na4//P2O7	8	8	Na//Cl/SO4/BO2	11	2
Средние погрешности				5	6

 

Результаты сравнения поверхностей ликвидусов для 30 эвтектических систем из справочника [19] и монографии [20] показали, что принципиального, качественного расхождения в картинах расположения изотерм не наблюдается. Средняя погрешность по концентрациям составила 6 %, по температуре 5 %. Для 30 исследованных авторами поверхностей ликвидусов погрешности по температуре не превысили 11 %, среди погрешностей по концентрациям для одной из систем она составила 14 %.

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

About the authors

Olga S. Afanaseva

Samara State Technical University

Email: afa_samara@inbox.ru

PhD. Sci. (Techn.), Associate Professor

Russian Federation, 244, Molodogvardeyskaya st., Samara, 443100

Elena A. Afanaseva

Samara State Technical University

Email: afa_samara@inbox.ru

Postgraduate Student

Russian Federation, 244, Molodogvardeyskaya st., Samara, 443100

Galina F. Egorova

Samara State Technical University

Author for correspondence.
Email: afa_samara@inbox.ru

PhD. Sci. (Techn.), Associate Professor

Russian Federation, 244, Molodogvardeyskaya st., Samara, 443100

References

Supplementary files