Журнал физической химии, 2021, T. 95, № 4, стр. 655-660

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

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

Обсуждается подход с позиции молекулярно-кинетической теории, описывающей объемные состояния трехагрегатных систем, а также их границ, на базе дискретной модели решеточного газа (МРГ) [13–16]. Для простоты принято, что компоненты трех агрегатных состояний имеют соизмеримые по размеру компоненты, и теория учитывает межмолекулярные взаимодействия в квазихимическом приближении (КХП), описывающем прямые корреляций. Теория является единственной, которая гарантирует равноточное описание как самих трех агрегатных состояний, так и их границ раздела фаз.

Ранее получено [15], что совместное действие поверхностного потенциала и межмолекулярных взаимодействий формирует подобласти общей пористой системы, в которых реализуются фазовые переходы первого рода [13, 14, 17]. Поверхностный потенциал стенок стремится стратифицировать молекулы адсорбата по двумерным подобластям общего пространства, чтобы сформировать фазы примерно одинаковой плотности в этих двумерных подобластях с одинаковым потенциалом. Диаграммы расслаивания в микронеоднородных системах характеризуются множественностью фазовых переходов первого рода пар–жидкость [15, 18]. Каждой кривой расслаивания отвечает своя критическая точка. Вид и число кривых расслаивания, а также параметры критических точек адсорбционных систем определяются поверхностным потенциалом, зависящим от состава и строения неоднородных поверхностей адсорбента [15, 18]. Это относится к адсорбции, как на открытых неоднородных поверхностях, так и внутри пористых тел.

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

Первый этап. В работе [19] была предложена двухуровневая структурная модель для пористых систем: надмолекулярный уровень отражал тип пор и их распределение по размерам, включая учет неоднородности стенок пор в сочленениях пор разного типа, а также их связанность через парные структурные ФР. Масштаб L₂ относится к надмолекулярному уровню (λ = L₁ ≪ L₂ ≪ L₃, где λ – размер адсорбата, L₃ – размер образца). Надмолекулярный уровень включает в себя участки пористого тела с характерным размером Н определенной геометрии (щелевидные, цилиндрические и сферические), где Н – ширина для щелей или диаметр для сфер и цилиндров, либо элементарные объемы структуры твердого тела в случае глобулярных систем (тогда геометрия пор задается через характеристики глобул). Размеры моделей щелей или цилиндров не могут быть бесконечными, и требуется явное задание конечной протяженности по разным осям. В общем случае H ≤ L₂. При H = L пора занимает весь рассматриваемый участок. При H = 0 пора отсутствует, что позволяет отразить наличие тупиковых пор, соседствующих с данным участком q.

Зададим надмолекулярную структуру ФР: F_q – характеризует долю участков типа q, и F_qp – характеризует вероятность нахождения рядом с участком типа q участка типа p, $\sum\nolimits_p^{} {{{F}_{{qp}}}} $ = F_q, 1 ≤ q, p ≤ t, где t – число типов рассматриваемых участков пористого тела. Также введем функцию Н_qp – условную вероятность нахождения участка поры типа p около участка поры типа q (в некотором выбранном направлении), F_qp = F_qH_qp, причем $\sum\nolimits_p^{} {{{H}_{{qp}}}} $ = 1. Аналогично для молекулярного уровня можно составить балансовые связи для внутренних типов узлов  f и g, где ${{F}_{{q,f}}}$ – доля узлов типа  f для поры типа q, $\sum\nolimits_f^{t(q)} {{{F}_{{q,f}}}} = {{F}_{q}}$, где t(q) – число типов узлов внутри участка q; ${{F}_{{q,f;p,g}}}$ – парная ФР для узлов q, f  и p, g. Их отношение определяет условные вероятности ${{H}_{{q,f;p,g}}}_{~} = {{F}_{{q,f;p,g}}}{\text{/}}{{F}_{{q,f}}}$ нахождения узла p, g рядом с узлом q, f, ${{\sum\nolimits_{g = 1}^{t(p)} {{{H}_{{q,f;p,g}}}} }_{~}} = {{H}_{{qp}}}$.

