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

Журнал физической химии, 2022, T. 96, № 10, стр. 1411-1420

Численный анализ термодинамического определения поверхностного натяжения парожидкостной системы в модели решеточного газа

Е. С. Зайцева a, Ю. К. Товбин a, *

a Институт общей и неорганической химии им. Н.С. Курнакова РАН
119991 Москва, Россия

* E-mail: tovbinyk@mail.ru

Поступила в редакцию 26.02.2022
После доработки 26.02.2022
Принята к публикации 28.02.2022

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

Аннотация

Проведен численный анализ термодинамического определения поверхностного натяжения (ПН) парожидкостной системы по Гиббсу как избыточной величины свободной энергии ΔF двухфазной системы с учетом и без учета наличия границы раздела фаз. Расчет проведен в простейшем варианте модели решеточного газа (МРГ) при учете взаимодействия ближайших соседей в квазихимическом приближении. Сопоставлены разные способы расчета ПН, которое выражается через разные парциальные вклады $M_{q}^{i}$ в избыточную свободную энергию ΔF (здесь i = A – молекула А и i = V – вакансии, 1 ≤ q ≤ κ, q – номер монослоя внутри границы, κ – ее ширина). Получена неоднозначность значений ПН в зависимости от вида функций $M_{q}^{i}$. Эти различия в величинах ПН демонстрируются на примере температурной зависимости ПН для плоской границы и для зависимости ПН от размера капли при фиксированной температуре. Результаты расчета ПН по термодинамическому определению сопоставлены с аналогичным расчетом по варианту определения ПН, который учитывает специфику вакансий в МРГ как механическую характеристику системы.

Ключевые слова: поверхностное натяжение, система пар–жидкость, модель решеточного газа, квазихимическое приближение, малые системы

Поверхностное натяжение (ПН) по Гиббсу определяется как избыточная свободная энергия ΔF двухфазной системы с учетом и без учета наличия границы раздела фаз [1–6]. Стартуя с общих термодинамических определений для свободной энергии Гельмгольца (F) и Гиббса (G): F = G – PVsys, и $G = \sum\nolimits_{i = 1}^{{{s}_{c}}} {{{\mu }_{i}}{{N}_{i}}} $, где P – давление и Vsys – объем системы, μi и Ni – химический потенциал и число молекул компонента i смеси, состоящей из sc сортов молекул [7–9], запишем

(1)
$F = \sum\limits_{i = 1}^{{{s}_{c}}} {{{\mu }_{i}}} {{N}_{i}} - P{{V}_{{sys}}}.$
Тогда избыток свободной энергии на границе определяется как $\Delta F = F - {{F}_{\alpha }} - {{F}_{\beta }} = \sum\nolimits_{i = 1}^{{{s}_{c}}} {{{\mu }_{i}}\Delta {{N}_{i}}} + $ $ + \;\sigma A$, где ${{F}_{\alpha }}$ и ${{F}_{\beta }}$ – свободные энергии однородных областей жидкой (α) и паровой (β) фаз, $\Delta {{N}_{i}}$ – избыток компонента i на границе, σ и А – ПН и площадь границы раздела фаз. Выражения для ${{F}_{\alpha }}$ и ${{F}_{\beta }}$ имеют вид аналогичный (1), с характеристиками, относящимися к фазе пара и жидкости. Условие $\sum\nolimits_{i = 1}^{{{s}_{c}}} {{{\mu }_{i}}} \Delta {{n}_{i}} = 0$ определяет положение эквимолекулярной разделяющей поверхности, к которой относится величина ПН

(2)
$\sigma = \Delta F{\text{/}}A.$

До последнего времени на основе изложенного подхода Гиббса [1] существовало четыре способа расчета ПН [2–6], что фактически свидетельствовало об отсутствии однозначного способа расчета данной величины. Эти термодинамические построения проводились при разных предположениях о характере распределения напряжений для локального тензора давления внутри переходной области границы раздела фаз метастабильных капель. Хотя физические основы термодинамики не позволяют ее использовать внутри границы фаз [10], тем не менее, эти подходы позже были трансформированы в соответствующие молекулярно-статистические теории [2, 6, 10].

В работах [11, 12] был предложен строго равновесный способ расчета ПН на основе определения Гиббса, который исключает появление метастабильных капель. Этот способ качественно отличается от классического термодинамического определения ПН [1–10] тем, что на переходную область переносится требование выполнения соотношения для времен релаксаций процессов переноса импульса и массы, которые отсутствуют в классической термодинамике (оперирующей уравнением Лапласа для искривленных границ) и в статистической термодинамике (интегральные уравнения и метод МД) [2, 6, 10]. Эта специфика приводит при расчете ПН к использованию средних величин локальных химпотенциалов и давлений внутри локальных областей границы, а не их тензорных компонентов.

Данный подход был сформулирован на основе, так называемой, модели решеточного газа (МРГ) [13–15]. МРГ является наиболее распространенной при исследовании фазовых состояний веществ и с ее помощью получены наиболее важные результаты по теории фазовых переходов, включая критические области парожидкостной системы [13–21]. Данная модель давно и активно применяется при исследовании плоских границ раздела фаз [15, 22–27]. Позже в рамках МРГ были разработаны подходы для описания искривленных поверхностей (сферические и цилиндрические капли) [10, 28–31], а также для описания искривленных границ раздела пар–жидкость, имеющих сложную геометрию, в трехагрегатных системах [32–34].

