1. На главную
  2. Электронные версии
  3. Журнал физической химии
  4. T. 93, № 6, 2019

Журнал физической химии, 2019, T. 93, № 6, стр. 810-817

Об определении фазового состава сложных термодинамических систем

Г. В. Белов a, b, c, *

a Московский государственный университет имени М.В. Ломоносова, Химический факультет
Москва, Россия

b Российская академия наук, Объединенный институт высоких температур
Москва, Россия

c Московский государственный технический университет имени Н.Э. Баумана
Москва, Россия

* E-mail: gbelov@yandex.ru

Поступила в редакцию 30.08.2018
После доработки 23.10.2018
Принята к публикации 11.11.2018

Полный текст (PDF)

Аннотация

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

Ключевые слова: термодинамическое равновесие, фазовый состав, гетерогенные системы

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

Методы равновесной термодинамики широко применяются для решения научных и технических задач, к числу которых относится, в частности, расчет равновесного состава сложных термодинамических систем и построение фазовых диаграмм [1–3]. В последнее время получила распространение практика встраивания процедур расчета равновесного состава в более сложные программные комплексы, предназначенные для моделирования гидро- и газодинамических процессов (см., например, [4]).

Современные программы расчета равновесного состава сложных термодинамических систем, как правило, сопряжены с базами данных по термодинамическим свойствам индивидуальных веществ. Информация из этих баз данных обычно выбирается автоматически на основании списка химических элементов, включенных в термодинамическую систему. При этом число возможных фаз в системе чаще всего превышает число независимых компонентов в ней (обычно число независимых компонентов совпадает с числом химических элементов). Поэтому одной из главных задач расчета термодинамического равновесия является определение фазового состава термодинамической системы. Иными словами, прежде чем рассчитать равновесный состав необходимо определить, какие фазы системы являются стабильными при заданных условиях равновесия. Задача определения перечня стабильных фаз является основной и при построении фазовых диаграмм [5, 6]. Однако число компонентов при построении фазовых диаграмм, как правило, относительно невелико, поэтому методы, применяемые для их построения в общем случае непригодны для анализа фазового состава сложных термодинамических систем, число независимых компонентов к которых может составлять несколько десятков, а число продуктов химических реакций может достигать нескольких сотен. С расчетами такого рода приходится сталкиваться, например, при моделировании процессов в облученном топливе ядерных реакторов [7].

Поскольку фазовый состав термодинамической системы априори не задан, неизвестен и вид фундаментального уравнения состояния (характеристических функций). В своей работе “Метод геометрического представления термодинамических свойств” [8] Дж.В. Гиббс пишет, что поверхность, описывающая взаимосвязь термодинамических параметров равновесной многофазной термодинамической системы (термодинамическая поверхность), в общем случае состоит из двух частей – первичной и производной. Первичная поверхность соответствует гомогенным фазам, а производная – гетерогенным. При этом первичная поверхность – свойство фазы, а производная – свойство системы. Объединение двух поверхностей (первичной и производной) Гиббс назвал поверхностью рассеянной энергии. По сути дела, эта поверхность и является геометрическим изображением характеристической функции термодинамической системы.

В цитированной выше работе Гиббса указано, что производную поверхность легко построить, если известен вид первичных поверхностей. И это действительно так, если речь идет о геометрическом построении поверхности системы малой размерности. Однако аналитический способ ее построения не так прост, даже для бинарных систем. В работах Г.Ф. Воронина и др. [5, 9, 10] показано, что эту поверхность рассеянной энергии можно рассматривать как выпуклую оболочку множества точек, которые характеризуют значения одной из характеристических функций в соответствующей системе координат. В частности, для построения фазовых диаграмм целесообразно анализировать расположение точек, характеризующих значение энергии Гиббса в координатах GTP(x), где x – вектор состава.

Вопрос о выпуклости характеристических функций термодинамической системы неоднократно рассматривался в прошлом (см. например [9, 11–14]). Установлено, что из выпуклости энтропии и внутренней энергии системы следует выпуклость остальных характеристических функций. В работе [15] показано, что термодинамические потенциалы являются выпуклыми (обращенными выпуклостью вниз) функциями экстенсивных и вогнутыми (обращенными выпуклостью вверх) интенсивных переменных.

Поскольку далее будут использоваться понятия “выпуклая область” и “выпуклая функция”, приведем соответствующие определения. Область является выпуклой, если отрезок прямой, соединяющий две любые точки области, принадлежит этой области. Функция $f{\kern 1pt} \left( x \right)$ выпукла на выпуклой области X, если для любых двух точек ${{x}_{1}},{{x}_{2}} \in X$выполняется соотношение:

$\begin{gathered} f\left[ {\theta {{x}_{2}} + \left( {1 - \theta } \right){\kern 1pt} {{x}_{1}}} \right] \leqslant \theta f{\kern 1pt} \left( {{{x}_{2}}} \right) + \left( {1 - \theta } \right){\kern 1pt} f{\kern 1pt} \left( {{{x}_{1}}} \right), \\ 0 < \theta < 1. \\ \end{gathered} $

Для вогнутой функции, определенной на выпуклом множестве, знак неравенства следует изменить [16]. Таким образом, выпуклая функция имеет минимум, а вогнутая – максимум, причем если поменять знак выпуклой функции, получим вогнутую функцию.

ОБЩАЯ ФОРМУЛИРОВКА ЗАДАЧИ РАСЧЕТА РАВНОВЕСИЯ

В свое время Дж.В. Гиббс предложил следующие формулировки условий равновесия: для равновесия изолированной термодинамической системы, рабочие координаты которой зафиксированы, необходимо и достаточно, чтобы ${{\left( {\delta S} \right)}_{U}} \leqslant 0$ либо ${{\left( {\delta U} \right)}_{S}} \geqslant 0$. Символ δ здесь означает виртуальное смещение (вариацию) в смысле теоретической механики. Более подробно физический смысл этих формулировок раскрыт позже рядом авторов (см., например, [14]).

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

$S{\kern 1pt} \left( {U,V,\vec {n}} \right) \to {\text{max,}}$
$U = {\text{const}},\quad V = {\text{const,}}$
$\mathop \sum \limits_{i \in {{I}_{T}}} \,{{a}_{{ki}}}{{n}_{i}} = {{b}_{k}},\quad k \in {{I}_{E}},$
${{n}_{i}} \geqslant 0,\quad i \in {{I}_{T}},$
IT – множество индексов веществ в системе, IE – множество индексов химических элементов в системе, ${{a}_{{ki}}}$– количество атомов элемента k в молекуле вещества i, ni – количество вещества i, ${{b}_{k}}$ – содержание элемента k в системе, $\vec {n}$ – вектор состава, компонентами которого являются значения ni.

Функция Лагранжа для расчета равновесного состава химически реагирующей системы имеет вид

$\Lambda = S + \mathop \sum \limits_{k \in {{I}_{E}}} \,{{\lambda }_{k}}\left( {\mathop \sum \limits_{i \in {{I}_{T}}} \,{{a}_{{ki}}}{{n}_{i}} - {{b}_{k}}} \right),$
λk – неопределенные множители Лагранжа. Неизвестными являются значения λk и ni.

Координаты условного максимума энтропии определяются следующей системой уравнений:

(1)
${{\left( {\frac{{\partial S}}{{\partial {{n}_{i}}}}} \right)}_{{U,V,{{n}_{{j \ne i}}}}}} + \mathop \sum \limits_{k \in {{I}_{E}}} \,{{a}_{{ki}}}{{\lambda }_{k}} = 0,\quad i \in {{I}_{T}},$
(2)
$\mathop \sum \limits_{i \in {{I}_{T}}} \,{{a}_{{ki}}}{{n}_{i}} = {{b}_{k}},\quad k \in {{I}_{E}}.$

Используя фундаментальное уравнение Гиббса, уравнение (1) можно преобразовать к виду

(3)
$ - {{\mu }_{i}} + \mathop \sum \limits_{k \in {{I}_{E}}} \,{{a}_{{ki}}}{{\lambda }_{k}} = 0,\quad i \in {{I}_{T}}.$

Для веществ в состоянии идеального газа

${{\mu }_{i}} = \mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right) + RT\ln {{\bar {p}}_{i}},$
для компонентов идеального раствора
${{\mu }_{i}} = \mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right) + RT{\text{ln}}{{x}_{i}},$
а для конденсированных веществ, образующих однокомпонентные фазы
${{\mu }_{i}} = \mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right),$
$\quad\mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right) = H_{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right) - TS_{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right),$
$H_{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right)$, $S_{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right)$ – энтальпия и энтропия вещества i при температуре T и стандартном давлении.

Умножив соотношение (3) на ni (здесь ni – число молей вещества i) и просуммировав его по всем веществам, можно получить равенство вида

