Моделирование адсорбции молекул на ступенчатой поверхности в случае квадратной решетки Текст научной статьи по специальности «Математика»

ХИМИЧЕСКИЕ НАУКИ

УДК 541124/128 С. И. ЕВСЕЕВА

А. В. МЫШЛЯВЦЕВ М. Д. МЫШЛЯВЦЕВА

Омский государственный технический университет Институт проблем переработки углеводородов СО РАН

МОДЕЛИРОВАНИЕ АДСОРБЦИИ МОЛЕКУЛ НА СТУПЕНЧАТОЙ ПОВЕРХНОСТИ В СЛУЧАЕ КВАДРАТНОЙ РЕШЕТКИ

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

1. Модель и метод ные взаимодействия. Моделирование неоднородных

Большинство реальных поверхностей являются поверхностей, в том числе и ступенчатых поверхностей,

неоднородными [1-3]. Причиной неоднородности, в значительно сложнее, чем моделирование однородных

частности, могут быть точечные дефекты. Столь же часто систем [5-7].

встречаются и линейные дефекты, среди которых особое В работе авторов [8] была построена и исследо-

место занимают ступеньки [1]. Адсорбция простых моле- вана простейшая модель ступенчатой поверхности, кул на неоднородных поверхностях изучалась, начиная с основанной на стандартной модели решеточного газа

30-х годов ХХ века, и тогда же были построены простые [9,10] на квадратной решетке, бесконечной в обоих теоретические модели, приводящие при различных направлениях. Поверхность, соответствующая этой

предположениях к классическим изотермам Темкина, модели, схематически показана на рис. 1. Ширину Фрейндлиха-Зельдовича и т.д. [4]. Во всех этих моделях террасы обозначали через Ь (Ь > 2) и считали, что

рассматривались только случайно распределенные узлы первого, (Ь + 1)-го, (2Ь + 1)-го и следующих

точечные неоднородности и не учитывались латераль- рядов вдоль Х направления обладают свойствами, от-

личными от свойств остальных узлов. Рассматривали мономолекулярную адсорбцию молекул одного сорта. Предполагали, что узлы каждого первого ряда (адсорбционные центры, расположенные вблизи ступенек) обладают потенциальной ямой глубины А (сильные адсорбционные центры). Энергия адсорбента для частиц, адсорбированных на невыделенные узлы, равнялась нулю (слабые адсорбционные центры). Считали, что взаимодействуют только ближайшие адсорбированные частицы и энергия взаимодействия е одинакова для всех возможных расположений взаимодействующей пары. Моделирование проводилось с использованием метода трансфер-матрицы [10,11]. Были построены фазовые диаграммы и локальные изотермы построенной модели. Было показано, что в рамках исследуемой модели влияние ступеньки распространяется в глубь террасы не более чем на 2^3 атомных расстояния. Показано, что при 2 < А/е < 3 и L > 3 заполнение первого выделенного ряда зависит от химического потенциала немонотонно, что может оказать сильное влияние на кинетику химических реакций, протекающих на ступенчатых поверхностях. Установлено отсутствие упорядоченных структур для данной модели при глубине потенциальной ямы, более чем в четыре раза превосходящей величину энергии латерального отталкивания.

В данной работе исследуется модель, отличная от предыдущей лишь тем, что узлы каждого первого ряда и последнего ряда (рис. 1) обладают потенциальной ямой глубины А1 (А1 > 0) и высоты А (А < 0) соответственно. В данной модели ширина террасы L > 3, так как случай L = 2 аналогичен исследованному случаю [8]. Моделирование проводилось с помощью метода трансфер-матрицы (МТМ) [10,11] и метода Монте-Карло (МК) [12].

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

(1)

<т> з' к к

где ц - химический потенциал, <пп> означает суммирование по всем парам ближайших соседей. Во втором члене суммирование выполняется по всем узлам решетки, в третьем члене суммирование распространяется на сильные адсорбционные центры, в четвертом члене по выделенным последним рядам (рис. 1).

Метод трансфер-матрицы.

Построенная модель неоднородна, но транс-ляционно-инвариантна, что позволяет использовать метод трансфер-матрицы [10,11]. Для изучаемой модели можем записать следующее соотношение:

Jr(K) (2)