В МРГ выражение для свободной энергии (1) в объемной фазе записывается в виде F = $ = \sum\nolimits_{i = 1}^s {{{\mu }_{i}}{{N}_{i}}} $, где свободные ячейки (вакансии) являются частицами сорта i = s, (${{\mu }_{s}} = - P{{v}_{0}}$, т.к. Vs = = ${{{v}}_{0}}$N, ${{{v}}_{0}}$ – объем ячейки) и они отражают объем системы, не блокированный реальными молекулами сорта 1 ≤ i ≤ sc = s – 1; $N = \sum\nolimits_{i = 1}^s {{{N}_{i}}} $. Уравнение (1) перепишется (в нормированном виде на один узел системы) как

(3)
$F{\text{/}}N = \sum\limits_{i = 1}^s {{{\mu }_{i}}{{\theta }_{i}}} ,$
здесь ${{\theta }_{i}}$ – мольная доля частиц компонента i решеточной системы в однородной фазе. По аналогии с выражением (1) выражение (3) отвечает химическому потенциалу компонента решеточной структуры, где введен химпотенциал вакансий μs.

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

(4)
$F{\text{/}}N = \sum\limits_{q = 1}^\kappa {\sum\limits_{i = 1}^s {M_{q}^{i}(k)\theta _{q}^{i}} } ,$
где число узлов N относится к переходной области, состоящей из κ монослоев, 1 ≤ q ≤ κ. Величины Мi(k) в (4) характеризуют вклады компонентов i в свободную энергию объемной фазы и этих же компонентов в локально неоднородных областях q границы, через которые идет расчет ПН.

В зависимости от способа построения выражения (4) возможно три вида функции $M_{q}^{i}(k)$: k = 1 – перегруппировкой слагаемых выражения для F с парным потенциалом, k = 2 – дифференцированием по мольной доли частиц сорта i при фиксированном числе узлов типа q, или k = 3 – переменном числе узлов типа q [36]. Все виды функции $M_{q}^{i}(k)$ в (4) связаны с химическими потенциалами компонентов системы в (1) и (3), и анализ этих связей позволяет выявить связь функций $M_{q}^{i}(k)$ с ПН на основе термодинамического определения ПН [11, 12]. Для анализа сути этого термодинамического определения ПН достаточно ограничиться случаем бинарной смеси решеточной системы, в которой компонентами являются молекулы А и вакансии V, отвечающие чистому флюиду.

В данной работе проводится численный анализ полученных уравнений для ПН [11, 12, 35, 36] для плоской и сферической границ при учете парного потенциала взаимодействия между частицами с целью выяснения влияния вида функций $M_{q}^{i}$ на рассчитываемые значения ПН по формуле (2), где избыток свободной энергии ΔF находится по микроскопической модели как разность выражений (4) и (3).

Модель

В расчетах используется простейший вариант МРГ с учетом взаимодействия ближайших соседей в квазихимическом приближении (КХП) на жесткой решеточной структуре с числом соседей z [13–15]. Рассмотрим систему, состоящую из капли радиусом R с границей раздела пар–жидкость и окружающим паром как аналог равновесной двухфазной системы при данной температуре Т (плоская граница раздела фаз представлена предельным случаем R → ∞). Переходная область разделяется на мономолекулярные слои шириной λ однородные по своим свойствам (λ – среднее расстояние между молекулами в жидкой фазе). Эти слои нумеруются индексом q, где q – номер узла, относящийся к рассматриваемому монослою, 1 ≤ q ≤ κ, здесь κ – ширина переходной области плюс по одному монослою от объемных фаз (q = 1 отвечает жидкости и q = κ отвечает пару).

Структуру флюида в объемной фазе будем характеризовать набором величин $z_{{qp}}^{*}$, обозначающих числа ближайших соседних узлов слоя p вокруг узлов слоя q; $\sum\nolimits_{p = q - 1}^{q + 1} {z_{{qp}}^{*}} $ = z. Общий баланс узлов связей между соседними молекулами запишется в виде $\sum\nolimits_{p = q - 1}^{q + 1} {{{z}_{{qp}}}} $(R) = z. Для сферических капель в термодинамической версии модели структурные числа для искривленной решетки zqp(R) выражаются через аналогичные числа для плоской решетки $z_{{qp}}^{*}$ в виде поправок, зависящих от радиуса монослоя в переходной области [12, 28]

(5)
$\begin{gathered} {{z}_{{q < p}}}(R) = z_{{qp}}^{*},\quad {{z}_{{q > p}}}(R) = z_{{qp}}^{*}\left[ {1 + {\text{ }}1{\text{/}}(R + q - 1)} \right], \\ {{z}_{{q = p}}}(R) = z - {{z}_{{q < p}}}(R) - {{z}_{{q > p}}}(R). \\ \end{gathered} $

В асимптотическом пределе больших капель все значения zqp(R) стремятся к своим пределам $z_{{qp}}^{*}$ для плоской границы раздела.

Молекулярные распределения частиц сорта А (и, соответственно, вакансий V $\theta _{q}^{V} = 1 - \theta _{q}^{А}$) задаются плотностями $\theta _{q}^{А}$ частиц А в слое q, 1 ≤ q ≤ κ, которые описываются в КХП следующей системой уравнений