$\mathop \sum \limits_{i \in {{I}_{T}}} \,{{n}_{i}}{{\mu }_{i}} = \mathop \sum \limits_{k \in {{I}_{E}}} \mathop \sum \limits_{i \in {{I}_{T}}} \,{{n}_{i}}{{a}_{{ki}}}{{\lambda }_{k}} = \mathop \sum \limits_{k \in {{I}_{E}}} \,{{b}_{k}}{{\lambda }_{k}} = G,$
из которого следует, что λk имеет смысл химического потенциала k-го химического элемента термодинамической системы, поскольку ${{\lambda }_{k}} = {{\left( {\partial G{\text{/}}\partial {{b}_{k}}} \right)}_{{p,T,{{b}_{{j \ne k}}}}}}.$

Соотношения (2) можно представить в виде

$\mathop \sum \limits_{j \in {{I}_{P}}} \,{{n}_{j}}\mathop \sum \limits_{i \in {{I}_{j}}} \,{{a}_{{ki}}}{{x}_{{ij}}} = {{b}_{k}},\quad k \in {{I}_{E}},$
в котором nj – количество вещества фазы j, IP – множество индексов фаз, включенных в модель термодинамической системы, Ij – множество индексов веществ, образующих фазу j, xij – мольная доля вещества i в фазе j. Такая форма записи облегчает понимание сложности расчета равновесия в термодинамической системе с неопределенным составом фаз.

В [17] предложено вычислять мольную долю вещества i стабильной (с ненулевым количеством) фазы по формуле

(4)
${{x}_{i}} = {\text{exp}}\left[ {\left( { - \mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right) + \mathop \sum \limits_{k \in {{I}_{E}}} \,{{a}_{{ki}}}{{\lambda }_{k}}} \right){\text{/}}RT} \right],$
здесь R – газовая постоянная. В более общем случае справедливо неравенство [18]
${{x}_{i}} \geqslant {\text{exp}}\left[ {\left( { - \mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right) + \mathop \sum \limits_{k \in {{I}_{E}}} \,{{a}_{{ki}}}{{\lambda }_{k}}} \right){\text{/}}RT} \right],$
которое можно интерпретировать так: набор значений λk считается допустимым, если вклад атомов в энергию молекулы не превышает ее химического потенциала
$\mathop \sum \limits_{k \in {{I}_{E}}} \,{{a}_{{ki}}}{{\lambda }_{k}} \leqslant \mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right) + RT\ln {{x}_{i}}.$
Из этого неравенства следует, что сумма
$\mathop \sum \limits_{i \in {{I}_{j}}} \exp \left[ {\left( { - \mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right) + \mathop \sum \limits_{k \in {{I}_{E}}} \,{{a}_{{ki}}}{{\lambda }_{k}}} \right){\text{/}}RT} \right]$
не должна превышать 1, ее можно рассматривать как характеристику возможности образования фазы j.

Равенство (3) должно выполняться для всех веществ фазы, присутствующей в равновесной системе, таким образом, оно выражает условие термодинамической стабильности фазы. В соответствующей системе координат соотношение $\sum\nolimits_{k \in {{I}_{E}}}^{} {{{a}_{{ki}}}{{\lambda }_{k}}} = {{\mu }_{i}}$ является уравнением гиперплоскости, касательной к поверхности энергии Гиббса стабильных фаз, [19]. Для того, чтобы термодинамическая система находилась в равновесии, необходимо, чтобы химические потенциалы компонентов стабильных растворов и стехиометрических фаз принадлежали этой гиперплоскости [20]. Если мольные энергии Гиббса всех возможных фаз лежат выше гиперплоскости $\sum\nolimits_{k \in {{I}_{E}}}^{} {{{a}_{{ki}}}{{\lambda }_{k}}} \quad$ или принадлежат ей, это – достаточное условие равновесия [19]. При этом предполагается, что выполняются условия материального баланса и неотрицательности чисел молей веществ, а сумма мольных долей веществ каждой стабильной фазы равна 1. Иными словами, разность $\sum\nolimits_{i \in {{I}_{j}}}^{} {{{x}_{{ij}}}} \left[ {\mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right) - \sum\nolimits_{k \in {{I}_{E}}}^{} {{{a}_{{ki}}}{{\lambda }_{k}}} } \right]$ характеризует расстояние от мольной энергии Гиббса фазы j до гиперплоскости, и это расстояние в равновесной системе не может быть меньше нуля [21, 22]. Таким образом, выполнение условия