где ü v - большая статистическая сумма для решетки из N рядов с циклическими граничными условиями, Ti — взятие следа. Здесь

Рис. 1. Схематический вил ступенчатой поверхности при L = 3. Числами показаны номера рядов. Кружочками показаны адсорбированные молекулы

В, - Т, Т,, Т,3Т31,

(3)

где ТгТг,~3Т,;Т31. - стандартно определенные трансфер-матрицы, соответствующие переходу от выделенного первого ряда к невыделенному ряду, переходу между двумя невыделенными рядами, переходу от невыделенного ряда к выделенному последнему ряду и переходу от выделенного последнего ряда к выделенному первому ряду. Матрица В, соответствует переходу между двумя ближайшими выделенными первыми рядами, т.е. от первого ряда к (I + 1)-му ряду. Отметим, что для рассматриваемого случая матрица В(, определенная формулой (3), является несимметричной.

Выражение (2) может быть записано для любой матрицы:

(4)

полученной циклической перестановкой матриц в произведении (3). Матрица В соответствует переходу от г -го ряда к (I + г)-му ряду и позволяет вычислить вероятности различных конфигураций адсорбированных частиц в г -х рядах.

Всего будет I различных матриц В , соответствующих переходам от первого ряда к (I + 1)-му ряду, от второго ряда к (I + 2)-му ряду и т.д. Все эти матрицы имеют тождественный спектр, но, вообще говоря, различные собственные векторы. Собственные векторы матриц В , соответствующие наибольшему собственному значению, описывают распределение вероятностей расположения адсорбированных частиц в г -м ряду. Исходя из сказанного, мы использовали следующие расчетные формулы:

Г, RT, 7

П = --¡n 1 ,

LM

(5)

где О - большой термодинамический потенциал, приходящийся на один узел решетки, М - число узлов в кольце вдоль X направления, , - наибольшее собственное значение матрицы В , Т — абсолютная температура в Кельвинах, Д — универсальная газовая постоянная. Формула для локальной степени покрытия г-го ряда имеет вид:

(6)

где - к -я компонента нормированного

Рис. 2. Химический потенциал (Кдж/Моль)

левого/правого собственного вектора, соответствующего наибольшему собственному значению матрицы [} , ¡^ - число занятых узлов в кольце, находящемся в к — ом состоянии. Здесь евклидовое скалярное произведение (yfi :.)равно единице. Полная степень покрытия может быть определена как среднее арифметическое значение локальных степеней покрытия. Вычисления проводились при М = 10.

Алгоритм моделирования методом Монте-Карло.

Моделирование процессов на поверхности методом Монте-Карло имеет ряд преимуществ, в частности метод Монте Карло позволяет вводить в модель по существу любые представления о реальном процессе. Основой метода является так называемое основное уравнение (уравнение Колмогорова, master equation):

ИР

Z Г" ; , 7 , (7)

где t - время, а, и в — конфигурации состояния поверхности, Ра, и Рр - их вероятности, Wa|3 и Wpa - так называемые вероятности перехода в единицу времени, которые определяют скорость, с которой изменяется состояние поверхности. Основное уравнение — наиболее общее уравнение изменения вероятности состояния системы (в данном случае поверхности). Из уравнения (7) видно, что полная вероятность сохраняется

ЯР

. у . -г " , ■(-. ■ о. (8)

Значения вероятностей перехода могут быть вычислены при помощи квантово-химических методов или из экспериментальных данных. Существует много различных алгоритмов реализующих метод Монте-Карло: одни выводятся из Основного Уравнения и приводят к результатам, которые являются статистически идентичными, другие не могут быть получены из Основного Уравнения и от них нужно отказаться.

При моделировании нами использовался упрощенный вариант кинетического метода Монте-Карло.

-тъ dt

Проблема реального времени в алгоритмической формулировке метода Монте-Карло была решена Fichthorn и Weinberg. Их метод называют Кинетическим методом Монте-Карло (KMC). Основное предположение, лежащее в основе метода — интервалы времени между последовательными событиями распределены в соответствии с распределением Пуассона, а переходные вероятности определяются из распределения Больцмана.

Блок алгоритма кинетического метода Монте-Карло состоит в следующем:

1. Инициализация. Генерируется начальная конфигурация а. Устанавливается начальный момент времени t. Выбирается условие, при котором нужно остановить моделирование.

2. Время реакции. Генерируется интервал времени At, в которое не происходит реакция. Время изменяется от t к t + At.

3. Реакция. Изменяется конфигурация, в зависимости от вероятности перехода, то есть происходит реакция.

4. Проверка условий. Если условия остановки выполнены, тогда происходит остановка. Если нет, то повторение со второго шага.

Так как нас интересует только равновесное состояние поверхности мы упрощаем алгоритм, по существу исключая время, но сохраняем выбор вероятностей перехода. Такой упрощенный алгоритм оказывается, за исключением выбора вероятностей перехода, тождественным алгоритму Метрополиса [13].

Подсчет вероятности изменения состояния для пустой ячейки производился по формуле:

(9)

,

1-ехр ---

Ч RT )