Будем для конкретности рассматривать ограниченные поры с выделенными длинными осями: прямоугольные параллелепипеды (сокращено “щель”) и цилиндры. Обозначим общую длину поры через L_q, где q – символ типа поры. Пусть центральная ось поры обозначается символом ξ. Координаты узлов вдоль оси ξ меняются от 1 до L_q, 1 ≤ ξ ≤ L_q , и символ ξ отражает сечение поры. Символ узла  f является векторным (кроме символа ξ): узел имеет две координаты для “щели” (f = h и ψ, где h – номер монослоя в одном из направлений сечения поры, 1 ≤ h ≤ H_q, перпендикулярном длинной центральной оси ξ, ψ – номер монослоя в другом направлении сечения поры, 1 ≤  ≤ ψ ≤ ψ_q) и для цилиндра (f = R и u, где R – радиус цилиндрического монослоя, u – номер оси параллельной главной оси поры, находящейся в этом монослое, 0 ≤ R ≤ R_q). Для глобулярных структур символ f отражает специфику этих геометрий [15].

Второй этап – описание молекулярных распределений по неоднородным узлам общей системы в МРГ отражает собственный объем молекул, исключая двойное заполнение узла разными молекулами, и взаимодействия между ними в КХП. Для каждого узла типа  f, 1 ≤ f ≤ t(q), участка q имеется своя вероятность заполнения ${{\theta }_{{q,f}}}$. Узел q, f  характеризуется локальной константой Генри ${{a}_{q}}_{{,f}}$ = $a_{{q,f}}^{0}\exp \beta {{Q}_{{q,f}}})$, β = (${{k}_{{\text{B}}}}$T)^–1, ${{Q}_{{q,f}}}$ – энергия связи молекулы со стенками поры, она выражается через параметр взаимодействия адсорбат–адсорбент, ${{k}_{{\text{B}}}}$ – константа Больцмана. Расчет локальных заполнений узлов разных групп ${{\theta }_{{q,f}}}$ с учетом энергетической неоднородности узлов решетки и взаимодействия между ближайшими соседями проводится по формулам [15, 18]:

где Р_А – давление адсорбтива, функция ${{\Lambda }_{{q,f}}}$ учитывает неидеальность адсорбционной системы в КХП. Индекс p, g пробегает всех соседей ${{z}_{{q,f}}}$ узла q, f внутри поры, ${{d}_{{q,f;p,g}}} = {{z}_{{q,f}}}{{H}_{{q,f;p,g}}}$ при реализации слоевых распределения узлов разного типа. Функции $t_{{q,f;p,g}}^{{AA}}$ – условная вероятность нахождения частицы А в узле и p, g рядом с частицей А в узле q, f, рассчитанная в КХП. Дополнительно учитываются следующие нормировочные соотношения: $\theta _{{qp,fg}}^{{AA}}$ + $q_{{qp,fg}}^{{AV}}$ = ${{\theta }_{{q,f}}}$ ≡ $\theta _{{q,f}}^{A}$ и $\theta _{{qp,fg}}^{{VA}}$ + $\theta _{{qp,fg}}^{{VV}}$ = $\theta _{{q,f}}^{V}$ = = 1 – ${{\theta }_{{q,f}}}$ (индекс V означает вакансию – свободный узел); а также $\theta _{{q,f}}^{A}$ + $\theta _{{q,f}}^{V}$ = 1. Равновесное распределение частиц по узлам разного типа ${{\theta }_{{q,f}}}$ находилось из системы уравнений (1), (2) задавая θ итерационным методом Ньютона (точность решения этой системы не менее чем 0.1%).

Нахождение функции ${{\theta }_{{q,f}}}$ позволяет рассчитать свободную энергию системы:

где ${{F}_{{q,f}}}({v})$ – локальная свободная энергия ячейки ${{{v}}_{{q,f}}}$, ее слагаемые описывают решеточные вклады F_lat и вклады колебательного F_vib и поступательного F_tr движения компонентов в каждом из слоев неоднородной области q, f; $F_{{q,f}}^{{{\text{lat}}}} = ~\,\,E_{{q,f}}^{{{\text{pot}}}} - TS_{{q,f}}^{{{\text{lat}}}}~$, где формулы для потенциальной энергии и энтропии неоднородных слоев, а также $F_{{q,f}}^{{{\text{vib}}}}$ и $F_{{q,f}}^{{{\text{tr}}}}$ приведены в работах [20, 21].