(5)
$\mathop \sum \limits_{i \in {{I}_{j}}} \,{{x}_{{ij}}}\left[ {\mu _{i}^{ \circ }{\kern 1pt} \left( T \right) - \mathop \sum \limits_{k \in {{I}_{E}}} \,{{a}_{{ki}}}{{\lambda }_{k}}} \right] \geqslant 0,\quad j \in {{I}_{P}}$
при соблюдении ограничений $\sum\nolimits_{i \in {{I}_{j}}}^{} {{{x}_{{ij}}}} = 1$, $j \in {{I}_{{AP}}}$, где IAP – множество индексов стабильных фаз, необходимо и достаточно для существования термодинамического равновесия в системе [19, 23].

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

$\begin{gathered} \Lambda = S + \mathop \sum \limits_{k \in {{I}_{E}}} \,{{\lambda }_{k}}\left( {\mathop \sum \limits_{j \in {{I}_{{AP}}}} \,{{n}_{j}}\mathop {\sum \,}\limits_{i \in {{I}_{j}}} {{a}_{{ki}}}{{x}_{{ij}}} - {{b}_{k}}} \right) + \\ + \;\mathop \sum \limits_{j \in {{I}_{{AP}}}} {{\lambda }_{j}}\left( {\mathop \sum \limits_{i \in {{I}_{j}}} \,{{x}_{{ij}}} - 1} \right). \\ \end{gathered} $

АНАЛИЗ УСТОЙЧИВОСТИ ФАЗ

Двухступенчатый расчет равновесия

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

Подробное описание процедуры анализа устойчивости фаз, используемой в одной из программ расчета равновесия, приводится в работе [20]. Отмечается, что при анализе фазового состава необходимо принимать во внимание три главных критерия:

1) баланс масс химических элементов (уравнение (2)),

2) условия термодинамической стабильности фаз вида (3),

3) правило фаз Гиббса.

Задача расчета равновесия разбивается на две подзадачи:

1) определение фазового состава,

2) расчет равновесного состава для данной комбинации стабильных фаз.

Практичную процедуру расчета равновесия можно описать так. Исходя из некоторых априорных соображений, задается набор стабильных фаз. В частности, если проводится серия расчетов, можно заимствовать перечень стабильных фаз из предыдущего расчета. С этим набором осуществляется несколько итераций вычисления равновесного состава и значений неопределенных множителей Лагранжа λk (мольных энергий химических элементов), проверяется выполнение условия (5) для компонентов всех фаз системы. В случае необходимости перечень стабильных фаз корректируется, расчет равновесного состава с измененным фазовым составом продолжается. Если для данного набора стабильных фаз удается рассчитать равновесный состав и значения λk и при этом для всех компонентов стабильных фаз выполняются условия (5), решение считается найденным.

Подобный подход использован в целом ряде работ (см., например [24–28]). В работе [25] возможность появления конденсированной фазы вещества i оценивалась путем сравнения химических потенциалов этого вещества в газообразном (g) и конденсированном (c) состояниях. Фаза добавлялась в систему при выполнении неравенства

${{[\mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right)]}_{c}} < {{[\mu _{i}^{{^{ \circ }}}{\kern 1pt} \left( T \right) + RT\ln {{\bar {p}}_{i}}]}_{g}}.$

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

В работе [27] на первом шаге определяется равновесный состав газовой фазы и значения λk. Далее для каждого конденсированного вещества вычисляется разность

${{\delta }_{i}} = {{\mu }_{i}} - \sum\nolimits_{k \in {{I}_{E}}}^{} {{{a}_{{ki}}}{{\lambda }_{k}}} $. Если для некоторого конденсированного вещества значение ${{\delta }_{i}}$ отрицательно, оно добавляется в систему, и производится новый расчет равновесного состава. Если условие отрицательности ${{\delta }_{i}}{\;}$ выполняется для нескольких конденсированных веществ, в систему добавляется вещество с наименьшим значением ${{\delta }_{i}}$. Процедура повторяется до тех пор, пока не останется конденсированных веществ с отрицательными значениями ${{\delta }_{i}}$. Процедура анализа фазового состава предполагает и удаление конденсированного вещества из системы, если его концентрация становится отрицательной.

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

Расчет равновесного состава с данной комбинацией стабильных фаз предполагает необходимость решения системы нелинейных уравнений. В цитируемых выше работах для этой цели используется метод Ньютона, который предполагает применение матрицы Гессе H:

${\mathbf{Hx}} = {\mathbf{A}},$
где x – вектор неизвестных, A – вектор констант. При изменении перечня фаз возможна ситуация, когда две или более строк матрицы матрица Гессе станут линейно зависимыми, матрица станет вырожденной. Поскольку в этом случае решение найти нельзя, необходимо избегать подобных комбинаций фаз.