(6)
$\beta {{v}_{0}}{{a}_{q}}P = \theta _{q}^{A}{\text{/}}(1 - \theta _{q}^{A})\prod\limits_{p = q - 1}^{q + 1} {{{{\left[ {1 + t_{{qp}}^{{AA}}x_{{АА}}^{{}}} \right]}}^{{{{z}_{{qp}}}(R)}}}} ,$
где $t_{{qp}}^{{AA}}$ – условная вероятность нахождения частицы сорта А в ячейке слоя p рядом с другой частицей А в ячейке слоя q: $t_{{qp}}^{{AA}} = \theta _{{qp}}^{{AA}}{\text{/}}\theta _{q}^{A} = $ $ = 2\theta _{p}^{A}{\text{/}}[\delta _{{qp}}^{{AA}} + b_{{qp}}^{{AA}}]$, $\delta _{{qp}}^{{AA}} = 1 + x_{{AA}}^{{}}(1 - \theta _{q}^{A} - \theta _{p}^{A})$ $b_{{qp}}^{{AA}} = $ $ = {{\{ {{[\delta _{{qp}}^{{AA}}]}^{2}} + 4x_{{AA}}^{{}}\theta _{q}^{A}\theta _{p}^{A}\} }^{{1{\text{/}}2}}}$, $\theta _{{qp}}^{{AA}}$ – вероятность нахождения пары частиц АА на соседних ячейках монослоев q и p соответственно; P – давление в системе; ${{x}_{{AA}}} = \exp \{ - \beta {{\varepsilon }_{{AA}}}\} - 1$, $\beta = {{({{R}_{B}}T)}^{{ - 1}}}$, ${{R}_{B}}$ – газовая постоянная, ${{\varepsilon }_{{AA}}}$ – энергия взаимодействия пары частиц АА, описываемая потенциальной функцией Леннард-Джонса. Взаимодействия с вакансиями равны нулю ${{\varepsilon }_{{AV}}} = {{\varepsilon }_{{VA}}} = 0$.

Система уравнений (6) относительно локальных плотностей $\theta _{q}^{А}$ строится из условия равенства химического потенциала μА частиц А во всех слоях 1 ≤ q ≤ κ.

В МРГ вводятся одночастичные вклады $\nu _{q}^{i}$ в свободную энергию F компонента i на узлах типа q неоднородной системы с границей раздела фаз. Разность этих вкладов $\nu _{q}^{A} - \nu _{q}^{V} = {{\beta }^{{ - 1}}}\ln (\beta {{v}_{0}}{{a}_{q}}P)$ включает в себя статсуммы внутренних движений компонентов i, влияние внешних полей и химические потенциалы [15], где ${{a}_{q}} = {{F}_{q}}{\text{/}}{{F}_{0}}$ – отношение статсумм молекулы в решеточной структуре (Fq) и в объемной фазе (F0), формально для вакансий следует положить $\nu _{q}^{V} = 0$. Или величина aqP фиксирует значение химического потенциала вещества в разных слоях q (в уравнениях (6)).

Зная функцию $\theta _{{qp}}^{{AA}}$, через нормировочные соотношения $\sum\nolimits_{j = 1}^s {\theta _{{qp}}^{{ij}}} = \theta _{q}^{i}$ (s = 2) находятся остальные парные функции $\theta _{{qp}}^{{ij}}$. Размерность система уравнений (6) относительно локальных плотностей $\theta _{q}^{i}$ равна числу слоев (κ – 2) переходной области между паром и жидкостью. Она решается итерационным методом Ньютона при заданных значениях плотности пара для q = 1 и жидкости для q = κ. Точность решения этой системы не менее чем 0.1%. Плотности сосуществующих газовой и жидкой фаз в объеме и равновесное давление в системе P определялись с помощью построения Максвелла [12–15].

По концентрационному профилю (6) рассчитывается ПН по трем вариантам функции $M_{q}^{i}(k)$, k = 1–3. Согласно уравнениям [11, 12] термодинамическое определение ПН (2) выражается через указанные функции $M_{q}^{i}(k)$ как

(7)
$\begin{gathered} A\sigma = \frac{1}{{{{F}_{\rho }}}}\sum\limits_{i = 1}^s {\left( {\sum\limits_{q \leqslant q{\kern 1pt} *} {{{F}_{q}}(M_{q}^{i}(k) - M_{1}^{i}(k))\theta _{q}^{i}} + } \right.} \\ \, + \left. {\sum\limits_{q > q*} {{{F}_{q}}(M_{q}^{i}(k) - M_{\kappa }^{i}(k))\theta _{q}^{i}} } \right), \\ \end{gathered} $
где функции $M_{q}^{i}(k)$ определяются формулами [35], А – площадь поверхности ячейки решеточного газа,

(8)
$\begin{gathered} 1)\,\,M_{q}^{i}(1) = \nu _{q}^{i} + kT\ln \theta _{q}^{i} + \frac{{kT}}{2}\sum\limits_{p = q - 1}^{q + 1} {{{z}_{{qp}}}\ln \frac{{\hat {\theta }_{{qp}}^{{ii}}\hat {\theta }_{{qp}}^{{ik}}}}{{{{{(\theta _{q}^{i})}}^{2}}\theta _{{qp}}^{{ki}}}}} , \\ \hat {\theta }_{{qp}}^{{ij}} = \theta _{{qp}}^{{ij}}\exp ( - \beta {{\varepsilon }_{{AA}}}), \\ \end{gathered} $
(9)
$2)\,\,M_{q}^{i}(2) = \nu _{q}^{i} + kT\ln \theta _{q}^{i} + kT\sum\limits_{p = q - 1}^{q + 1} {{{z}_{{qp}}}\ln \frac{{\hat {\theta }_{{qp}}^{{iV}}}}{{\theta _{q}^{i}}}} ,$
(10)
$3)\,\,M_{q}^{i}(3) = \nu _{q}^{i} + kT\ln \theta _{q}^{i} + \frac{{kT}}{2}\sum\limits_{p = q - 1}^{q + 1} {{{z}_{{qp}}}\ln \frac{{\hat {\theta }_{{qp}}^{{ii}}}}{{\theta _{q}^{i}\theta _{p}^{i}}}} .$

