Моделирование неидеального адсорбционного слоя для каталитических реакций с несколькими интермедиатами на примере механизма Ленгмюра - Хиншельвуда Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 19, № 5, 2014

Моделирование неидеального адсорбционного слоя для каталитических реакций с несколькими интермедиатами на примере механизма Ленгмюра — Хиншельвуда

В. И. Быков1, Е.А. Мамаш2, О. Э. Лешаков3

Для модели неидеального адсорбционного слоя с двумя типами адсорбированных веществ с применением метода трансфер-матрицы построены изотермы адсорбции при различных значениях энергий латеральных взаимодействий. Для неидеальной модели механизма Ленгмюра — Хиншельвуда проведён параметрический анализ.

Ключевые слова: каталитические реакции, неидеальные модели, трансфер-матрица, параметрический анализ.

К настоящему времени накоплен значительный опыт по построению и анализу нелинейных и нестационарных кинетических моделей [1]. Традиционно моделирование сложных реакций, в том числе процессов горения, проводится в рамках закона действующих масс (для гетерогенных реакций — действующих поверхностей), т. е. в рамках идеальных моделей химической кинетики [1-3]. Однако для описания реальных систем, поведение которых существенно отличается от "идеального", возникает необходимость использования адекватных моделей, учитывающих влияние осложняющих факторов, прежде всего химической неидеальности (кинетика Марселена Де-Донде [4]) и неидеальности адсорбционного слоя катализатора. Кинетика Марселена Де-Донде и её обобщение на случай неизотермических реакций были проанализированы авторами [5, 6], где показана важность учёта термодинамических ограничений и корректных поправок на неидеальность при моделировании такого рода систем. Отдельно возникает вопрос о возможности построения моделей, сочетающих химическую неидеальность с неидеальностью поверхности катализатора, однако он не является задачей настоящего исследования.

В данной работе рассматривается кинетическая модель, построенная в рамках неидеальной модели адсорбционного слоя на примере классического адсорбционного механизма Ленгмюра — Хиншельвуда. В [7, 8] уже исследовалась возможность получения математических моделей некоторых базовых каталитических процессов на основе знания энергий латеральных взаимодействий адсорбированных веществ. Эти параметры, характеризующие микроуровень, могут быть определены на основе экспериментальных методов изучения каталитических систем, например, с использованием термореакционных и термодесорбционных спектров. С другой стороны, изотермы адсорбции можно получить и расчётным путём методом трансфер-матрицы (МТМ) [9].

1 Российский химико-технологический университет им. Д.И. Менделеева, Москва

2Институт вычислительных технологий СО РАН, Новосибирск, Россия

3 Тувинский институт комплексного освоения природных ресурсов СО РАН, Кызыл, Россия Контактный e-mail: elenamamash@gmail.com

В представленном исследовании на примере механизма Ленгмюра — Хиншельвуда с двумя промежуточными веществами приводится математическая модель, связывающая макроскопические параметры (степени покрытия поверхности адсорбированными веществами) с микроскопическими (энергии латеральных взаимодействий адсорбированных частиц). Обсуждаются результаты аналитического и численного изучения этой модели на основе использования МТМ. Ранее результаты применения МТМ к анализу неидеальной модели механизма Ленгмюра — Хиншельвуда обсуждались в работе [10], однако там были приведены результаты численного решения модели без сопоставления с аналитическим описанием.

Модель совместной адсорбции нескольких частиц имеет не только теоретическое, но и прикладное значение. Самым простым и практически важным случаем является наличие двух типов адсорбированных частиц, когда каждый узел решетки может находиться в трёх состояниях: узел пуст 0, узел занят частицей типа А, узел занят частицей типа В. Для наглядности все возможные комбинации состояний двух соседних узлов можно представить в следующем виде:

00 0А 0В

А0 АА АВ

В0 ВА ВВ

поэтому трансфер-матрица однородной модели решеточного газа (МРГ) на линейной цепочке с взаимодействием только ближайших соседей имеет размерность 3 х 3 и её элемент можно записать как