Изменение фазового состава не должно нарушать правила фаз Гиббса. В соответствии с теоремой Дюгема [29], равновесное состояние термодинамической системы, исходные массы которой известны, однозначно характеризуется заданными значениями двух термодинамических параметров. С учетом этого обстоятельства правило фаз Гиббса ограничивает число фаз в системе – оно не должно превышать числа независимых компонентов.

Удаление фазы из системы (набора стабильных фаз) производится в том случае, если в процессе вычислений число молей фазы принимает нулевое или отрицательное значение. Фаза j может быть добавлена в систему, если вычисленное по формуле ${{\pi }_{j}} = \sum\nolimits_{i \in {{I}_{j}}}^{} {{{x}_{{ij}}}} \left[ {\mu _{i}^{{^{ \circ }}}\left( T \right) - \sum\nolimits_{k \in {{I}_{E}}}^{} {{{a}_{{ki}}}{{\lambda }_{k}}} } \right]$ значение πj отрицательно [20]. Нетрудно видеть, что в этом случае πj характеризует расстояние между значением мольной энергии Гиббса фазы и касательной гиперплоскостью. Если добавление новой фазы приводит к нарушению правила фаз Гиббса, то после добавления одной фазы другую следует удалять из перечня.

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

Эффективность описываемого алгоритма во многом зависит от частоты обновления перечня стабильных фаз. Если обновлять его слишком часто, есть вероятность получить неверное решение. Если же обновлять его слишком редко, время вычислений возрастает. Наконец, возможна ситуация, когда рассмотренным выше критериям замены отвечают сразу несколько фаз. Алгоритм программы [4] предполагает возможность добавления и (или) удаления только одной фазы.

Метод “больших молекул”

В третьем томе справочника [30] для анализа фазового состава предложено использовать метод “больших молекул”. В соответствии с этим подходом, частицы вещества в конденсированном состоянии рассматриваются как большие молекулы, состоящие из достаточно большого количества обычных молекул. Предполагается, что вещество, состоящее из больших молекул, образует отдельную газообразную фазу. Реакцию диссоциации большой молекулы $M_{i}^{x}$ на атомы можно записать в виде

$M_{i}^{x} = {{n}^{x}}\mathop \sum \limits_{k \in {{I}_{E}}} \,{{a}_{{ki}}}{{A}_{k}},$
где nx – число нормальных молекул, образующих большую молекулу, aki – число атомов химического элемента k в нормальной молекуле. В [30] показано, что для практических целей величину nx целесообразно выбирать из диапазона 1000–10000. Например, для конденсированного Al2O3 при nx = 1000 химическая формула большой молекулы имеет вид Al2000O3000 .

Если известны свойства вещества в конденсированном состоянии $H_{i}^{{^{ \circ }}},\;S_{i}^{{^{ \circ }}}$, свойства вещества, образованного большими молекулами, можно рассчитать по формулам

$H_{i}^{{^{ \circ }x}} = {{n}^{x}}H_{i}^{{^{ \circ }}},\quad\quad S_{i}^{{^{ \circ }x}} = {{n}^{x}}S_{i}^{{^{ \circ }}}.$

Формула для расчета константы равновесия реакции диссоциации большой молекулы на атомы имеет вид

${{K}^{x}}\, = \,\exp \left( {\frac{{{{n}^{x}}\sum\limits_{k \in {{I}_{E}}} {{{a}_{k}}S_{k}^{{^{ \circ }}}\, - \,{{S}^{{^{ \circ }x}}}} }}{R}\, - \,\frac{{{{n}^{x}}\sum\limits_{k \in {{I}_{E}}} {{{a}_{k}}H_{k}^{{^{ \circ }}}\, - \,{{H}^{{^{ \circ }x}}}} }}{{RT}}} \right).$

Как отмечается в [1], модель реагирующей смеси может предполагать наличие как больших, так и обычных молекул одного вещества в разных агрегатных состояниях (газ и конденсат).

Дальнейшее развитие метод больших молекул получил в работах Б.Г. Трусова [31–33]. Условие равновесия для отдельных конденсированных фаз было предложено записать в виде