Дополнительно результаты расчета ПН по термодинамическому определению сопоставляются с определением ПН, указанному в [10, 28], в котором учитывается, что вакансии отражают механические свойства системы, т.е. добавлен четвертый вариант определения ПН.

(11)
$4)\,\,A\sigma = \frac{1}{{{{F}_{\rho }}}}\left( {\sum\limits_{q \leqslant q*} {{{F}_{q}}(\mu _{q}^{V} - \mu _{1}^{V})} + \sum\limits_{q > q{\kern 1pt} *} {{{F}_{q}}(\mu _{q}^{V} - \mu _{\kappa }^{V})} } \right),$
здесь функции $\mu _{q}^{V} = M_{q}^{V}(3)$ определены в (10) с $\nu _{q}^{V} = 0$ для любого q.

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

(12)
$\sum\limits_{q \leqslant q*} {{{F}_{q}}(\theta _{q}^{A} - \theta _{1}^{A})} + \sum\limits_{q > q*} {{{F}_{q}}(\theta _{q}^{A} - \theta _{\kappa }^{A})} = 0.$

При q ≤ q* находятся слои с повышенной плотностью, при q > q* – слои с пониженной плотностью. Вклад каждого монослоя выражается через весовые функции Fq = Nq/N, N = $\sum\nolimits_{q = 2}^{k - 1} {{{N}_{q}}} $, 2 ≤ q ≤ ≤ κ – 1.

Условия расчета

Для моделирования объемной фазы флюида использовалась кубическая примитивная решетка с числом соседей в первой к.с. z = 6. Учитывались взаимодействия только между ближайшими соседями с потенциалом Леннарда-Джонса использовались параметры взаимодействия, дающие εАА = 238 кал/моль.

Изотермы для объемных свойств и профили поверхностных свойств построены для температур τ = 0.68 и 0.85 (τ = T/Tcr, где Tcr – критическая температура флюида в объемной фазе). Для простоты расчетов принято ${{a}_{q}} = 1$, означающее отсутствие внешнего потенциала.

Ширина границы раздела фаз (κ) и радиус равновесных капель (R) измеряется в единицах параметра решеточной структуры λ (λ = 1.12σ, σ – параметр потенциала Леннарда-Джонса) – среднее расстояние между молекулами в плотной фазе, или в числах монослоев. То есть κ и R являются безразмерными величинами.

Значения ПН выведены как произведение на площадь ячейки решеточного газа А в тех же единицах, что задается ${{\varepsilon }_{{AA}}}$, кал/моль. Функции $M_{q}^{i}(k)$ также выведены в единицах кал/моль.

Рассмотрены температурные зависимости ПН на плоской границе фаз в диапазоне температур от τ = 0.6 до 1 (до критической точки), а также размерные зависимости ПН от величины минимального размера R = R0 капли до R = 800 монослоев при τ = 0.68 и 0.85.

Результаты расчета

Объемная фаза. На рис. 1а показаны изотермические зависимости плотности компонента А, θA, и химического потенциала ln(a0P) в системе при τ = 0.68 (кривая 1) и 0.85 (кривая 2), здесь a0= β${{{v}}_{0}}$. Пунктиром показа петля Ван-дер-Ваальса.

Рис. 1.

Изотермические зависимости (а) химического потенциала от плотности А; давление a0P в фазах и температурные зависимости свойств сосуществующих фаз (б).

Реальная изотерма (сплошные кривые на рис. 1а) распадается на две ветви пара (слева при малых плотностях) и жидкости (справа при больших плотностях), отделенные друг от друга горизонтальным участком, изображающим процесс изотермического расслоения двухфазной системы при равновесном химическом потенциале ln(a0P). При этом плотности жидкости и газа остаются неизменными и равными значениям $\theta _{{{\text{liq}}}}^{{\text{A}}}$ и $\theta _{{{\text{vap}}}}^{{\text{A}}}$ соответственно. Плотности $\theta _{{{\text{liq}}}}^{{{\text{A(1)}}}}$ и $\theta _{{{\text{vap}}}}^{{{\text{A(1)}}}}$ отвечают кривой 1, а плотности $\theta _{{{\text{liq}}}}^{{{\text{A(2)}}}}$ и $\theta _{{{\text{vap}}}}^{{{\text{A(2)}}}}$ – кривой 2. С ростом температуры, согласно рис. 1а, увеличивается равновесный химический потенциал ln(a0P) и плотности сосуществующих фаз $\theta _{{{\text{liq}}}}^{{\text{A}}}$ и $\theta _{{{\text{vap}}}}^{{\text{A}}}$ сближаются. Участок петли, показанной пунктиром на рис. 1а, с отрицательной сжимаемостью не отвечает термодинамической устойчивости, т.е. точки, лежащие на таких участках, соответствуют неустойчивым состояниям вещества и такие состояния вещества не могут быть реализованы. Эти участки отражают суть правила Максвелла [12–15].