Третий этап – расчет полной изотермы адсорбции путем взвешивания локальных изотерм с помощью двухуровневых структурных ФР узлов разного типа, введенных на этапе 1. Этот этап необходим для сравнения с экспериментально наблюдаемой изотермой адсорбции. Изотерма адсорбции θ(P_А) запишется как

Четвертый этап связан с выявлением областей конденсации адсорбата внутри разных участков пор и нахождением сосуществующих плотностей пара и жидкости при данной температуре. Анализируются поры всех типов q, f – в каждой из них на локальных изотермах проверяется наличие петли типа ван-дер-Ваальса, и для нее строится своя секущая [10–12]. Оценки минимальных размеров пор некоторых простых геометрий, в которых реализуется расслаивание адсорбата, рассмотрены в [22]. Здесь также опущены вопросы о введении калибровочных функций двух видов (для двух- или трехмерных решеток и размерных эффектов) для повышения точности описания фазового состояния систем [15].

Средние плотности сосуществующих паровой $\theta _{{q,f}}^{{{\text{vap}}}}$ и жидкой $\theta _{{q,f}}^{{{\text{liq}}}}$ фаз адсорбата в порах q, f определяются с помощью построения Максвелла [13–15]. Среднему значению плотности пара соответствуют множество величин локальных плотностей $\theta _{{q,f}}^{{{\text{vap}}}}(h,\psi )$. Для конкретности рассматриваем “щель”, в которой выделяется ее площадь перпендикулярная оси, содержащая число узлов M(q, h, ψ) в данном сечении поры ξ (аналогично для жидкости). Этим выделяются двухфазные локальные области $\theta _{{q,f}}^{{{\text{vap}}}} \leqslant \theta _{{q,f}}^{{}} \leqslant \theta _{{q,f}}^{{{\text{liq}}}}$. Из общего множества пор отделяются поры, в которых реализуются мениски пар–жидкость (это формирует множества q, f*), или внутри пор будут однофазные твердое–пар или твердое–жидкость системы. Сами стенки пор могут быть сложными: неоднородными по составу поверхности и ее структуре (шероховатость). Если области по θ соизмеримы с условиями заполнения пор разного размера, то купола кривых расслаивания будут искривлены [15].

На пятом этапе для каждой петли ван-дер-Ваальса из множества q, f* проводится расчет переходных областей между выявленными фазами пар–жидкость внутри данного типа пор $\theta _{{q,f}}^{{{\text{vap}}}} \leqslant \theta _{{q,f}}^{{}} \leqslant \theta _{{q,f}}^{{{\text{liq}}}}$. Этот этап требует расширения на единицу размерности структурных ФР этапа 1, чтобы отразить переходную область мениска между двумя мобильными фазами внутри рассматриваемой поры q, f вдоль длинной оси ξ, т.е. выделим сечения переходной области, нумеруемые от 1 до L. Решение системы уравнений для изотерм второго этапа (1), (2) на расширенном множестве узлов внутри переходной области осуществляется при фиксированных локальных плотностях сосуществующих фаз {$\theta _{{q,f}}^{{{\text{vap}}}}(1)$} и {$\theta _{{q,f}}^{{{\text{liq}}}}(L)$} из пункта 4. Это решение формирует поле плотностей в промежуточной области $\theta _{{q,f}}^{{}}(\xi ) \equiv {{\theta }_{{q,h,\psi }}}(\xi )$, где символ оси ξ нумерует сечения поры вдоль ее стенок от жидкости ξ = 1 к пару ξ = L. Величина L находится итерационным путем последовательного ее увеличения при условии ΔF/ΔL → 0, где ΔF – приращение свободной энергии (3) во всей переходной области жидкость–пар, и ΔL = 1 (так как L дискретная величина). Уравнения для изотерм (1) и (2) обеспечивают постоянство химического потенциала во всей системе, исключающие какие-либо ее метастабильные состояния.

Шестой этап включает в себя расчет локального ПН менисков пар–жидкость в каждой точке ($q,h,\psi (\xi )$) переходной области мениска по найденным величинам концентрационных полей $\theta _{{q,f}}^{{}}(\xi )$ пункта 5 [23–26]. Локальное ПН определяется через локальные избыточные значения свободной энергии (3), как в макросистемах [27],