(6)
$\left( { - \mu _{i}^{{^{ \circ }}} + \mathop \sum \limits_{k \in {{I}_{E}}} \,{{a}_{{ki}}}{{\lambda }_{k}}} \right){{n}_{i}} = 0,\quad i \in {{I}_{{SP}}},$
ISP – множество индексов веществ, которые могут образовать однокомпонентные конденсированные фазы. Как нетрудно видеть, данное равенство выполняется либо в случае, если ni > 0: $ - \mu _{i}^{{^{ \circ }}} + \sum\nolimits_{k \in {{I}_{E}}}^{} {{{a}_{{ki}}}{{\lambda }_{k}}} = 0$, либо если ni = 0: $ - \mu _{i}^{{^{ \circ }}} + $ $ + \;\sum\nolimits_{k \in {{I}_{E}}}^{} {{{a}_{{ki}}}{{\lambda }_{k}}} < 0$.

При разработке алгоритма программы АСТРА Б.Г. Трусов предложил заменить на этапе определения фазового состава соотношение (6) его приближенным аналогом вида ${{n}_{i}} = $ $ = \;{\text{exp}}\left[ {A\left( { - \mu _{i}^{{^{ \circ }}} + \sum\nolimits_{k \in {{I}_{E}}}^{} {{{a}_{{ki}}}{{\lambda }_{k}}} } \right)} \right]$, где A – некоторое большое число (≈1000). Данное выражение используется в расчете до тех пор, пока фазовый состав термодинамической системы не установлен. Фазовый состав считается установленным после того, как относительная погрешность результатов двух последовательных итераций становится меньше 10–4. Условием присутствия i-й однокомпонентной фазы в системе является выполнение неравенств:

$ - \mu _{i}^{{^{ \circ }}} + \mathop \sum \limits_{k \in {{I}_{E}}} \,{{a}_{{ki}}}{{\lambda }_{k}} \leqslant 0,\quad\quad {{n}_{i}} > 0.$
При этом число фаз в системе не должно нарушать правила фаз Гиббса.

Возможность присутствия в системе конденсированного раствора проверяется аналогичным образом. До определения фазового состава уравнение, определяющее возможность образования конденсированного раствора с индексом j, имеет вид:

$\mathop \sum \limits_{i \in {{I}_{j}}} \,{{x}_{{ij}}} - \frac{{\ln {{n}_{j}}}}{A} = 1,$
xij – мольная доля вещества i в фазе j, nj – число молей фазы j, A = 1000. После того, как фазовый состав термодинамической системы установлен, в дальнейших расчетах применяется одно из уравнений: $\sum\nolimits_{i \in {{I}_{j}}}^{} {{{x}_{{ij}}}} = 1$, если фаза присутствует в системе, либо nj = 0. Мольная доля компонента раствора вычисляется по формуле (4).

Использование процедуры линейного программирования

Анализ литературы показывает, что на практике чаще применяется формулировка условия равновесия термодинамической системы с использованием энергии Гиббса

$G{\kern 1pt} \left( {T,p,\vec {n}} \right) \to {\text{min,}}$
$T = {\text{const}},\quad p = {\text{const,}}$
$\mathop \sum \limits_{i \in {{I}_{T}}} \,{{a}_{{ki}}}{{n}_{i}} = {{b}_{k}},\quad k \in {{I}_{E}},$
${{n}_{i}} \geqslant 0,\quad i \in {{I}_{T}}.$
Если система содержит только однокомпонентные конденсированные фазы, то целевая функция в задаче условной минимизации имеет вид
$\mathop \sum \limits_j \,{{n}_{j}}G_{j}^{{^{ \circ }}} \to {\text{min,}}$
т.е. является линейной. Ограничения, выражающие закон сохранения вещества также линейны. Поэтому фазовый и равновесный состав такой термодинамической системы можно легко рассчитать с использованием процедуры линейного программирования.

Однако процедуру линейного программирования можно использовать и для поиска координат условного экстремума нелинейной функции. Соответствующий алгоритм расчета равновесного состава смеси веществ был предложен в [34]. Смысл подхода, предложенного в этой работе, заключается в замене нелинейных слагаемых вида $\alpha {\text{ln}}\alpha $ кусочно-линейной функцией. Это можно сделать, положив ${{\alpha }_{i}} = {{x}_{{i\quad}}}$ и $\quad{{\beta }_{i}} = {{\alpha }_{i}}{\text{ln}}{{\alpha }_{i}}$. В частности, для однофазной системы такой прием позволяет определить химический состав термодинамической системы путем итерационного решения следующей задачи:

$\mathop \sum \limits_{i \in {{I}_{T}}} \,{{n}_{i}}G_{i}^{{^{ \circ }}} + \mathop \sum \limits_{i \in {{I}_{T}}} \,\mathop \sum \limits_{j \in J} \,{{\beta }_{{ij}}}{{y}_{{ij}}} \to {\text{min,}}$
при условиях (2),
$n = \mathop \sum \limits_{i \in {{I}_{T}}} \,{{n}_{i}},$
${{n}_{i}} = \mathop \sum \limits_{j \in J} \,{{\alpha }_{{ij}}}{{y}_{{ij}}},\quad i \in {{I}_{T}},$
${{y}_{{ij}}} \geqslant 0,\quad i \in {{I}_{T}},\quad j \in J.$
Здесь J – множество индексов точек, в которых отрезки прямых сопрягаются с аппроксимируемой кривой, yij – вспомогательные переменные. Для практических целей достаточно выбрать два отрезка, тогда J = 1, 2, 3.

Нетрудно видеть, что в такой формулировке мы имеем дело с задачей линейного программирования, и для ее решения можно использовать одну из отработанных процедур. Данный подход был использован для расчета равновесного состава многокомпонентных гетерогенных термодинамических систем с несколькими растворами [35, 36]. Было установлено, что получить окончательное решение с помощью только процедуры линейного программирования удается не всегда. В частности, линейная аппроксимация плохо работает в области малых концентраций компонентов смеси. Однако, фазовый состав можно определить с ее помощью довольно надежно. Поэтому целесообразно дополнить вышеописанный алгоритм расчета равновесного состава процедурой нелинейной оптимизации.

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Большинство рассмотренных выше способов определения фазового состава сложной термодинамической системы предназначены для случая, когда температура задана. Однако в практике термодинамических вычислений довольно часто встречаются ситуации, когда равновесная температура неизвестна и ее значение необходимо найти (адиабатное горение, изоэнтропное расширение и т.д.). Задача расчета равновесного состава в этом случае усложняется, поскольку, не зная температуру, нельзя вычислить значение химического потенциала. Решать ее можно либо путем многократного расчета равновесия с заданными значениями температуры и давления, как предлагается в [37], либо за счет расширения расчетной системы уравнений, как сделано в программе АСТРА [31]. Трудоемкость вычислений в первом случае существенно возрастает, поскольку вычислять равновесие приходится многократно, до тех пор, пока не будет найдена такая комбинация температуры и давления, при которой будут выполнены заданные условия равновесия. Однако надежность определения фазового состава в этом случае выше, чем при использовании расширенной системы уравнений.

ЗАКЛЮЧЕНИЕ

Основные трудности расчета равновесного состава сложной термодинамической системы обусловлены тем, что необходимо найти условный экстремум функции, вид которой априори неизвестен. Как показывает проведенный анализ литературы, для реализации алгоритмов расчета равновесного состава многокомпонентной гетерогенной системы обычно используются эвристические правила, цель которых – обеспечить по возможности быстрое и надежное определение фазового состава системы. При этом, по всей видимости, ни один из алгоритмов не гарантирует надежного и эффективного определения перечня стабильных фаз. На наш взгляд, применение процедуры линейного программирования позволяет повысить надежность расчета фазового и равновесного состава сложных термодинамических систем. Работа выполнена в рамках темы “Химическая термодинамика” ( АААА-А16-116061750195-2).