На рис. 1б показаны температурные зависимости плотностей сосуществующих фаз $\theta _{{{\text{liq}}}}^{{\text{A}}}$ (кривая 1) и $\theta _{{{\text{vap}}}}^{{\text{A}}}$ (2), значения которых отложены на оси ординат слева, и равновесного химического потенциала ln(a0P) (3), значения которого отложены на оси ординат справа. С ростом температуры химический потенциал в системе (3) растет, а плотности сосуществующих фаз (1 и 2) сближаются.

На рис. 2а показаны концентрационные профили компонента А в монослоях переходной области, 1 ≤ q ≤ κ, где q = 1 и κ – монослои от фаз жидкости и пара соответственно, на плоской границе раздела фаз (кривые 1–2) и для равновесной капли радиусом R = 200 (кривые 3, 4) при τ = 0.68 (1, 3) и 0.85 (2, 4). На рис. 2б показана температурная зависимость ширины переходной области κ с плоской границей раздела фаз.

Рис. 2.

Концентрационные зависимости (а) плотности А между сосуществующими фазами (пояснения в тексте); температурная зависимость ширины переходной области между фазами (б).

Кривые на рис. 2а показывают, что с увеличением температуры ширина переходной области плоской границы увеличивается: кривая 1 переходит в кривую 2, а кривая 3 в кривую 4, а плотности сосуществующих фаз сближаются.

Переход от кривой 1 к кривой 3 на рис. 2а (аналогично переход от кривой 2 к кривой 4) связан с переходом от плоской границы к сферической равновесной капле радиуса R = 200 – здесь наблюдается уменьшение ширины переходной области при неизменных плотностях сосуществующих фаз. Одновременно с уменьшением радиуса равновесной капли происходит смещение концентрационного профиля к жидкости, т.е. эквимолекулярная поверхность сдвигает от центра переходной области, в которой она располагается при плоской границе фаз, к капле.

Температурная зависимость ширины переходной области между фазами показана на рис. 2б. Ширина переходной области κ увеличивается с ростом температуры. Изменение ширины κ происходит дискретно по количеству монослоев, что отражает дискретную природу вещества на молекулярном уровне.

Функции Мi(k) в объемной фазе. Перед расчетом ПН сопоставим представленные функции (8)–(10) в объемной фазе. Кроме трех вариантов функций Мi(k), k = 1–3, представим еще один вид функций Мi(k = 2) [36], который обозначим как Мi(k = 2*). Он выражается как

(9*)
$2*)\,\,M_{q}^{i}(2) = \nu _{q}^{i} + kT\ln \theta _{q}^{i} + kT\sum\limits_{p = q - 1}^{q + 1} {{{z}_{{qp}}}\ln \frac{{\hat {\theta }_{{qp}}^{{iV}}}}{{\theta _{q}^{i}\theta _{q}^{V}}}} ,$

Его отличие от функций Мi(k = 2) связано с делением на плотность вакансий в логарифмическом слагаемом (при p = q). Вид данной функции становится близким по виду функций Мi(k = 3) (10). Для вакансий в функции МV(k = 2*) логарифмическое слагаемое совпадает по виду с аналогичным логарифмическим слагаемым в (10), исключением коэффициента ½ перед знаком логарифма.

На рис. 3 показаны температурные зависимости функций Mi(k) в жидкости $M_{{{\text{liq}}}}^{i}$(k) (кривые 1–4) и паре $M_{{{\text{vap}}}}^{i}$(k) (кривые 5–8) компонентов i = A (рис. 3а) и V (рис. 3б). Вариант k = 1 показан кривыми 1, 5, вариант k = 2 – кривыми 2, 6, вариант k = 2* – кривыми 3 и 7 и вариант k = 3 – кривыми 4, 8. На оси ординат слева отложены значения на кривых 1, 3–5 и 7, 8 для вариантов k = 1, 2* и 3, на оси ординат справа – на кривых 2 и 6 для варианта k = 2. Кривые, относящиеся к разным осям, заметно отличаются друг от друга.

Рис. 3.

Температурные зависимости функций Mi(k) компонентов i = A (а) и V (б).

Для вариантов k = 1 и 3 получаем, что значения функций Mi(k) совпадают в жидкости и в паре, $M_{{{\text{liq}}}}^{i}$(k) = $M_{{{\text{vap}}}}^{i}$(k), а также между данными вариантами, как для i = A (рис. 3а) и V (рис. 3б): кривые 1, 4, 5 и 8 совпадают на обоих полях рис. 3. Для вариантов k = 2 и 2* значения функций Mi(k) отличаются в жидкости и в паре, и с ростом температуры их значения сближаются: кривые 2 и 6 для k = 2 и кривые 3 и 7 для k = 2* сближаются с ростом температуры.

Химический потенциал. Разность введенных в (4) функций Мi(k) для i = A и V определяет химический потенциал в объемной фазе. Учитывая выписанные выше связи и условия расчета, имеем следующие выражения

(13)
$\begin{gathered} M_{q}^{A}(k)--M_{q}^{V}(k) = {{M}_{A}}(k)--{{M}_{V}}(k) = \\ = {{\beta }^{{ - 1}}}\ln (\beta {{v}_{0}}P) = {{\nu }_{A}} - {{\nu }_{V}}. \\ \end{gathered} $

Для всех четырех вариантов функций k = 1 (8), 2 (9), 3 (10) и 2* (9*). То есть все рассматриваемые версии разностей функций Мi(k) в объемной фазе являются эквивалентными, так как отражают общее состояние объемной фазы, никак не связанное с наличием границы между фазами. Этот факт имеет важное значение для сопоставления разных способов расчета ПН при полностью идентичных условиях в объемных фазах.