где п - энергия взаимодействия ближайших соседей.

Подсчет вероятности изменения состояния для занятой ячейки производился по формуле:

0,4 0,6 0,а 1,0 0,0 0.2 0.4 Рис. 3. Полная степень покрытия

1— ехр

I Пб-Ц.V

(10)

Ч *т }

Здесь под химическим потенциалом и понимаем: для средних рядов террасы — заданное ¡; для сильных адсорбционных центров — и = 1 - Ди; для слабых адсорбционных центров — ¡ = ¡1 + Ди, при условии равенства Д1 и Д2.

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

Д. Д, „

нами принято, что —!_ _ ], _1 — _|. Под локальной

е е

изотермой мы будем понимать зависимость степени покрытия конкретного ряда от химического потенциала. При построении изотерм при каждом значении химического потенциала производилось от 105 до 106 монтекарловских шагов (мы пользуемся стандартным определением шага как количества попыток, равного количеству ячеек на решетке). Усреднение проводилось по второй половине шагов. Размер решетки варьировался от 100х100 до 400х400.

2. Результаты и обсуждение.

Численное моделирование проводилось с использованием МТМ и метода МК. Вычислительный эксперимент показал совпадение результатов, полученных используемыми методами. Вычисления проводились при различных значениях параметров Т (Т е [100; 700]

К), Ь (Ь = 3.....15), Д1 (Д1 е [2; 50] кДж/моль), Д2 (Д2 е

[-50; -2] кДж/моль) и при £ = 10 кДж/моль. Рассмотрим полученные результаты.

На рис. 2 и 3 показаны зависимости локальных степеней покрытия от химического потенциала и от полной степени покрытия соответственно при различных значениях Т, Ь, Д2, которые приведены на рисунках, и при Д1/£ = 2. Сплошная линия показывает предельные зависимости при Д1 = Д2 = 0. Числами показаны степени заполнения соответствующих рядов

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

При фиксированной температуре Т с увеличением Ь степени покрытия для первого, второго и последнего рядов практически совпадают. Степень покрытия (Ь/2 + 1) — го ряда быстрее стремится к предельной зависимости. При возрастании Д1 степень покрытия для первых рядов отдаляется налево от предельной зависимости, а при возрастании Д2 степени покрытия для вторых и последних выделенных рядов приближаются к предельной зависимости справа.

При фиксированной Ь с повышением температуры Т степени покрытия для первых, вторых и последних рядов приближаются к предельной зависимости. Степени покрытия для остальных рядов при Ь > 4 находятся вблизи предельной зависимости.

Рассмотрим рис. 3. В случае Ь = 3 предельные зависимости в области параметров Д1/£ < 1 качественно подобны предельным зависимостям при Ь = 2 [8]. То же самое можно сказать и про область параметров Д1/£ > 3. Рассмотрим область параметров 1 < Д1/£ < 3. Наиболее интересная качественная особенность, появляющаяся в этой области - немонотонный характер зависимости степени покрытия первого выделенного ряда от полной степени покрытия поверхности при 81Ы < 0,5. При Д1/£ = 1, Д2/£ = -1 немонотонный характер имеет и зависимость степени покрытия третьего выделенного ряда от полной степени покрытия 81Ы > 0,5. Отметим, что такая же особенность будет наблюдаться и в зависимости локальной степени покрытия первого ряда от химического потенциала. Вначале происходит заполнение первого ряда, при этом остальные ряды остаются пустыми. Далее начинается заполнение второго ряда, а затем и третьего ряда. Затем происходит уменьшение степени покрытия первого ряда до 81Ы = = 0,5. После этого вновь начинается заполнение первого ряда при постоянном заполнении второго и медленном заполнении третьего рядов. И наконец, первый