Список литературы

  1. Алемасов В.Е., Дрегалин А.Ф., Крюков В.Г. и др. Математическое моделирование высокотемпературных процессов в энергосиловых установках. М.: Наука, 1989. 256 с.

  2. The SGTE Casebook, Second Edition: Thermodynamics at Work / Ed. by K. Hack. CRC Press, 2008. 454 p.

  3. Чудненко К.В. Термодинамическое моделирование в геохимии: теория, алгоритмы, программное обеспечение, приложения. Новосибирск: Академическое издательство “Гео”, 2010. 287 с.

  4. Piro M.H.A., Simunovic S., Besmann T.M. et al. // Computational Materials Science. 2013. V. 67. P. 266. http://dx.doi.org/https://doi.org/10.1016/j.commatsci.2012.09.011

  5. Воронин Г.Ф. // Журн. физ. химии. 2003. Т. 77. № 10. С. 1874.

  6. Emelianenko M., Zi-Kui Liu, Qiang Du // Computational Materials Science. 2006. V. 35. P. 61. https://doi.org/10.1016/j.commatsci.2005.03.004

  7. Besmann T.M., McMurray J.W., Simunovic S. // CALPHAD: Computer Coupling of Phase Diagrams and Thermochemistry. 2016. V. 55. P. 47. http://dx.doi.org/https://doi.org/10.1016/j.calphad.2016.04.004

  8. Гиббс Дж. В. Термодинамика, статистическая механика. М.: Наука, 1982. 584 с.

  9. Воронин Г.Ф. // Журн. физ. химии. 2005. Т. 79. № 12. С. 2126.

  10. Воронин Г.Ф., Восков А.Л. // Вестн. Моск. ун-та. Сер. 2. Химия. 2013. Т. 54. № 1. С. 3.

  11. Galgani L., Scotti A. // Physica. 1968. V. 40. P. 150.

  12. Galgani L., Scotti A. // Ibid. 1969. V. 42. P. 242.

  13. ter Horst H.J. // Ann. Phys. 1987. V. 176. P. 183.

  14. Callen H.B. Thermodynamics and an Introduction to Thermostatics. Wiley, 1985. 495 p.

  15. Galgani L., Scotti A. // Pure and Appl. Chem. 1970. V. 22. № 3–4. P. 229.

  16. Банди Б. Методы оптимизации. Вводный курс. М.: Радио и связь, 1988. 128 с.

  17. White W.B., Johnson S.M., Dantzig G.B. // J. Chem. Phys. 1958. V. 28. № 5. P. 751.

  18. Dorn W.S. // Ibid. 1960. V. 32. № 5. P. 1490.

  19. Smith J.V., Missen R.W., Smith W.R. // AIChE Journal. 1993. V. 39. № 4. P. 707. https://doi.org/10.1002/aic.690390421

  20. Piro M.H.A. // CALPHAD: Computer Coupling of Phase Diagrams and Thermochemistry. 2017. V. 58. P. 115. https://doi.org/https://doi.org/10.1016/j.calphad.2017.06.002

  21. Michelsen M.L. // Fluid Phase Equilibria. 1982. V. 9. P. 1.

  22. Zhang H., Bonilla-Petriciolet A., Rangaiah G. // The Open Thermodynamics Journal. 2011. V. 5. P. 71.

  23. Baker L.E., Pierce A.C., Luks K.D. // Society of Petroleum Engineers Journal. 1982. V. 22. P. 731.

  24. Boll R.H. // J. Chem. Phys. 1961. V. 34. P. 1108. https://doi.org/10.1063/1.1731708

  25. Zeleznik F.J., Gordon S. A General IBM 704 or 7090 Computer Program for Computation of Chemical Equilibrium Compositions, Rocket Performance, and Chapman-Jouguet dentonations. NASA TN D-1454. 1962. 150 p.

  26. Eriksson G. // Acta Chemica Scandinavica. 1970. V. 25. P. 2651.

  27. Gordon S., McBride B.J. Computer Program for Calculation Complex Chemical Equilibrium Compositions and Applications. I. Analysis. RP-1311. 1994. 58 p.

  28. Shobu K., Tabaru T. // Materials Transactions. 2005. V. 46. № 6. P. 1175. https://doi.org/https://doi.org/10.2320/matertrans.46.1175

  29. Пригожин И., Дефэй Р. Химическая термодинамика. Новосибирск: Наука, 1966. 502с.

  30. Термодинамические и теплофизические свойства продуктов сгорания, Том III / Под ред. акад. В.П. Глушко. М.: ВИНИТИ, 1973. 624 с.

  31. Синярев Г.Б., Ватолин Н.А., Трусов Б.Г. и др. Применение ЭВМ для термодинамических расчетов металлургических процессов. М.: Наука, 1982. 263 с.

  32. Трусов Б.Г. Термодинамический метод анализа высокотемпературных состояний и процессов и его практическая реализация: Дис. … докт. техн. наук. М.: МВТУ им. Н.Э. Баумана, 1984. 210 с.

  33. Ватолин Н.А., Моисеев Г.К., Трусов Б.Г. Термодинамическое моделирование в высокотемпературных неорганических системах. М.: Металлургия, 1994. 352 с.

  34. White W.B., Johnson S.M., Dantzig G.B.// Management Science. 1958. V. 5. № 1. P. 38.

  35. Белов Г.В. // Вычислительные методы и программирование. 2009. Т. 10. С. 56.

  36. Belov G. // J. Mathem. Chem. 2010. V. 47. Issue 1. P. 446.

  37. Smith W.R. // Theoretical Chemistry. Advances and Perspectives. 1980. V. 5. P. 185. https://doi.org/10.1021/i160073a001

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал физической химии

Архивы выпусков Информация о журнале Отправить рукопись в журнал