Уравнения на концентрационный профиль (6) имеют эти фазы в качестве граничных условий, и, следовательно, найденные решения, продемонстрированные на рис. 2, точно также являются идентичными для всех последующих вариантов расчета ПН. (Ниже, однако, отсутствуют расчеты ПН для варианта 2*, так как для него все величины оказываются отрицательными, что противоречит физическому смыслу ПН для процесса расслаивания.) Обсуждение варианта 2* обусловлено тем фактом, что для него, как и для всех других вариантов функций Mi(k) выполняется выражение (13). Это указывает на множественность разных вариантов построений функций Mi(k), которым отвечает постоянство химического потенциала в объемных фазах и постоянство концентрационного профиля.

Функции $M_{q}^{i}$(k) границ раздела фаз. Значения функций $M_{q}^{i}$(k) для разных монослоев q переходных областей границ раздела фаз показывают их относительные вклады в локальные значения свободной энергии. Как следует из вышеизложенного сам профиль (6) зафиксирован, а для расчета ПН, необходимо знание разностей этих функций для компонентов i = A и V. На рис. 4 показаны профили вкладов ($M_{q}^{i}$(k) – Mi(k)) компонентов i = A (кривые 1–3) и V (кривые 4–6) в значение ПН в переходной области на плоской решетке (рис. 4а, б) и для капли радиусом R = 200 (рис. 4в, г) при τ = 0.68 (рис. 4а, в) и 0.85 (рис. 4б, г). Вариант k = 1 показан кривыми 1, 4 с символами-квадратами, вариант k = 2 – кривыми 2, 5 с символами-кругами и вариант k = 3 – кривыми 3, 6 с символами-треугольниками. На оси ординат слева отложены значения на кривых 1, 3, 4 и 6 для вариантов k = 1 и 3, на оси ординат справа – на кривых 2 и 5 для варианта k = 2. Пунктиром показан уровень нулевых значений ($M_{q}^{i}$(k) – Mi(k)) на осях ординат слева и справа.

Рис. 4.

Профили вкладов в ПН. Обозначения см. текст.

На всех полях рис. 4 кривые 1 и 2, построенные для компонента А, совпадают с кривыми 4 и 5, построенными для V. Это является следствием выполнения тождества $M_{q}^{i}$(k) – $M_{q}^{j}$(k) = Mi(k) – ‒ Mj(k) = β–1ln(a0P) для вариантов k = 1 и 2. Для варианта k = 3 подобное тождество, как показано выше, имеет место только для объемных фаз Mi(3) – ‒ Mj(3) = β–1ln(a0P), но не для переходной области $M_{q}^{i}$(3) – $M_{q}^{j}$(3) ≠ β–1ln(a0P).

Кривые 1 и 4 по варианту k = 1, кривые 2, 5 по варианту k = 2 и кривая 6 по варианту k = 3 дают положительные вклады ($M_{q}^{i}$(k) – Mi(k)) со стороны жидкости и отрицательные – со стороны пара. Кривая 3 по варианту k = 3 единственная дает положительные вклады ($M_{q}^{i}$(k) – Mi(k)) со стороны пара, а отрицательные – со стороны жидкости.

По варианту k = 3 кривые 3 и 6 на плоской границе раздела фаз (рис. 4а и б) имеют симметричный вид относительно θ = 0.5, который также не сильно нарушается для капли с радиусом R = 200 (рис. 4в и г). По другим вариантам k = 1 и 2 кривые 1, 4 и 2, 5 соответственно не обладают той же симметрией вкладов со стороны жидкости и со стороны пара. Величины ПН (7) по определениям 1–3 рассчитываются как сумма положительных и отрицательных вкладов ($M_{q}^{i}$(k) – Mi(k)) от компонентов А (кривые 1–3) и V (кривые 4–6), взвешенные на локальные плотности компонента А и V соответственно. Отметим, что в силу тождеств $M_{q}^{A}$(k) – $M_{q}^{V}$(k) = MA(k) – MV(k) для вариантов k = 1 и 2, величины ПН могут быть также рассчитаны как сумма вкладов ($M_{q}^{i}$(k) – Mi(k)) только по одному из компонентов А или V (без взвешивания на плотность).

Аналогичное утверждение для варианта k = 3 несправедливо в силу неравенства $M_{q}^{A}$(3) – $M_{q}^{V}$(3) ≠ ≠ MA(3) – MV(3). Поэтому определение (7) через функции (10), являясь средневзвешенным от плотностей молекул А и вакансий, отличается от ПН через функции (11), содержащие только вклады вакансий, и это приводит к разным величинам ПН.

Плоская граница. На рис. 5 показаны температурные зависимости ПН по вариантам 1–4 (номера кривых соответствуют номерам вариантов). На оси ординат слева отложены значения на кривых 1, 3 и 4 для соответствующих определений ПН, на оси ординат справа – на кривой 2 для второго определения ПН.

Рис. 5.

Температурные зависимости ПН в рассматриваемых вариантах определений ПН.

На рис. 5 для всех определений получено уменьшение ПН с ростом T до нуля в критической точке (τ = 1). Кривые 1 и 4 по соответствующим определениям ПН имеют близкие значения, так что практически совпадают.