Т (иг, иг+1) = ехр [^ К + <0 + ^ (п? + -

-£АЛПАПА+! - £??П?П?+1 - £л? (<4?1 + п?И^)] , (2)

где £лл, £??, £л? = £?л — энергии латеральных взаимодействий ближайших соседних

А?

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

Согласно методу трансфер-матрицы степени покрытия поверхности веществами А и В могут быть вычислены по формулам

вл = в? = VI, (3)

где вл, в? — степени покрытия поверхности катализатора веществами А и В, —

координаты нормированного собственного вектора трансфер-матрицы (2), соответствующего наибольшему собственному значению. Из (2) следует, что есть функции, зависящие от значений потенциалов ^л и . Пример изотермы ва = в л (^л,^?) приведён на рис. 1, из которого в отличие от изотермы адсорбции для идеального слоя видно наличие характерной ступеньки. На рис. 2 приведены кривые, полученные в сечении изотермы, изображённой на рис. 1, при фиксированных значениях химических потенциалов ^л и .

Хорошим приближением при описании адсорбированных слоёв на поверхности твёрдых тел могут служить двумерные решёточные модели. Термодинамический гамильтониан такой МРГ имеет вид [9-11]

Ней = £лл ^ ПлПл + £?? ^ п? п? + £л? ^ плп? - ^л ^ пл - ^ п?. (4)

(гага) (гага) (гага) г г

Рис. 1. Изотерма вл {цл, Цв ) для линейной решёточной модели, построенная при елл = —1,

£вв = 2 елв = —2

Рис. 2. Кривые, полученные в результате сечения изотермы, изображённой на рис. 1, плоскостями; ца = —20 (1), 0 (2), 10 (3), 20 (4)

Здесь символ (пп) означает, что суммирование распространяется на все пары ближайших соседних узлов квадратной решетки. Фактически мы имеем решетку на бесконечном цилиндре с периметром М, где М — число узлов решетки. Тогда элемент трансфер-матрицы будет иметь вид

Тц = ехр -

м

м

^л ^ (пА + — -лл (плкпАк+1 + пЛк ^+1) +

к=1

к=1

м

м

+Мв (п»В + п?к) — еВВ ^ (пВкпВк+1 + ПИпВк+0

к=1

м

к=1

£ЛВ/ у (пАк пВк+1 + пВк пАк+1 + пАкпВк+1 + пВкпАк+1)

м

/ у \'Нк' к=1

,Л„Л

м

м

—£лл ^ (пЛк4) — *вв£ (пВкп,Вк) — -лв ^ (пЛкп% + пВкпЛк)

к=1

к=1

к=1

Чтобы сократить размерность полученной трансфер-матрицы, был использован приём, основанный на трансляционной инвариантности системы, позволяющий более чем в М раз уменьшить её размеры [9]. Выборочные расчёты показали, что при увеличении параметра М (М > 6) меняются только количественные характеристики системы, тогда как качественное поведение остается тем же. Изотермы адсорбции для МРГ на квадратной решетке рассчитываются по формулам

1 зм м

«а = М £ V? £ пАк, »В

i=1 к=1

1 3м м

= М ^ ^ ^ пВк i=1 к=1

где щ — ¿-я компонента нормированного собственного вектора трансфер-матрицы, отвечающего наибольшему собственному значению. График изотермы 9 а = 0а (^а,^б) приведён на рис. 3. Заметим, что по сравнению с изотермой для МРГ на линейной цепочке (см. рис. 1) ступенька здесь более выражена. Это в большей степени соответствует реальным изотермам. Как и выше, на рис. 4 представлены сечения полученной поверхности при фиксированных значениях химических потенциалов. Форма изотерм полностью определяется значениями параметров £аа, £бб , £аб. Так, при увеличении энергии £аб площадь пологой ступеньки изотермы значительно уменьшается, что приводит к снижению вероятности возникновения критических явлений.

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

А2 + 2Х о 2АХ, В + X о ВХ,

А^ + ВХ ^ АВ + 2Х. (7)

Соответствующие (7) балансовые нестационарные уравнения имеют вид

^ = 2кЛРл (1 - 9а - 9б)2 - 2кЛ9А - кг9а9б,

аъ

Мг,

кБрб (1 - 9а - 9б) - кВ9б - кг9а9б- (8)

аъ

Здесь Ъ — время, к^А, к^, кБ, кБ — константы скоростей адсорбции и десорбции для веществ А и В соответственно, кг — константа скорости реакции, Рл, Рб — парциальные

Рис. 3. Изотерма в а {ца, Цб ) для квадратной решётки при в ал = -1, вбб = 2, еаб = -2

Рис. 4. Кривые, полученные в результате сечения изотермы, изображённой на рис. 3, плоскостями; ца = -20 (1), 0 (2), 10 (3), 20 (4)

давления веществ в газовой фазе, которые считаются постоянными. В рамках теории переходного состояния и МРГ без учёта латеральных взаимодействий активированных комплексов эти константы имеют вид [12]

кЛ = кЛоРоо/ (1 — »л — »в)? , кЛ = кАо ехр (2дл/ДТ) Роо/»Л,

кв = кВ ко. = ка,о,

к| = кВо ехр (дв/ДТ) (1 — «Л — »в) /»в,

кг = кг,о ехр ((иа + Дв) /ДТ) Роо/»л»в, (9)

где &Ло, кЛо, кВо, кВо, кг,о — константы элементарных процессов при малых степенях покрытия; ил и дв — химические потенциалы веществ А и В соответственно; Роо — вероятность того, что пара ближайших соседних мест, зависящая от потенциалов Да и Дв, пуста; Д — универсальная газовая постоянная; Т — абсолютная температура. Подставляя выражения (9) в систему (8), получаем:

Л»Л 2кЛо Роо Ра — 2кЛоРоо ехр (2дл/ДТ) — кг,оРоо ехр ((да + Дв) /ДТ),

л»

_»£ = (1 — »Л — »В) (кВоРв — кВо ехр (дв/ДТ)) — кг,оРоо ехр ((да + Дв) /ДТ). (10)

Химические потенциалы да и Дв есть некоторые функции от »а, »в и энергий -ал, -вв, -лв. Поэтому правая часть системы (10) зависит от параметров микроуровня, тогда как в левой части находятся производные от степеней покрытия — макроскопических параметров. Численно решать систему (10) относительно переменных »л и »в может оказаться достаточно сложно, поскольку возникает задача аппроксимации значений Да (»л, »в), Дв (»л, »в) по заранее построенной изотерме на каждом шаге вычислений с заданной точностью. Альтернативный приём здесь состоит в переходе от совокупности переменных »а, »в к совокупности переменных да, Дв. Вычислив, согласно МТМ, решения да (¿), Дв (¿), можно получить кинетические зависимости »л (¿), »в (¿). Таким образом, в процессе реализации вычислительного алгоритма все операции проводятся относительно заданных параметров модели и переменных да, Д в, а затем на основании МТМ находятся »л (дл, ДВ) и »В (дл, ДВ) и строятся искомые зависимости.

Кроме получения кинетических зависимостей, параметрический анализ динамической системы включает в себя определение стационарных состояний системы, анализ их зависимости от параметров, построение линий локальных бифуркаций: кривой кратности Ьд (границы области множественности) и кривой нейтральности (границы устойчивости), построение параметрических и фазовых портретов.

Уравнение стационарности системы (10) можно получить, приравняв её правые части к нулю. Например, относительно да оно имеет вид

(1 — »л (да, дт 1п v (йл)) — »в (дл, дт 1п v (йл))) х . 1о — О^лЙ — кг,оРоо

где

X (к^Во — к£оР (Дл)) — кг,оРоо (дл, ДТ 1п V (Дл)) V (Дл) = 0, (11)

2 (кЛо — кЛое2М

V (Да) = 1",ок -^, Да = Дл/ДТ. (12)

кг

Учитывая, что функция V (уа) = V {ул,кАА0,кАА0,кАА0), удобно выписать обратные параметрические зависимости для параметров кВ0 и кв0 (здесь и далее для краткости опустим, что функции Роо, 9а, 9б и их производные вычисляются в стационарном состоянии, т.е. в точке (уа, ЯТ 1пV (уа))):

кБ

ка , 0

2ил

Ро

оо

1 - 9а - 9

+

1Б

кБ

0

кг0в^А

13)

кБ

0

кг0в^А

кБ

ка ,0

2иА

Р0

00

1 - 9а - 9

'Б

14)

Здесь ил = кА0 - кАре2»А.

Вычислив производную по у а правой части выражения (13) для параметра кБ0 и приравняв её к нулю, из полученного уравнения выразим кБ0. В результате получим аналитическое выражение для кривой кратности Ьд {кБ0, кБ0

кБ

, 0

дР

ЯТ-Р^ил - 2Р00кА0в2^А дул ,

- , -9а -9б , -

кг^А Р00илЯТ\ — + т— кгпе^А

+

дул дул

(1 - 9а - 9Б) (кА0 + кА^А) (1 - 9А - 9Б)2 (кА^ + А^а)

к?0 = 2ил

Р0

00

+

кБ

1 - 9а - 9Б кт^А

Линия нейтральности Ьа для системы (10) имеет вид

кБ = -Р кБ = кБ (кБ \

ка,0 = ка,0 = ка,0 \к4,0]

где

Р

2илР00А

1 - 9 а - 9

-дР00тт 4 ,„А О „2дд ттА

Ь1

+ иА ^ - иБ

8,1 иА - ЯТкА0Р00^А - иА

-9а

89а

Б

ду

Б

дул

д9

Б

дув

2 8Р°0 ил - ив

ду

Б

д9

Б

Я

дул' 89а

-А(. к^ТА - V Ц - ЯтР (уА)(1 - 9А - 9Б) ду А ■

А

двл двл -9в -9в

и В = и,

д"А д"В

дР00

_ Р00 ду Б + ЯТ

иА = и (— + —

и ил ду А + ЯТ

+ 100 , и, = кг^А V (уа) .

Численный анализ показал, что как в рамках линейной МРГ, так и для МРГ на квадратной решетке при построении линии нейтральности Ьа параметры можно подобрать таким образом, что кривая Ьа в их плоскости будет иметь петлю. Это указывает на

возможность возникновения в системе автоколебаний. Пример параметрического портрета для линейной МРГ представлен на рис. 5. Здесь области 1, 2 соответствуют одному стационарному состоянию типа "узел". Для области 1 оно устойчиво, для области 2 — неустойчиво. Области 3-6 отвечают трём стационарным состояниям, два из которых устойчивы и являются "узлом", или "фокусом", а их области притяжения делятся сепаратрисами третьего неустойчивого стационарного состояния типа "седла". Подробный анализ механизма (7), содержащего два интермедиата, приведён в работе [13].

Если процесс адсорбции происходит в присутствии трёх, четырёх и более интер-медиатов, то основная схема построения соответствующих изотерм при помощи МТМ принципиально не будет отличаться от схемы, приведённой выше для двух типов частиц. Работа с двумерными МРГ может быть связана с рядом трудностей технического характера, возникающих из-за резко возрастающей размерности трансфер-матрицы с ростом количества частиц. Однако для трёх и четырёх частиц эта проблема вполне решаема, а задачи, в которых рассматриваются каталитические механизмы с числом интермедиатов более пяти, встречаются нечасто и могут быть решены на основе подбора оптимальных алгоритмов МТМ и привлечения дополнительных вычислительных ресурсов. Результаты численного и качественного анализа динамической системы для химических потенциалов д на основе изотермы адсорбции могут быть представлены в терминах степеней покрытия в. В многочисленных экспериментальных работах последних десятилетий было показано, что адсорбционный слой на поверхности катализатора может быть существенно неидеальным [14]. В этом случае обычно используемая модель поверхности, не учитывающая латеральные взаимодействия адсорбированных частиц, недостаточно адекватна для ряда реальных систем. Учёт энергий латеральных взаимодействий адсорбированных частиц не только может приводить к возникновению различных критических эффектов [15, 16], но и позволяет проанализировать влияние микроуровня на ход реакции на макроуровне [7-10]. В настоящей работе тема неидеальности рассмотрена в значительно более общей постановке — для каталитических систем с несколькими интермедиатами как в рамках линейной решеточной модели, так и для

Рис. 5. Параметрический портрет системы (12) для одномерной МРГ при елл = -2, £вв = 10, елв = -2, 0 = 15, &Л0 = 2, кг,0 = 10, ЕТ = 3. Сплошная линия — кривая нейтральности, штриховая — кривая кратности

классической двумерной МРГ на квадратной решетке. Установить связь между макроскопическими и микроскопическими параметрами позволяет метод трансфер-матрицы решёточных систем [9], занимающий промежуточное положение между детерминированным подходом и имитационным моделированием (метод Монте-Карло [17]) и позволяющий описывать гетерогенно-каталитические системы с учётом упорядоченного характера расположения адсорбированных частиц на поверхности катализатора.

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

[1] Быков В.И., Цыбенова С.Б. Нелинейные модели химической кинетики. М.: КРА-САНД, 2011. 400 с.

[2] Димитров В.И. Простая кинетика. Новосибирск: Наука, 1982. 382 с.

[3] Азатян В.В., Андрианова З.С., Иванова А.Н. Моделирование ингибирования распространения пламени в водородо-воздушной среде // Кинетика и катализ. 2010. Т. 51, № 4. С. 483-491.

[4] Feinberg M. On chemical kinetics of a certain class // Arch. Rat. Mech. Anal. 1972. Vol. 46, No. 1. P. 1-41.

[5] Bykov V.I., Gorban A.N., Yablonskii G.S. Description of nonisothermal reactions in terms of Marcelin-de Donder kinetics and its generalizations // React. Kinet. Catal. Lett. 1982. Vol. 20, No. 34. P. 261-256.

[6] Быков В.И., Иванова А.Н. Химическая неидеальность как причина критических явлений // Кинетика и катализ. 1986. Т. 27, вып. 1. С. 73-79.

[7] Myshlyavtsev A.V., Samdanohap R.T. Criterion of multiplicity of steady states and the explicit expression for bifurcation set in the simplest model of a catalytic reaction with an arbitrary isotherm // React. Kinet. Catal. Lett. 1994. Vol. 53, No. 2. P. 269-276.

[8] Быков В.И., Мамаш Е.А. Нестационарные кинетические модели некоторых типовых каталитических механизмов // Кинетика и катализ. 2012. Т. 53, № 5. С. 556-562.

[9] Мышлявцев А.В., Мышлявцева М.Д. Латеральные взаимодействия в адсорбционном слое и критические явления в реакции, протекающей по помеханизму Ленгмюра — Хиншельвуда // Там же. 2007. Т. 48, № 4. C. 576-585.

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

[11] Товбин Ю.К., Черкасов А.В. Влияние неидеальности адсорбционной системы на число стационарных решений простейшего каталитического процесса // Теор. и эксперим. химия. 1984. Т. 20, № 4. С. 507-512.

[12] Zhdanov Vol. P. Elementary Physicochemical Processes on Solid Surfaces. N.-Y.: Plenum Press, 1991. 324 p.

[13] Быков В.И., Мамаш Е.А. Автоколебания в модели адсорбционного каталитического механизма // Вестник КрасГУ. 2004. Вып. 3. C. 4-9.

[14] Imbihl R., Ertl G. Oscillatory kinetics in heterogeneous catalysis // Chem. Rev. 1995. Vol. 95, No. 3. P. 697-733.

[15] Hilderbrand M. Self-organized nanostructures in surface chemical reaction mechanisms and mesoscopic modeling // Chaos. 2002. Vol. 12, No. 1. P. 144-156.

[16] Oertzen A., Rotermund H.H., Mikhailov A.S., Ertl G. Standing wave patterns in the CO oxidation on a Pt(110) surface: Experiments and modeling //J. Phys. Chem. B. 2000. Vol. 104. P. 3155-3178.

[17] Matveev A.V., Latkin E.I., Elokhin V.I., Gorodetskii V.V. Turbulent and stripes wave patterns caused by limited CO diffusion during CO oxidation over Pd(110) surface: Kinetic Monte-Carlo studies // Chem. Eng. J. 2005. Vol. 107. P. 181-189.

Поступила в 'редакцию 4 февраля 2014 г., с доработки — 23 июня 2014 г.

Modelling of a non-ideal adsorbed overlayer for catalytic reactions with several intermediations on the example of a Langmuir — Hinshelwood mechanism

Bykov Valerii I.1, Mamash Elena A.2, Leshakov Oleg E.3

Traditionally, modeling of complex reactions is implemented in accordance with the mass action law (active surfaces for heterogeneous reactions), i.e. within frameworks of ideal models of chemical kinetics. However, the description of real systems, whose behavior is significantly different from the "ideal", requires implementation of more adequate models. These models need to take into account the influence of complicating factors, first of all, the chemical non-ideality and non-ideality of the adsorbed overlayer of the catalyst. Purpose. In this paper we consider a kinetic model formulated under assumptions of a non-ideal model of the adsorbed overlayer based on the classic example of the Langmuir — Hinshelwood adsorption mechanism with two intermediates.

Methodology. The mathematical model of the mechanism linking the macroscopic parameters (surface coverage) with microscopic (lateral interaction energies of adsorbed particles) is presented. We use both 1D and 2D lattice gas model on a square lattice as the model of the adsorbed overlayer. The adsorption isotherms are calculated using the transfer matrix method well known in statistical physics.

Findings. It is shown that the profile of the isotherms is completely determined by the values of energies of lateral interactions. The qualitative analysis of the model is carried out; exact analytical expressions for the bifurcation curves are obtained. Originality. The numerical analysis of the model, conducted through the use of the transfer matrix method showed the possibility of self-oscillations both within linear lattice gas model and the square lattice gas model.

Keywords: catalytic reactions, non-ideal models, transfer-matrix, parametric analysis.

Received 4 February 2014 Received in revised form 23 June 2014

1 Emanuel Institute of Biochemical Physics RAS, Moscow, Russia

2Institute of Computational Technologies SB RAS, 630090, Novosibirsk, Russia

3 Tuvinian Institute for Exploration of Natural Resources SB RAS, 667000, Kyzyl, Tuva Republic, Russia

Corresponding author: e-mail: elenamamash@gmail.com