Аннотация
Теория функционала плотности (DFT) активно используется и развивается
в последнее время. DFT является важным инструментом для описания широкого круга явления на наномасштабах: смачивание, капиллярная конденсация, адсорбция и многое другое. В данной работе рассматривается подход, позволяющий использовать теорию функционала плотности для расчета термодинамических свойств флюида в нанопорах без вычисления вариации свободной энергии
(Variation Free Density Functional Theory — VF-DFT). Такой подход позволяет исследовать флюиды со сложным типом взаимодействия и ускорить вычисление для
простых флюидов. В безвариационном методе равновесная плотность флюида в
поре представляется в виде разложения по ограниченному набору базисных функций. Коэффициенты разложения функции плотности ищутся с помощью стохастических методов оптимизации (генетический алгоритм — GA, метод роя частиц
— PSO) так, чтобы минимизировать свободную энергию системы. В работе рассматриваются два флюида: азот при температуре 77.4 K и аргон при температуре
87.3 K в поре 3.6 нм и проводится сравнение алгоритмов оптимизации. Также
рассматривается гибридный подход (Hybrid Density Functional Theory — H-DFT)
на основе стохастических методов оптимизации и метода простой итерации для
поиска равновесной плотности флюида в поре. Такое объединение методов помогает значительно сократить время поиска решения без потери качества. Также
рассматривается задача восстановления распределения пор по размерам на основе
экспериментальных данных по низкотемпературной адсорбции методами теории
функционала плотности. В работе приводятся результаты по решению модельных
задач восстановления распределения пор по размерам, рассматриваются случаи
унимодальных и бимодальных гауссовских распределения с различными значениями параметров дисперсии и математического ожидания.
2
Содержание
1 Введение
4
2 Теория функционала плотности
2.1 Короткодействующее отталкивание
2.2 Диполь-дипольное взаимодействие
2.3 Метод простой итерации (МПИ) . .
2.4 Геометрия поры . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
8
9
13
15
16
3 Безвариационный подход в теории функционала плотности (VF-DFT)
3.1 Построение базисных функций . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Стохастическая оптимизация . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Генетический алгоритм (ГА) . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Метод роя частиц (МРЧ) . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Гибридный подход (H-DFT) . . . . . . . . . . . . . . . . . . . . . . . . . .
18
19
22
23
24
25
4 Результаты работы безвариационного алгоритма
4.1 Азот . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Высокое давление . . . . . . . . . . . . . . .
4.1.2 Низкое давление . . . . . . . . . . . . . . . .
4.2 Аргон . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Высокое давление . . . . . . . . . . . . . . .
4.2.2 Среднее давление . . . . . . . . . . . . . . .
.
.
.
.
.
.
25
26
26
29
31
31
32
.
.
.
.
.
.
34
34
38
38
41
44
46
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5 Восстановление распределения пор по размерам
5.1 Введение в задачу восстановления распределения пор по размерам . .
5.2 Результаты восстановления распределения пор по размерам . . . . . .
5.2.1 Влияние параметров распределения на вид изотермы . . . . .
5.2.2 Зависимость решения от параметр регуляризации . . . . . . .
5.2.3 Восстановление сдвинутого и бимодального распределения . .
5.2.4 Зависимость решения от качества экспериментальных данных
.
.
.
.
.
.
.
.
.
.
.
.
6 Заключение
48
Список Литературы
49
3
1
Введение
Задача исследования флюидов на масштабах от микро- до нанометров становится все более актуальной как для нефтегазовой индустрии (сланцевые углеводороды),
так и для других технологических отраслей, в которых необходимо производить сепарацию флюидов, улавливание и хранение углекислого газа, исследовать различного
рода течения флюидов в природных или искусственных наноразмерных структурах
[1—4]. В частности, для нефтегазовой индустрии, в нетрадиционных коллекторах большая часть порового пространства представляет собой структуры размера нескольких
десятков нанометров. Однако, на таких масштабах классические подходы имеют ограничения, становится необходимым более точно учитывать силы межмолекулярного взаимодействия, чтобы корректно описывать и учитывать такие физические явления, как:
адсорбция [5; 6], образование двойного электрического слоя [4; 7], капиллярная конденсация [4], смачивание [4; 8], исследовать гетерогенность стенок коллектора [9—12].
Корректное описание этих явлений позволяет изучать структуру материалов на наномасштабе, поведение флюидов и их взаимодействие как между собой, так и с твердой
средой.
С открытием высоко упорядоченный мезопористых материалов начали развиваться методы описания поведения флюида в нанопорах и характеризации порового
пространства. Традиционные методы, такие как BET [13], BJH [14] не подошли для таких задач так, как они не способны учесть различную морфологию порового пространства, учесть влияние микропористости и предсказать распределение пор по размерам,
которые могли бы быть независимо определены с помощью рентгеновской дифракции
(XRD) и просвечивающей электронной микроскопии (TEM) с точностью, недоступной
ранее. Новые нанопористые материалы и новые экспериментальные подходы высокого
разрешения требуют новых адекватных теоретических методов для анализа экспериментальных данных.
Одним из наиболее распространенных инструментов для описания флюидов и
их взаимодействий с поверхностями на таких масштабах является Теория Функционала Плотности (DFT). DFT предоставляет собой компромисс между классическими полуэмпирическими методами и молекулярным моделированием. С одной стороны, DFT
способна учитывать микроскопическую структуру макроскопической системы при относительно низких вычислительных затратах, а с другой, DFT — более строгая теория,
чем классические феноменологические подходы. Несмотря на то, что DFT, как инструмент исследования термодинамических свойств флюидов в классических системах, был
успешно применен в 1976 году [15], теория функционала плотности только в последнее
время стала активно использоваться и развивается для изучение равновесия и кинетики фазовых переходов, свойств полимерных материалов, тонких пленок, биологических
систем.
4
В 1989 году, теория функционала плотности была впервые применена для расчета распределения пор по размерам из данных по низкотемпературной изотерме адсорбции [16], и вскоре было признано, что DFT обеспечивает более обоснованный и
универсальный подход к расчету параметров структуры пор по сравнению с традиционными методами, основанными на уравнении Кельвина [17].
С математической точки зрения, молекулярная теория функционала плотности
напоминает DFT из квантовой механики (оба имеют один и тот же акроним) за исключением того, что в первом случае функционал плотности относится к строению
атомов или элементам полимерной молекулы, тогда как во втором случае идет речь
об электронах. Молекулярная теория функционала плотности опирается на теорему о
том, что система с определенной температурой 𝑇 , объемом 𝑉 , химическим потенциалом каждого из компонент флюида 𝜇𝑖 , внешним потенциалом для каждой компоненты
𝑉𝑒𝑥𝑡 (𝑟) однозначно определяется равновесной плотностью распределения молекул 𝜌 (𝑟)
[18]. Следствием этой теоремы является то, что для энергия Гельмгольца для системы
может быть записана как функционал плотности 𝐹 [𝜌 (𝑟)] и в равновесии свободная энерия минимальна. Математическую основу для DFT в квантово-механических терминах
сформулировали Хоэнберг и Кон [18].
Таким образом, для исследования сложных систем методами теории функционала плотности необходимо две вещи. Во-первых, корректно задать энергию Гельмгольца
для системы, именно это определяет физику системы. Во-вторых, найти распределение плотности флюида, которое бы минимизировало свободную энергию системы. На
основе минимальности свободной энергии делается вывод о том, что вариация свободной энергии по плотности должна равняться нулю. На этом утверждении строится вся
дальнейшая теория вычисления равновесной плотности [19; 20]. К сожалению, для реальных флюидов свободная энергия Гельмгольца имеет сложную структуру, из-за чего
вычисление ее вариации по плотности является трудоемкой задачей. Поэтому приходится строить различного рода приближения, которые бы позволяли с определенной
точностью описывать физику системы. В последнее время разрабатываются подходы,
позволяющие учитывать гетерогенность поверхности пор [9—11], кулоновскую природу частиц флюида. В работе [9] наглядно продемонстрировано, что усложнение физической модели и учет гетерогенности стенок поры через потенциал взаимодействия
молекул флюида с поверхностью может значительно изменить вид свободной энергии
и ее вариации, наиболее значительное изменение претерпело слагаемое с вариацией от
диполь-дипольного взаимодействия. Кроме изменения физики системы на вид вариации
влияет геометрия задачи. В зависимости от того, какой материал исследуется, вариацию свободной энергии приходится считать в разных системах координат [19], потому
что для одних материалов поры можно рассматривать как цилиндрические трубки, а
5
для других материалов стенки поры представляются как две параллельные стенки, и
задача в таком случае имеет одно выделенное направление.
В данной работе представлен новый подход вычисления равновесной плотности,
который не требует вычисления вариации свободной энергии Гельмгольца (Variation
Free Density Functional Theory). На начальном этапе один раз для всех последующих
задач рассчитывается датасет, состоящий из функций, поведение которых максимально близко к поведению исследуемой функции плотности флюида. Конкретно в данной
работе, в качестве функций для датасета были взяты равновесные плотности азота
при различных относительных давлениях в поре 3.6 нм (эти функции были рассчитаны с помощью классического подхода с методом простой итерации [19]). Также были
добавлены плотности флюидов с параметрами взаимодействия азота, но с отличными радиусами молекул. Добавление таких «искусственных» функций позволяет внести
разнообразие в данные и применять алгоритм для флюида, отличного от азота.
На следующем этапе с помощью метода главных компонент (principle component
analysis — PCA) из датасета выделяются наиболее значимые паттерны, а неизвестная
функция плотности флюида представляется в виде разложения по рассчитанным с помощью PCA функциям, которые содержат эти наиболее значимые паттерны. Таким
образом, задача поиска равновесной плотности сводится к задаче оптимизации коэффициентов разложения по базису. Благодаря тому, что используется не весь базис, а
только его наиболее значимая часть, удается значительно понизить размерность оптимизационной задачи без заметной потери качества решения. На основе энергетического
критерия удалось сократить количество базисных функций со 140 до 10. Полученные
10 векторов содержат 95% информации о исходном базисе.
В безвариационном алгоритме поиск равновесной плотности осуществлялся с помощью стохастических методов оптимизации, в отличие от классических вариационных
подходов, на основе метода простой итерации или ньютоновского метода [19—22]. Среди
широкого спектра стохастических оптимизационных алгоритмов были выбраны и сравнены два наиболее распространенных и популярных генетический алгоритм (genetic
algorithm — GA), метод роя частиц (particle swarm optimization — PSO). Оптимизационные алгоритмы ищут коэффициенты разложения плотности флюида по базису так,
чтобы минимизировать свободную энергию системы.
Также был разработан гибридный подход (Hybrid Density Functional Theory), который основывается на безвариационном методе и методе простой итерации. Работа гибридного алгоритма H-DFT отличается от VF-DFT только тем, что решение, которое
выдает стохастический алгоритм оптимизации подается как начальное приближение
для метода простой итерации. Такая комбинация позволяет сократить время поиска
6
равновесной плотности без потери качества решения.
Важно отметить, что несмотря на то, что базис был построен на основе флюидов
с параметрами азота, удается вычислить равновесную плотность флюида, информации
о котором в базисе не было, VF-DFT хорошо справился с задачей вычисления равновесной плотности аргона. Из преимуществ VF-DFT то, что для его работы не нужно
вычислять вариацию свободной энергии, а время поиска решения значительно меньше,
чем у DFT с методом простой итерации, но качество решения получается незначительно хуже того, что дает классический DFT с методом простой итерации и зависит от
выбора базисных функций. Чтобы сохранить преимущество в скорости и улучшить качество получаемого решения, в работе рассматривается гибридный алгоритм H-DFT на
основе стохастических методов оптимизации и метода простой итерации. Комбинирование подходов позволило ускорить в несколько раз вычисления равновесной плотности
без потери качества.
Разработанные подходы были применены к задаче нахождения равновесной плотности флюида в поре 3.6 нм. Рассматривались флюиды азот и аргон, стенки поры
— углерод, геометрия задачи — планарная, то есть стенки поры представляются как
две параллельные пластины. В качестве потенциала взаимодействия флюида со стенкой был взят потенциал Стилла 10-4-3. Для описания флюид-флюид взаимодействия
использовался подход фундаментальной теории меры (Fundamental Measure Theory —
FMT) для короткодействующего отталкивания, разработанный Розенфельдом [19; 23],
и схема WCA [24] для учета диполь-дипольного взаимодействия. Исследуемый флюид
и поверхность, азот, аргон и карбон, достаточно хорошо изучены как экспериментально, так и с помощью классического DFT. Выбор именно такого флюида и поверхности
обусловлен простотой их моделирования и наличием большого количества данных для
валидации результатов разработанных методов.
В разделе 2 будет приведена справка по теории функционала плотности и методе простой итерации, в разделе 3 будет приведено описание метода анализа главных
компонент, методов оптимизации, которые использовались в данной работе, и подробнее рассмотрены преимущества подходов VF-DFT и H-DFT. Затем, в разделе 4 будут
рассмотрены задачи поиска равновесной плотности азота и аргона и валидация полученных решений с классическим подходом DFT с методом простой итерации. В разделе
5 будет рассмотрен один из классических подходов для восстановления распределения
пор по размерам, а в разделе 5.2 буду приведены результаты работы этого подхода для
решения нескольких модельных задач.
7
2
Теория функционала плотности
Будем рассматривать систему с заданными температурой 𝑇 , объемом 𝑉 , химическим потенциалом 𝜇. Согласно второму началу термодинамики, для системы с внешним
потенциалом 𝑉𝑒𝑥𝑡 (𝑟) термодинамический потенциал большого канонического ансамбля
Ω минимален в равновесном состоянии. Так как Ω потенциал связан со свободной энергией Гельмгольца 𝐹 [𝜌 (𝑟)], то Ω потенциал также является функционалом плотности
[4; 9; 19]:
∫︁
Ω [𝜌 (𝑟)] = 𝐹 [𝜌 (𝑟)] +
𝑑3 𝑟𝜌 (𝑟) (𝑉𝑒𝑥𝑡 (𝑟) − 𝜇) ,
(2.1)
где 𝐹 [𝜌 (𝑟)] — энергия Гельмгольца, 𝜌 (𝑟) — плотность флюида в точке 𝑟, 𝑉𝑒𝑥𝑡 — внешний потенциал (потенциал стенки поры в нашем случае), 𝜇 — химический потенциал
флюида; температура, объем системы и химический потенциал считаются фиксированными.
Из условия экстремальности функционала следует, что его вариация должна
обращаться в ноль в равновесии [4]
𝛿𝐹 [𝜌 (𝑟)]
𝛿Ω [𝜌 (𝑟)]
=
+ 𝑉𝑒𝑥𝑡 (𝑟) − 𝜇 = 0.
𝛿𝜌 (𝑟)
𝛿𝜌 (𝑟)
(2.2)
Решением уравнения (2.2) является равновесная плотность, чтобы решить это
уравнение необходимо знать вид энергии Гельмгольца. Несмотря на то, что формализм
свободной энергии является точным, свободная энергия неизвестна для большинства
систем. Задача вычисления свободной энергии Гельмгольца равносильна задаче вычисления статистической сумм системы, поэтому для решения прикладных задач часто
разбивают свободную энергию на идеальную часть и поправки к ней, которые учитывают различного рода взаимодействия [4]:
𝐹 [𝜌] = 𝐹 𝑖𝑑 [𝜌] + 𝐹 𝑒𝑥 [𝜌]
∫︁
(︀
)︀
𝑖𝑑
𝐹 [𝜌] = 𝑘𝐵 𝑇 𝑑𝑟 𝜌 (𝑟) ln (Λ3 𝜌 (𝑟)) − 1 .
(2.3)
(2.4)
Выражение (2.4) учитывает свободную энергию флюида, как если бы он был
ℎ
—
идеальным газом, где 𝑘𝐵 — константа Больцмана, 𝑇 — температура, Λ = √2𝜋𝑚𝑇
тепловая длина волны де-Бройля, ℎ — постоянная Планка, 𝑚 — масса молекулы газа.
Слагаемое 𝐹 𝑒𝑥 содержит в себе информацию о взаимодействии молекул и именно
оно определяет физику молекул флюида. В общем виде его можно представить через
функции Майера через разложение по плотности следующим образом:
𝐹
𝑒𝑥
𝑘𝐵 𝑇
[𝜌 (𝑟)] = −
2
∫︁
∫︁
𝑑𝑟1
(︀ )︀
𝑑𝑟2 𝜌 (𝑟1 ) 𝜌 (𝑟2 ) 𝑓 (𝑟12 ) + 𝒪 𝜌3 ,
(2.5)
здесь 𝑟12 = |𝑟1 − 𝑟2 |, 𝑓 (𝑟) — функция Майера, которая выражается через потенциал
8
парного межмолекулярного взаимодействия 𝑈 (𝑟1 , 𝑟2 )
(︂
𝑈 (𝑟12 )
𝑓 (𝑟12 ) = exp −
𝑘𝐵 𝑇
)︂
(2.6)
− 1.
Парным потенциалом взаимодействия простых молекул является потенциал Леннарда–
Джонса, но так как этот потенциал обладает особенностью, его проблематично напрямую учесть в выражении для свободной энергии. Поэтому в работе в дальнейшем будет
идти разделение добавочной части свободной энергии на два слагаемых, сумма которых обеспечит корректный учет взаимодействия простых молекул через потенциал
Леннарда–Джонса:
𝐹 𝑒𝑥 = 𝐹 𝐻𝑆 + 𝐹 𝑎𝑡𝑡 .
(2.7)
Первое слагаемое отвечает за короткодействующее отталкивание, а второе за дальнодействующее притяжение. Каждое из слагаемых будет подробнее рассмотрено ниже.
2.1
Короткодействующее отталкивание
Для учета взаимодействия двух атомов на расстояниях порядка радиуса в теории
функционала плотности атомы представляют как абсолютно жесткие сферы, таким
образом каждая частица имеет объем, недоступный для других атомов [25]. Корректное
описание короткодействующего взаимодействия позволяет исследовать такое сложное
явление, как фазовые переходы.
В теории функционала плотности для описания короткодействующего отталкивания популярен подход, основанный на теории меры (Fundamental measure theory —
FMT), который разработал Розенфельд [23]. Подход Розенфельда обладает преимущетвом в универсальности модели для описания различных явлений и позволяет рассматривать систему твердых сфер: в свободном объеме (3D); флюид жестких дисков,
зажатый между двумя плоскостями (2D); флюид жестких стержней в цилиндрической
поре (1D), и в кавернах (0D). Эта замкнутая, математически строгая теория применима
как к однокомпонентному флюиду, так и к полидисперсной смеси флюидов в разных
фазах, состоящей из несферических частиц. В данной работе будет рассматриваться модель жестких сфер для однокомпонентного флюида. В данном случае функция Майера
имеет наглядную геометрическую интерпретацию (𝑅 — здесь и в дальнейшем радиус
частицы, если это не оговорено):
𝑈 (𝑟12 ) =
{︃
∞ 𝑟12 6 2𝑅
0
иначе,
{︃
=⇒
𝑓 (𝑟12 ) =
− 1 𝑟12 6 2𝑅
0
иначе,
(2.8)
Таким образом, функции Майера позволяет учитывать объем, недоступный для
центров различных частиц. Так как выражение (2.5) с учетом вида функции Майера для жестких сфер (2.8) является сложным для прямого аналитического решения,
функцию Майера представляют в виде разложение на более простые слагаемые.
9
Объем двух тел 𝑖 и 𝑗 можно представить следующим образом:
𝑉𝑖+𝑗 = 𝑉𝑖 + 𝑅𝑖 𝑆𝑗 + 𝑅𝑗 𝑆𝑖 + 𝑉𝑗 ,
(2.9)
где 𝑉𝑘 , 𝑆𝑘 , 𝑅𝑘 — объем, площадь поверхности и радиус кривизны соответственно. С
учетом этого, Розенфельд [23] отметил, что функцию Майера можно представить в
виде следующей суммы:
− 𝑓 (𝑟𝑖𝑗 ) = 𝜔3𝑖 ⊗ 𝜔0𝑗 + 𝜔0𝑖 ⊗ 𝜔3𝑗 + 𝜔2𝑖 ⊗ 𝜔1𝑗 + 𝜔1𝑖 ⊗ 𝜔2𝑗 − 𝜔2𝑖 ⊗ 𝜔1𝑗 − 𝜔1𝑖 ⊗ 𝜔2𝑗
(2.10)
с весовыми функциями
𝜔3 (𝑟) = Θ (𝑅 − 𝑟) ,
𝜔2 (𝑟) = 𝛿 (𝑅 − 𝑟) ,
𝑟
𝜔2 (𝑟) = 𝛿 (𝑅 − 𝑟) ,
𝑟
(2.11)
(2.12)
(2.13)
и 𝜔1 (𝑟) = 𝜔2 (𝑟) / (4𝜋𝑅), 𝜔0 (𝑟) = 𝜔2 (𝑟) / (4𝜋𝑅2 ), 𝜔1 (𝑟) = 𝜔2 (𝑟) / (4𝜋𝑅). Здесь Θ (𝑟)
— функция Хэвисайда, а 𝛿 (𝑟) — функция Дирака, символом ⊗ обозначена свертка
функций (𝛼, 𝛽 = 1,2,3)
𝜔𝛼𝑖
⊗
𝜔𝛽𝑗
(𝑟 = 𝑟𝑖 − 𝑟𝑗 ) =
∫︁
𝑑𝑟 ′ 𝜔𝛼𝑖 (𝑟 ′ − 𝑟𝑖 ) 𝜔𝛼𝑗 (𝑟 ′ − 𝑟𝑗 ) .
(2.14)
Связь между разложением функции Майера (2.10) и выражением для объема
двух тел (2.9) становится очевидной, если проинтегрировать весовые функции 𝜔𝛼𝑖 (𝑟),
для 𝛼 = 3 получается объем частицы 𝑉𝑖 , 𝑆𝑖 для 𝛼 = 2, 𝑅𝑖 для 𝛼 = 1 и 1 для 𝛼 = 0, что
является основными геометрическими мерами для сферической частицы в трехмерном
пространстве (Fundamental Measures) [19; 23].
С помощью весовых функций 𝜔𝛼 можно перейти к взвешенным плотностям 𝑛𝛼 ,
через которые в дальнейшем будет выражена свободная энергия Гельмгольца
∫︁
𝑛𝛼 (𝑟) =
𝑑3 𝑟′ 𝜌 (𝑟 ′ ) 𝜔𝛼 (𝑟 − 𝑟 ′ ) .
(2.15)
Если рассмотреть предел в выражении для свободной энергии для короткодействующего отталкивания (2.5) при плотности стремящейся к нулю, то с учетом (2.15)
можно получить точное выражение, в дальнейшем это будет учтено для явной записи
свободной энергии взаимодействия жестких сфер
lim 𝐹
𝜌→0
𝐻𝑆
∫︁
[𝜌 (𝑟)] = 𝑘𝐵 𝑇
𝑑𝑟 {𝑛0 (𝑟) 𝑛3 (𝑟) + 𝑛1 (𝑟) 𝑛2 (𝑟) − 𝑛1 (𝑟) · 𝑛2 (𝑟)}
(2.16)
Получив предельное выражение (2.16) Розенфельд предположил, что и для высоких плотностей вид свободной энергии должен зависеть только от взвешенных плот-
10
ностей [19; 23]
𝐹
𝐻𝑆
∫︁
[𝜌 (𝑟)] = 𝑘𝐵 𝑇
𝑑𝑟 Φ𝑅𝐹 [𝑛𝛼 (𝑟)]
(2.17)
Из анализа размерности подинтегральной функции можно предположить, что ее
вид должен быть следующим:
Φ𝑅𝐹 = 𝑔1 (𝑛3 ) 𝑛0 + 𝑔2 (𝑛3 ) 𝑛1 𝑛2 + 𝑔3 (𝑛3 ) 𝑛1 · 𝑛2 + 𝑔4 (𝑛3 ) 𝑛32 + 𝑔5 (𝑛3 ) 𝑛2 𝑛2 · 𝑛2
(2.18)
Каждое слагаемое в выражении (2.18) имеет размерность (длина−3 ), функции
𝑔1 , . . . ,𝑔5 являются неизвестными. Однако с учетом выражения (2.16) можно сделать
вывод о поведении функций 𝑔1 , . . . ,𝑔5 при низки плотностях 𝑔1 = 𝑛3 + 𝑛23 /2 + 𝒪 (𝑛33 ),
𝑔2 = 1 + 𝑛3 + 𝒪 (𝑛23 ), 𝑔3 = −1 − 𝑛3 + 𝒪 (𝑛23 ), 𝑔4 = 1/24𝜋 + 𝒪 (𝑛3 ), 𝑔5 = −3/24𝜋 + 𝒪 (𝑛3 ).
Отсюда можно сделать вывод о взаимосвязи функций:
𝑔5 = −3𝑔4 ,
𝑔2 = −𝑔3 .
С учетом этой взаимосвязи выражение для плотности свободной энергии будет
иметь вид
(︀
)︀
Φ𝑅𝐹 = 𝑔1 (𝑛3 ) 𝑛0 + 𝑔2 ((𝑛3 ) 𝑛1 𝑛2 − (𝑛3 ) 𝑛1 · 𝑛2 ) + 𝑔4 (𝑛3 ) 𝑛32 − 3 (𝑛3 ) 𝑛2 𝑛2 · 𝑛2 .
(2.19)
Функции 𝑔1 , 𝑔2 , 𝑔4 могут быть определены из требований выполнения термодинамических соотношений для функционала Φ𝑅𝐹 . В работе [23] Розенфельд использовал
соотношение из теории масштабированных частиц (Scaled Particle Theory — SPT):
𝑝
𝜇𝐻𝑆
=
𝑅→∞ 𝑘𝐵 𝑇 𝑉
𝑘𝐵 𝑇
lim
(2.20)
где 𝑉 — объем сферической частицы с радиусом 𝑅, 𝑝 — давление, 𝜇𝐻𝑆 — химический
потенциал от взаимодействия жестких сфер. Химический потенциал от взаимодействия
жестких сфер может быть найден варьированием функционала Φ𝑅𝐹 по плотности (из
определения химического потенциала по свободной энергии):
𝜇𝐻𝑆
𝛿𝐹 𝑒𝑥 [𝜌 (𝑟)] ∑︁ 𝜕Φ𝑅𝐹 𝜕𝑛𝛼
=
=
.
𝑘𝐵 𝑇
𝛿𝜌 (𝑟)
𝜕𝑛𝛼 𝜕𝜌
𝛼
(2.21)
∫︀
Учитывая геометрический смысл весовых функций 𝜕𝑛3 /𝜕𝜌 = 𝑑𝑟 𝜔3 (𝑟) ≡ 𝑉 , 𝜕𝑛2 /𝜕𝜌 ≡
𝑆, 𝜕𝑛1 /𝜕𝜌 ≡ 𝑅, 𝜕𝑛0 /𝜕𝜌 ≡ 1 (производные с векторными взвешенными плотностями
равны нулю), предельное выражение (2.20) получается следующим:
𝜇𝐻𝑆
𝜕Φ𝑅𝐹
𝑝
=
.
=
𝑅→∞ 𝑘𝐵 𝑇 𝑉
𝜕𝜌
𝑘𝐵 𝑇
lim
11
(2.22)
С другой стороны можно воспользоваться тождеством для большого канонического ансамбля для флюида в свободном объеме (bulk) Ω𝑏𝑢𝑙𝑘 = −𝑝𝑉 , а плотность
свободной энергии большого канонического ансамбля Ω𝑏𝑢𝑙𝑘 /𝑉 = Φ + 𝑓𝑖𝑑 + 𝜇𝜌𝑏𝑢𝑙𝑘 . Таким
образом, получается система уравнений:
∑︁ 𝜕Φ𝑅𝐹
𝜕Φ𝑅𝐹
= −Φ𝑅𝐹 +
𝑛𝛼 + 𝑛0 .
𝜕𝑛3
𝜕𝑛
𝛼
𝛼
(2.23)
Собирая все члены, содержащие 𝑛0 получается дифференциальное уравнение для 𝑔1
(’ обозначает производную по 𝑛3 )
𝑔1′ (𝑛3 ) (1 − 𝑛3 ) = 1
𝑔1 (𝑛3 ) = 𝑐𝑜𝑛𝑠𝑡1 − ln (1 − 𝑛3 )
=⇒
(2.24)
Аналогично и для других неизвестных функций 𝑔2 , 𝑔4
𝑔2′ (𝑛3 ) (1 − 𝑛3 ) = 𝑔2 (𝑛3 )
=⇒
𝑔2 (𝑛3 ) =
𝑐𝑜𝑛𝑠𝑡2
1 − 𝑛3
(2.25)
𝑔4′ (𝑛3 ) (1 − 𝑛3 ) = 2𝑔4 (𝑛3 )
=⇒
𝑔4 (𝑛3 ) =
𝑐𝑜𝑛𝑠𝑡3
(1 − 𝑛3 )
(2.26)
Константы интегрирования можно подобрать, если подставить полученные выражения в выражение для предела свободной энергии от взаимодействия жестких сфер при
плотности, стремящейся к нулю (2.16). Таким образом, функционал Розенфельда [19;
23] окончательно можо записать в следующем виде:
Φ𝑅𝐹 = −𝑛0 ln (1 − 𝑛3 ) +
𝑛1 𝑛2 − 𝑛1 · 𝑛2 𝑛32 − 3𝑛2 𝑛2 · 𝑛2
+
.
1 − 𝑛3
24𝜋 (1 − 𝑛3 )2
(2.27)
Важно отметить, что функционал Розенфельда согласуется с уравнением состояние Перкус–Йевика и прямой корреляционной функцией Перкус–Йевика. На рис. 1 приведен профиль плотности флюида жестких сфер в поре с абсолютно жесткими стенками. Результаты расчета свалидированы с результатами из [19] при значении плотности
упаковки флюида 𝜂 = 4𝜋𝑅3 𝜌𝑏𝑢𝑙𝑘 /3 = 0.4783. Внешний потенциал в такой поре задается
следующим образом:
{︃
∞ 𝑟 > 𝐻, 𝑟 6 0
𝑉𝑒𝑥𝑡 (𝑟) =
(2.28)
0 иначе,
Таким образом, в данной работе учитывается короткодействующее отталкивание
между молекулами флюида в свободной энергии Гельмгольца через потенциал Розенфельда.
12
12
FMT расчет
Данные из Roth 2010
10
𝜌𝜎 3
8
6
4
2
0
0.6 0.8
1
1.2 1.4 1.6 1.8
𝑧/𝜎
2
2.2 2.4 2.6 2.8
3
Рисунок 1: Распределения профиля плотности жестких сфер в поре из абсолютно жестких стенок при 𝜂 = 0.4783 [19]
2.2
Диполь-дипольное взаимодействие
Большинство версий теорий функционала плотности учитывают диполь-дипольное
взаимодействие между атомами и молекулами через приближение среднего поля. Несмотря на то, что использование приближения среднего поля накладывает ряд ограничений
на исследуемую систему и обладает меньшей общностью, существуют модификации
функционала Гельмгольца, которые помогают обойти это препятствие. В данной работе
диполь-дипольное взаимодействие рассматривается как возмущение над потенциалом
короткодействующего отталкивания по схеме, которую разработали Weeks, Chandler,
Anderson (WCA) [24]
𝑈 (𝑟) = 𝑈0 (𝑟) + 𝑈𝑎𝑡𝑡 (𝑟)
(2.29)
где 𝑈0 (𝑟) — потенциал взаимодействия жестких сфер, 𝑈𝑎𝑡𝑡 (𝑟) — возмущение. Из работы
[24] потенциал возмущения представляется следующим образом:
⎧
⎪
𝑟<𝜆
⎨−𝜀𝑓 𝑓
𝑈𝑎𝑡𝑡 (𝑟) = 𝑈𝐿𝐽 𝜆 < 𝑟 < 𝑟𝑐𝑢𝑡
⎪
⎩
0
𝑟 > 𝑟𝑐𝑢𝑡
𝑈𝐿𝐽 = 4𝜀𝑓 𝑓
(︂(︁
𝜎 )︁12
𝑓𝑓
𝑟
13
−
(︁ 𝜎 )︁6 )︂
𝑓𝑓
𝑟
(2.30)
(2.31)
𝜆 = 21/6 𝜎𝑓 𝑓 — положение минимума потенциала Леннарда–Джонса, 𝜎𝑓 𝑓 — диаметр
частицы флюида, 𝜀𝑓 𝑓 — значение потенциала Леннарда–Джонса в минимуме, 𝑟𝑐𝑢𝑡 —
расстояние действия потенциала. На рис. 2 схематично проиллюстрирован потенциал
взаимодействия двух молекул [24].
Рисунок 2: Схематичная иллюстрация модифицированного потенциала Леннарда–
Джонса по схеме WCA из [24] (на рисунке 𝑈 (𝑟) = 𝑈𝑎𝑡𝑡 (𝑟))
В итоге, в приближении среднего поля свободная энергия диполь-дипольного
взаимодействия записывается следующим образом:
𝐹
𝑎𝑡𝑡
𝑘𝐵 𝑇
=
2
∫︁ ∫︁
𝑑3 𝑟1 𝑑3 𝑟2 𝜌 (𝑟1 ) 𝜌 (𝑟2 ) 𝑈𝑎𝑡𝑡 (𝑟) .
(2.32)
На рис. 3 приведен график распределения плотности азота в поре, который был посчитан с помощью приведенной в данной работе модели DFT и свалидирован с результатами из статьи [6].
14
0.18
DFT расчет
Данные из Ravikovitch 2001
𝜌, Å−3
0.14
0.1
0.06
0.02
2
3
4
5
6
7
8
9
10
11
12
13
𝑧, Å
Рисунок 3: Распределение профиля плотности азота в поре, шириной 𝐻𝑐𝑐 = 10𝜎𝑓 𝑓 при
относительном давлении в свободном объеме 𝑃/𝑃0 = 0.7 [6]
2.3
Метод простой итерации (МПИ)
Чтобы найти равновесную плотность флюида, которая бы минимизировала свобоную энергию большого канонического ансамбля (2.1), методом просто итерации, необходимо знать вариацию омега-потенциала [19]. С учетом выражения (2.3) и (2.4) равенство нулю вариации свободной энергии большого канонического ансамбля запишется в
виде:
𝛿𝐹 𝐻𝑆 [𝜌 (𝑟)] 𝛿𝐹 𝑎𝑡𝑡 [𝜌 (𝑟)]
+
+ 𝑉𝑒𝑥𝑡 (𝑟) − 𝜇 = 0.
𝑘𝐵 𝑇 ln Λ 𝜌 (𝑟) +
𝛿𝜌 (𝑟)
𝛿𝜌 (𝑟)
3
(2.33)
Аналогично энергии Гельмгольца, химический потенциал флюида можно также условно представить как сумму трех слагаемых
𝜇 = 𝜇𝑖𝑑 + 𝜇𝐻𝑆 + 𝜇𝑎𝑡𝑡 ,
(︀
)︀
𝜇𝑖𝑑 (𝜌) = 𝑘𝐵 𝑇 ln Λ3 𝜌𝑏𝑢𝑙𝑘 = 𝑐𝑜𝑛𝑠𝑡,
(︂∑︁
)︂
𝜕Φ 𝜕𝑛𝛼
𝐻𝑆
𝜇 (𝜌) = 𝑘𝐵 𝑇
= 𝑐𝑜𝑛𝑠𝑡,
𝜕𝑛𝛼 𝜕𝜌
∫︁
𝑎𝑡𝑡
𝑏𝑢𝑙𝑘
𝜇 (𝜌) = 𝜌
𝑑3 𝑟 𝑈𝑎𝑡𝑡 (𝑟) = 𝑐𝑜𝑛𝑠𝑡.
15
(2.34)
(2.35)
(2.36)
(2.37)
Под 𝜌𝑏𝑢𝑙𝑘 = Λ−3 exp (𝜇𝑖𝑑 /𝑘𝐵 𝑇 ) подразумевается плотность флюида в свободном объеме.
В уравнении (2.33) плотность входит явно в первом слагаемом и ее можно выразить
𝜌ˆ (𝑟) = 𝜌
𝑏𝑢𝑙𝑘
{︂
exp
1
𝑘𝐵 𝑇
(︂
)︂}︂
𝛿𝐹 𝐻𝑆 [𝜌 (𝑟)] 𝛿𝐹 𝑎𝑡𝑡 [𝜌 (𝑟)]
𝐻𝑆
𝑎𝑡𝑡
−
−
− 𝑉𝑒𝑥𝑡 (𝑟) + 𝜇 + 𝜇
. (2.38)
𝛿𝜌 (𝑟)
𝛿𝜌 (𝑟)
Таким образом, получается выражение для равновесной плотности флюида, однако плотность неявно входит в правую часть (2.38), через вариацию энергии Гельмгольца. Для того, чтобы решить это нелинейное уравнение обычно используют метод
простой итерации с релаксацией. В качестве начального приближения выбирается плотность флюида в свободном объеме 𝜌𝑗=0 = 𝜌𝑏𝑢𝑙𝑘 (от выбора начального приближения
зависит устойчивость и скорость сходимости метода). Решение на следующем шаге итерации определяется как сумма решений на предыдущей итерации и на текущей
𝜌𝑗+1 = (1 − 𝛾) 𝜌𝑗 + 𝛾 𝜌ˆ𝑗
(2.39)
Параметр 𝛾 ∈ [0,1] определяет степень доверия новому решению на последующей итерации, чем он больше, тем быстрее сходится метод, но при слишком больших
значениях решение будет неустойчиво и сходимости не будет. Метод простой итерации
— весьма надежен и с его помощью, рано или поздно решение будет найдено, вот только
чаще всего работает он достаточно долго. Кроме того, для работы этого метода необходимо вычислять вариацию свободной энергии Гельмгольца и это накладывает некоторые ограничения на вид энергии Гельмгольца. Усложнение модели, например учет
гетерогенности поверхности [9—11], или учет кулоновского взаимодействия, приводит
к изменению вида энергии Гельмгольца, что влечет изменение соответствующих вариаций. Становится необходимым изменение практически каждого слагаемого в (2.38) и
добавление новых. В работе [9] наглядно видно то, как учет неоднородности поверхности может сильно повлиять на выражения для вариаций, и какие сложности могут при
этом возникнуть в вычислении равновесной плотности, несмотря на то, что рассматриваемый флюид является простым для моделирования.
2.4
Геометрия поры
В работе рассматривается планарная геометрия поры (стенки поры — параллельные пластины). В такой постановке, задача становится симметричной по двум направлениям (𝑥, 𝑦 условно) и плотность будет зависеть только от одной координаты
𝜌 (𝑟) = 𝜌 (𝑧). С учетом такой геометрии можно проинтегрировать по переменным 𝑥, 𝑦
выражения для взвешенных плотностей:
∫︁
𝑛𝛼 (𝑧) =
𝑑3 𝑧 ′ 𝜌 (𝑧 ′ ) 𝜔𝛼 (𝑧 − 𝑧 ′ ) ,
16
(2.40)
где к весовым функциям после интегрирования: 𝜔3 = 𝜋 (𝑅2 − 𝑧 2 ) Θ (𝑅 − |𝑧|), 𝜔2 =
2𝜋𝑅Θ (𝑅 − |𝑧|), 𝜔2 = 2𝜋𝑧Θ (𝑅 − |𝑧|), 𝜔1 = 𝜔2 /4𝜋𝑅, 𝜔1 = 𝜔2 /4𝜋𝑅, 𝜔0 = 𝜔2 /4𝜋𝑅2 . Также
можно выписать вариацию для свободной энергии взаимодействия жестких сфер
∑︁ ∫︁
𝛿𝐹 𝐻𝑆 [𝜌 (𝑧)]
𝜕Φ𝑅𝐹 {𝑛𝛼 (𝑧 ′ )} 𝛿𝑛𝛼 (𝑧 ′ )
= 𝑘𝐵 𝑇
,
𝑑𝑧 ′
𝛿𝜌 (𝑧)
𝜕𝑛𝛼 (𝑧 ′ )
𝛿𝜌 (𝑧 ′ )
𝛼
𝛿𝑛𝛼 (𝑧 ′ )
𝛿
=
′
𝛿𝜌 (𝑧 )
𝛿𝜌 (𝑧)
∫︁
𝑑𝑧 ′′ 𝜌 (𝑧 ′′ ) 𝜔𝛼 (𝑧 ′ − 𝑧 ′′ ) = 𝜔𝛼 (𝑧 ′ − 𝑧) .
(2.41)
(2.42)
Также можно упростить выражения для вариации энергии от диполь-дипольного
взаимодействия
𝐹
𝑎𝑡𝑡
𝑘𝐵 𝑇
=
2
∫︁
∫︁
𝑑𝑧1 𝑑𝑧2 𝜌 (𝑧1 ) 𝜌 (𝑧2 )
𝑑𝑥1 𝑑𝑥2 𝑑𝑦1 𝑑𝑦2 𝑈𝑎𝑡𝑡 (𝑟) =
∫︁ ∫︁
𝑘𝐵 𝑇
=
𝑑𝑧1 𝑑𝑧2 𝜌 (𝑧1 ) 𝜌 (𝑧2 ) 𝐺 (𝑧1 ,𝑧2 ) , (2.43)
2
где с учетом (2.30) при |∆𝑧| 6 𝜆
⎫
⎧
√2
√
2
𝑟
−Δ𝑧
2
2
⎪
𝑐𝑢𝑡
⎪
⎪
[︂
]︂⎪
∫︁
⎬
⎨ 𝜆∫︁−Δ𝑧
𝜎𝑓6𝑓
𝜎𝑓12𝑓
−
=
−𝜀𝑓 𝑓 𝑟𝑑𝑟 +
𝑑𝑟 4𝜀𝑓 𝑓 𝑟
𝐺 (𝑧1 ,𝑧2 ) = 2𝜋
2 + 𝑟 2 )6
2 + 𝑟 2 )3 ⎪
⎪
(∆𝑧
(∆𝑧
⎪
⎪
√
⎭
⎩ 0
𝜆2 −Δ𝑧 2
(︂
)︂
)︂
(︂
(︀ 2
)︀ 4
1
1
1
1
2
12
6
= −𝜋𝜀𝑓 𝑓 𝜆 − ∆𝑧 + 𝜋𝜀𝑓 𝑓 𝜎𝑓 𝑓
− 10 − 2𝜋𝜀𝑓 𝑓 𝜎𝑓 𝑓
− 4
, (2.44)
5
𝜆10
𝑟𝑐𝑢𝑡
𝜆4
𝑟𝑐𝑢𝑡
при |∆𝑧| > 𝜆 и |∆𝑧| 6 𝑟𝑐𝑢𝑡
√
2 −Δ𝑧 2
𝑟𝑐𝑢𝑡
∫︁
[︂
0
𝜎𝑓6𝑓
]︂
−
=
(∆𝑧 2 + 𝑟2 )6
(∆𝑧 2 + 𝑟2 )3
(︂
)︂
(︂
)︂
1
1
4
1
1
12
6
= 𝜋𝜀𝑓 𝑓 𝜎𝑓 𝑓
− 10 − 2𝜋𝜀𝑓 𝑓 𝜎𝑓 𝑓
− 4
, (2.45)
5
𝑟𝑐𝑢𝑡
𝑟𝑐𝑢𝑡
∆𝑧 10
∆𝑧 4
𝑑𝑟 4𝜀𝑓 𝑓 𝑟
𝐺 (𝑧1 ,𝑧2 ) = 2𝜋
𝜎𝑓12𝑓
а при |∆𝑧| > 𝑟𝑐𝑢𝑡 , 𝐺 (𝑧1 ,𝑧2 ) = 0. Таким образом, вариация энергии Гельмгольца от
диполь-дипольного взаимодействия
𝛿𝐹 𝑎𝑡𝑡 [𝜌 (𝑧)]
𝑘𝐵 𝑇
=
𝛿𝜌 (𝑧)
2
∫︁
𝑑𝑧 ′ 𝜌 (𝑧 ′ ) 𝐺 (𝑧 − 𝑧 ′ )
(2.46)
Стенки поры в работе рассматриваются углеродные и потенциал взаимодействия
флюида со стенкой описывается потенциалом Стилла 10-4-3 [26]
𝑉𝑠𝑓 (𝑧) =
3
2𝜋𝜀𝑠𝑓 𝜌𝑉 𝜎𝑠𝑓
∆
(︂ (︁
)︂
4
𝜎𝑠𝑓
2 𝜎𝑠𝑓 )︁10 (︁ 𝜎𝑠𝑓 )︁4
.
−
−
5 𝑧
𝑧
3∆ (0.61∆ + 𝑧)3
17
(2.47)
Полный потенциал в поре от двух стенок:
𝑉𝑒𝑥𝑡 (𝑧) = 𝑉𝑠𝑓 (𝑧) + 𝑉𝑠𝑓 (𝐻𝑐𝑐 − 𝑧) ,
(2.48)
−3
𝜌𝑉 — плотность материала стенки 0.114 Å в нашем случае для углерода, ∆ — расстояние между слоями атомов углерода 3.35 Å, 𝐻𝑐𝑐 — диаметр поры рассчитанный между центрами атомов углерода противоположных стенок, ширина поры доступная для
флюида 𝐻 = 𝐻𝑐𝑐 − 2𝑧0 + 𝜎𝑓 𝑓 , где 𝑧0 — глубина проникновения молекулы флюида в
поверхность (значение 𝑧, при котором потенциал стенки обращается в ноль), в нашей
работе мы использовали приближение 𝑧0 = 0.9𝜎𝑠𝑓 .
3
Безвариационный подход в теории функционала плотности (VF-DFT)
В отличие от классического DFT, который использует метод простой итерации
для поиска равновесного решения, безвариационный подход в теории функционала
плотности (Variation Free Density Functional Theory VF-DFT) использует стохастические методы оптимизации. На рис. 4 схематично представлен принцип работы разработанных алгоритмов в сравнении с классическим методом. Применение безградиентных стохастических методов оптимизации, таких как генетический алгоритм (Genetic
Algorithm — GA) и метод роя частиц (Particle Swarm Optimization — PSO), позволяет найти равновесную плотность 𝜌* (𝑧) без вычисления вариаций и сократить время
поиска решения. Отсутствие необходимости считать вариацию свободной энергии является главным преимуществом VF-DFT. Для разных систем вид свободной энергии
Гельмгольца различен, следовательно, отличается и ее вариация. Для некоторых систем
вычисление вариации не представляется возможным из-за сложной структуры свободной энергии системы. Например, в работе [9] от стандартной модели адсорбции газа
на идеально гладкой поверхности перешли к рассмотрению адсорбции на гетерогенной
поверхности, это повлекло сильные изменения в вариации энергии Гельмгольца. Для
того, чтобы учитывать более сложную структуру молекул флюида, требуется усложнять и менять вид свободной энергии Гельмгольца. Для описания сложных молекул
активно развивается в настоящее время статистическая теория ассоциативных флюидов (Statistical Associating Fluid Theory — SAFT) [27; 28], из-за учета цепочечных типов
связей в молекуле в выражении для свободной энергии появляется дополнительное слагаемое. Помимо того, что необходимо выбирать физическую модель и соответствующий
ей вид свободной энергии, исследуемые задачи часто наделяют определенными геометрическими свойствами, из-за которых диктуется тип системы координат и упрощается
вычисление вариаций [19], переход от одной системы координат к другой может также
вызывать трудности с точки зрения пересчета вариаций. Для некоторых физических
моделей, из-за сложной структуры энергии Гельмгольца не всегда возможно посчитать
18
вариацию по плотности не прибегая к каким-либо ухищрениям и упрощениям модели,
что является проблемой для исследования этих систем. Использование VF-DFT, дает
возможность обойти эти преграды и получить достаточно точное решение, зная лишь
вид свободной энергии Гельмгольца.
Модель флюида 𝐹 [𝜌 (𝑧)]
VF-DFT
Classical DFT
0.3
0.2
0.1
Вычисление 𝛿𝐹 [𝜌]
вариаций
𝛿𝜌 (𝑧)
Набор функций
PCA
Стохастическая
оптимизация
*
𝑎𝑖 )
𝜌𝑠𝑜 = arg min Ω (̃︀
̃︀
𝑎𝑖
Метод простой
итерации
H-DFT
Равновесная
плотность 𝜌* (𝑧)
Рисунок 4: Схематическое описания принципа работы классического подхода в DFT с
методом простой итерации (черная линия), подхода Variation Free DFT (красная линия),
гибридного подхода Hybrid DFT (красные и синии линии)
3.1
Построение базисных функций
Задача оптимизации свободной энергии системы для поиска равновесной плотности флюида обладает слишком большим пространством поиска. Чтобы уменьшить
размерность оптимизационной задачи можно использовать информацию о характерном
поведении искомой функции, вычленить основные закономерности и искать плотность
флюида в виде линейной комбинации функций, которые содержат в себе эту информацию. В таком представлении искомой функции, задача поиска равновесной плотности
19
сводится к поиску оптимальных коэффициентов разложения 𝑎𝑖 , которые бы минимизировали Ω потенциал Ω [𝜌 (𝑧)] = Ω [𝑎𝑖 ]. На качество и скорость поиска решения будет
влиять количество и вид базисных функций.
Для того, чтобы базисные функции наиболее точно описывали поведение плотности произвольного флюида в поре, было решено составить датасет 𝑋 из функций
плотности азота при температуре 77.4 К и разных относительных давлениях (от 10−6
до 0.99). Эти функции плотности были посчитаны классическим DFT с методом простой итерации один раз для всех рассматриваемых кейсов в разделе 4. Также, для того,
чтобы добавить в базис как можно больше информации о поведении искомой функции,
были взяты равновесные плотности флюидов с радиусами молекул от 0.85 Å, до 2 Å при
температуре 77.4 K, остальные параметры взаимодействия соответствуют параметрам
азота и были зафиксированы. Такой выбор базисных функций никак не ограничивает
область применимости алгоритма, в разделе 4 будет продемонстрировано, что VF-DFT
применим для флюида, отличного от азота, при других термодинамических условиях.
В итоге, получился датасет 𝑋 из 140 векторов рис. 5. Чтобы выделить характерные
0.4
Матрица 𝑋
0.3
0.3
𝜌 (𝑧), Å−3
50
0.2
0.2
100
0.1
0.1
0
0
5
10
15
𝑧, Å
Рисунок 5: Три различных элемента 𝜌 (𝑧) из матрицы 𝑋. На вставленном изображении
представлено изображение датасета 𝑋, где цветом отмечено значение соответствующей
функции при заданной координате
паттерны для функции плотности и наиболее точно описать поведение плотности произвольного флюида в поре, можно применить метод главных компонент [29] (Principle
Component Analysis — PCA) к датасету 𝑋 с разнообразными решениями. Пусть векторстолбец 𝑥 — одно из распределений плотности (в нашей работе их 140). Отнормируем
матрицу 𝑋(𝑁 𝑥 × 𝑁 𝑠), 𝑁 𝑥 — размерность по 𝜌, 𝑁 𝑠 — число функций в датасете 𝑋,
следующим образом:
20
0.2
0.1
Энергия компонент
0
−0.1
−0.2
−0.3
−0.4
1
0.8
0.6
0.4
0.2
0
101
102
Число компонент
8
10
12
14
−0.5
0
2
4
6
𝑧, Å
Рисунок 6: Пять первых главных компонент матрицы 𝑋. Относительная ошибка реконструкции полного датасета составляет 4.84%, график распределения ошибки показан
темно-зеленой линией
𝑋 (𝑖) = 𝑥 − 𝑥𝑚𝑒𝑎𝑛 ,
(3.1)
где 𝑥𝑚𝑒𝑎𝑛 (𝑁 𝑥 × 1) — среднее распределение плотности по датасету 𝑋.
Согласно [30], разложим центрированную матрицу 𝑋 с помощью сингулярного
разложения (Singular Value Decomposition — SVD):
𝑋 = 𝑈 · Σ · 𝑉 T,
(3.2)
где 𝑈 (𝑁 𝑥 × 𝑁 𝑠) — матрица левых сингулярных векторов (ортонормированный базис),
Σ (𝑁 𝑠 × 𝑁 𝑠) — квадратная диагональная матрица, в которой содержатся сингулярные
числа, отсортированные по убыванию, 𝑉 (𝑁 𝑠 × 𝑁 𝑠) — матрица правых сингулярных
векторов (тоже ортонормированный базис). Как будет показано далее, левые сингулярные вектора как раз и несут информацию об искомых паттернах, а соответствующие им сингулярные числа показывают, насколько часто встречается данный паттерн в нашем датасете. Для облегчения вычислений в РСА и оптимизации применим
усеченный (truncated) SVD, выбрав число компонент 𝐾 из энергетического критерия
на Σ. Квадрат сингулярного числа равен дисперсии по соответствующему направлению в пространстве признаков, проекция на такое направление будет минимальна. Под
энергией элемента базиса понимается квадрат соответствующего сингулярного числа,
нормированный на сумму квадратов всех сингулярных чисел, соответственно чем больше энергии у выбранных направлений, тем меньше будет проекция у признаков и тем
21
лучше выбранные компоненты будут описывать датасет [31]
𝐼 = 𝜎2/
∑︁
𝜎2.
(3.3)
̃︁ = 𝑈
̃︀ · Σ
̃︀ · 𝑉̃︀ T . Для датасета 𝑋
̃︁ можно вычислить сэмплоТаким образом, 𝑋 ≈ 𝑋
̃︁ · 𝑋
̃︁T /(𝑁 𝑠 − 1). Обратим внимание,
вую (оценочную) ковариационную матрицу 𝑄 = 𝑋
̃︁ · 𝑋
̃︁T = 𝑈
̃︀ · Σ
̃︀ · 𝑉̃︀ T · 𝑉̃︀ · Σ
̃︀ T · 𝑈
̃︀ T = 𝑈
̃︀ · Σ
̃︀ 2 · 𝑈
̃︀ T , то есть квадрат сингулярных
что 𝑋
̃︁ равен собственным числам 𝑆 матрицы 𝑄 с точностью до поправчисел матрицы 𝑋
̃︀ равны собственным векторам
ки Бесселя 1/(𝑁 𝑠 − 1), а левые сингулярные вектора 𝑈
̃︀ и 𝑆, можно воспользоваться подходом из [32],
ковариационной матрицы 𝑄. Зная 𝑈
в котором новая реализация генерируется с использованием ковариационной матрицы
следующим образом:
̃︀ · 𝑆 · 𝑎 + 𝑥𝑚𝑒𝑎𝑛 .
𝑥𝑛𝑒𝑤 = 𝑈
(3.4)
Здесь 𝑎 — вектор 𝐾 × 1 из распределения Гаусса 𝒩 (0,1), его можно трактовать как
вектор оптимизируемых параметров (компоненты этого вектора по сути являются коэффициентами разложения искомой плотности по базисным функциям, которые содержат информацию о характерном поведении плотности флюида вблизи стенок), изменяя
который мы будем подбирать 𝑥𝑛𝑒𝑤 , минимизирующий Омега-потенциал.
В результате, после анализа главных компонент на основе матрицы 𝑋 был вычислен новый базис размером 10 векторов (вид первых пяти главных компонент базиса
на рис. 6), которые охватывают 95% информации по энергетическому критерию. Это
значит, что вместо того, чтобы искать 140 неизвестных чисел 𝑎, придется искать всего 10, причем качество такого подхода практически не будет отличаться от решения с
поиском всех 140 коэффициентов.
3.2
Стохастическая оптимизация
Как уже было показано, рассматриваемая задача оптимизации характеризуется тем, что вычислить вариацию свободной энергии достаточно сложно аналитически,
а изменение вида потенциала Гельмгольца ведет к необходимости пересчета каждого
слагаемого в (2.38). В работе была применена стохастическая оптимизация, в которой исследуемая функция используется как черный ящик, без вычисления вариаций.
Искомая функция представляется в виде (3.4), а алгоритмы оптимизации подбирают
элементы вектора 𝑎 так, чтобы минимизировать функционал свободной энергии Ω [𝑎].
Размерность оптимизационной задачи равняется размерности вектора 𝑎, в данной работе эта размернасть равна 10×1. Были рассмотрены два итерационных эвристических
подхода — ГА и МРЧ, широко применяющихся в науке и технике и зарекомендовавшие
себя во многих отраслях [33—36]. В этом разделе будут представлены краткие описания
принципов работы алгоритмов, а также рекомендации по настройке их параметров.
22
3.2.1
Генетический алгоритм (ГА)
Создание генетический алгоритм было вдохновлено идеей Дарвиновской теории
эволюции. В основе алгоритма лежит идея, что наиболее приспособленные особи, которые эволюционируют, наследуют признаки у своих родителей, мутируют, именно они и
выживают («выживает сильнейший»). ГА был успешно применен для решения самых
разнообразных задач из различных дисциплин: отбор признаков в машинном обучении,
обработки изображений [37—39]. Описать все модификации и методы разработанные на
основе классического генетического алгоритма не представляется возможным, наиболее
полное изложение можно найти в статьях [34; 35].
Инициализация: будем называть вектор 𝑎 хромосомой, набор хромосом — поколение. Рекомендуем брать 𝑁 𝑐 хромосом равным число оптимизируемых параметров.
Сначала они все случайные в заданных границах, для каждой вычисляется целевая
функция (фитнес функция) — в задаче минимизации чем она меньше, тем более хромосома «приспособлена». В данном случае роль целевой функции играет омега потенциал Ω [𝑎]. В генетическом алгоритме оптимизации каждая итерация представляет
собой процесс формирования дочернего поколения из родительского. Этот процесс происходит с помощью генетических операторов. В работе использовались самые базовые
операторы, распишем их подробнее.
Селекция: одна хромосома с лучшим значением целевой функции переносится в
дочернее поколение из родительского без изменений. Это позволяет как минимум не потерять лучшее найденное к данному моменту решение, что обеспечивает немонотонную
сходимость алгоритма.
Скрещивание: каждая новая хромосома (𝑁 𝑐−1 штук) получается путем скрещивания пары родительских. Вероятность выбора хромосомы из родительского поколения
для скрещивания может быть как заданной, так и зависеть от фитнес-функции (roulette
wheel [35]). Использование того или иного подхода изменяет свойства алгоритма, в данной работе используется комбинированный подход. Наиболее распространенные типы
скрещивания — через точку, когда у родительских хромосом выбирается точка 𝑐 и все
гены с индексами меньше чем 𝑐 у родительских хромосом меняются между друг другом,
или через индексы, когда выбирается несколько точек в родительских хромосомах и обмен происходит блоками генов. В данной работе был использован подход скрещивания
через индексы.
Мутация: новая хромосома, полученная на предыдущем этапе из пары родительских, с заданной вероятностью подвергается мутации, т.е. случайному изменению
одной или нескольких ее компонент. Высокий уровень мутации замедляет сходимость
алгоритма, но позволяет более успешно находить глобальный оптимум. Уровень мутации в работе был определен на уровне 0.5 — это значит, что мутации подвержено
половина хромосом популяции.
После формирования дочернего поколения для всех хромосом вычисляется фитнесфункция (если она не была известна). Алгоритм останавливается по достижении отве23
денного числа вызовов фитнес-функции.
3.2.2
Метод роя частиц (МРЧ)
Генетические алгоритмы на сегодняшний день являются одними из самых популярных инструментов решения задач оптимизации. Однако их использование не лишено недостатков. Решение сходится к оптимуму лишь после большого числа вычислений, при этом его глобальность не гарантируется. Нет четкого критерия останова – при
недостаточной мутации решение может очень долго находиться в локальном минимуме,
монотонная сходимость не гарантируется. Все это приводит к появлению новых, более
развитых ГА, использование ГА в связке с другими алгоритмами, а также развитие
новых стохастических подходов.
Одним из таких подходов является метод роя частиц (в англоязычной литературе
PSO — Particle Swarm Optimization). Впервые он был разработан в 1995 году Кеннеди,
Эберхартом и Ши [40] для исследования поведения роя насекомых или косяка рыб,
передвигающихся по пространству в поисках пищи. В дальнейшем он был упрощен
и применен для решения различных задач оптимизации [41]. Рассмотрим этот метод
подробнее, пользуясь терминологией из [36].
Инициализация: вектор 𝑎 называем частицей, набор частиц — рой (по аналогии с
хромосомой и поколением). Как и в ГА, каждая частица хранит свое текущее положение
𝑥𝑖 (𝑘) и соответствующее значение целевой функции 𝐹 𝑜𝑏𝑗 (𝑥), а также вектор скорости
𝑣𝑖 (𝑘) (того же размера, что и вектор положений). Помимо этого, для каждой частицы
хранится лучшее значение целевого функционала, на какой-либо итерации достигнутого ей и соответствующее положение. Целевая функция (Ω-потенциал) чем меньше, тем
лучше.
Число частиц в работе задается равным числу оптимизируемых коэффициентов,
а начальная скорость равна нулю. Начальное положение частиц задается из равномерного распределения.
МРЧ — итерационный алгоритм, в котором новое положение частицы 𝑖 на итерации 𝑘 + 1, 𝑥𝑖 (𝑘 + 1), определяется путем прибавления скорости 𝑣𝑖 (𝑘 + 1) к 𝑥𝑖 (𝑘) :
𝑥𝑖 (𝑘 + 1) = x𝑖 (𝑘) + 𝑣𝑖 (𝑘 + 1). Компоненты 𝑣𝑖 (𝑘 + 1) вычисляются следующим образом:
𝑣𝑖 (𝑘 + 1) = 𝜔𝑣𝑖 (𝑘) + 𝑐1 𝑟1 (𝑦ˆ𝑖 (𝑘) − 𝑥𝑖 (𝑘)) + 𝑐2 𝑟2 (𝑦 * (𝑘) − 𝑥𝑖 (𝑘)) ,
(3.5)
где 𝜔, 𝑐1 , 𝑐2 — некоторые веса, обозначающие степень привлекательности того или иного
направления, 𝑟 — случайные величины в интервале от 0 до 1.
Из формулы (3.5) видно, что скорость складывается из трех компонент, называемых инерциальной, когнитивной и социальной соответственно. Инерциальная компонента 𝜔𝑣𝑖 (𝑘) отвечает за движение частицы в направлении, в котором она двигалась на
предыдущей итерации 𝑘. Когнитивное слагаемое 𝑐1 𝑟1 (𝑦ˆ𝑖 (𝑘) − 𝑥𝑖 (𝑘)) отвечает за движение в зависимости от предыдущего лучшего положения для частицы 𝑖. Социальная
компонента 𝑐2 𝑟2 (𝑦 * (𝑘) − 𝑥𝑖 (𝑘)) включает в себя информацию о лучшем положении для
24
всех частиц и отвечает за движение по направлению к нему. На каждой итерации каждая частица в рое перемещается на новое место в пространстве поиска, и для каждого
положения вычисляется целевая функция, определяющая, насколько хорошим является данное решение.
Алгоритм останавливается по достижении отведенного числа вызовов целевой
функции, а лучшее положение частицы за все итерации считается ответом.
В работе коэффициенты 𝜔 = 0.7298, 𝑐1 = 1.4962, 𝑐2 = 1.4962. Эти значения
являются наиболее универсальными и подходят для многих задач, именно с них стоит
начинать тестировать PSO, стоит отметить, что эти три коэффициента не могут быть
произвольными и связаны между собой соотношением из [42]. В ходе работы изменение
этих параметров не дало лучшего решения.
3.3
Гибридный подход (H-DFT)
Подход VF-DFT не требует вычисление вариации энергии Гельмгольца и позволяет быстрее получать решение, чем классический DFT с методом простой итерации.
Однако, классический DFT дает более точное решение. Чтобы сохранить преимущества
в скорости вычислений и качестве решения от обоих подходов было решено совместить
их в один гибридный алгоритм Hybrid Density Functional Theory (H-DFT) рис. 7, в котором на начальном этапе приближенное решение ищется с помощью VF-DFT и стохастических методов оптимизации. Затем, то решение, которое выдает VF-DFT подается
в качестве начального приближения классическому DFT с методом простой итерации и
решение уточняется. Такая комбинация позволила получить значительный выигрыш в
скорости поиска равновесной плотности в сравнении с классическим DFT (результаты
по скорости и качеству приведены в разделе 4), при этом качество решения осталось
на том же уровне, что и у классического алгоритма. Из-за того, что решение VF-DFT
близко к равновесному, для метода простой итерации можно устанавливать значение
𝛾 достаточно большим. Это влияет на скорость метода простой итерации и на число
необходимых итераций для алгоритма.
Свойства
флюида
𝐹 [𝜌 (𝑟)]
VF-DFT
+
Стохастическая 𝜌0 = 𝜌*
𝑠𝑜
оптимизация
Классический
DFT
+
Метод простой итерации
Равновесная
плотность
𝜌* (𝑧)
Рисунок 7: Схема работы гибридного алгоритма для теории функционала плотности.
4
Результаты работы безвариационного алгоритма
В этом разделе представлены результаты тестирования метода VF-DFT и гибридного метода H-DFT, которые используют ГА и МРЧ для поиска равновесной плот25
ности флюида в поре, с классическим DFT, который использует метод простой итерации. Было проведено сравнение методов по времени поиска решения и значению
Ω потенциала для найденного решения. Методы тестируются на задачах поиска равновесного распределения плотности азота и аргона в поре 3.6 нм. Все необходимые
параметры взаимодействия флюид-флюид и поверхность-флюид приведены в таблице
1. Параметры для взаимодействия поверхности с флюидом были посчитаны по правилу
Лоренца-Бертло (4.1), где 𝜀𝑠𝑠 = 28 K, 𝜎𝑠𝑠 = 3.4 Å
𝜀𝑠𝑓 =
√
𝜀𝑠𝑠 𝜀𝑓 𝑓 ,
𝜎𝑠𝑓 =
1
(𝜎𝑠𝑠 + 𝜎𝑓 𝑓 )
2
(4.1)
Таблица 1: Параметры взаимодействия флюид-флюид и поверхность-флюид. Для потенциала ЛД 𝑟𝑐𝑢𝑡 равен 5𝜎𝑓 𝑓 . Плотность поверхности стенки поры 𝜌𝑉 = 0.114 Å−3 ,
расстояние между слоями атомов углерода ∆ = 3.35 Å.
Флюид-Флюид
Флюид
𝑁2
𝐴𝑟
4.1
Поверхность-Флюид
𝑇, K
𝜀𝑓 𝑓 /𝑘𝐵 , K
𝜎𝑓 𝑓 , Å
𝜀𝑠𝑓 /𝑘𝐵 , K
𝜎𝑠𝑓 , Å
77.4
87.3
94.45
111.95
3.575
3.358
51.43
55.99
3.487
3.379
Азот
В самом простом приближении стоит тестировать алгоритм с определения равновесной плотности для флюида, информации о котором больше содержится в базисе.
Как уже говорилось, базис строился на основе равновесных плотностей азота в поре
3.6 нм при фиксированной температуре 77.4 К и разных относительных давлениях.
Также в базис вошли равновесные плотности «искусственных» флюидов с разными
радиусами молекул (от 0.85 Å, до 2 Å) с фиксированными остальными параметрами
взаимодействия. Во всех задачах ниже используется один и тот же базис.
4.1.1
Высокое давление
При расчете равновесной плотности на больших значениях относительного давления, как правило, классические алгоритмы долго ищут равновесное решение. Именно
в таких задачах безвариационный подход дает наибольший выигрыш по времени расчета. Несмотря на то, что метод простой итерации весьма надежен, для обеспечения
сходимости в таких задачах параметр 𝛾 приходится выбирать очень маленьким, что
сильно влияет на скорость сходимости. На рис. 8 показана плотность азота в поре 3.6
нм при температуре 77.4 K и относительном давлении для флюида в свободном объеме
𝑃/𝑃0 = 0.6924 (𝑃0 — давление насыщения). Важно отметить, что в датасет 𝑋 не входит функция при тех относительных давлениях, при которых производятся расчеты в
26
задачах, рассматриваемых ниже. Разница во времени расчета весьма значительна, при
этом решение достаточно близко по значению Ω потенциала к классическому DFT.
0.2
Время расчета
𝜌, Å−3
МРЧ
ГА
МПИ
590 с
600
400
0.1
200
16 с
МРЧ
0
0
2
4
6
8
10
17 с
ГА МПИ
12
14
𝑧, Å
Рисунок 8: График равновесной плотности при 𝑃/𝑃0 = 0.6924 для азота. Красная
сплошная линия — VF-DFT с МРЧ ΩМРЧ = −0.2630, синяя прерывистая линия —
VF-DFT с ГА ΩГА = −0.2628, желтыми кругами обозначен классический DFT с МПИ
ΩМПИ = −0.2634
Значение Омега потенциала для метода простой итерации убывает значительно
медленнее в сравнении с VF-DFT как видно на рис. 10. Также можно заметить, что
сходимость стохастических алгоритмов не является строго монотонной. Это связано
с принципом работы самих алгоритмов. Для стохастических алгоритмов нет точных
критериев сходимости, а это значит, что решение, которое выдает алгоритм, не обязано
быть лучшим. В таком случае, приходится либо продолжить расчет и надеяться, что
решение улучшится, либо запустить алгоритм заново. От того, как настроен алгоритм,
зависит насколько устойчив он будет и насколько быстро он будет находить решение.
Алгоритмы на основе ГА, МРЧ справились с задачей на порядок быстрее, но
качество полученных решений немного уступает классическому DFT на основе метода простой итерации. Однако стоит отметить, что структурную особенность поведения
флюида VF-DFT удалось воспроизвести хорошо. VF-DFT на основе МРЧ выдал значение Ω потенциала для посчитанной равновесной плотности Ω [𝜌*МРЧ ] = −0.2630, VF-DFT
на основе ГА Ω [𝜌*ГА ] = −0.2628, оба оказались очень близки к значению, которое выдал
классический DFT Ω [𝜌*МПИ ] = −0.2634.
На рис. 9 показан результат работы гибридных методов в сравнении с классическим DFT с методом простой итерации. Время работы немного отличается от VFDFT, но конечное значение Ω потенциала у всех трех профилей плотности совпало
27
0.2
Время расчета
590 с
𝜌, Å−3
Г-МРЧ 600
Г-ГА
МПИ
400
0.1
200
34 с
0
0
2
4
6
37 с
Г-МРЧ Г-ГА
8
10
12
МПИ
14
𝑧, Å
Рисунок 9: График равновесной плотности при 𝑃/𝑃0 = 0.6924 для азота. Красная
сплошная линия — H-DFT с МРЧ, синяя прерывистая линия — H-DFT с ГА, желтыми
кругами обозначен классический DFT с МПИ. Графики полностью совпадают, но время
работы гибридных методов на порядок меньше ΩГ-МРЧ = ΩГ-ГА = ΩМПИ = −0.2634
Ω [𝜌*Г-МРЧ ] = Ω [𝜌*Г-ГА ] = Ω [𝜌*Г-МПИ ] = −0.2634.
Безвариационный подход позволяет получать решение в несколько раз быстрее,
причем объединение VF-DFT с классическим DFT обеспечивает, кроме высокой скорости вычислений, приемлемое качество решения.
28
0.5
Метод простой итерации
1
0.4
0.3
МРЧ
ГА
0.6
0.2
0.2
0.1
−0.2
0
−0.1
1
1.5
2 ×105
−0.2
−0.3
0
3000
6000
9000
Число итераций стохастических алгоритмов
Рисунок 10: Значение Ω потенциала на итерациях для классического DFT с методом
простой итерации (желтая сплошная линия), значение Ω потенциала для VF-DFT с
методом роя частиц (красная сплошная линия) и VF-DFT с генетическим алгоритмом
(синяя пунктирная линия)
4.1.2
Низкое давление
На рис. 11 видно, что при низких относительных давлениях 𝑃/𝑃0 = 0.0044 выигрыш в скорости вычислений у VF-DFT не так значителен, как при высоких давлениях. При низких относительных давлениях классические алгоритмы работают быстрее,
параметр 𝛾 больше. Если посмотреть на полученный профиль плотности, качество решения кажется хуже, чем у классического DFT. Для гибридных алгоритмов рис. 12
оказалось, что генетический алгоритм работает медленнее, чем метод простой итерации, МРЧ по прежнему быстрее обоих (однако стоит учитывать., что в гибридном
подходе тратится время на имплементацию оптимизационных алгоритмов), а значение
Ω потенциала для всех трех одинаково Ω [𝜌*Г-МРЧ ] = Ω [𝜌*Г-ГА ] = Ω [𝜌*МПИ ] = −0.0308.
29
0.2
Время расчета
6
𝜌, Å−3
МРЧ
ГА
МПИ
3.7 с
4
0.1
1.9 с
2
0
0
2
4
6
4.3 с
МРЧ
8
10
ГА МПИ
12
14
𝑧, Å
Рисунок 11: График равновесной плотности при 𝑃/𝑃0 = 0.0044 для азота. Красная
сплошная линия — VF-DFT с МРЧ ΩМРЧ = −0.0305, синяя прерывистая линия —
VF-DFT с ГА ΩГА = −0.0286, желтыми кругами обозначен классический DFT с МПИ
ΩМПИ = −0.0308
0.2
Время расчета
𝜌, Å−3
Г-МРЧ
Г-ГА
МПИ
0.1
6
4.6 с 4.3 с
4
3.1 с
2
0
0
2
4
6
Г-МРЧ Г-ГА
8
10
12
МПИ
14
𝑧, Å
Рисунок 12: График равновесной плотности при 𝑃/𝑃0 = 0.0044 для азота. Красная
сплошная линия — H-DFT с МРЧ, синяя прерывистая линия — H-DFT с ГА, желтыми
кругами обозначен классический DFT с МПИ. Графики полностью совпадают, но время
работы гибридных методов на порядок меньше ΩГ-МРЧ = ΩГ-ГА = ΩМПИ = −0.0308
30
4.2
Аргон
В этом разделе представлены наиболее ценные результаты работы нового алгоритма VF-DFT. Чтобы убедиться в том, что безвариационный подход можно применять
не только для расчета равновесной плотности азота при температуре 77.4 K (на основе азота считался базис), но и для других флюидов, было решено применить VF-DFT
для аргона при температуре 87.3 K. Результаты данного раздела подтверждают идею
применимости безвариационного метода для расчета плотности флюида со сложной
структурой свободной энергии Гельмгольца без вычисления вариаций.
4.2.1
Высокое давление
На рис. 13 видно, что положение пиков совпало не точно, но стоит отметить,
что количество и их амплитуда VF-DFT смог восстановить. Кроме того, VF-DFT по
прежнему справляется с расчетами в разы быстрее, чем классический DFT с методом
простой итерации.
0.2
Время расчета
𝜌, Å−3
МРЧ
ГА
МПИ
45 с
40
20
0.1
12.7 с
0
0
2
4
6
МРЧ
8
10
16.2 с
ГА МПИ
12
14
𝑧, Å
Рисунок 13: График равновесной плотности при 𝑃/𝑃0 = 0.9932 для аргона. Красная
сплошная линия — VF-DFT с МРЧ ΩМРЧ = −0.1674, синяя прерывистая линия —
VF-DFT с ГА ΩГА = −0.1626, желтыми кругами обозначен классический DFT с МПИ
ΩМПИ = −0.1712
При использовании гибридного алгоритма рис. 14 решения полностью совпали
при меньшем времени расчета.
31
0.2
Время расчета
Г-МРЧ
Г-ГА
МПИ
45 с
40
𝜌, Å−3
26.2 с
20
0.1
12.1 с
0
0
2
4
6
Г-МРЧ Г-ГА
8
10
12
МПИ
14
𝑧, Å
Рисунок 14: График равновесной плотности при 𝑃/𝑃0 = 0.9932 для аргона. Красная
сплошная линия — H-DFT с МРЧ, синяя прерывистая линия — H-DFT с ГА, желтыми
кругами обозначен классический DFT с МПИ. Графики полностью совпадают, но время
работы гибридных методов на порядок меньше ΩГ-МРЧ = ΩГ-ГА = ΩМПИ = −0.1712
4.2.2
Среднее давление
С уменьшением относительного давления безвариационный подход дает меньший
выигрыш в скорости рис. 15. В этом кейсе алгоритмам удалось уловить положение
пиков, но величину пиков смог восстановить только VF-DFT с МРЧ.
Гибридный подход в случае с аргоном при средних давлениях рис. 16 не дал
ожидаемого выигрыша в скорости расчета, H-DFT с генетическим алгоритмом оказался даже медленнее классического DFT. Для таких давлений в методе простой итерации
для аргона можно использовать достаточно большое значение 𝛾 и время, которое тратится в H-DFT на получение приближенного решение, оказалось очень близким ко
времени работы классического DFT.
32
0.2
Время расчета
𝜌, Å−3
МРЧ
ГА
МПИ
6
5.5 с
4
2.8 с
0.1
3.7 с
2
0
0
2
4
6
МРЧ
8
ГА МПИ
10
12
14
𝑧, Å
Рисунок 15: График равновесной плотности при 𝑃/𝑃0 = 0.4739 для аргона. Красная
сплошная линия — VF-DFT с МРЧ ΩМРЧ = −0.1049, синяя прерывистая линия —
VF-DFT с ГА ΩГА = −0.1019, желтыми кругами обозначен классический DFT с МПИ
ΩМПИ = −0.1077
0.2
Время расчета
𝜌, Å−3
Г-МРЧ
Г-ГА
МПИ
10
8с
8
6
5.5 с
5с
4
0.1
2
0
0
2
4
6
Г-МРЧ Г-ГА
8
10
12
МПИ
14
𝑧, Å
Рисунок 16: График равновесной плотности при 𝑃/𝑃0 = 0.4739 для аргона. Красная
сплошная линия — H-DFT с МРЧ, синяя прерывистая линия — H-DFT с ГА, желтыми
кругами обозначен классический DFT с МПИ. Графики полностью совпадают, но время
работы гибридных методов на порядок меньше ΩГ-МРЧ = ΩГ-ГА = ΩМПИ = −0.1077
33
Выводы Таким образом, наглядно видно, что безвариационный и гибридный подходы превосходят в скорости расчета классические методы. Также было продемонстрировано, что VF-DFT применим для описания поведения не только азота, но и других
флюидов. Наибольший выигрыш в скорости у разработанных методов наблюдается для
случаев с высоким относительным давлением.
5
5.1
Восстановление распределения пор по размерам
Введение в задачу восстановления распределения пор по
размерам
Сланцевый газ является перспективным нетрадиционным источником углеводородов, месторождения которого найдены по всему миру [43]. В последнее время наблюдается значительный прогресс в оценке и разведке сланцевого газа в Китае [44]. Многочисленные исследования показали, что в поровой системе богатых органикой сланцев
преобладают нанопоры [45]. Нанопоры обычно меньше 100 нм и включают микропоры
(<2 нм), мезопоры (2-50 нм) и макропоры (50-100 нм) рис. 17 [46]. Богатые органи-
Рисунок 17: Изображение образца породы с месторождения сланцевого газа методом
электронной микроскопии (FE-SEM) [46]
кой сланцы имеют сложные и гетерогенные поровые системы с широкой вариабельностью типа пор, размера, формы и механизма течения [47]. Таким образом, исследование
структуры нанопор может помочь лучше понять пути хранения и миграции газа в сланцах.
Для характеризации структуры порвого пространства в сланцах используют различные методы, например, адсорбция газов, ртутная порометрия, полевая эмиссионная электронная микроскопия или трансмиссионная электронная микроскопия (field
emission scanning electron microscopy or transmission electron microscopy — FE-SEM/TEM),
ионная сканирующая электронная микроскопия (focused ion beam scanning electron
microscopy — FIB-SEM) и малоугловое или ультра-малоуглового рассеяния нейтронов
34
(small-angle or ultra-small-angle neutron scattering — SANS/USANS). Каждый метод или
техника применимы к определенному диапазону для характеристики размера пор. Низкотемпературная адсорбция газ стала наиболее популярным методом для характеристики нанопористой структуры сланцев, так как она позволяет оценить полный диапазон
размеров нанопор [43].
В работе [16] представлена методика определения распределения пор по размерам для пористых углеродсодержащих материалов, основанная на применении классической теории функционала плотности к решению интегрального уравнения адсорбции.
В этом подходе измеренное количество зондирующего газа, адсорбированного в пористом материале, вычисляется из свертки между распределением пор по размерам и ядром адсорбции. Эта процедура показала лучшие результаты для описания адсорбции
газа в нанопористых материалах и восстановления распределения пор по размерам,
по сравнению с методами, основанными на уравнении Кельвина, такими как метод
Барретта-Джойнера-Халенды (Barrett-Joyner-Halend — BJH).
Большинство моделей адсорбции газа в нанопорах представляют структуру пор
как совокупность щелевидных пор со стенками из атомов углерода. Это представление
аппроксимирует типичную пластинчатую геометрию пор, наблюдаемую в углеродистых
материалах, геометрическая неоднородность моделируется изменением размера пор.
Таким образом, экспериментальная изотерма описывается как комбинация изотерм в
отдельных щелевидных порах с использованием интегрального уравнения адсорбции
[48]
𝐻
∫︁𝑚𝑎𝑥
𝑁𝑒𝑥𝑝 (𝑃 ) =
𝑁𝑆 (𝑃,𝐻) 𝜒 (𝐻) 𝑑𝐻,
(5.1)
𝐻𝑚𝑖𝑛
где 𝑁𝑒𝑥𝑝 — экспериментальная изотерма адсорбции, 𝐻 — ширина поры, 𝑃 — давление,
𝜒 (𝐻) — распределение пор по размерам в материале, 𝑁𝑆 (𝑃,𝐻) — изотерма адсорбции
в одиночной поре фиксированной ширины. Величина 𝑁𝑆 (𝑃,𝐻) может быть рассчитана
с помощью теории функционала плотности следующим образом:
1
𝑁𝑆 (𝑃,𝐻) =
2𝑁𝐴
∫︁𝐻
(𝜌 (𝑧) − 𝜌𝑏𝑢𝑙𝑘 (𝑃 )) 𝑑𝑧,
(5.2)
0
где 𝑁𝐴 — число Авогадро, 𝜌 (𝑧) — равновесный профиль плотности флюида в поре,
𝜌𝑏𝑢𝑙𝑘 — плотность флюида в свободном объеме при давлении 𝑃 . На рис. 18 представлен характерный вид изотермы адсорбции в одиночной плоской поре, одна точка на
такой кривой по сути является проинтегрированным профилем плотности типа рис. 3
при фиксированном давлении. Рассчет изотермы производился классическим DFT с
методом простой итерации, а результаты были свалидированы с [6].
Обращение уравнения (5.1) для решения обратной задачи и нахождения распределения пор по размерам 𝜒 (𝐻) как известно является плохо поставленной или некорректной задачей. Можно переписать интегральное уравнение (5.1) в матричном виде
35
0.08
𝑁𝑆 (𝑃,𝐻), ммоль/м2
Данные из Ravikovitch 2001
DFT расчет
0.06
0.04
0.02
0
0.1
0.2
0.3
0.4
0.5
𝑃/𝑃0
0.6
0.7
0.8
0.9
1
Рисунок 18: Изотерма адсорбции азота в поре 𝐻𝑐𝑐 = 71.5 Å, рассчитанная с помощью
классического DFT и свалидированная с [6]
[48]
𝐴𝑥 = 𝑏
(5.3)
где 𝐴 (𝑚 × 𝑛) — матрица изотерм для одиночных пор, которая считается с помощью
теории функционала плотности (𝑚 — количество точек по давлению 𝑃 , при которых
считается изотерма для одной поры; 𝑛 — количество одиночных пор, для которых были рассчитаны изотермы), 𝑏 — вектор с экспериментальными значениями изотермы
адсорбции размера (𝑚 × 1), 𝑥 — неизвестный вектор (𝑛 × 1) распределения пор по размерам. На рис. 19 и рис. 20 изображены изотермы адсорбции азота в различных порах,
всего в данной работе 𝑛 = 20 изотерм и 𝑚 = 286 точек по давлению для каждой изотермы. Набор таких изотерм адсорбции для одиночных пор называют ядрами адсорбции
и каждый столбец матрицы 𝐴 является изотермой для одиночной поры.
В данной работе используется метод регуляризации Тихонова, в котором решение
вычисляется путем минимизации следующей квадратичной формы
ℒ = ‖𝐴𝑥 − 𝑏‖22 + 𝜆2 ‖𝐿𝑥‖22 → min
(5.4)
где 𝐿 — линейный оператор (матрица 𝑝 × 𝑛, 𝑝 6 𝑛), который накладывает ограничения
на решение, например, на гладкость решения, а 𝜆 — параметр регуляризации. Выбор
метода регуляризации не уникален и для разных задач может быть своя матрица 𝐿. В
данной работе полагается, что 𝐿 = 𝐼, единичная матрица. Продифференцируем выра-
36
𝑁𝑆 (𝑃,𝐻), ммоль/м2
0.08
0.06
0.04
0.02
0
0.1
0.2
0.3
0.4
0.5
𝑃/𝑃0
0.6
0.7
0.8
0.9
1
3.35 Å
6.20 Å
9.04 Å
11.89 Å
14.73 Å
17.57 Å
20.42 Å
23.26 Å
26.07 Å
28.94 Å
31.75 Å
34.63 Å
37.43 Å
40.28 Å
43.12 Å
45.96 Å
48.80 Å
51.65 Å
54.49 Å
57.33 Å
Рисунок 19: Набор двадцати изотерм адсорбции азота для пор шириной от 𝐻𝑖𝑛 = 3.35 Å
до 𝐻𝑖𝑛 = 57.33 Å
жение (5.4), чтобы выразить искомое распределение пор по размерам 𝑥
𝜕ℒ
= 0 = 𝐴T (𝐴𝑥 − 𝑏) + 𝜆2 𝑥
𝜕𝑥
(5.5)
(︀
)︀−1
𝑥 = 𝐴T 𝑏 𝐴T 𝐴 + 𝜆2 𝐼
.
(5.6)
C помощью сингулярного разложения матрицы 𝐴 можно избежать обращения матриц.
𝐴 = 𝑈 𝑆𝑉
T
=
𝑛
∑︁
𝑢𝑖 𝑠𝑖 𝑣𝑖T
(5.7)
𝑖=1
где 𝑈 = (𝑢1 , . . . , 𝑢𝑛 ) и 𝑉 = (𝑣1 , . . . , 𝑣𝑛 ) ортогональные матрицы, а 𝑆 — диагональная
матрица сингулярных чисел матрицы 𝐴, которые расположены в порядке убывания.
Таким образом, решение (5.4) дается выражением [48]
𝑥=
𝑛
∑︁
𝑖=1
𝑠2𝑖 𝑢T
𝑖 𝑏
𝑣𝑖 .
2
2
𝑠𝑖 + 𝜆 𝑠𝑖
37
(5.8)
Матрица 𝐴
𝑁𝑆 (𝑃,𝐻), ммоль/м2
0.08
0.06
0.08
0.06
100
0.04
200
0.02
0.04
5
10
15
20
0
0.02
10−7
10−6
10−5
10−4
10−3
𝑃/𝑃0
10−2
10−1
100
3.35 Å
6.20 Å
9.04 Å
11.89 Å
14.73 Å
17.57 Å
20.42 Å
23.26 Å
26.07 Å
28.94 Å
31.75 Å
34.63 Å
37.43 Å
40.28 Å
43.12 Å
45.96 Å
48.80 Å
51.65 Å
54.49 Å
57.33 Å
Рисунок 20: Набор двадцати изотерм адсорбции азота для пор шириной от 𝐻𝑖𝑛 = 3.35 Å
до 𝐻𝑖𝑛 = 57.33 Å в логарифмическом масштабе. На вставленном изображении продемонстрирована матрица ядер адсорбции 𝐴, где цветом отмечено значение изотермы
для заданной поры и относительного давления
5.2
Результаты восстановления распределения пор по размерам
В качестве иллюстрации работы метода восстановления пор по размерам, рассматривается несколько модельных задач. В этих задачах будем задавать распределение пор по размерам, и с помощью (5.1) считать «экспериментальную» изотерму адсорбции. Затем, на основании (5.8) будет рассчитываться искомое распределение, в такой
модельной задаче есть несколько параметров, которые будут влиять на ответ, изменяя
их можно увидеть сильные и слабые стороны данного подхода. По своей природе, распределение пор по размерам подчиняется гауссовскому распределению, поэтому именно
такие распределения и будут исследоваться ниже
(︃
)︃
(𝑥 − 𝜇)2
1
exp −
.
𝑓 (𝑥) = √
2𝜎 2
2𝜋𝜎
5.2.1
(5.1)
Влияние параметров распределения на вид изотермы
В данном разделе приводится иллюстрация того, насколько будет меняться вид
рассчитанной изотермы адсорбции, при изменении параметров распределения пор по
38
Распределение пор по размерам pdf
размерам (по сути это отражает чувствительность приведенного в работе метода восстановления распределения пор по размерам к начальным данным). На рис. 21 показаны
два распределения с одинаковой дисперсией 𝜎1 = 𝜎2 = 7, но с различным математическим ожиданием 𝜇1 = 25, 𝜇2 = 40. При этом, после свертки данных распределений
с ядрами адсорбции из матрицы 𝐴 получаются две различные изотермы. Из рис. 22
видно, что расчетная изотерма адсорбции чувствительна к значению математического
ожидания распределения пор по размерам.
𝜇1 = 25
𝜇2 = 40
0.05
0.03
0.01
0
5
10
15
20
25
30
35
40
45
50
55
60
Ширина пор, Å
Рисунок 21: Два различны распределения пор по размерам с одинаковыми дисперсиями
𝜎1 = 𝜎2 = 7, но с различными математическим ожиданием 𝜇1 = 25, 𝜇2 = 40
39
0.02
𝜇1 = 25
𝜇2 = 40
𝑁𝑒𝑥𝑝 , ммоль/м2
0.016
0.012
0.008
0.004
0
0.1
0.2
0.3
0.4
0.5
P/P0
0.6
0.7
0.8
0.9
1
Рисунок 22: Расчетные изотермы адсорбции для распределений пор по размерам с одинаковой дисперсией, но различными математическими ожиданиями 𝜇1 = 25, 𝜇2 = 40
Дисперсия распределения также влияет на вид расчетной изотермы. На рис. 23
показаны два распределения с одинаковыми математическими ожиданиями 𝜇1 = 𝜇2 ,
но с различными дисперсиями 𝜎1 = 8, 𝜎2 = 4. После свертки таких распределений с
ядрами адсорбции рис. 19 получились две различные изотермы рис. 24.
40
Распределение пор по размерам pdf
𝜎1 = 8
𝜎2 = 4
0.09
0.07
0.05
0.03
0.01
0
5
10
15
20
25
30
35
40
45
50
55
60
Ширина пор, Å
Рисунок 23: Два различны распределения пор по размерам с одинаковым математическим ожиданием 𝜇1 = 𝜇2 = 33, но с различными дисперсиями 𝜎1 = 8, 𝜎2 = 4
𝜎1 = 8
𝜎2 = 4
𝑁𝑒𝑥𝑝 , ммоль/м2
0.016
0.012
0.008
0.004
0
0.1
0.2
0.3
0.4
0.5
P/P0
0.6
0.7
0.8
0.9
1
Рисунок 24: Расчетные изотермы адсорбции для распределений пор по размерам с одинаковой математическими ожиданиями, но различными дисперсиями 𝜎1 = 8, 𝜎2 = 4
5.2.2
Зависимость решения от параметр регуляризации
Отсутствие регуляризации 𝜆 = 0. В самом простом приближении стоит рассмотреть задачу, когда регуляризационный параметр отсутствует и 𝜆 = 0. В таком случае
41
восстановленное распределение должно полностью совпасть с изначально заданным
распределением рис. 26. Исходное распределение было смоделировано по гауссовскому
закону распределения с параметрами 𝜇𝑖𝑛𝑖𝑡 = 33, 𝜎𝑖𝑛𝑖𝑡 = 7, у восстановленного распределения получились параметры 𝜇𝑟𝑒𝑐 = 32.99, 𝜎𝑟𝑒𝑐 = 6.99. На рис. 25 представлена
«экспериментальная» изотерма адсорбции в образце и восстановленная, которая рассчитывалась как свертка восстановленного распределения с ядрами адсорбции, видно
хорошее соответствие между двумя изотермами.
Исходная
Восстановленная
𝑁𝑒𝑥𝑝 , ммоль/м2
0.016
0.012
0.008
0.004
0
0.1
0.2
0.3
0.4
0.5
P/P0
0.6
0.7
0.8
0.9
1
Рисунок 25: Сравнение между восстановленной изотермой адсорбцией в образце (синяя
линия) и исходной «экспериментальной» изотермой для случая 𝜆 = 0
42
Распределение пор по размерам pdf
Восстановленное
Начальное
0.05
0.03
0.01
0
5
10
15
20
25
30
35
40
45
50
55
60
Ширина пор, Å
Рисунок 26: Сравнение восстановленного (синий цвет) и исходного (красный цвет) распределения пор по размером при 𝜆 = 0
Решение при регуляризации 𝜆 = 0.07. При изменении параметра регуляризации
исходное и восстановленное распределение получаются отличными друг от друга, причем при увеличении 𝜆 у восстановленного распределения увеличивается дисперсия,
а восстановленная изотерма адсорбции становится более гладкой. У исходного распределения были параметры 𝜇𝑖𝑛𝑖𝑡 = 33, 𝜎𝑖𝑛𝑖𝑡 = 7, а у восстановленного получились
𝜇𝑟𝑒𝑐 = 32.76, 𝜎𝑟𝑒𝑐 = 8.17. На рис. 27 продемонстрирована восстановленная изотерма
адсорбции и первоначальная изотерма, видно, что соответствие не такое хорошее как
на рис. 25 для случая 𝜆 = 0. На рис. 28 показаны распределения пор по размерам для
случая 𝜆 = 0.07.
43
Исходная
Восстановленная
𝑁𝑒𝑥𝑝 , ммоль/м2
0.016
0.012
0.008
0.004
0
0.1
0.2
0.3
0.4
0.5
P/P0
0.6
0.7
0.8
0.9
1
Распределение пор по размерам pdf
Рисунок 27: Сравнение между восстановленной изотермой адсорбцией в образце (синяя
линия) и исходной «экспериментальной» изотермой для случая 𝜆 = 0.07
Восстановленное
Начальное
0.05
0.03
0.01
0
5
10
15
20
25
30
35
40
45
50
55
60
Ширина пор, Å
Рисунок 28: Сравнение восстановленного (синий цвет) и исходного (красный цвет) распределения пор по размером при 𝜆 = 0.07
5.2.3
Восстановление сдвинутого и бимодального распределения
Для случая, когда распределение не центрировано, метод из раздела выше достаточно хорошо справляется с восстановлением распределения пор по размерам. На
44
Распределение пор по размерам pdf
рис. 29 приведены результаты работы метода при параметрах 𝜆 = 0.07, 𝜇𝑖𝑛𝑖𝑡 = 45,
𝜎𝑖𝑛𝑖𝑡 = 4, В результате получилось распределение с параметрами 𝜇𝑟𝑒𝑐 = 45.05, 𝜎𝑟𝑒𝑐 = 5.1.
При этом стоит отметить, что как и для самого первого случая с 𝜆 = 0, при выборе
такого значения параметра регуляризации начальное и восстановленное распределения
полностью совпадают.
Восстановленное
Начальное
0.1
0.08
0.06
0.04
0.02
0
5
10
15
20
25
30
35
40
45
50
55
60
Ширина пор, Å
Рисунок 29: Сравнение восстановленного (синий цвет) и исходного (красный цвет) распределения пор по размером при 𝜆 = 0.07, 𝜇𝑖𝑛𝑖𝑡 = 45, 𝜎𝑖𝑛𝑖𝑡 = 4, 𝜇𝑟𝑒𝑐 = 45.05, 𝜎𝑟𝑒𝑐 = 5.1
Чаще всего, реальное распределение пор по размерам является нормальным бимодальным. Для того, чтобы проверить работу алгоритма восстановления распределения пор по размерам было задано исходное распределение следующим образом:
(︃
2
1
(𝑥 − 𝜇1 )
𝑓 (𝑥) = √
exp −
2𝜎12
2𝜋𝜎1
)︃
(︃
+√
2
5
(𝑥 − 𝜇2 )
exp −
2𝜎22
2𝜋𝜎2
)︃
,
(5.2)
1
2
где параметры распределения 𝜇1𝑖𝑛𝑖𝑡 = 20, 𝜎𝑖𝑛𝑖𝑡
= 4, 𝜇2𝑖𝑛𝑖𝑡 = 40, 𝜎𝑖𝑛𝑖𝑡
= 5, а параметр
регуляризации в алгоритме был выбран 𝜆 = 0.07. На рис. 30 представлен результат
работы алгоритма, в сравнении с начальным распределением.
45
Распределение пор по размерам pdf
0.4
Восстановленное
Начальное
0.3
0.2
0.1
0
5
10
15
20
25
30
35
40
45
50
55
60
Ширина пор, Å
Рисунок 30: Сравнение восстановленного (синий цвет) и исходного (красный цвет) распределения пор по размером для случая с бимодальным распределением при 𝜆 = 0.07,
1
2
𝜇1𝑖𝑛𝑖𝑡 = 20, 𝜎𝑖𝑛𝑖𝑡
= 4, 𝜇2𝑖𝑛𝑖𝑡 = 40, 𝜎𝑖𝑛𝑖𝑡
=5
5.2.4
Зависимость решения от качества экспериментальных данных
Данная модельная задача демонстрирует как важна точность экспериментальных данных. С помощью заданного нормального распределения с параметрами 𝜇𝑖𝑛𝑖𝑡 =
33, 𝜎𝑖𝑛𝑖𝑡 = 7 была рассчитана «экспериментальная» изотерма адсорбции. Затем, к
изотерме была добавлена случайна величина, равномерно распределенная на отрезке
[10−5 , 10−4 ]. В результате получилась изотерма, которая показана на рис. 32 красным
цветом. С помощью зашумленной изотермы было восстановлено распределение пор по
размерам рис. 31. Для алгоритма восстановления распределения пор по размерам был
выбран параметр 𝜆 = 0.05. При таком значении регуляризационного параметра посчитанная изотерма адсорбции лучше всего совпадает с исходной изотермой. В такой
задаче, с зашумленными данными у восстановленного распределения появляются артефакты на хвостах распределения и регуляризационный параметр 𝜆 требует более
тонкой настройки.
46
Распределение пор по размерам pdf
Восстановленное
Начальное
0.06
0.04
0.02
0
5
10
15
20
25
30
35
40
45
50
55
60
Ширина пор, Å
Рисунок 31: Сравнение восстановленного (синий цвет) и исходного (красный цвет)
распределения пор по размером для случая с зашумленной изотермой при 𝜆 = 0.05,
𝜇𝑖𝑛𝑖𝑡 = 33, 𝜎𝑖𝑛𝑖𝑡 = 7, 𝜇𝑟𝑒𝑐 = 32.5, 𝜎𝑟𝑒𝑐 = 8.2
Исходная с шумом
Восстановленная
𝑁𝑒𝑥𝑝 , ммоль/м2
0.016
0.012
0.008
0.004
0
0.1
0.2
0.3
0.4
0.5
P/P0
0.6
0.7
0.8
0.9
1
Рисунок 32: Сравнение между восстановленной изотермой адсорбцией в образце (синяя
линия) и исходной зашумленной изотермой для случая 𝜆 = 0.05
Выводы Данный раздел не содержит в себе научной новизны, но наглядно демонстрирует применение теории функционала плотности для исследования структуры нанопористых материалов. Методы восстановление распределения пор по размерам опи47
раются на результаты DFT по адсорбции флюида в изолированной поре фиксированной
ширины. Усовершенствование моделей DFT, принципов расчета равновесной плотности
в поре и скорость этих расчетов становятся существенными для быстрого и качественного описания порового пространства. Разработанные в данной работе методы VF-DFT
и H-DFT позволят исследовать сложные флюиды, учитывать различную геометрию
порового пространства и значительно сократить время расчета изотермы адсорбции
флюида в поре заданной ширины.
6
Заключение
Являясь универсальным методом статистической механики, классическая теория
функционала плотности предлагает мощную альтернативу множеству традиционных
теоретических методов и молекулярному моделированию для связывания микроскопических свойств физико-химических систем со структурными и термодинамическими
свойствами. Практическая ценность DFT проявляется не только в ее универсальности
для описания термодинамических свойств классических систем, но и в универсальности для решения задач, которые не могут быть достигнуты традиционными теориями. Разработанный в данной работе безвариационный подход для поиска равновесной
плотности может позволить в дальнейшем описывать сложные взаимодействия флюидов, для которых представляется трудоемким вычисления вариации свободной энергии
Гельмгольца. Комбинирование безвариационного подхода с методом простой итерации
позволило в разы ускорить расчет равновесной плотности при высоких относительных
давлениях, но с незначительной потерей качества решения (∼ 5%). Наибольший эффект
от ускорения расчетов наблюдается для случаев с высоким относительным давлением.
В дальнейшем, комбинированный алгоритм можно использовать для ускорения расчета изотерм адсорбции. При этом стоит отметить, что скорость и качество решения,
которые показывает безвариационный подход VF-DFT, напрямую зависет от базиса.
Размер базиса и вид базисных функций можно постоянно совершенствовать и дорабатывать, чтобы получать решения более высокого качества. Кроме того, на скорость
работы алгоритма влияет настройка параметров оптимизационных алгоритмов. Было
проведено исследование по зависимости скорости работы стохастических алгоритмов
в данных задачах и предложены рекомендации по первоначальной настройке. В дальнейшем предлагается исследовать безвариационный подход в комбинации с другими
стохастическими методами оптимизации или с любыми другими методами оптимизации нулевого порядка.
Как одно из наиболее важных применений теории функционала плотности, в
данной работе была рассмотрена задача восстановления распределения пор по размерам на основании экспериментальных данных по низкотемпературной адсорбции. Рассмотренный подход для решению данной проблемы показал себя весьма действенным
и наглядным для описания структуры порового пространства.
48
Список литературы
1. New insights into application of nanoparticles for water-based enhanced oil recovery in
carbonate reservoirs / S. Rezaei Gomari [и др.] // Colloids and Surfaces A: Physicochemical
and Engineering Aspects. — 2019. — т. 568, January. — с. 164—172. — ISSN 18734359. —
DOI: 10.1016/j.colsurfa.2019.01.037. — URL: https://doi.org/10.1016/j.
colsurfa.2019.01.037.
2. Carbon dioxide capture in metal-organic frameworks / K. Sumida [и др.] // Chemical
Reviews. — 2012. — т. 112, № 2. — с. 724—781. — ISSN 00092665. — DOI: 10.1021/
cr2003272.
3. Atomistic insights into the nanofluid transport through an ultra-confined capillary /
X. Wang [и др.] // Physical Chemistry Chemical Physics. — 2018. — т. 20, № 7. —
с. 4831—4839. — ISSN 14639076. — DOI: 10.1039/c7cp08140e.
4. Wu J. Density functional theory for chemical engineering: From capillarity to soft
materials // AIChE Journal. — 2006. — март. — т. 52, № 3. — с. 1169—1193. —
ISSN 00011541. — DOI: 10.1002/aic.10713.
5. Pore size analysis of MCM-41 type adsorbents by means of nitrogen and argon adsorption /
A. V. Neimark [и др.] // Journal of Colloid and Interface Science. — 1998. — т. 207,
№ 1. — с. 159—169. — ISSN 00219797. — DOI: 10.1006/jcis.1998.5748.
6. Ravikovitch P. I., Vishnyakov A., Neimark A. V. Density functional theories and
molecular simulations of adsorption and phase transitions in nanopores // Physical
Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics. —
2001. — т. 64, № 1. — с. 20. — ISSN 1063651X. — DOI: 10.1103/PhysRevE.64.011602.
7. Yu Y. X., Wu J., Gao G. H. Density-functional theory of spherical electric double
layers and 𝜁 potentials of colloidal particles in restricted-primitive-model electrolyte
solutions // Journal of Chemical Physics. — 2004. — т. 120, № 15. — с. 7223—7233. —
ISSN 00219606. — DOI: 10.1063/1.1676121.
8. Berim G. O., Ruckenstein E. Nanodrop on a nanorough solid surface: Density functional
theory considerations // Journal of Chemical Physics. — 2008. — т. 129, № 1. — ISSN
00219606. — DOI: 10.1063/1.2951453.
9. Aslyamov T., Khlyupin A. Density functional theory formulation for fluid adsorption
on correlated random surfaces // Journal of Chemical Physics. — 2017. — т. 147, №
15. — с. 1—12. — ISSN 00219606. — DOI: 10.1063/1.4997001.
10. Quenched solid density functional theory and pore size analysis of micro-mesoporous
carbons / A. V. Neimark [и др.] // Carbon. — 2009. — июнь. — т. 47, № 7. — с. 1617—
1628. — ISSN 00086223. — DOI: 10.1016/j.carbon.2009.01.050.
49
11. Jagiello J., Olivier J. P. 2D-NLDFT adsorption models for carbon slit-shaped pores
with surface energetical heterogeneity and geometrical corrugation // Carbon. — 2013. —
т. 55, № 2. — с. 70—80. — ISSN 00086223. — DOI: 10.1016/j.carbon.2012.12.011. —
URL: http://dx.doi.org/10.1016/j.carbon.2012.12.011.
12. Solvation forces between molecularly rough surfaces / K. Yang [и др.] // Journal of
Colloid and Interface Science. — 2011. — т. 362, № 2. — с. 382—388. — ISSN 00219797. —
DOI: 10.1016/j.jcis.2011.06.056.
13. Brunauer S., Emmett P. H., Teller E. Adsorption of gases in multimolecular layers //
Journal of the American chemical society. — 1938. — т. 60, № 2. — с. 309—319.
14. Barrett E. P., Joyner L. G., Halenda P. P. The determination of pore volume and
area distributions in porous substances. I. Computations from nitrogen isotherms //
Journal of the American Chemical society. — 1951. — т. 73, № 1. — с. 373—380.
15. Ebner C., Saam W., Stroud D. Density-functional theory of simple classical fluids. I.
Surfaces // Physical Review A. — 1976. — т. 14, № 6. — с. 2264.
16. A new analysis method for the determination of the pore size distribution of porous
carbons from nitrogen adsorption measurements / N. Seaton, J. Walton [и др.] //
Carbon. — 1989. — т. 27, № 6. — с. 853—861.
17. Cohan L. H. Sorption hysteresis and the vapor pressure of concave surfaces // Journal
of the American Chemical Society. — 1938. — т. 60, № 2. — с. 433—435.
18. Gas I. E. P. Hohenberg and W. Kohn // Phys. Rev. B. — 1964. — т. 136. — с. 864.
19. Roth R. Fundamental measure theory for hard-sphere mixtures: A review. — 2010. —
DOI: 10.1088/0953-8984/22/6/063102.
20. Sears M. P., Frink L. J. A new efficient method for density functional theory calculations
of inhomogeneous fluids // Journal of Computational Physics. — 2003. — т. 190, № 1. —
с. 184—200. — ISSN 00219991. — DOI: 10.1016/S0021-9991(03)00270-5.
21. Wu J. Variational Methods in Molecular Modeling / под ред. J. Wu. — Singapore :
Springer Singapore, 2017. — с. 324. — (Molecular Modeling and Simulation). — ISBN
978-981-10-2500-6. — DOI: 10 . 1007 / 978 - 981 - 10 - 2502 - 0. — URL: https : / /
www . amazon . com / Variational - Methods - Molecular - Modeling - Simulation / dp /
9811025002 / ref = sr _ 1 _ 4 ? s = english - books & ie = UTF8 & qid = 1475122294 & sr =
1 - 4 & keywords = phase + field % 5C % 5Cnhttp : / / www . springer . com / la / book /
9789811025006%20http://link.springer.com/10.1007/978-981-.
22. Edelmann M., Roth R. A numerical efficient way to minimize classical density functional
theory // Journal of Chemical Physics. — 2016. — т. 144, № 7. — ISSN 00219606. —
DOI: 10.1063/1.4942020. — URL: http://dx.doi.org/10.1063/1.4942020.
50
23. Rosenfeld Y. Free-energy model for the inhomogeneous hard-sphere fluid mixture and
density-functional theory of freezing // Physical Review Letters. — Beer-Sheva, 1989. —
авг. — т. 63, № 9. — с. 980—983. — ISSN 0031-9007. — DOI: 10.1103/PhysRevLett.
63.980. — URL: https://link.aps.org/doi/10.1103/PhysRevLett.63.980.
24. Weeks J. D., Chandler D., Andersen H. C. Role of repulsive forces in determining the
equilibrium structure of simple liquids // The Journal of Chemical Physics. — 1971. —
т. 54, № 12. — с. 5237—5247. — ISSN 00219606. — DOI: 10.1063/1.1674820.
25. Chandler D., Weeks J. D., Andersen H. C. Van der Waals picture of liquids, solids,
and phase transformations // Science. — 1983. — т. 220, № 4599. — с. 787—794.
26. Steele W. A. The interaction of gases with solid surfaces. т. 3. — Pergamon, 1974.
27. Application of the SAFT-𝛾 Mie group contribution equation of state to fluids of relevance
to the oil and gas industry / V. Papaioannou [и др.] // Fluid Phase Equilibria. —
2016. — т. 416. — с. 104—119. — ISSN 03783812. — DOI: 10.1016/j.fluid.2015.12.
041. — URL: http://dx.doi.org/10.1016/j.fluid.2015.12.041.
28. Accurate statistical associating fluid theory for chain molecules formed from Mie segments /
T. Lafitte [и др.] // Journal of Chemical Physics. — 2013. — т. 139, № 15. — ISSN
00219606. — DOI: 10.1063/1.4819786.
29. Elizarev M., Mukhin A., Khlyupin A. Objective-Sensitive Principal Component Analysis
for High-Dimensional Inverse Problems // arXiv preprint arXiv:2006.04527. — 2020.
30. Abdi H., Williams L. J. Principal component analysis. wiley interdisciplinary reviews:
computational statistics // Wiley Interdisplinary Reviews: Computational Statistics. —
2010. — с. 1—47.
31. Sun W., Durlofsky L. J. A New Data-Space Inversion Procedure for Efficient Uncertainty
Quantification in Subsurface Flow Problems // Mathematical Geosciences. — 2017. —
т. 49, № 6. — с. 679—715. — ISSN 18748953. — DOI: 10.1007/s11004-016-9672-8.
32. Efficient real-time reservoir management using adjoint-based optimal control and model
updating / P. Sarma [и др.] // Computational Geosciences. — 2006. — т. 10, № 1. —
с. 3—36. — ISSN 14200597. — DOI: 10.1007/s10596-005-9009-z.
33. Eberhart, Yuhui Shi. Particle swarm optimization: developments, applications and resources //
Proceedings of the 2001 Congress on Evolutionary Computation (IEEE Cat. No.01TH8546).
т. 1. — IEEE, 2001. — с. 81—86. — ISBN 0-7803-6657-3. — DOI: 10.1109/CEC.2001.
934374. — URL: http://ieeexplore.ieee.org/document/934374/.
34. Montes G., Bartolome P., Udias A. The Use of Genetic Algorithms in Well Placement
optimization // Proceedings of SPE Latin American and Caribbean Petroleum Engineering
Conference. — Society of Petroleum Engineers, 03.2001. — с. 1—10. — DOI: 10.2523/
69439 - MS. — URL: http : / / www . spe . org / elibrary / servlet / spepreview ? id =
00069439.
51
35. Genetic algorithm-based pore network extraction from micro-computed tomography
images / A. Nejad Ebrahimi [и др.] // Chemical Engineering Science. — 2013. — т.
92. — с. 157—166. — ISSN 00092509. — DOI: 10.1016/j.ces.2013.01.045. — URL:
http://dx.doi.org/10.1016/j.ces.2013.01.045.
36. Onwunalu J. E., Durlofsky L. J. Development and application of a new well pattern
optimization algorithm for optimizing large-scale field development // Proceedings SPE Annual Technical Conference and Exhibition. — 2009. — т. 3. — с. 1926—1943. —
DOI: 10.2118/124364-ms.
37. Huang J., Cai Y., Xu X. A hybrid genetic algorithm for feature selection wrapper
based on mutual information // Pattern Recognition Letters. — 2007. — т. 28, № 13. —
с. 1825—1844. — ISSN 01678655. — DOI: 10.1016/j.patrec.2007.05.011.
38. Ouellette R., Browne M., Hirasawa K. Genetic algorithm optimization of a convolutional
neural network for autonomous crack detection // Proceedings of the 2004 Congress
on Evolutionary Computation, CEC2004. — 2004. — т. 1. — с. 516—521. — DOI:
10.1109/cec.2004.1330900.
39. Sharma S., Singh H., Balint-Kurti G. G. Genetic algorithm optimization of laser pulses
for molecular quantum state excitation // Journal of Chemical Physics. — 2010. — т.
132, № 6. — ISSN 00219606. — DOI: 10.1063/1.3314223.
40. Kennedy J., Eberhart R. Particle swarm optimization // Proceedings of ICNN’95 International Conference on Neural Networks. — IEEE, 1995. — ISBN 0-7803-27683. — DOI: 10.1109/ICNN.1995.488968. — URL: http://ieeexplore.ieee.org/
document/488968/.
41. Banks A., Vincent J., Anyakoha C. A review of particle swarm optimization. Part I:
Background and development // Natural Computing. — 2007. — т. 6, № 4. — с. 467—
484. — ISSN 15677818. — DOI: 10.1007/s11047-007-9049-5.
42. Trelea I. C. The particle swarm optimization algorithm: Convergence analysis and
parameter selection // Information Processing Letters. — 2003. — т. 85, № 6. — с. 317—
325. — ISSN 00200190. — DOI: 10.1016/S0020-0190(02)00447-7.
43. DFT modeling of CO2 and Ar low-pressure adsorption for accurate nanopore structure
characterization in organic-rich shales / L. Zhang [и др.] // Fuel. — 2017. — т. 204. —
с. 1—11.
44. Geological characteristics and resource potential of shale gas in China / C. Zou [и
др.] // Petroleum exploration and development. — 2010. — т. 37, № 6. — с. 641—653.
45. Morphology, genesis, and distribution of nanometer-scale pores in siliceous mudstones
of the Mississippian Barnett Shale / R. G. Loucks [и др.] // Journal of sedimentary
research. — 2009. — т. 79, № 12. — с. 848—861.
52
46. Chalmers G. R., Bustin R. M., Power I. M. Characterization of gas shale pore systems
by porosimetry, pycnometry, surface area, and field emission scanning electron microscopy/
transmission electron microscopy image analyses: Examples from the Barnett, Woodford,
Haynesville, Marcellus, and Doig unitsCharacterization of Gas Shale Pore Systems //
AAPG bulletin. — 2012. — т. 96, № 6. — с. 1099—1119.
47. Pore structure characterization of North American shale gas reservoirs using USANS/SANS,
gas adsorption, and mercury intrusion / C. R. Clarkson [и др.] // Fuel. — 2013. — т.
103. — с. 606—616.
48. Unified approach to pore size characterization of microporous carbonaceous materials
from N2, Ar, and CO2 adsorption isotherms / P. I. Ravikovitch [и др.] // Langmuir. —
2000. — март. — т. 16, № 5. — с. 2311—2320. — ISSN 07437463. — DOI: 10.1021/
la991011c.
53
Отзывы:
Авторизуйтесь, чтобы оставить отзыви хорошего настроения
удачи
успехов в конкурсе
Наверное было затрачено много времени и труда на работу
Продолжай свое исследование
Админам респект
И продвижения статьи в топы?
Как на счет взаимных комментариев под работами?)
Красиво написанная работа
Так держать
Молодец
Интересная работа!