При температурах τ от 0.6 до 0.8 все кривые 1–4 имеют приблизительно линейный вид. Вблизи критической точки (τ = 1) определения 1, 3 и 4 (кривые 1, 3 и 4) дают положительное значение второй производной ПН от температуры T, а кривая 2 по определению 2 имеет отрицательный изгиб в этой области.

Граница капель. На рис. 6 показаны размерные зависимости ПН σ, нормированные на значение для плоской границы раздела фаз σbulk, по вариантам 1–4 (номера кривых соответствуют номерам вариантов) при τ = 0.85 (поле а), 0.68 (поле б) и 0.54 (поле в). Пунктиром на графиках показан уровень σ/σbulk = 1.

Рис. 6.

Размерные зависимости ПН в рассматриваемых вариантах определений ПН для трех температур: (а) τ = 0.85, (б) τ = 0.68, (в) τ = 0.54.

На поле 6а имеет место монотонное убывание ПН с уменьшением радиуса равновесной капли по всем 4 вариантам (кривые 1–4) от значения ПН σbulk на плоской границе раздела фаз до нуля при некотором минимальном радиусе R0 капли, отвечающей ее состоянию как фазы. Наибольшее R0 наблюдается для первого определения (кривая 1), затем идет R0 по варианту 4 (кривая 4), а наименьшее R0 дают варианты 2 и 3 (кривые 2 и 3 соответственно) – эти соотношения сохраняются для всех рассмотренных температур на рис. 6а–в.

При пониженных температурах (рис. 6 б,в) наблюдается скачкообразное изменение ПН с уменьшением R. Это является следствием дискретного изменения ширины переходной области κ, т.е. дискретной природы вещества, которая явно проявляется при малых размерах фазы для низких температур. В вариантах ПН по определениям 1, 2 и 4 сохраняется монотонное уменьшение ПН с уменьшением R (кривые 1, 2 и 4 соответственно). Вариант 3 при низких температурах (кривая 3 на рис. 6б) выше некого R* дает рост ПН с уменьшением R от значения ПН σbulk на плоской границе раздела фаз и выше. Ниже R* на кривой 3 на рис. 6б в результате скачкообразного уменьшения ПН приобретает значение ниже σbulk, после чего начинает монотонно убывать с уменьшением R.

При температурах вблизи тройной точки (рис. 6в) относительное расположение кривых меняется и уже кривая 1 выше некого R* дает рост ПН с уменьшением R (его максимальное значение превышает линию σ/σbulk = 1 примерно на 5%) и ниже R* на кривой 1 на рис. 6в в результате скачкообразного уменьшения ПН опускается до нуля. Остальные кривые 2–4 дают монотонное уменьшение ПН с уменьшением R.

Таким образом, все рассмотренные определения ПН представляют собой сумму локальных избыточных величин $M_{q}^{i}$(k) в монослоях переходной области жидкость–пар по сравнению со значениями в фазах, которые синусоидальным образом изменяются от жидкости к пару, образуя минимум в области отрицательных значений вкладов в ПН и максимум – в области положительных.

Расчеты по термодинамическому определению ПН (7), основанные на трех функциях $M_{q}^{i}$(k), k = 1–3, показали разное поведение как для плоской границы, так и для искривленной границы капли. Это указывает на неоднозначность термодинамического определения ПН (7), так как в объемной фазе все виды функций Mi(k) однозначно связаны с величинами химического потенциала молекул (13). Сравнения температурных зависимостей ПН и размерных зависимостей ПН при τ = const для разных вариантов определений ПН проведены на основе общего концентрационного профиля молекул между фазами.

Для плоской границы (рис. 5) кривые 1, 2 и 3 существенно отличаются друг от друга, хотя для всех определений получено уменьшение ПН с ростом T до нуля в критической точке (τ = 1). Сопоставление этих кривых с определением ПН (11) на основе представлений о дополнительном учете специфики вакансий в МРГ как о механической характеристике системы показывает, что кривые 1 и 4 практически совпадают, тогда как кривые 2 и 3 сильно отличаются от кривой 4. Причем кривая 2 дает качественное отличие в температурной зависимости от других вариантов, состоящее в отрицательном изгибе (вторая производная по температуре меньше нуля) кривой зависимости вблизи критической точки.

Размерные зависимости рассмотренных определений ПН (рис. 6) расходятся между собой для всех вариантов, хотя их общей тенденцией является монотонное убывание нормированной величины ПН от общего объемного значения с убыванием радиуса капли R до минимального размера фазы R0 за исключением кривой 3 при τ = 0.68 и кривой 1 для τ = 0.54. Малые превышения линии σ(R)/σbulk = 1 значениями ПН на этих кривых (до нескольких процентов) связаны с резким изменением профиля при малых радиусах капель в рамках МРГ. Учет мягкости решеточной структуры уменьшает максимальные значения σ(R)/σbulk в окрестности линии σ(R)/σbulk = 1 [37]. Сопоставление зависимостей для трех видов ПН, рассчитанных по термодинамическому определению, с кривой 4 (формула (11)) указывает на сильную зависимость относительного положения кривой 4 и размерных кривых σ(R)/σbul с разными функциями $M_{q}^{i}$(k), k = 1–3, в зависимости от температуры.