где А₁ – площадь внутри ячейки, к которой относится рассчитываемая величина локального ПН [23], ρ_е – положение эквимолекулярной поверхности в сечении ξ для оси ψ параллельной центральной оси поры q, f. Она находится из уравнения материального баланса вдоль оси ψ: где функция ${{F}_{{q,h,\psi }}}(\xi )$ является обобщением функции этапа 1 для описания распределенной модели ячеек вдоль оси ψ по узлам сечений ξ для группы узлов f поры q.

Величина локального давления ${{\pi }_{{q,h,\psi }}}(\xi )$ в (5) для сечения ξ точки $q,h,\psi (\xi )$ выражается как

где символ g пробегает всех ближайших соседей узла h, ψ из сечения ξ в поре q по сечениям ξ и ξ'; $\theta _{{q,h,\psi ;p,g,\psi '}}^{{VV}}(\xi \xi {\kern 1pt} ')$ – парная функция для узлов q, h, ψ, ξ и p, g, ψ', ξ'.

Локальные ПН (5), зависящие от символов $q,h,\psi (\xi )$, построены на общем множестве концентрационных плотностей $\theta _{{q,f}}^{{}}(\xi )$. Из них можно составить выражения для традиционных простых геометрий границ капель или площади плоской границы (см. [23–26]).

Сочленения. Положение мениска внутри пор зависит от размера пор и свойств их стенок [15]. Выше рассмотрен случай мениска вдали от их обоих концов длинных ограниченных по размеру пор. Наиболее распространенным является случай, когда мениск существует на границе между соседними порами разной ширины в окрестности их сочленения. Этот случай является характерным в системе пор разного размера: величина давления насыщенного пара зависит от ширины поры [1, 5, 10, 11, 15] – сначала заполняются наиболее узкие поры, а затем широкие. Для участков сочленений вводятся свои области (число разных типов сочленений обозначим через t(qp)). В них следует формировать свои типы неоднородных типов узлов, как например, между цилиндрами и сферами в системе сфероцилиндров [15]. Однако из-за их малого размера в них не может существовать “собственное” расслаивание фаз, поэтому мениск в них является следствием условий конденсаций между соседними порами.

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

Также возможно наличие мениска между изгибами длинных пор примерно одинакового диаметра из-за изменения структуры стенок (и их коэффициентов Генри) около мест изгибов [28–30]: специфику сочленений отражают области изгибов пор.

Обсуждение. Задача расчета ПН для менисков пар–жидкость состоит из двух частей: расчет локальных распределений адсорбата (этапы 1–3) и собственно расчет локальных равновесных величин ПН (этапы 4–6). Предложенная методика дает способ расчета локальных ПН σ_loc для менисков пар–жидкость в пористых системах на базе молекулярных распределений адсорбата, получаемых из микроскопической теории (они ищутся на базе молекулярных потенциалов из решения уравнений (1)). Эти распределения формируют на микроуровне границу раздела фаз. Информация о ПН дает представление о пространственном распределении фаз в пористых системах. Методика позволяет получить среднюю величину ПН σ, которая нужна для сравнения с экспериментом, и ее локальные изменения σ_loc в зависимости от физических свойств системы. Отличие теоретической информации от термодинамической заключается в ее более детальном характере, отражающем размерные и энергетические особенности системы.

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

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

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

Помимо указанных трех видов ПН на открытых поверхностях активно обсуждается вопрос о расчете линейного натяжения между тремя агрегатными состояниями [33–36], чтобы уточнить описание свойств трехагрегатных систем. Такая же постановка задачи, естественно, возникает и для менисков в поре. Вопрос о путях выделения линейного натяжения для менисков в порах требует своей разработки.

В заключение отметим пути повышения точности моделирования: 1) Приближенный учет корреляций в КХП (2) завышает значения критической температуры. Для ее уточнения в ограниченных по размерам порах необходимо использовать калибровочные функции с учетом размерных эффектов [15]. 2) Система уравнений (1) и (2) при использовании мягкой структуры МРГ должна быть дополнена уравнениями на изменения длин связей, отражающими локальные механические равновесия [24, 25, 37]. Это меняет решение этой системы, но сохраняет данную методику на остальных этапах.

Работа выполнена при поддержке Российского фонда фундаментальных исследований (код проекта № 18-03-00030а).