ряд полностью заполняется, начинается заполнение второго ряда при заполнении третьего ряда (рис. 3а).

Рассмотрим случай I = 4 в области параметров 1 < Д1/е < 3. Вначале происходит заполнение первого ряда, при этом остальные ряды остаются пустыми. Далее начинается заполнение третьего ряда, а затем второго и четвертого ряда. Затем происходит уменьшение степени покрытия первого ряда до 81Ы = 0,5. После этого вновь начинается заполнение первого ряда при постоянном заполнении третьего, второго и медленном заполнении четвертого рядов. И, наконец, первый ряд полностью заполняется, начинается заполнение третьего и второго рядов при медленном заполнении четвертого ряда (рис. 3б). При 2,5 < Д1/е < 4 зависимость степени покрытия третьего ряда от полной степени покрытия немонотонна.

При возрастании температуры рассматриваемые зависимости монотонны (рис. 2в, 3в).

При I = 5 каких-либо качественных изменений не происходит, немонотонность зависимости степеней покрытия первого и пятого рядов от полной степени покрытия появляется при 1 < Д1/е < 3, пятого последнего ряда при 2,5 < Д1/е < 5.

При I > 7 качественных изменений не происходит. В качестве примера для I = 7 такие зависимости приведены на рис. 3г.

С возрастанием Д1 заполнение первого ряда происходит при малых степенях покрытия. Чем ниже Д2, тем дольше пустует последний выделенный ряд.

3. Выводы

Анализ полученных результатов позволяет сделать следующие выводы:

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

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

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

Библиографический список

1. Zangwill A. Physics at Surfaces, Cambridge University Press, 1988.

2. Rudzinsky W., Everett D.H. Adsorption of gases on heterogeneous surfaces, Academic press, London, 1991.

3. Rudzinsky W., Steele W.A., Zgrablich G. (Eds.). Equilibria and Dynamics of Gas Adsorption on Heterogeneous Solid Surfaces, Elsevier, Amsterdam, 1997.

4. Tamaru K. Dynamic heterogeneous catalysis, Academic Press, New York, 1978.

5. Uebing С. Phys. Rev. B. 49 (1994) 13913.

6. Bulnes F., Ramirez-Pastor A.J., Zgrablich G. Phys. Rev. E. 65 (2002) 31603.

7. Tovbin Yu.K., Rabinovich A.B. Langmuir, 20 (2004) 6041.

8. Myshlyavtsev A.V., Myshlyavtseva M.D. Modeling of adsorption and phase diagrams for stepped surfaces: Transfer matrix approach. // Applied surface science, 2007. — V. 253, № 13. - P. 5591 -5595.

9. Шулепов Ю.В., Аксененко Е.В. Решеточный газ. — Киев: Наукова думка. — 1981. — 267 с.

10. Бэкстер P. Точно решаемые модели в статистической механике - М.: Мир, 1985 - 486 с.

11. Мышлявцев А.В., Мышлявцева М.Д. Вычислительные аспекты метода трансфер-матрицы. — Кызыл: ТувИКОПР СО РАН. — 2000. — 101 с.

12. Биндер К. Методы Монте Карло в статистической физике - М.: Мир, 1982 — 398 c.

13. Metropolis N., Rosenbluth A.W., Rosenbluth M.N., Teller A.W., Teller E. Equation of state calculations by fast computing machines // J. Chem. Phys. — 1953. — V.21. - P. 1087.

ЕВСЕЕВА Светлана Ивановна, ассистент кафедры высшей математики.

МЫШЛЯВЦЕВ Александр Владимирович, доктор химических наук, проректор по учебной работе Омского государственного технического университета, профессор кафедры «Химическая технология переработки углеводородов».

МЫШЛЯВЦЕВА Марта Доржукаевна, кандидат физико-математических наук, доцент кафедры высшей математики Омского государственного технического университета.

Дата поступления 12.12.2007 г.

©Евсеева С.И., Мышлявцев А.В., Мышлявцева М.Д.