Таким образом, проведенный анализ разных определений ПН, обусловленных использованием разных функции $M_{q}^{i}$(k), в МРГ показал неоднозначность определений ПН (k = 1, 2, 2* и 3) по термодинамическому определению ПН при одинаковых фазовых состояниях сосуществующих фаз – пара и жидкости, и одинаковом концентрационном профиле между этими фазами. Идентичным условиям состояния системы отвечают разные функции $M_{q}^{i}$(k), k = 1, 2, 2* и 3, отражающие локальные парциальные вклады компонентов решеточной системы i на узлах типа q, в выражении для избытка свободной энергии ΔF на базе МРГ. При этом для функций $M_{q}^{i}$(k), k = 1, 2 и 2*, выполняются равенства $\mu _{q}^{A}$(k) = $M_{q}^{A}$(k) – $M_{q}^{V}$(k) = = μА, а для функций $M_{q}^{i}$(3) это равенство не выполняется (выше отмечено, что для варианта 2* все величины ПН имеют нефизические отрицательные значения, хотя $\mu _{q}^{A}$(2*) = $M_{q}^{A}$(2*) – ‒ $M_{q}^{V}$(2*) = μА).

Указанная неоднозначность определений ПН означает, что ПН не является чисто равновесной термодинамической функцией, для которой задание внешних условий однозначно определяет внутренние характеристики системы.

Работа выполнена в рамках государственного задания ИОНХ РАН в области фундаментальных научных исследований (№ 44.2).

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

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

  2. Оно С., Кондо С. Молекулярная теория поверхностного натяжения. М.: Изд-во иностр. лит., 1963. 292 с.

  3. Русанов А.И. Фазовые paвновесия и поверхностные явления. Л.: Химия, 1967. 387 с.

  4. Адамсон А. Физическая химия поверхностей. М.: Мир, 1979. 568 с.

  5. Джейкок М., Парфит Дж. Химия поверхностей раздела фаз. М.: Мир, 1984. 269 с.

  6. Роулинсон Дж., Уидом Б. Молекулярная теория капиллярности. М.: Мир, 1986. [J. Rowlinson and B. Widom, Molecular Theory of Capillarity (Oxford Univ., Oxford, UK, 1978).]

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

  8. Базаров И.П. Термодинамика. М.: Высш. школа, 1991. 376 с.

  9. Воронин Г.Ф. Основы термодинамики. М.: изд-во МГУ, 1987. 192 с.

  10. Товбин Ю.К. Малые системы и основы термодинамики. М.: Физматлит, 2018. 408 с.

  11. Товбин Ю.К. // Журн. физ. химии. 2018 Т. 92. № 12. С. 1902.

  12. Товбин Ю.К. // Там же. 2019. Т. 93. № 9. С. 1311.

  13. Хилл Т. Статистическая механика. М.: Изд-во иностр. лит., 1960. 486 с. (Hill T.L. Statistical Mechanics. Principles and Selected Applications. N.Y.: McGraw–Hill Book Comp.Inc., 1956.)

  14. Хуанг К. Статистическая механика. М.: Мир, 1966. 520 с.

  15. Товбин Ю.К. Теория физико-химических процессов на границе газ–твердое тело, М.: Наука, 1990. 288 с.

  16. Onsager L. // Phys Rev. 1944. V. 65. P. 117.

  17. Domb C. // Adv. Phys. 1960. V. 9. P. 149.

  18. Стенли Г. Фазовые переходы и критические явления. М.: Мир. 1973. 400 с.

  19. Паташинский А.З., Покровский В.П. Флуктуационная теория фазовых переходов. М.: Наука, 1975. 256 с.

  20. Ма Ш. Современная теория критических явлений. М.: Мир, 1980. с.

  21. Глазов В.М., Павлова Л.М. Химическая термодинамика и фазовые равновесия. М.: Металлургия, 1981. 336 с.

  22. Ono S. // Mem. Fac. Eng. Kyusgu Univ. 1947. V. 10. P. 195.

  23. Lane J.E. // Austr. J. Chem. 1968. V. 21. P. 827.

  24. Пиотровская Е.М., Смирнова Н.А. // Коллоидн. журн. 1979. Т. 41. С. 1134.

  25. Товбин Ю.К. // Там же. 1983. Т. 45. № 3. С. 707.

  26. Окунев Б.Н., Каминский В.А., Товбин Ю.К. // Там же. 1985. Т. 47. № 6. С. 1110.

  27. Смирнова Н.А. Молекулярные растворы. Л.: Химия, 1987. 334 с.

  28. Товбин Ю.К., Рабинович А.Б. // Изв. АН. Сер. химическая. 2010. № 4. С. 663.

  29. Товбин Ю.К., Рабинович А. Б. // Там же. 2009. № 11. С. 2127.

  30. Товбин Ю.К., Зайцева Е.С., Рабинович А.Б. // Журн. физ. химии. 2017. Т. 91. № 10. С. 1730.

  31. Товбин Ю.К., Зайцева Е.С. // Теплофизика высоких температур. 2018. Т. 56. № 3. С. 381 с.

  32. Товбин Ю.К., Рабинович А.Б. // Докл. АН. Сер. физическая химия. 2008. Т. 422. № 1. С. 59.

  33. Товбин Ю.К. Молекулярная теория адсорбции в пористых телах. М.: Физматлит, 2012. 624 с.

  34. Зайцева Е.С., Товбин Ю.К. // Журн. физ. химии. 2020, Т. 94. № 12. С. 1889.

  35. Товбин Ю.К. // Там же. 1992. Т. 66. № 5. С. 1395.

  36. Товбин Ю.К. // Там же. 2016. Т. 90. № 7. С. 1059.

  37. Товбин Ю.К., Зайцева Е.С. // Там же. 2019. Т. 93. № 9. С. 1437.

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

Инструменты

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

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

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