Нелинейная динамика молекулярных процессов в гетерогенных системах Песков Николай Владимирович

Содержание к диссертации

Введение

2 Модели колебаний в реакции окисления СО на поверхности палладия 39

2.1 Механизмы колебаний в реакции окисления СО на палладии 39

2.1.1 Окисление-восстановление поверхности 39

2.1.2 Растворение и сегрегация кислорода 42

2.1.3 Модель регулярных колебаний в Pd-цеолитном катализаторе . 44

2.2 Автоколебания в точечной модели реакции 48

2.2.1 Быстрые и медленные переменные 49

2.2.2 "Квазистационарное"приближение 53

2.2.3 Циклы-утки в точечной модели 55

2.3 Основные результаты главы 59

3 Модели колебаний в реакции окисления СО в зерне палладиевого цеолитного катализатора 60

3.1 Распределенные модели гетерогенного катализа 60

3.2 Строение цеолитного катализатора 64

3.3 Двухслойная модель реакции в зерне цеолита 67

3.4 Динамика двухслойной модели 70

3.4.1 Динамика реакции в одном слое зерна 70

3.4.2 Фазовая диаграмма двухслойной модели 74

3.4.3 Бифуркационные диаграммы колебаний 76

3.5 Хаотические колебания в двухслойной модели 84

3.5.1 Переход к хаосу через бифуркации удвоения периода 87

3.5.2 Перемежаемость 94

3.6 Непрерывная модель реакции 98

3.6.1 Фаза и частота сложных колебаний 101

3.6.2 Динамика непрерывной модели при Ds = 0 103

3.6.3 Динамика непрерывной модели при Ds > 0 115

3.7 Основные результаты главы 124

4 Модели колебаний в реакции окисления СО на наночастицах палладия 126

4.1 Моделирование влияния размера частицы катализатора на скорость реакции 126

4.1.1 Точечная модель реакции с учетом окисления объема катализатора 126

4.1.2 Учет влияния размера частицы катализатора 132

4.2 Стохастическая модель реакции на частице 137

4.2.1 Формулировка стохастической модели реакции 137

4.2.2 Уравнения для средних концентраций 140

4.2.3 Алгоритм Монте-Карло . 143

4.3 Динамика стохастической модели 144

4.3.1 Численные результаты 144

4.3.2 Фазовые траектории стохастической модели 150

4.4 Модель реакции на цепочке частиц 153

4.4.1 Модель цеолитного катализатора 153

4.4.2 Модель сорбции и диффузии газа в цеолите 156

4.4.3 Модель поверхностной реакции 160

4.4.4 Колебания скорости реакции в модели 163

4.5 Основные результаты главы 172

5 Анализ экспериментальных данных 174

5.1 Некоторые алгоритмы анализа временных серий 174

5.1.1 Подавление шума 175

5.1.2 Реконструкция фазовых переменных 176

5.1.3 Определение размерности вложения 182

5.1.4 Вычисление показателей Ляпунова 185

5.2 Комплекс программ CDP 191

5.2.1 Модуль ввода данных 191

5.2.2 Модуль определения времени задержки 193

5.2.3 Модуль изображения аттрактора 195

5.2.4 Модуль определения размерности вложения 196

5.2.5 Модуль оценки старшего показателя Ляпунова 198

5.2.6 Модуль отображений 200

5.3 Анализ временных серий скорости реакции окисления СО в Pd цеолитном катализаторе 201

5.3.1 Экспериментальные данные 201

5.3.2 Оценка величины временной задержки 203

5.3.3 Оценка размерности вложения 207

5.3.4 Оценка старшего показателя Ляпунова 211

5.4 Основные результаты главы 213

6 Моделирование нелинейной динамики на микроуровне. Стоха стическая модель МЛЭ 215

6.1 Введение 215

6.1.1 Общая характеристика метода МЛЭ 215

6.1.2 Механизм роста соединений типа III-V в МЛЭ 220

6.2 Стохастическая модель МЛЭ 222

6.2.1 Модель кристаллической решетки соединения III-V 223

6.2.2 Модели элементарных актов 226

6.2.3 Алгоритм Монте-Карло 232

6.2.4 Комплекс компьютерных программ 235

6.3 Результаты моделирования 237

6.3.1 Послойный рост GaAs(OOl) . 237

6.3.2 Образование объемных вакансий 242

6.3.3 Роль миграции Ga в МЛЭ GaAs 245

6.3.4 Рост слоев AlGaAs на ступенчатых гранях 251

6.4 Основные результаты главы 258

7 Выводы 260

Литература 266

• Модель регулярных колебаний в Pd-цеолитном катализаторе
• Переход к хаосу через бифуркации удвоения периода
• Формулировка стохастической модели реакции
• Оценка величины временной задержки

Введение к работе

1.1 Нелинейная динамика в гетерогенном катализе

Теория динамических систем развивалась с конца XIX века в трудах таких выдающихся ученых как А. Пуанкаре, А. М. Ляпунов, Л. И. Мандельштам, А. А. Андронов, Н. Н. Боголюбов, М. Морс, X. Уитни, В. И. Арнольд и др. По-видимому, наиболее простым определением динамической системы является следующее: система называется динамической, если задание начальных условий полностью определяет ее поведение в последующие моменты времени [1]. Во многих случаях динамическая система задается системой обыкновенных дифференциальных уравнений (ОДУ)

ix(*) = F(X(t),p), (1.1.1)

где X,F К.", р - скалярный или векторный параметр, t > О - время; с начальными условиями Х(0) = Х, которые выделяют конкретную фазовую траекторию. Особый интерес представляют нелинейные динамические системы и, в частности, класс нелинейных диссипативных систем, для которых divF < 0. Такие системы используются в качестве моделей различных явлений и процессов в открытых неравновесных системах в физике, химии, биологии и т.д.

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

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

В настоящее время методы нелинейной динамики широко используются в самых разных областях фундаментальной и прикладной науки. Они позволяют рассматривать с единой точки зрения самые разные явления в природе, в технике, в медицине, в экономике и в других областях. Общим для всех этих явления является то, что они происходят в открытых неравновесных системах и для них характерны элементы самоорганизации. Одной из областей науки, где в последнее время активно и успешно применяются методы нелинейной динамики, является гетерогенный катализ [2, 3]. Особенно ярко нелинейные эффекты проявляются в модельных реакциях окисления СО и редукции N0 на гранях с малыми индексами Миллера монокристаллов металлов платиновой группы, которые протекают в сверхвысоком вакууме [4]. Современные экспериментальные методы позволяют in situ наблюдать пространственно-временные структуры в адсорбционном слое на поверхности монокристалла.

Нелинейные явления в гетерогенных каталитических реакциях наблюдаются на разных пространственных масштабах. На атомарном микроскопическом масштабе (10~⁹ м, нм) нелинейные взаимодействия между атомами в адсорбционном слое (латеральные взаимодействия) могут приводить к образованию сверхструктур различного типа, т.е. регулярных решеток на поверхности катализатора, в узлах которого располагаются адатомы одного типа; к формированию стационарных и нестационарных островковых и доменных структур в адсорбционном слое. В процессе роста кристаллов в неравновесных условиях нелинейные взаимодействия выражаются в различных режимах роста, - послойном, трехмерном, фрактальном. На мезоскопическом масштабе (Ю^-6 м, мкм) в адсорбционном слое наблюдаются пространственно-временные диссипативные структуры, бегущие и стоячие волны

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

Активное изучение нелинейной динамики гетерогенных каталитических реакций началось в начале семидесятых годов с открытием автоколебаний в реакции окисления СО на нанесенном платиновом катализаторе [5] и в реакции окисления Н₂ на поликристаллическом никеле [6]. За тридцать лет после этого открытия нелинейные явления были обнаружены во многих каталитических системах. Особенно широко изучались реакции экологического катализа, в частности, реакции окисления СО, СО + 0₂ -> С0₂, и редукции NO, NO + СО -» N₂ -f С0₂, на поверхности металлов платиновой группы. Эти реакции играют важную роль в процессе обезвреживания выхлопных газов двигателей внутреннего сгорания. В ходе экспериментального изучения этих реакций на поверхности катализатора были обнаружены различные стационарные и нестационарные диссипативные структуры: стоячие и бегущие волны скорости реакции, спиральные волны, мишеневидные структуры и другие пространственно-временные структуры [7, 8, 4, 9]. Экспериментальное открытие нелинейных явлений в гетерогенном катализе стимулировало многочисленные работы по математическому моделированию нелинейной динамики этих явлений [10, 11, 12].

Нелинейная динамика каталитических реакций является основанием для разработки теории и практики катализа [2]. Изучение нелинейной динамики гетерогенных реакций на всех пространственных масштабах важно для понимания механизмов реакций, для разработки новых промышленных технологий и для определения оптимальных условий протекания реакций. Поведение реакций в динамическом режиме в гораздо большей степени зависит от межмолекулярных взаимодействий, чем в равновесном или стационарном режимах. Поэтому изучение динамики реакций дает больше информации о физико-химических процессах на атомарном уровне, чем изучение тех же реакций в стационарном режиме. Огромную роль в разработке современной теории гетерогенного катализа играет математическое моделирование и вычислительный эксперимент, концепция которого была разработана

академиком А. А. Самарским [13].

1.2 Математическое моделирование реакций

на различных пространственных масштабах

Гетерогенные каталитические реакции и процессы обладают многоуровневой упорядоченной структурой. Можно выделить три главных уровня: микроуровнь (0.1 -100 нм), мезоуровнь (1 - 100 мкм) и макроуровнь (> 100 мкм) [14]. Уровень должен выделяться так, чтобы он мог быть изучен экспериментально. Наблюдаемые закономерности протекания процесса на данном уровне не должны зависеть от его масштаба. Влияние масштаба должно определяться граничными и начальными условиями. Основные взаимодействия происходят между соседними уровнями, однако не исключается также и влияние более высоких уровней. Важное значение имеют координация, согласование и интеграция протекающих процессов на разных уровнях. Особое значение при решении проблем катализа имеет атомно-молекулярный уровень. Процессы, протекающие на атомно-молекулярном уровне, определяют активность катализаторов и селективность реакций [15].

Математические модели процессов на разных уровнях каталитических реакций обладают спецификой, характерной для данного уровня. Наиболее полное и детальное описание процессов достигается в моделях микроуровня. В принципе, модели мезо- и макроуровня должны получаться из моделей микроуровня с помощью процедур типа процедуры усреднения по статистическому ансамблю при некоторых дополнительных предположениях. Однако, не смотря на значительные усилия в этом направлении, вопрос о выводе моделей мезо- и макроуровня из модели микроуровня до сих пор остается открытым. Поэтому многие модели гетерогенного катализа носят феноменологический характер [16, 17]. Ниже кратко перечисляются характерные особенности моделей разных уровней.

Макроскопический уровень. Модели макроскопического уровня гетерогенных каталитических реакций описывают эволюцию усредненных по поверхности катализатора характеристик реакции. Чаще всего такими характеристиками являются

средние покрытия (поверхностные концентрации) поверхности катализатора молекулами реагентов и физико-химические характеристики поверхности, которые могут изменяться в ходе реакции. В моделях макроуровня нет зависимости переменных от пространства, поэтому их называют сосредоточенными или точечными моделями.

С поверхностью катализатора можно связать регулярную решетку адсорбционных центров. В этом случае поверхностная концентрация 0,- молекул 4-го типа представляет собой среднюю по решетке вероятность нахождения молекулы г'-го сорта в узле решетки. Точечные модели, в которые входят одноузельные вероятности, называют моделями среднего поля. Для более адекватного описания процесса в некоторых случаях в качестве переменных модели рассматривают также средние вероятности заполнения пары соседних узлов. Для вывода уравнений для парных вероятностей используют различные приближения, в частности, квазихимическое приближение.

В точечной модели в приближении среднего поля уравнение для эволюции переменной в і имеет вид:

где Rf_tk - скорость элементарной стадии реакции с участием реагента г-го типа, которая изменяет величину #,-. Скорости элементарных стадий могут зависеть от переменных модели #!, ...,0jv и от внешних параметров, например, температуры реактора, парциальных давлений реагентов и т.п. Обычно внешние параметры считаются постоянными в ходе реакции, поэтому точечная модель реакции представляет собой автономную систему обыкновенных дифференциальных уравнений. Реакция может протекать неограниченно долго только при условии постоянного притока реагентов в реактор, поэтому реакционная система является открытой неравновесной системой, а ее математическая модель является диссипативной системой уравнений. Уравнения системы, как правило, нелинейны и для изучения модели широко применяется качественный анализ систем ОДУ [18]. Скорости элементарных стадий могут различаться на несколько порядков величины, поэтому уравнения модели являются, как поавило. жесткими и лля их численного лептения слелует ппименять

специальные методы.

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

В настоящей работе представлены макроскопические модели реакции окисления моноксида углерода СО на поверхности монокристалла палладия при атмосферном давлении (Глава 2).

Мезоскопический уровень. К моделям гетерогенного катализа мезоскопического уровня обычно относят модели, представленные системами уравнений в частных производных типа "реакция-диффузия". Переменные модели обычно те же, что и в макроскопических моделях, однако теперь они зависят не только от времени, но и от пространственных координат, поэтому такие модели называют распределенными моделями. С помощью моделей мезоскопического уровня описывают нелинейные явления типа автоволн разных типов, автосолитонов, стационарных и нестационарных диссипативных структур [1, 10, 19, 20, 21].

Модели гетерогенных реакций мезоскопического уровня обычно имеют вид системы дифференциальных уравнений:

— = D(u)u + R(u),

для вектор-функции u(x,i) в некоторой области пространства с заданными начальными и граничными условиями. Динамическими переменными модели, - компоненты вектор функции и, - в общем случае являются концентрации реагентов на поверхности катализатора и в газовой фазе, концентрации промежуточных веществ, а также определенные физико-химические характеристики поверхности. Оператор

диффузии D(u) описывает диффузию решеточного газа на поверхности катализатора и диффузию реагентов в газовой фазе. Члены уравнений системы, обозначенные символом R, описывают химические реакции на поверхности катализатора и обмен веществом между газовой фазой и поверхностью.

При выводе уравнений типа "реакция-диффузия" обычно рассматриваются соотношения баланса массы и энергии в физически бесконечно малом объеме пространства [22]. (По этой причине модель называется мезоскопической.) В качестве модели реакции R часто берут макроскопическую модель. Т.е. пространственно однородное решение u(x, t) = u(f), на котором оператор диффузии равен нулю, должно удовлетворят уравнениям макроскопической модели реакции.

Изучение моделей типа "реакция-диффузия" имеет достаточно большую историю. В качестве первых образцов можно упомянуть работы Колмогорова А. Н., Петровского И. Г. и Пискунова Н. С. [23], Тюринга А. М. [24]. Не смотря на значительные успехи в изучении отдельных моделей, связанные с анализом автомодельных решений, спектральными методами и др., достаточно развитой качественной теории диссипативных систем нелинейных дифференциальных уравнений в частных производных до сих пор нет. Методы исследования моделей связаны главным образом с численным решением начально-краевых задач.

В настоящей работе рассматривается модель типа "реакция-диффузия", которая описывает реакцию окисления СО в кристаллите цеолита диаметром « 5 мкм с внедренными в него наночастицами палладия (Глава 3).

В последнее время в литературе рассматривается новый класс моделей мезоско-пического уровня. Эти модели предназначены для описания реакции в реакторах малого объема, например, на поверхности наночастиц или на металлических остриях размером от нескольких единиц до нескольких десятков нанометров. В реакции, протекающей в ограниченном объеме, принимает участие относительно небольшое (~ 10³ -г 10⁴) число частиц реагентов. Для описания процессов в такой системе используется марковская модель с конечным числом состояний и с непрерывным временем. Состояние модели представляет собой набор чисел Л/* = {Ni,N₂,...,N_n}, где N{ - число молекул г'-го реагента. Допустимые переходы между состояниями

соответствуют элементарным стадиям реакции. Марковская модель естественным образом учитывает влияние статистических флуктуации и корреляций числа частиц реагентов на динамику реакции. Вероятности переходов зависят от текущего состояния системы. Для газофазных реакций такие модели были предложены в работах Гиллеспи [25] (см. также [26]).

В настоящей работе мезоскопическая стохастическая модель используется для моделирования реакции окисления СО на наночастице палладия (Глава 4).

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

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

частиц.

Адсорбционный слой на поверхности твердого тела можно рассматривать как решеточный газ, молекулы которого располагаются в узлах регулярной плоской или пространственной (многослойная адсорбция) решетки. Модель решеточного газа является одной из основных моделей статистической механики. В настоящее время эту модель широко используют для моделирования поверхностных процессов на микроскопическом уровне. Основные параметры модели - энергии латеральных взаимодействий частиц, расположенных в различных узлах решетки. Каталитический процесс состоит из совокупности элементарных стадий адсорбции, десорбции, миграции и реакции. Для моделирования эволюции решеточного газа широко применяется стохастический подход, основанный на динамическом методе Монте-Карло [27, 28], однако известны также детерминистические решеточные модели, построенные с помощью приближения среднего поля и квазихимического приближения [29, 30].

В настоящей работе на микроскопическом уровне моделируется рост полупроводниковой кристаллической пленки соединения типа А³ В⁵ в процессе молекулярно-лучевой эпитаксии (см. Главу 6).

1.3 Содержание работы и основные результаты 1.3.1 Основные задачи работы

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

Ряд задач, представленных в главах 2-4, связаны с моделированием реакции окисления СО на палладиевом цеолитном катализаторе. Большое количество экспериментальных данных по этой реакции, на основании которых разрабатывались математические модели разных уровней, были получены М. М. Слинько, Н. Егером и др. Было установлено, что при давлении реакционной смеси на входе реактора близком к атмосферному и температуре ~ 500 К, реакция протекает в колебательном режиме если величина парциального давления СО, Рсо? лежит в интервале 1 < Рсо < 10 Торр. Причем регулярные колебания наблюдались в относительно

узком интервале Рссъ вне этого интервала колебания были нерегулярными. Также были получены данные, касающиеся влияния размера частиц палладия и степени их предварительного окисления на ход реакции и характер колебаний скорости реакции. Основная задача моделирования реакции состояла в том, чтобы на основе известных представлениях о возможных механизмах реакции разработать математическую модель, динамическое поведение которой качественно совпадает с динамическом поведением реакции в заданной области значений внешних параметров. В ходе исследований эта задача распалась на несколько отдельных задач.

1. Разработать точечную математическую модель автоколебаний в реакции окисления СО на палладии при атмосферном давлении, которая описывает регулярные колебания скорости реакции.

2. Разработать распределенную модель реакции окисления СО в палладиевом цеолитном катализаторе, которая описывает сложные и хаотические колебания.

3. Разработать точечную модель реакции на поверхности наночастицы палладия с учетом влияния размера частицы и степени окисления катализатора.

4. Разработать модель реакции на наночастицах палладия, встроенных в матрицу цеолита.

В ходе работы над этими задачами были получены новые экспериментальные данные в объеме, позволяющем провести анализ временных серий. В связи с эти возникла следующая задача

5) Разработать комплекс программ для нелинейного анализа временных серий и
выполнить обработку экспериментальных данных.

Работы по математическому моделированию реакции окисления СО на палладиевом цеолитном катализаторе проводились под общим руководством члена-корреспондента РАН М. Г. Слинько.

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

6. Разработать стохастическую решеточную модель роста соединения типа А³В⁵ в молекулярно-лучевой эпитаксии, которая качественно верно воспроизводит экспериментально установленные механизмы роста.

7. Исследовать кинетику образования объемных вакансий в процессе роста GaAs в МЛЭ.

8. Исследовать кинетику роста слоев трехкомпонентного соединения А^В^С⁵ на ступенчатой грани А³С⁵.

1.3.2 Содержание работы

Работа состоит из шести глав и выводов.

Во второй главе представлена новая точечная модель колебаний в реакции окисления СО на поверхности палладия при атмосферном давлении. Модель построена на основе двух известных моделей: модели Sales-Turner-Maple (STM) [31] и модели Bassett-Imbihl (IB) [32, 33]. Первая модель описывает колебания в реакции на поликристаллическом палладии при атмосферном давлении, вторая - на грани (110) монокристалла палладия в сверхвысоком вакууме. Для объяснения колебаний в модели STM используется концепция периодического окисления-восстановления поверхности, а в модели IB - концепция подповерхностного кислорода. В настоящее время концепция подповерхностного кислорода считается более предпочтительной, поэтому она используется в новой модели, однако значения констант скоростей в новой модели близки к константам модели STM.

Механизм колебаний в точечной модели связан с периодическим изменением концентрации кислорода в первом подповерхностном слое в кристаллической решетке палладия. За один период колебаний происходят следующие процессы. В

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

Точечная модель реакции представляет собой систему трех нелинейных ОДУ первого порядка, которые описывают эволюцию поверхностных покрытий СО и кислорода и концентрации подповерхностного кислорода. Колебания в системе возникают в результате сверхкритической бифуркации Андронова-Хопфа. Параметры модели подбирались так, чтобы колебания в системе существовали в заданном диапазоне Р_Со и при этом качественно воспроизводились особенности изменения характеристик колебаний (периода и формы колебаний) при изменении давления

Рсо-

Из трех переменных точечной модели скорость изменения одной, - концентрации подповерхностного кислорода, много меньше скоростей изменения двух других переменных. Для систем с быстрыми и медленными переменными характерны траектории-утки, которые существуют, когда значения управляющего параметра лежат в достаточно малой окрестности точки бифуркации Хопфа. Особенностью траекторий-уток, в частности, циклов-уток, является чрезвычайно сильная зависимость от управляющего параметра. Период и размер предельного цикла может изменяться на несколько порядков при исчезающе малом изменении параметра. Циклы-утки в точечной модели реакции рассмотрены во второй части главы 2.

Третья глава посвящена результатам моделирования реакции окисления СО в зерне цеолитного катализатора на мезоуровне с помощью распределенной модели типа "реакция-диффузия". Предполагается, что зерно цеолитного катализатора

представляет собой кристаллит цеолита сферической формы. Диаметр цеолитной сферы порядка 10~⁶-10^-5 м. С помощью специальной технологии внутри кристаллита цеолита выращиваются наночастицы палладия с весьма узким распределением по размерам и равномерным распределением частиц по объему цеолита. В зависимости от размера наночастиц и общей массы палладия внутри сферы может располагаться ~ 10⁵-10⁷ частиц. Молекулы СО и Ог проникают внутрь зерна по порам цеолита, адсорбируются на поверхность Pd и реагируют между собой.

В простейшем варианте распределенной модели, который называется двухслойной моделью реакции, зерно цеолита делиться на две части равного объема, - внутреннюю (центральную) и внешнюю (приповерхностную). Концентрации реагентов в газовой фазе в каждой части считаются постоянными, т.е. независящими от пространственных координат. Реакция в каждой части зерна описывается точечной моделью. При расчете массообмена между поверхностью металла и газовой фазой учитывается суммарная площадь поверхности металлических частиц и объем зерна цеолита. Внешняя и внутренняя части зерна связаны диффузией реагентов в газовой фазе в порах цеолита. Когда реакция протекает в автоколебательном режиме модель представляет собой систему двух связанных химических осцилляторов. Особенность модели состоит в том, что осцилляторы связаны не через динамические переменные, как обычно бывает при диффузионной связи, а через концентрацию СО в газовой фазе, которая является параметром точечной модели. Таким образом, модель представляет собой систему параметрически связанных осцилляторов. Хорошо известно, что в системе связанных нелинейных осцилляторов возможны сложные динамические режимы, включая детерминированный хаос [34, 35, 36]. В работе построена фазовая диаграмма двухслойной модели на плоскости параметров (Рсо, D), где Рсо - парциальное давление СО вне зерна, D = D0/R2, А) ^_ коэффициент диффузии СО в матрице цеолита, R - радиус зерна. На диаграмме выделены области периодических, квазипериодических и хаотических колебаний. Показано, что периодические колебания переходят в квазипериодические через вторичную бифуркацию Хопфа. Подробно исследованы хаотические колебания, которые существуют в относительно узком интервале значений давления СО вблизи правого

конца интервала колебаний.

Тонкая структура колебаний в реакции в зерне впервые обнаружена и описана с помощью непрерывной модели, которая описывает реакцию в пористом шаре с помощью системы уравнений типа "реакция-диффузия". В этой модели предполагается, что активная поверхность катализатора равномерно распределена по объему пористого шара. В силу сферической симметрии задачи динамические переменные системы зависят только от одной пространственной координаты - радиуса г, и от времени. В ходе реакции в модели формируется неоднородное распределение концентрации СО вдоль радиуса шара: концентрация СО уменьшается от поверхности к центру. Поэтому локальные химические осцилляторы, расположенные на разных расстояниях от центра, оказываются в неодинаковых условиях, что приводит к сложному динамическому поведению модели. В работе представлен численный анализ модели в интервале значений внешнего давления Рсо? (Ръ Рг)₅ на котором существует колебательное решение, при некотором фиксированном значении отношения D₀/R². Показано, что на большей части интервала автоколебаний в системе существуют сферические волны скорости реакции, распространяющиеся от поверхности зерна к его центру, причем скорость волны зависит от г. Анализ средних частот локальных осцилляторов показал, что в системе формируется неоднородное и нелинейное распределение средних частот w_a(r), такое что ш_а(г) постоянно при г < _с, где _с, 0 < _с < R, некоторая точка внутри шара. Таким образом, все осцилляторы в центральной области г < _с зерна образуют один кластер. Вид зависимости w_a(r) при _с < г < R зависит от значения Рсо? а также от наличия или отсутствия диффузии адсорбированного СО по поверхности катализатора, однако всегда частота осцилляторов на поверхности зерна больше частоты в центре.

В отсутствии диффузии СО по поверхности металла характерная особенность пространственно-временного поведения реакции связана с критической точкой _с в распределении средних частот локальных осцилляторов. В окрестности этой точки скорость волн резко замедляется и справа от точки _с на интервале _с < г < R число волн растет линейно со временем. С течением времени отдельные импульсы вблизи точки _с становятся настолько узкими, что дальнейшее численное решение системы

оказывается невозможным без сгущения сетки в окрестности _с. При приближении давления Fco к правому концу Р₂ интервала колебаний пространственно-временное поведение системы усложняется. Справа от критической точки _с между центральным кластером и приповерхностной областью, в которой накапливаются волны, образуется новый кластер осцилляторов, причем колебания в этом кластере происходят в противофазе с колебаниями в центральном кластере. В некоторой малой окрестности Pi колебания еще более усложняются и становятся нерегулярными.

При наличии поверхностной диффузии СО волны скорости реакции не могут скапливаться вблизи точки _с, так как диффузия разрушает достаточно узкие импульсы. В этом случае распределение частот w_a(r) локальных осцилляторов имеет вид кусочно-постоянной функции. На интервалах постоянства распределения oj_a(r) локальные осцилляторы образуют кластер. В левой части интервала колебаний (Pi, Р2), при низких давлениях, имеется два кластера, центральный и приповерхностный, причем частота колебаний внешнего кластера больше чем внутреннего. Через равные интервалы времени на границе между кластерами происходит потеря одного импульса. При увеличении давления между внешним и внутренним кластерами появляется промежуточный кластер с промежуточной частотой. При дальнейшем росте давления кластеров становится четыре и т.д. Также как и в отсутствии поверхностной диффузии, вблизи Р₂ колебания становятся нерегулярными.

В четвертой главе представлена новая точечная модель реакции окисления СО на поверхности частицы палладия с учетом влияния на скорость реакции степени окисления палладия и размера частицы. Экспериментально установлено, что реакция происходит по разному на предварительно окисленном катализаторе и на восстановленном катализаторе. Восстановленный катализатор более активен. Вид колебаний на окисленном и восстановленном катализаторах также различный. На восстановленном катализаторе большую часть периода реакция протекает в активном режиме с высокой скоростью реакции, а на окисленном - в пассивном режиме. Если реакция на предварительно окисленном катализаторе протекает достаточно

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

Во второй части главы 4 представлена стохастическая модель реакции на нано-частице. Кристаллическую наночастицу металла с гранецентрированной кристаллической решеткой можно представить в виде октаэдра, каждая грань которого совпадает с кристаллографической плоскостью (111). Для частиц такой формы с диаметром не более 10 нм число поверхностных атомов не превышает нескольких тысяч. Следовательно, число атомов реагентов, адсорбированных на поверхности частицы, также относительно невелико (~ 1 — 10³). Для столь малых систем описание реакции в приближении среднего поля может быть слишком грубым. В настоящее время наиболее адекватными моделями таких систем являются микроскопические стохастические модели, основанные на динамическом методе Монте-Карло [37]. Однако, в случае цеолитного катализатора, который состоит из огромного числа наночастиц, важнейшую роль в динамике реакции играет диффузионное взаимодействие между процессами на отдельных наночастицах. Поскольку в данной работе главное внимание сосредоточено на изучении роли взаимодействия в динамике реакции, то детали протекания реакции на отдельных частицах представляются не столь важными. Для изучения реакции на больших ансамблях частиц достаточно иметь простую и грубую модель реакции на одной частице. В качестве такой

модели предложена мезоскопическая стохастическая модель, которая представляет собой марковскую модель с конечным числом состояний. Каждое состояние модели характеризуется набором чисел {іУі,і\Гг, ...,iV„}, где iV; - число молекул г'-го сорта на поверхности частицы, 0 < ]Г)»-ЭД < N_tot. Переходы в модели соответствуют элементарным стадиям реакции, скорости переходов согласуются со скоростями элементарных стадий в приближении среднего поля, так что при N_tot —) со динамика средних концентраций в марковской модели описывается уравнениями среднего поля. Марковская модель естественным образом учитывает влияние флуктуации концентраций реагентов на динамику реакции и позволяет определить границы применимости приближения среднего поля для описания реакции на наночастице.

В третьей части главы сделана попытка учесть реальную структуру цеолитного катализатора и основные физико-химические процессы. Модель описывает динамику реакции на поверхности N наночастиц Pd октаэдральной формы, расположенных вдоль прямой линии в матрице цеолита, т.е. на одномерной цепочке наночастиц. Реагенты поступают внутри цепочки через крайние звенья, на внешних границах которых поддерживаются постоянные концентрации реагентов. Реакция на каждой частице описывается точечной моделью. Для описания диффузии смеси газов в порах цеолита используется теория Максвелла-Стефана [38]. Сорбция газов матрицей цеолита описывается моделью адсорбции Лэнгмюра.

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

СО колебания скорости реакции возникают сначала на крайних частицах и являются периодическими. При небольшом увеличении давления колебания наблюдаются на двух крайних частицах, причем с разными частотами так что суммарные колебания являются квазипериодическими. При дальнейшем возрастании давления наблюдается достаточно узкий интервал давлений, на котором центральная часть цепочки частиц ведет себя как возбудимая среда. Колебания на крайних частицах в цепочке действуют как возмущение на центральную часть и, если это возмущение достаточно велико, реакция на все цепочке спонтанно совершает один цикл колебаний. Такие колебания происходят в случайные моменты времени и суммарные колебания скорости реакции являются хаотическими. При дальнейшем увеличении давления СО характер колебаний изменяется. Частицы в центральной части цепочки образуют большой центральный кластер, на котором наблюдаются регулярные синхронные колебания, в то время как на концах цепочки реакция проявляет нерегулярные колебания. Результирующие колебания выглядят как сумма регулярных колебаний и высокочастотного шума с относительно малой амплитудой. Таким образом эта модель наиболее близко описывает экспериментальные данные, в которых колебания имеют нерегулярный вид в широком диапазоне давлений СО, а регулярные колебания наблюдаются в узком диапазоне давлений СО.

В пятой главе представлен комплекс компьютерных программ для анализа нелинейных колебаний и результаты анализа экспериментальных данных по реакции окисления СО на палладиевом цеолитном катализаторе с помощью этих программ. В комплекс включены наиболее известные алгоритмы нелинейного анализа временных серий [39], а именно, алгоритмы реконструкции фазового пространства, оценки корреляционной размерности и размерности вложения, оценки максимального показателя Ляпунова. Кроме того в комплекс входят стандартные программы анализа Фурье, фильтрации шума, построения одномерных и двумерных отображений и некоторые другие алгоритмы.

В работе приведены результаты анализа временных серий скорости реакции окисления СО на палладиевом цеолитном катализаторе. Оказалось, что имеющиеся

данные весьма трудно поддаются обработке и результаты анализа неоднозначны. Известные алгоритмы определения величины времени задержки для реконструкции фазовых переменных из скалярной временной серии не дают определенного ответа. Размерность вложения d^ также плохо оценивается. Из результатов анализа следует, что dg f» 10. Алгоритмы оценки максимального показателя Ляпунова для временных серий дают положительную величину, лежащую в интервале (0.05, 0.15). Отметим, что в двухслойной модели колебаний значения старшего показателя Ляпунова хаотических колебаний также лежат в этом интервале. Кроме того, результаты анализа временных серий концентрации С0₂ в модели колебаний в реакции на цепочке частиц весьма близки к результатам анализа экспериментальных данных. Таким образом, можно сделать вывод, что представленные в работе распределенные модели реакции качественно верно описывают динамику реакции.

Шестая глава посвящена моделированию динамики поверхностной каталитической реакции на микроскопическом уровне. Здесь приведено описание модели и результаты моделирования роста кристалла полупроводникового соединения типа А³ В⁵ (наиболее важным представителем соединений типа А³ В⁵ является GaAs) в процессе молекулярно-лучевой эпитаксии (МЛЭ).

С помощью МЛЭ выращивают сверхтонкие монокристаллические пленки на ориентированной грани монокристалла-подложки. Рост пленки происходит в сверхвысоком вакууме. Атомы и/или молекулы вещества растущей пленки доставляются на поверхность подложки посредством молекулярных пучков, испускаемых специальными эффузионными ячейками. Рост пленки происходит в термодинамически неравновесных условиях; интенсивность падающих на поверхность молекулярных потоков намного больше интенсивности потоков десорбированных частиц. Для получения ровной поверхности пленки подложку нагревают до температуры ~ 800 К. Интенсивности потоков образующих элементов таковы, что средняя скорость роста пленки порядка одного атомного монослоя в секунду. Метод МЛЭ позволяет использовать различные средства контроля качества поверхности и химического состава пленки в процессе роста и предоставляет широкие возможности управления

процессом роста.

Процесс роста кристалла в МЛЭ является весьма сложным неравновесным многофакторным процессом. В настоящее время нет аналитической модели которая описывала бы все стадии роста соединения А³В⁵ в достаточно широком диапазоне внешних условий. Поэтому наиболее адекватное моделирование МЛЭ получается с помощью стохастических решеточных моделей, основанных на динамическом методе Монте-Карло. В шестой главе описана стохастическая решеточная модель роста в МЛЭ трехкомпонентного соединения А^В^С⁵ как на плоской, так и на вици-нальных (ступенчатых) гранях (001). Отличительной особенностью модели является возможность образования объемных вакансий в процессе роста. С помощью стохастической модели было выполнено изучение кинетики формирования объемных вакансий в МЛЭ GaAs и кинетика формирования неоднородного латерального распределения компонент А и В в процессе МЛЭ соединения на ступенчатых гранях.

1.3.3 Основные результаты работы

Результаты работы можно разделить на три группы. Первая группа результатов относится к математическому моделированию каталитической реакции окисления СО на палладиевом цеолитном катализаторе.

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

2. Разработана и исследована дискретная (в частности, двухслойная) распределенная модель реакции в зерне цеолитного катализатора, которая представляет собой систему параметрически связанных нелинейных химических осцилляторов. Модель описывает периодические, квазипериодические и хаотические колебания скорости реакции.

3. Ралпаботана и исслелована модель пеакпии в зепне типа "пеактшя-лисЬсЪузия"

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

4. Разработана и исследована детерминистическая точечная модель реакции на поверхности наночастицы палладия с учетом влияния на скорость реакции степени окисленности палладия и размера частицы.

5. Разработана и исследована мезоскопическая стохастическая марковская модель реакции на поверхности наночастицы. Изучена роль флуктуации числа адатомов на поверхности в динамике реакции. Описаны колебания, индуцированные шумом.

6. Разработана и исследована модель реакции на линейной цепочке наночастиц катализатора, встроенных в матрицу цеолита. Модель наиболее полно учитывает характерные особенности строения цеолитного катализатора. В модели обнаружено два качественно различных динамических режима, которые реализуются при достаточно малом и достаточно большом отношении к скорости сорбции к скорости диффузии смеси газов в цеолите. При ас < 1 колебания скорости реакции на частицах катализатора происходят синхронно. При к > 1 наблюдаются сложное динамическое поведение реакции, которое качественно соответствует экспериментальным данным.

Вторая группа результатов связана с анализом экспериментальных временных серий.

7) Разработан комплекс программ для нелинейного анализа временных серий.

В комплекс включены наиболее известные алгоритмы восстановления фазовых переменных из скалярной серии, вычисления корреляционной размерности и минимальной размерности вложения, оценки максимального показателя Ляпунова. В комплекс включены также традиционные алгоритмы линейного анализа, в частности, алгоритмы фильтрации шума, анализа Фурье, вейвлет-анализа и др. Программы написаны на языке C++ и работают в среде Windows.

8) Проведен анализ экспериментальных временных серий скорости реакции окис
ления СО на палладиевом цеолитном катализаторе. Показано, что времен
ные серии являются хаотическими с весьма высокой размерностью вложения
(~ 10). Каталитическая система является достаточно сложной и результаты
анализа временных серий говорят о том, что она не может быть адекватно
описана простой конечномерной моделью.

И, наконец, третья группа результатов связана с моделированием роста кристаллов в молекулярно-лучевой эпитаксии.

9) Разработана стохастическая решеточная модель роста в МЛЭ полупроводни
кового соединения типа А³В⁵, которая качественно воспроизводит основные
характерные особенности процесса МЛЭ. Модель основана на динамическом
методе Монте-Карло и моделирует процесс роста в виде марковской после
довательности элементарных актов, соответствующих элементарным стадиям
каталитической реакции роста соединения А³В⁵, - адсорбции, поверхностной
миграции, десорбции, встраивания в кристаллическую решетку.

10) С помощью стохастической модели изучена кинетика формирования объемных вакансий в процессе МЛЭ GaAs. Показано, что при увеличении отношения интенсивностей потоков мышьяка и галлия увеличивается концентрация объемных вакансий, что приводит к ухудшению электрофизических свойств кристаллической пленки. Результаты моделирования объяснили данные экспериментов по выращиванию пленок при избыточных давлениях мышьяка.

И) С помощью стохастической модели описана динамическая неустойчивость процесса роста, которая возникает на начальной стадии роста слоев Ajj ₅Bq ₅С⁵ на ступенчатой поверхности А³С⁵ и приводит к латерально неоднородному распределению компонент А и В в нескольких слоях в перпендикулярном ступеням направлении. Неоднородное распределение может быть использовано для создания латеральных сверхрешеток.

Результаты опубликованы в 21 статье в отечественных и зарубежных изданиях и докладывались на 10 российских и международных конференциях.

1.4 Методы исследования

Для решения поставленных в работе задач автором были разработаны пакеты программ для численного решения и комплексного анализа решений систем нелинейных ОДУ, для нелинейного анализа временных серий и для моделирования МЛЭ соединений типа А³В⁵ методом Монте-Карло. Пакеты программ для анализа временных серий и моделирования роста кристаллов описаны, соответственно, в главах 4 и 5. В этом разделе описаны пакеты программ для бифуркационного анализа систем ОДУ и для численного анализа периодических и хаотических решений систем ОДУ, которые использовались для работы с точечными и распределенными моделями в главах 2 и 3.

1.4.1 Качественный анализ систем ОДУ

Изучение любой нелинейной модели начинается, как правило, с определения возможных типов асимптотических решений в заданной области значений параметров модели. Для этой цели применяются методы бифуркационного анализа [40] и прежде всего - метод продолжения по параметру стационарных решений системы. В комплексе программ использовался алгоритм продолжения по параметру, который состоит в следующем.

Пусть система эволюционных уравнений с достаточно гладкой правой частью

±X(t) = F(X(t),p), (1.4.2)

где t - время, X Ж" - вектор независимых переменных и F К" - вектор правых частей, зависящий от параметра р, имеет при р = р₀ стационарное решение Х₀(ро). Т.е. Хо = Хо(ро) является решением системы уравнений

F(X,p₀) = 0. (1.4.3)

Если вектор-функция F(X,p) дифференцируема в окрестности (Хо,ро), то для достаточно малого 5р можно написать

д д

F(X₀ + SX,_Po + Sp) и F(X₀,p₀) + g^SX + —Sp, (1.4.4)

где производные вычисляются в точке (Х₀,ро)- Для построения стационарного решения (1.4.2) дляр (ро—Sp,po+5p) используется метод продолжения по параметру с помощью алгоритма типа "предиктор - корректор".

Пусть Х₀ + SX - решение уравнения (1.4.3) при р = р₀ + Sp и определитель матрицы Якоби J(X,p) = dF/дХ не равен нулю в точке (Х₀,ро), то по теореме о неявной функции в некоторой окрестности ро существует непрерывная функция Х(р) такая, что F(X(p),p) = 0. Тогда из (1.4.4) имеем

6Х =

dF"¹"¹

В^Р. (1-4.5)

Принимая во внимание это равенство, алгоритм продолжения по параметру построен следующим образом.

Начальное значение ро параметра, по возможности, выбирают таким, чтобы решение уравнения (1.4.3) определялось достаточно легко. Для тех случаях, когда угадать хорошее приближение для стационарного решения не удается, в программе предусмотрено численное решение системы (1.4.2) до установления стационарного (стабильного) решения. На первом шаге при р = ро методом Ньютона решается уравнение (1.4.3) и находится стационарное решение Хо- Затем на стационарном решении вычисляется матрица производных правой части F системы (1.4.2) по независимым переменным X и по параметру р. Если определитель det(J(Xo,p)) Ф 0, задается приращение параметра Sp, вычисляется вектор приращений независимых переменных SX и полагается ХІ⁰⁾ = Хо + SX. Вектор Х[⁰⁾ (прогноз) используется

в качестве начального приближения для метода Ньютона при решении уравнения F(X,p) = 0 (коррекция). Если нужная точность решения не достигается за заданное число итераций, то приращение параметра 5р уменьшается в два раза и описанная процедура повторяется. В результате находится новое значение параметра р_г = р₀ + Sp и новое стационарное решение Xi. Описанные действия повторяются до тех пор, пока параметр не пробежит заданный интервал значений.

В ходе продолжения стационарного решения по параметру вычисляются собственные значения ReAi > ЫеАг > ... > ReA_n матрицы Якоби J(X,p)) и анализируется устойчивость стационарного решения в линейном приближении [41]. В этом приближении стационарное решение устойчиво, если ReAi < 0, т.е. все собственные значения матрицы Якоби лежат в левой полуплоскости комплексной плоскости. Свойства стационарного решения и/или число стационарных решений изменяются при тех значениях параметра, когда одно или больше собственных значений пересекают линию ReA = 0. Такие значения называются бифуркационными значениями параметра, в которых происходят бифуркации стационарного решения. Если в момент бифуркации мнимую ось пересекают вещественные собственные значения матрицы Якоби (при этом det(J) = 0), то бифуркация связана с изменением числа стационарных решений (точка ветвления стационарного решения); если мнимую ось пересекает пара комплексно-сопряженных собственных значений, то бифуркация связана с появлением или исчезновением периодического решения (бифуркация Андронова-Хопфа).

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

Точка поворота р* характеризуется тем, что в этой точке ранг матрицы Якоби J(X*,p*) равен п — 1 (для точки ветвления ранг J не больше п — 1). При движении

по стационарному решению в этой точке одно из вещественных собственных значений матрицы Якоби Л*(р) меняет знак. Если к — 1, т. е. максимальное собственное значение пересекает мнимую ось, то стационарное решение становится неустойчивым при перемене знака \\{р) с '+' на '—', либо устойчивым при обратной перемене знака. В этом случае при р = р* сливаются два стационарных решения: одно типа "устойчивый узел", другое типа "седло", поэтому такую бифуркацию называют также седло-узловой бифуркацией. В точке поворота число стационарных решений системы изменяется на 2; два стационарных решения либо появляются, либо исчезают при дальнейшем изменении параметра в том же направлении, что и до точки поворота.

При приближении к точке поворота р* прогноз стационарного решения по формуле (1.4.5) становится невозможным из-за вырождения матрицы J. Поэтому в алгоритме продолжения решения параметр р автоматически заменяется на одну из переменных системы Хк, к = 1,...,п, для которой матрица J*, полученная из матрицы J заменой к-то столбца dF/dxk на столбец dF/dp, не является вырожденной при р = р*. Таким образом, в окрестности точки поворота стационарное решение продолжается по переменной х^.

Для построения линии точек поворота на плоскости двух параметров (рі,Рг) в программе используется такой же алгоритм продолжения по параметру стационарного решения для системы п + 1 уравнений. Эта система состоит из п уравнений (1.4.3) и добавочного уравнения det(J(X,pi,p₂)) = 0. Процедура продолжения начинается с одной из точек поворота, найденных при продолжении стационарного решения по одному параметру.

Точки бифуркации Андронова-Хопфа при продолжении стационарного решения по параметру р определяются как такие точки р*, в которых пара комплексно-сопряженных собственных значений Лі(р), Л₂(р) = Ai(p), ImAi ф 0, пересекает мнимую ось. Программа отмечает такие точки при движении по стационарному решению, однако она не определяет характер бифуркации, которая может быть подкритической или сверхкритической. В простых случаях тип бифуркации можно установить по бифуркационной диаграмме предельного цикла, построенного с

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

Линия бифуркаций Андронова-Хопфа на плоскости двух параметров (рі,рг) строится с помощью алгоритма продолжения по параметру стационарного решения системы п + 1 уравнений. В этой системе к п уравнениям (1.4.3) добавлено уравнение Re(Ai + Аг) = 0, которое определяет переход пары собственных значений через мнимую ось. Заметим, что этому условию удовлетворяет также седло с вещественными Ai, А₂, такими, что Ai = —А2, однако этот случай легко отличить от бифуркации Андронова-Хопфа.

Все бифуркационные и фазовые диаграммы, которые приводятся в данной работе, построены с помощью представленной в этом разделе программы бифуркационного анализа систем ОДУ.

1.4.2 Численный анализ решений систем ОДУ

Системы обыкновенных дифференциальных уравнений, которые представляют модели химической кинетики, являются, как правило, достаточно жесткими. Характерные скорости быстрых и медленных стадий реакции могут различаться на много порядков величины/Поэтому для их решения используются специальные методы, предназначенные для численного решения жестких систем. В данной работе используется метод Гира с адаптивным выбором шага интегрирования и порядка метода и метод Розенброка четвертого порядка точности. Эти методы широко известны и описаны в стандартных учебниках по численным методам. Программы, реализующие алгоритмы обоих методов, взяты из [42].

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

В программу встроены некоторые алгоритмы анализа периодических и хаотических решений.

Анализ периодических решений

Периодические решения нелинейных диссипативных систем представляют собой предельные циклы. Предельные циклы являются асимптотическими решениями; решение системы с начальным условием из области притяжения предельного цикла асимптотически приближается к нему при t —> оо. По этой причине и в силу приближенного характера численных методов, численное решение можно рассматривать как периодическое только с определенной степенью точности.

Важную роль в исследовании периодических и апериодических осциллирующих решений играет функция последования Пуанкаре. Для построения функции по-следования используется поверхность в фазовом пространстве, трансверсальная к фазовым траекториям. В качестве такой поверхности часто используется секущая (гипер)плоскость. В программе секущая плоскость Р ставится следующим образом. По графикам динамических переменных системы, которые выдаются в ходе решения задачи, пользователь выбирает момент времени (или точку на фазовой траектории), когда нужно ставить плоскость. Решение задачи приостанавливается и на экран выдается диалоговое окно с фазовыми координатами и компонентами скорости фазовой точки в выбранный момент времени. По умолчанию секущая плоскость Р проходит через выбранную точку Х_р и перпендикулярна вектору скорости F(X_p) в этой точки, однако при необходимости положение и ориентацию плоскости в фазовом пространстве можно изменить.

Заданная таким образом плоскость Р используется для построения функции последования с помощью следующей процедуры [43]. В ходе численного решения задачи определяется ситуация, когда две последовательные точки X,_i и X,- фазовой траектории лежат по разные стороны от плоскости Р, причем переход через плоскость происходит в направлении вектора нормали (т.е. отклонение фазовой точки от секущей плоскости становится положительным). В этом случае решение возвращается в точку Х,_і, шаг интегрирования уменьшается в 10 раз и решение

продолжается с уменьшенным шагом до повторения описанной ситуации. Процесс уменьшения шага продолжается до тех пор, пока расстояние между точками, лежащими по разные стороны от секущей плоскости Р не станет меньше заданной величины. Затем фазовые координаты X,_i и X,- и скорости F(X,_i) и F(X,) используются для построения полиномов третьей степени Х{д) и Т(д), где д(Х) -отклонение точки X от плоскости Р. Полином Х(д) интерполирует фазовую траекторию между точками Х,_! и X,-, а Т{д) время как функцию д. Момент пересечения траекторией плоскости Р и координаты точки пересечения определяются как Т(0) и Д'(О), соответственно.

Пусть Х_р произвольная точка на плоскости Р и пусть фазовая траектория, проходящая через эту точку, в следующий раз пересекает плоскость Р, двигаясь в том же направлении, в точке Х_р. Соответствие Х_р —ї Х'_р задает функцию после-дования V, т.е. V{X_P) = Х_р. Пусть Хр , п = 1,2,3,..., - последовательные точки пересечения. Если эта последовательность имеет предел X*, то он удовлетворяет уравнению V(Xp) = X*. На практике одно-оборотный предельный цикл считается установившимся, если начиная с некоторого п выполняется неравенство

||_Х(")-^(Х_р"))||<є (1.4.6)

для заданного достаточно малого є. Для m-оборотного цикла должно выполняться неравенство

||X_pⁿ⁾ - X_p^n+m)(= V^m(X^))\\ < є. (1.6.5m)

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

Пусть Х_р - точка пересечения предельного цикла системы (1.4.2) с секущей плоскостью Р. Тогда V(X_P) = Хр и для близкой точки X' = Х_р + SX имеем

V(X') « Хр + (J^(X_P)) SX. (1.4.7)

_д_ ЭХ

П*о)

Предельный цикл устойчив, если при t —» со к нему сходятся все траектории динамической системы, которые начинаются в некоторой окрестности предельного цикла. Из (1.4.7) следует, что для устойчивости цикла Х₀() достаточно, чтобы выполнялось неравенство

Важной характеристикой предельного цикла являются его мультипликаторы (см., например. [41] стр. 187). Пусть Х₀() - периодическое решение системы (1.4.2) с периодом Т (Г-периодическое решение) и пусть ХДі) = Хо() + SX(t) - решение этой системы с начальным условием Xs(t₀) = Хо(о) + <^Хо. В линейном приближении Х() удовлетворяет системе линейных уравнений с Г-периодическими коэффициентами

Tt^{5X{t) =} (Jt^F(Xo{)))^sx{t)- ^(1А9)

Пусть Ф() - фундаментальная матрица решения линейной системы уравнений

X(i) = A(t)X{t) (1.4.10)

с Т-периодической матрицей коэффициентов А(), удовлетворяющая условию Ф(0) = Е. Тогда матрица Ф( + !Г) также является фундаментальной. Действительно, на основании тождества

Ф(*)=А(*)Ф(*)

имеем

^Ф(* + Т) = Ф(* + Г)^(* + Т) = A(t + Г)Ф(* + Т) = ЩЩ + Т).

Следовательно, Ф( + Т) также является фундаментальной.

Поскольку решение (1.4.10) с начальным условием Хо выражается через фундаментальную матрицу Х() = Ф()Х₀, для матрицы Ф( + Г) можно написать

Ф(* + Г) = Ф(і)С,

где С - постоянная не особенная матрица. Полагая здесь t = 0 и учитывая начальное условие для Ф, находим

Ф(* + Г) = Ф(*)Ф(Т). (1.4.11)

Матрица Ф(Т) называется матрицей монодромии.

Для системы в вариациях (1.4.9) матрица монодромии преобразует отклонение <Хо от предельного цикла в момент времени t₀ в отклонение SX-i — Ф(!Г)<5Хо в момент Ц +Т. Это преобразование аналогично действию функции последования (для достаточно малых отклонений) с той лишь разницей, что не фиксируется секущая плоскость. Об устойчивости цикла судят по собственным значениям матрицы монодромии, которые называются мультипликаторами цикла. Пусть ро, р_І5..., р„__х - собственные значения Ф(Т). Если взять отклонение <5Х₀ так, чтобы точка Х(о) + <5Хо лежала на предельном цикле, то тогда, очевидно, Ф(Т)Хо = &Х₀. Следовательно, один из мультипликаторов цикла всегда равен 1, поэтому положим р₀ = 1. Для устойчивости цикла необходимо, чтобы \pk\ < 1 при к > 1.

Бифуркации предельных циклов связаны с переходом мультипликаторов через единичную окружность на комплексной плоскости. Пусть мультипликаторы p₁,p2,...,p„_i упорядочены по убыванию модуля jpi| > ІР2І > ... > |р₀-і| и\рі\<1, т.е. цикл устойчивый. Переход мультипликатора р\ через единичную окружность в точке +1 связан с седлоузловой бифуркацией цикла, при которой устойчивый цикл сливается с неустойчивым и оба цикла исчезают. При переходе р\ через единичную окружность в точке —1 происходит бифуркация удвоения периода цикла. При этой бифуркации исходный цикл продолжает существовать, однако становится неустойчивым, и рождается цикл с периодом приблизительно равным удвоенному периоду исходного цикла. Если взять отклонение &Х₀ в направлении, соответствующему pi, то через период будем иметь SX-i = pi^Xo = -(JXo, а через два периода Хг = р?5Х₀ = SX-o, т.е. траектория возвращается в исходную точку. При переходе пары комплексно-сопряженных мультипликаторов pi и р₂ = Pi через единичную окружность происходит вторичная бифуркация Хопфа, при которой исходный цикл становится неустойчивым и рождается двумерный тор, соответствующий квазипериодическим колебаниям.

В программе для численного анализа решений систем ОДУ имеется блок вычисления мультипликаторов цикла. Для вычисления мультипликаторов задается секущая плоскость Р. Затем система (1.4.2) численно решается до выполнения условия (1.4.6), т.е. до выхода решения на предельный цикл с заданной точностью. После выхода решения на предельный цикл между последовательными пересечениями траекторией секущей плоскости вместе с системой (1.4.2) решается система в вариациях (1.4.9) с начальными условиями <Х^(_Р) = е^^к\ к = 1,2, ...,п, где eW = {0,...,0,1,0,...,0} - вектор-столбец с единицей в к-ой строке и нулями в остальных строках (е^ - орт в направлении k-ofi. оси декартовой системы координат в фазовом пространстве), t_p - момент пересечения траекторией секущей плоскости. Численная матрица монодромии получается в виде Ф(Г) = {SX.^(t_p + Г),5Х^(<р + Т),...,Х.(")(р + Т)}, а численные оценки мультипликаторов цикла определяются как собственные значения матрицы Ф(Т).

Отметим, что процедура вычисления мультипликаторов встроена в основной метод решения системы ОДУ и активизируется или прекращается при задании соответствующих команд в основном меню программы. Значения мультипликаторов выдаются на экран дисплея при каждом прохождении секущей плоскости и записываются в специальный файл. Процедура вычисления не значительно увеличивает общее время решения задачи, так как сводится к решению дополнительных систем линейных алгебраических уравнений с матрицей Якоби J, которая вычисляется при решении нелинейной системы (1.4.2).

Анализ хаотических решений

Хаотический режим является одним из (четырех) возможных асимптотических режимов в нелинейных диссипативных динамических системах. В отличие от "клас-сических"асимптотических режимов, к которым можно отнести стационарные решения (аттрактор - точка в фазовом пространстве), периодические решения (аттрактор - замкнутая траектория - цикл), квазипериодические решения (аттрактор - замкнутое fc-мерное многообразие - тор), хаотический режим был открыт относительно недавно (шестидесятые годы ХХ-го века [44]) и с тех пор подвергается

активному изучению.

Аттрактором хаотического режима является странный аттрактор, который обладает следующими свойствами [45].

а) Он является аттрактором, т.е. занимает ограниченную область фазового про
странства, к которой по истечение большого периода времени притягиваются
все траектории из области притяжения, которая может иметь очень сложную
структуру. Аттрактор состоит как бы из одной траектории, т.е. траектория с
течением времени должна пройти через каждую точку аттрактора.

б) "Странность"аттрактора выражается в высокой чувствительности к началь
ным условиям. Не смотря на сжатие в объеме, характерное для диссипативных
систем, не происходит сокращения длин во всех направлениях и расстояния
между первоначально сколь угодно близкими точками на аттракторе стано
вятся конечными через достаточно большое время.

в) Чтобы описывать физическую систему аттрактор должен быть структурно
устойчивым. Малые изменения управляющего параметра в 1.4.2 изменяют
структуру аттрактора непрерывным образом и множество параметров, для
которого 1.4.2 порождает странный аттрактор, не должно быть множеством
меры 0 - иначе аттрактор не является типичным и физически значимым.

Странный аттрактор в общем случае не является многообразием и представляет собой фрактальное множество с дробной хаусдорфовой размерностью. С точки зрения динамики хаотический режим характеризуется высокой чувствительностью к начальным данным. Это означает, что если Х_: и Х₂ две точки на хаотическом аттракторе, расстояние между которыми 5₀ = ||Хх — Х₂|| достаточно мало, и Х(; Xi) и Х(, Х₂) - выходящие из этих точек траектории динамической системы, то существует такой интервал времени и такое А > 0, что на этом интервале времени траектории экспоненциально расходятся, т.е. 8(t) = ||X(;Xi) — Х(;Х2)|| « S₀e^Xt. Скорость разбегания первоначально близких траекторий оценивается с помощью показателей Ляпунова.

Показатели Ляпунова можно определить следующим образом. Рассмотрим эволюцию в фазовом пространстве n-мерной диссипативной динамической системы сферы радиуса е. Под действием фазового потока системы сфера будет деформироваться в эллипсоид с п главными полуосями Pi(t). Показатель Ляпунова Л,- связан с длиной г-ой полуоси соотношением

Лі = lim lim log₂ ^-. (1.4.12)

Показатели Ляпунова характеризуют, таким образом, растяжение или сжатие фазового объема по разным направлениям. Однако эти направления не фиксированы, а изменяются сложным образом в процессе эволюции системы. Сумма показателей характеризует скорость изменения объема. Поскольку эволюция диссипативной системы уменьшает фазовый объем, то эта сумма должна быть отрицательной, однако по отдельным направлениям показатели могут быть положительными. Наличие положительных показателей является характерным признаком хаотического режима.

Каждому динамическому режиму соответствует определенный вид спектра показателей Ляпунова. Если система имеет только стационарные решения, то все показатели отрицательны и спектр имеет вид (—, —,..., —). Если есть периодическое решение, то старший показатель равен 0, а остальные отрицательны, спектр - (О,—,...,—). В случае квазипериодического решения число нулевых показателей равно размерности тора, спектр - (0,..., О, —,..., —). Для хаотического решения один или несколько показателей положительны, один показатель нулевой и остальные отрицательны. Например, для системы четырех уравнений возможные хаотические спектры: (+, 0, -, —) - хаос и (+, +, 0, —) - гиперхаос.

В программу решения систем ОДУ встроен модуль оценки показателей Ляпунова по методу Вольфа [46]. Исходя из определения показателей для их вычисления нужно следить за изменением инфинитезимального объема в фазовом пространстве в ходе эволюции системы. Однако на практике очень трудно определить критерий инфинитезимальности, поэтому удобнее работать в касательном пространстве с системой в вариациях (1.4.9).

Для оценки показателей Ляпунова одновременно с решением нелинейной системы (1.4.2) решается линейная система (1.4.9) для п векторов начальных данных, которые представляют собой произвольно ориентированный ортонормированный репер {vi,v₂,...,v_n} в фазовом пространстве. Для определенности будем считать, что г>^ = е(^к\ С течением времени длины векторов, в направлении которых показатели положительны, будут увеличиваться. Кроме того, каждый вектор стремится развернуться в направлении наибольшего роста длины. Поэтому со временем углы между векторами уменьшаются. Чтобы предотвратить "коллапс"репера к нему периодически применяют процедуру ортогонализации Грамма-Шмидта.

Пусть {vi,V2,...,v_n} - репер в некоторый момент времени t. Ортогонализация векторов Vk выполняется по формулам:

v[ = wx/IMI,..., vl = v_k - 5>ь4К-, v'_k = „Z/IKH,

i=i

где (,)- скалярное произведение в Е^п. При этом направление первого вектора v[ не изменяется, а каждый вектор v'_k поворачивается оставаясь в линейной оболочке векторов Vi,г»₂, , v_k так, чтобы стать перпендикулярным векторам v[, v'₂,..., v'_k_₁. Теоретически, частота выполнения ортогонализации не имеет значения, однако на практике рекомендуется выполнять ее не менее одного раза за средний период колебаний.

Предположим, что ортогонализация репера проводится в моменты времени t_n, п = 1,2,..., N. Тогда оценки показателей Ляпунова вычисляются по формулам

1 ^N

^{Afc~r^^l0g2(K(*fc)ll)' (1А13)}

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

Странный аттрактор в диссипативной системе может существовать в некоторой области значений внешних параметров. Предположим, что имеется только один управляющий параметр. При изменении параметра переход к хаосу может происходить только по вполне определенным сценариям. Таких сценариев перехода к хаосу известно три [47]. В литературе наиболее часто встречается переход к хаосу по сценарию Фейгенбаума [48]. В этом случае странный аттрактор образуется в результате бесконечной последовательности бифуркаций удвоения периода предельного цикла. Причем расстояния по параметру между последовательными бифуркациями убывают в геометрической прогресс, поэтому все бифуркации происходят на конечном интервале значений параметра. Вторым сценарием формирования странного аттрактора является сценарий Рюэля-Такенса. В этом сценарии процесс формирование странного аттрактора включает в себя несколько вторичных бифуркаций Хопфа с последующей трансформацией многомерного тора в странный аттрактор. Таким образом, первые два сценария связаны со сверхкритическими бифуркациями предельного цикла. Подкритические бифуркации связаны с третьим сценарием - перемежаемостью. Перемежаемость бывает нескольких типов, из которых наиболее известны типы I, II и III [49]. Эти типы перемежаемости связаны с подкри-тическими бифуркациями цикла, соответственно, седлоузловой, удвоения периода и вторичной бифуркацией Хопфа. При перемежаемости предельный цикл теряет устойчивость через одну из перечисленных бифуркаций и траектория динамической системы переходит на странный аттрактор. При перемежаемости странный аттрактор не формируется; он должен быть уже сформирован посредством одного из первых двух сценариев. В частности, перемежаемость часто встречается на границах интервалов периодичности, которые можно наблюдать в пространстве параметров в области хаоса.

Изменение структуры и типа аттрактора при изменении управляющего параметра можно наблюдать на бифуркационной диаграмме, которая строится следующим образом. При фиксированном значении параметра система (1.4.2) решается на достаточно большом интервале времени, так чтобы траектория приблизилась

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

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

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

максимума (минимума) вместе с диаграммой Ламерея заданной динамической переменной.

Для анализа сложных колебаний в частотной области широко используется анализ Фурье. Для построения амплитудных и энергетических спектров используются стандартные программы быстрого преобразования Фурье временных последовательностей и построения спектров (периодограмм) [42]. По виду спектров Фурье можно определить тип колебаний. Для хаотических колебаний характерен непрерывный спектр с быстрым (экспоненциальным) спадом энергии гармоник при увеличении частоты.

Наряду с анализом Фурье для изучения свойств системы связанных осцилляторов в работе применяется метод определения мгновенной частоты и фазы нелинейных колебаний, основанный на преобразовании Гильберта [50].

2 Модели колебаний в реакции окисления СО на поверхности палладия

Первым шагом в моделировании динамики реакции окисления СО в палладиевом цеолитном катализаторе была точечная модель реакции на поверхности палладия при атмосферном давлении [51]. Модель была построена на основе двух более ранних моделей, предложенных для описания автоколебаний скорости реакции на поликристаллическом палладии при атмосферном давлении - модели Sales-Turner-Maple (STM) [31], и на грани (НО) монокристалла палладия в сверхвысоком вакууме - модели Bassett-Imbihl (BI) [32, 33]. Эти модели различаются механизмами колебаний. В первой модели причиной колебаний является периодическое окисление и восстановление поверхности катализатора, во второй - периодический процесс проникновения адсорбированного кислорода в подповерхностный слой кристаллической решетки палладия (растворение кислорода) с последующим выходом кислорода на поверхность кристалла (сегрегация). В этой главе кратко описываются эти модели и подробно представлен наш вариант модели колебаний в реакции на поверхности монокристалла Pd при атмосферном давлении.

2.1 Механизмы колебаний в реакции окисления СО на палладии

2.1.1 Окисление-восстановление поверхности

Предложенная Sales, Turner и Maple в [31] модель автоколебаний в реакции окисления моноксида углерода 2СО + Ог —У 2СОг на поверхности металлов платиновой группы (Pt, Pd, Ir) основана на механизме окисления-восстановления поверхностного слоя кристалла. Несмотря на некоторые недостатки и выявленные впоследствии

противоречия с экспериментальными данными эта модель остается одной из основных моделей каталитической реакции окисления СО. В данной работе эта модель будет обозначаться аббревиатурой STM по первым буквам фамилий авторов.

В модели STM предполагается, что поверхностная реакция протекает по классическому механизму Лэнгмюра-Хиншельвуда, который содержит следующие элементарные стадии:

CO_ff + *A>CO, 0_2з + 2*-^20, СО + 0-^>С0_2!7 + 2*. (2.1.1)

В этих формулах индексом 'g' отмечены молекулы реагентов в газовой фазе, адсорбированные атомы кислорода и молекулы моноксида углерода не имеют индексов, звездочкой обозначены свободные адсорбционные центры на поверхности катализатора. Можно считать, что адсорбционные центры образуют регулярную поверхностную решетку, конкретный вид которой в данном случае не важен. Первая формула соответствует двум стадиям, - адсорбции и десорбции СО. Эта формула означает, что молекула СО из газовой фазы при адсорбции занимает один свободный узел поверхностной решетки, а при десорбции молекулы СО освобождается один узел; k±i обозначают константы скоростей адсорбции и десорбции СО. Вторая формула означает диссоциативную адсорбцию молекулы кислорода 0₂. Для адсорбции молекулы 0₂ нужны два соседних свободных узла в решетке, &₂ - константа скорости адсорбции кислорода. Третья формула описывает реакцию адсорбированных молекулы СО и атома О с образованием молекулы углекислого газа С0₂, которая немедленно десорбируется с поверхности катализатора, кз - константа скорости реакции образования С0₂.

Из перечисленных в [7] возможных механизмов автоколебаний скорости поверхностной реакции в модели STM используется "буфер". В данном случае буфером является оксидный слой на поверхности катализатора. Механизм колебаний в модели STM выглядит следующим образом.

1 Высокая скорость реакции. Поверхность преимущественно покрыта кислородом. Поверхностный слой окисляется. На окисленной поверхности нет адсорбции 0₂ и СО.

2. С увеличением площади окисленной поверхности уменьшается покрытие кислорода и растет покрытие СО. Скорость реакции падает.

3. Низкая скорость реакции. Поверхность покрыта преимущественно СО. Восстановление окисленной поверхности.

4. С уменьшением площади окисленной поверхности увеличивается поверхностная концентрация кислорода. Скорость реакции растет и система переходит в состояние 1.

В модели STM формирование слоя оксида происходит через реакцию адсорбированного кислорода с поверхностью, а восстановление - через реакцию окисленной поверхности с адсорбированным СО с образованием СОг. Поэтому к схеме (2.1.1) добавляются две стадии

ОЛО*, СО + О*- C0_2ff, (2.1.2)

где символ О* обозначает "окисленный атом "поверхности катализатора.

Дифференциальные уравнения, описывающие динамику реакции согласно схеме (2.1.1), (2.1.2), в модели STM имеют вид:

х — Pco&i(l — х — у — z) — k-ix - k_sxy — k₅xz,

У = P_0ik₂{l~x-y-z)²-k₃xy-k₄y(l-z), (2.1.3)

і = кіу(1 — z) — k₅xz]

где x(t), y(t) - зависящие от времени покрытия поверхности молекулами СО и атомами О, z(t) - относительная площадь окисленной поверхности; Рсо и Лэ₂ ~ Давление СО и кислорода в газовой фазе. Константы скоростей элементарных стадий зависят от температуры по закону Аррениуса

( RTJ'

к(Т) = к₀ ехр

где Т - температура катализатора, Е - энергия активации элементарной стадии, R - универсальная газовая постоянная (72=1.98589 кал/(К-моль)). В [31] приводятся следующие значения для параметров модели

При таких параметрах модель STM достаточно хорошо воспроизводит динамику реакции окисления СО на поверхности платины при атмосферном давлении.

Позже были получены некоторые экспериментальные подтверждения механизма окисления-восстановления поверхности катализаторов. В частности в работе [52] Turner и Maple измеряли скорость окисления и восстановления поверхности Pt, Ir и Pd. Оказалось, что измеренная скорость восстановления хорошо согласуется с теоретическим значением, а измеренная скорость окисления на порядок меньше используемой в модели. В экспериментальных работах других авторов (см. например, обзор [7]) также отмечаются противоречия с моделью STM. Однако, несмотря на это модель STM до сих пор является одной из основных моделей автоколебаний в гетерогенном катализе.

2.1.2 Растворение и сегрегация кислорода

Представленная в [32, 33] модель Bassett-Imbihl (BI), была разработана для описания колебаний в реакции 2СО + 0₂ —У 2СОг на грани (НО) монокристалла Pd, наблюдаемых при низких давлениях Ро₂ ^ 1-Ю^-3 Торр и Ро^/Рсо — Ю²-10³ и относительно низких температурах 300 К < Т < 450 К. Кинетические колебания скорости реакции при таких условиях экспериментально изучались в работах [53, 54, 55], где была построена подробная фазовая диаграмма реакции на плоскости параметров Рсо-Ро^ Однако в этих и других экспериментальных работах механизм окисления-восстановления поверхности не был обнаружен. Вместо этого было показано, что колебания скорости реакции при указанных условиях связаны с колебаниями концентрации подповерхностного кислорода. В колебательном цикле выделяются следующие стадии.

1' Высокая скорость реакции. Поверхность преимущественно покрыта кислородом. Адатомы кислорода медленно проникают внутрь кристалла палладия, -

"растворяются"в объеме катализатора.

2' Подповерхностный кислород блокирует адсорбцию кислорода. Концентрация кислорода на поверхности уменьшается, а концентрация СО возрастает, - скорость реакции падает.

3' Низкая скорости реакции. Поверхность покрыта преимущественно СО. Подповерхностный кислород медленно выходит на поверхность и вступает в реакцию с СО.

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

Таким образом, механизм растворения - выхода на поверхность кислорода действует аналогично механизму окисления-восстановления поверхности, однако с химической точки зрения это совершенно разные процессы. Кинетическая схема реакции с механизмом растворения-выхода на поверхность состоит из схемы Лэнгмюра-Хиншельвуда (2.1.1) и двух дополнительных стадий

О А О*, (2.1.4)

fc-4

которые описывают растворение и выход на поверхность кислорода с константами скоростей &4 и &_4, соответственно.

Однако одной замены механизма колебаний оказалось недостаточно. Первый вариант модели [32] не смог воспроизвести наблюдаемую фазовую диаграмму. Во втором варианте [33] были введены зависимости констант k-i и &₄ от плотности заполнения поверхности молекулами СО. В результате получилась гибридная модель, которая использовала принцип "буфера" (О*) и зависимость параметров от плотности заполнения. Эта модель хорошо описывала большинство наблюдаемых особенностей динамики реакции. Уравнения модели [33] имеют вид:

і - -PcoMl ~~ ^х) * k-ix — к₃ху,

у = PoMl-x-y)²e~"-kxV-kiV(l-z) + k-i{l-y)z, (2.1.5)

z = hy(l - z)- fc_₄(l - y)z;

k-₁ = K.₁exp\ -^7-2 1, fc₄ = ^expl _RT \. (2.1.6)

Переменные x иу имеют тот же смысл, что и в модели STM (2.1.3), a z обозначает концентрацию растворенного кислорода О* в первом приповерхностном атомном слое кристалла. Формулы для констант (2.1.6) означают, что в парах СО-СО и СО-О существует взаимодействие отталкивания. Скорости десорбции СО и растворения О увеличиваются с ростом х. Приведем значения параметров модели из [33]:

Дополнительные параметры имели следующие значения: а = 10, Е^' = 2.5, S_r^Ce?=5.0.

Отметим различия между уравнениями моделей STM и BI. В модели (2.1.5) не учитывается влияние кислорода (О и О*) на скорость адсорбции СО; влияние О* на скорость адсорбции О выражается экспонентой e~^az, а не квадратичным полиномом; стадия к_ъ заменена на к-±\ параметры к-\ и к± экспоненциально зависят от х. Таким образом, модель BI "более нелинейна", чем модель STM.

2.1.3 Модель регулярных колебаний в Pd-цеолитном катализаторе

Регулярные кинетические колебания в реакции окисления СО в палладиевом цео-литном катализаторе были обнаружены при следующих условиях: парциальное давление кислорода Ро₂ ~ 100 Торр, отношение парциальных давлений реагентов Ро₂/Рсо ~ Ю²-г 10³, температура газовой смеси в реакторе Т w 500К [56] (о структуре цеолитов и способах приготовления катализаторов см., например, обзор [57]). Авторы экспериментальных работ по изучению этой реакции (в частности, М.М. Слинько) полагали, что в основе наблюдаемых колебаний скорости реакции лежит механизм растворения и выхода на поверхность кислорода, т.е. механизм модели

BI (2.1.5). Однако рассматриваемая реакция происходит при высоких давлениях, где работает модель STM (2.1.3). Поэтому для описания регулярных колебаний в Pd-цеолитном катализаторе в [51] была предложена гибридная модель, имеющая некоторые общие черты с моделями (2.1.3) и (2.1.5). Ниже приводится подробное описание модели [51].

Предполагается, что реакция окисления СО на Pd происходит по кинетической схеме (2.1.1), (2.1.2). Однако стадии (2.1.2) означают не окисление-восстановление поверхностного слоя Pd, как в модели STM, а переход поверхностного кислорода в первый подповерхностный слой и реакцию подповерхностного кислорода О* с поверхностными молекулами СО. В отличие от модели BI здесь пропущена стадия выхода кислорода О* на поверхность. Так же, как в модели BI, предполагается, что кислород О* влияет на адсорбцию Ог (через е^_а2) и не влияет на адсорбцию СО. Однако допускается влияние поверхностного кислорода О на адсорбцию СО. Кроме того, так же, как в модели STM, константы скоростей не зависят от заполнений. В результате этих предположений получается следующая система уравнений, описывающая динамику реакции:

х = Х(х,у, z) = Pco^i(l -x — Sy) — к-хх - к_гху - k₅xz,

у = Y{x,y,z)=P₀₂k₂{l-x-y)²e~^az-hxy-hy(l-z), (2.1.7)

і = Z(x,y,z) = k₄y(l - z)- k₅xz;

где x(t), y(t) - поверхностные концентрации (покрытия) CO и О, соответственно; z{t) - концентрация растворенного кислорода О* в первом подповерхностном слое. Константы скоростей &, зависят от температуры Г по формуле Аррениуса. Используемые значения префакторов и энергий активации содержатся в следующей таблице.

Параметр а=10. Новый по сравнению с двумя предыдущими моделями подгоночный параметр S будет обсуждаться ниже; пока отметим, что его значение было выбрано равным 0.6. Давление кислорода считалось фиксированным и равным P_Qj = 160 Торр. Давление моноксида углерода Рсо рассматривалось как внешний бифуркационный параметр модели.

] Т 1—1'Г |"Г I I I ] і—I Г I ) I I I 1 | і—I Г 1 -J—I I I I \

0.0 0.5 1.0 1.5 2.0 2.5 3.0

0.0 0.5 1.0 1.5 2.0 2.5 3.0

Рсо Р"0РР]

Рис. 2.1: Вверху - стационарное решение х₀,2/о-^о системы (2.1.7) в зависимости от давления Рсо- Устойчивое решение показано сплошными линиями, неустойчивое -штриховыми. Пунктирными линиями отмечена амплитуда колебаний соответствующих переменных на устойчивом предельном цикле. Сверхкритическая бифуркация Андронова-Хопфа происходит в точках Рнх « 0.670 Торр и Р#2 % 2.623 Торр. Внизу - скорость реакции на стационарном решении.

На Рис.2.1 показана бифуркационная диаграмма стационарного решения системы (2.1.7), построенная методом продолжения по параметру. Из этой диаграммы

следует, что в области 0 < x,y,z < 1,. х+у < 1 фазового пространства, в которой фазовые переменные x,y,z имеют физический смысл, система (2.1.7) имеет в зависимости от Рсо либо одно устойчивое стационарное решение, либо одно неустойчивое стационарное решение вместе с устойчивым предельным циклом. На верхней части Рис. 2.1 компоненты стационарного решения Xo,yo,z₀ показаны разными цветами: х₀ (СО) - красным, уо (О) - синим и z₀ (О*) - зеленым цветом. Устойчивое решение показано сплошными линиями, неустойчивое - штриховыми. Стационарное решение теряет устойчивость при увеличении давления Рсо в точке Р_со = Р_Я1 « 0.670 Торр и при уменьшении давления Рсо в точке Р_Со = Ряг ^ 2.623 Торр. В этих точках при изменении давления в указанных направлениях происходит сверхкритическая бифуркация Андронова-Хопфа, при которой устойчивое стационарное решение становится неустойчивым и рождается устойчивый предельный цикл. Устойчивый предельный цикл существует на интервале Р_Н1 < Р_со < Ряг-Амплитуда колебаний фазовых переменных на предельном цикле показана на Рис. 2.1 пунктирными линиями. В нижней части Рис. 2.1 показана зависимость от Рсо скорости реакции R = к₃ху + k₅xz (в единицах монослой/'секупда) на стационарном решении Хо,yo,ZQ. Пунктирными линиями здесь также отмечена амплитуда колебаний скорости реакции на предельном цикле.

Подгоночный параметр модели 8 характеризует степень влияния адсорбированного кислорода на скорость адсорбции СО. При 5 = 0 выражение для скорости адсорбции СО получается таким же, как в модели BI (2.1.5). В нашей модели параметр 8 принимает некоторое промежуточное значение (8 = 0.6). В этой связи интересно проследить влияние параметра 8 на качественное поведение решений системы (2.1.7).

На Рис. 2.2 представлена фазовая диаграмма системы (2.1.7) на плоскости параметров (Рсо, 5) при 8 < 1.1, полученная при указанных выше значениях остальных параметров модели. На диаграмме серым цветом выделена область существования устойчивого предельного цикла. Граница области, - линия бифуркаций Андронова-Хопфа, построена с помощью программы продолжения по параметру бифуркаций коразмерности 1. Вне области колебаний система имеет единственное устойчивое

ИЗ

Р_со [Торр]

Рис. 2.2: Фазовая диаграмма системы (2.1.7) на плоскости параметров Рсо-& при 6 < 1.1. Серым цветом закрашена область автоколебаний.

стационарное решение в физической области фазового пространства. Из приведенной диаграммы следует существование минимального значения д~, ниже которого автоколебаний нет (при фиксированных значениях других параметров). С увеличением 6 область колебаний расширяется и изменяется форма колебаний.

2.2 Автоколебания в точечной модели реакции

В этом разделе рассматриваются особенности предельного цикла системы (2.1.7), связанные с релаксационным характером колебаний. Показано, что в системе можно выделить быстрые и медленные переменные. Число уравнений в системе можно уменьшить до двух с сохранением всех динамических свойств. При 8 « 1 вблизи правой точки бифуркации Андронова-Хопфа наблюдается так называемая "ложная бифуркация" предельного цикла. Изложение этого раздела следует в основном работе [58].

2.2.1 Быстрые и медленные переменные

В кинетической схеме рассматриваемой реакции предполагалось, что четвертая и пятая стадии (диффузия в объем и сегрегация кислорода) являются самыми медленными. Поэтому скорость изменения концентрации z(t) растворенного кислорода О* в модели (2.1.7) намного меньше скоростей изменения переменных x(t) и y(t).

2.5Торр

Р_го = 0.7Торр

2 4 6 8 10 Время [с]

2 4 6 8 10 Время [с]

2 4 6 8 10 Время [с]

Рис. 2.3: Установившиеся колебания в системе (2.1.7) при 8 = 0.6 и трех значений давления Рсо- (^а)>(^с)>(^е) Решение системы: красная линия - x(t), синяя - y(t), зеленая - z(t). (b),(d),(f) Скорости изменения переменных системы (правые части уравнений) x(t),y(t),z(t) показаны тем же цветом, что и сами переменные.

На Рис. 2.3а,с,е показаны графики зависимости от времени решений системы (2.1.7) в режиме установившихся колебаний (на предельном цикле) при S — 0.6 для трех значений давления Рсо'. (а) вблизи левой границы интервала автоколебаний при Рсо = 0.7, (с) в центре интервала автоколебаний при Рсо = 1-65, (е) вблизи правой границы при Рсо = 2.5. На Рис. 2.3b,d,f под графиками решений приведены графики зависимости от времени скоростей изменения переменных модели, т.е. правых частей системы (2.1.7). Скорость изменения каждой переменной показана

линией того же цвета, что и соответствующая переменная на верхнем графике. Из этих графиков видно, что пиковая скорость изменения z много меньше пиковых скоростей изменения х и у. Для приведенных на рисунке графиков максимальные значения абсолютных величин скоростей следующие: (Ь) 0.409, 0.245, 0.008; (d) 1.480,1.114, 0.029; (f) 0.698, 0.445, 0.018. Таким образом, во всех трех случаях (также как на всем интервале существования автоколебаний) пиковые значения скорости изменения z на порядок меньше пиковых скоростей изменения ж и у, которые имеют один и тот же порядок величины. Следовательно, для рассматриваемых значений параметров модели переменные х и у можно считать быстрыми, а переменную z -медленной.

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

X(x,y,z) = P_Co&i(l - х-5у) - k_iX - к₃ху - k_bxz = 0,

Y{x,y,z) = PoM~^az{^ -x-yf- k₃xy - hy{l - z) = 0. (2.2.8)

Эту линию легко можно построить с помощью программы продолжения по параметру стационарных решений нелинейных систем. Для этого в первых двух уравнениях системы (2.1.7) переменную z надо рассматривать, как внешний бифуркационный параметр. На Рис 2.4 показаны трехмерные проекции линии медленных движений S и траектории предельного цикла для тех же трех значений давления Рсо, нто и на Рис 2.3.

Линия медленных движений S имеет характерную форму кубической параболы с двумя локальными экстремумами. Как уже отмечалось, эту линию можно рассматривать как график стационарных решений системы двух уравнений для быстрых переменных х и у в зависимости от параметра z. В этом качестве (оливковый) участок линии S между экстремумами соответствует неустойчивым стационарным решениям, а две фиолетовых ветви - устойчивым стационарным решениям системы

(a)P_ro = 0.7 (b)P =1.65

(с) P_m = 2.5

Рис. 2.4: Трехмерные проекции линии медленных движений и траектории предельного цикла в системе (2.1.7) при S = 0.6 и трех значений давления Рсо- Черная точка показывает положение стационарного решения.

уравнения для быстрых переменных. Неустойчивое стационарное решение системы (2.1.7) показано на кривой черной точкой. Эта точка есть точка пересечения линии S и поверхности Z(x,y,z) = 0. Из Рис 2.4 видно, что при изменении Рсо форма кривой S качественно не изменяется.

Бифуркация Андронова-Хопфа рождения предельного цикла в системе (2.1.7) происходит тогда, когда при изменении Рсо стационарное решение проходит точку экстремума кривой S и переходит с устойчивой части S на неустойчивую. Родившийся предельный цикл в фазовом пространстве расположен вблизи точки экстремума S. При удалении стационарного решения от точки экстремума по неустойчивой части S диаметр предельного цикла постепенно увеличивается (Рис. 2.4а и 2.4с) и цикл принимает форму, характерную для релаксационных циклов (Рис. 2.4Ь). Как правило, релаксационный цикл имеет два быстрых и два медленных участка. Медленные участки траектории цикла расположены вдоль устойчивых ветвей кривой

S, где фазовая точка движется в сторону экстремума S. После прохождения экстремума происходит "срыв", - фазовая точка быстро падает в окрестность другой устойчивой ветви S.

Отметим, что в рассматриваемом случае при 6 = 0.6 увеличение диаметра предельного цикла при удалении стационарного решения от точки экстремума S происходит достаточно плавно. Однако в других случаях он может иметь "взрывной" характер. Ниже будет рассматриваться явление "взрывного" роста диаметра цикла при 6 = 1.

4.0-

3.5-

І 3.0-Q.

2.5-

2.0-1 г 1 1 1 1 1 1 1 1

0.5 1.0 1.5 2.0 2.5 3.0

PcoFopp]

Рис. 2.5: Зависимость периода автоколебаний в системе (2.1.7) от давления Рсо-

Сделаем также еще одно замечание, которое понадобится в следующей главе. В рассматриваемой системе в области автоколебаний при увеличении параметра Рсо длина траектории предельного цикла увеличивается на большей части интервала автоколебаний, за исключением небольшого участка вблизи Ряг> и вследствие этого увеличивается период предельного цикла. На Рис. 2.5 показан график зависимости периода колебаний от Рсо при S = 0.6 и указанных выше значениях других параметров.

2.2.2 "Квазистационарное'' приб лижение

В системе с быстрыми и медленными переменными при движении в ограниченной области фазовая точка не может надолго отклоняться от линии медленных движений S. При значительном отклонении от линии S фазовая точка с большой скоростью устремляется к той устойчивой ветви S, в области притяжения которой она находится. В рассматриваемом случае имеются две быстрых переменных и линия S есть линия пересечения поверхностей нулевой скорости быстрых переменных. Поскольку фазовая точка при движении по циклу большую часть времени находится вблизи S, можно ожидать, что динамическое поведение системы изменится незначительно, если считать, что фазовая точка всегда лежит на поверхности нулевой скорости одной из быстрых переменных. Т.е. можно ожидать, что динамическое поведение решения системы дифференциальных уравнений (2.1.7) не будет качественно отличаться от динамического поведения решения дифференциально-алгебраической системы, полученной из (2.1.7) заменой уравнения для одной из быстрых переменных на алгебраическое уравнение, приравнивающее нулю правую часть этого уравнения. В данном случае удобно в системе (2.1.7) заменить уравнение для у на алгебраическое, так как оно не содержит бифуркационный параметр Рсо- В результате получится система уравнений, которую можно назвать "квазистационарным" приближением для системы (2.1.7):

х = Х(х, у, z) — Pcoh (1 - х - 6у) - k-xx - k₃xy - k₅xz,
z = Z(x,y,z) = к^у(1 — z) — k₅xz, (2.2.9)

Y(x, y, z) = PoMl ~^x~ УТ*~^а* - ^к*^хУ - Ml - *) = 0.

Решая алгебраическое уравнение относительно у, находим

у = l-x + q± vV + 2g(l - х),

е^аг g=_op , (k₃x + k₄{l-z)).

В физической области фазового пространства лежит решение со знаком '—' перед

корнем. Поэтому положим

у = у (ж, z) = l-x + q- yjq² + 2q(l-x), (2.2.10)

и подставим это выражение в первые два уравнения системы (2.2.9). Получим систему из двух уравнений

х = X{x,y(x,z),z),

z = Z(x,y(x,z),z). (2.2.11)

(a) (b)

1 1-От

Р [Торр] Время [с]

Рис. 2.6: (а) Бифуркационная диаграмма стационарного решения системы (2.2.11) в зависимости от параметра Рсо- Сплошные линии - устойчивое решение, штриховые - неустойчивое. Черными значками отмечены точки сверхкритической бифуркации Андронова-Хопфа P#i ^w 0.671 и Р#2 ^w 2.624. Пунктирными линиями показана амплитуда колебаний переменных модели. (Ь) Установившиеся колебания в системе (2.2.11) при Рсо = 1-65 Торр. Красная линия - x(t), зеленая - z(t). Пунктирная синяя линия - y(t), вычисленное по формуле (2.2.10).

Динамическое поведение решений двумерной системы (2.2.11) полностью аналогично поведению решений трехмерной системы (2.1.7). На Рис. 2.6а показана бифуркационная диаграмма стационарного решения (xo,z₀) системы (2.2.11) (красная линия - ж₀(Рсо)₅ зеленая - z₀(P_Co)), которая и качественно и количественно

практически не отличается от бифуркационной диаграммы стационарного решения системы (2.1.7), показанной на Рис. 2.1а. В частности, положение точек бифуркации Андронова-Хопфа в обеих системах практически совпадают: P_Hi = 0.670, Р_Н2 - 2.623 в (2.1.7) и Р_Н1 = 0.671 и Р_Н2 = 2.624 в (2.2.11). Сравнение Рис. 2.6Ь и Рис. 2.3с показывает, что и графики установившихся колебаний переменных в обеих системах также очень похожи. Однако при фиксированном значении Рсо период колебаний в двумерной системе несколько меньше периода колебаний в трехмерной системе.

Таким образом, регулярные колебания в реакции окисления СО на поверхности Pd можно моделировать с помощью системы из двух уравнений (2.2.11), которая является "квазистационарным" приближением системы (2.1.7). Отметим, что с вычислительной точки зрения выигрыш от замены исходной системы (2.1.7) на систему меньшей размерности (2.2.11) не велик. Вычислительные затраты на решение обеих систем примерно одинаковы. Однако переход к системе (2.2.11) облегчает качественный анализ модели и может дать существенный выигрыш при решении пространственно распределенной модели (см. Главу 3).

2.2.3 Циклы-утки в точечной модели

Динамические системы с быстрыми и медленными переменными по своим свойствам аналогичны сингулярно возмущенным системам обыкновенных дифференциальных уравнений. Сингулярно возмущенная (автономная) система из двух уравнений имеет вид:

ех = /(ж, у), У = 9(х,у);

где є - малый параметр. При малых є, х можно считать быстрой переменной, а у - медленной. Особенность сингулярного возмущения состоит в том, что при є -* 0 система превращается в дифференциально-алгебраическую систему. Именно в таких системах (в частности, в системе уравнений Ван дер Поля) впервые были обнаружены и изучены так называемые траектории-утки и циклы-утки.

На Рис 2.7 показаны предельные циклы в системе Ван дер Поля

ех = у — (х³/3 — х),
у = а-х; (2.2.12)

при є = 0.05 и разных значений бифуркационного параметра а. В этой системе автоколебания существуют при — 1 < а < 1. В точках а = ±1 происходит надкритическая бифуркация Андронова-Хопфа. "Нормальный" релаксационный цикл (при а = 0) в системе показан на рисунке малиновыми кружками. Медленные участки цикла расположены вдоль устойчивых ветвей кривой медленных движений S, быстрые участки соответствуют переходам с одной ветви на другую. Так устроены все циклы в системе (2.2.12), если а достаточно далеко от точек ±1.

Рис. 2.7: Предельные циклы в системе Ван дер Поля при є = 0.05. Показан релаксационный цикл при а = 0 и циклы-утки при а « 1.

Однако вблизи точек бифуркации траектории предельных циклов иные. На Рис 2.7 показаны также циклы при а близком к 1. После рождения в точке а = 1 при дальнейшем уменьшении параметра а предельный цикл располагается в окрестности минимума кривой S (точка (1,-2/3)). Характерная особенность эти циклов состоит в том, что на траектории цикла имеются два следующих друг за другом участка, из которых первый расположен вдоль устойчивой ветви S, а второй

вдоль неустойчивой части S. Траектории, имеющие такую особенность, называют траекториями-утками, в частности, циклами-утками. Циклы-утки своей формой напоминают летящую утку. Подробные сведения по теории решений-уток сингулярно возмущенных систем представлен в обзоре [59] и в книге [60].

В соответствии с принятой в теории циклов-уток терминологией, участки траектории цикла, лежащие вблизи точки пересечения цикла с неустойчивой и устойчивой ветвями S, называют клювом и хвостом утки. При движении по параметру внутрь области автоколебаний клюв поднимается по неустойчивой части S и диаметр цикла увеличивается. После прохождения клювом точки второго экстремума S форма цикла изменяется. Появляется второй медленный участок, в этом случае говорят, что у утки вырастает голова. Как правило, все фазы изменения цикла от малого (гармонического) цикла, лежащего в окрестности экстремума S, до большого (релаксационного) цикла происходят в очень узком интервале значений бифуркационного параметра. В частности, на Рис. 2.7 циклу d соответствует а = 0.994491, а циклу cj - а = 0.9944909294, т.е. разница в значениях параметра меньше 10^_б. При уменьшении є циклы-утки теснее прижимаются к кривой медленных движений S и интервал значений параметра а, на котором существуют циклы-утки уменьшается.

Траектории циклов уток проходят внутри очень узкой полосы вдоль кривой S, при уменьшении є ширина этой полосы уменьшается и, соответственно, возрастают требования к точности решения системы уравнений для построения циклов-уток. В частности, для системы (2.2.12) уже при є = 0.01 ошибки округления, обусловленные конечной длиной представления чисел в компьютере, не позволяют построить все циклы-утки. В [59] необходимое число десятичных знаков в представлении вещественных чисел для построения всех циклов-уток в системе (2.2.12) при фиксированном є оценивается величиной е¹^. Без применения специальных программ, увеличивающих длину представления чисел, можно построить только циклы-утки малого диаметра (почти гармонические) и циклы большого диаметра (почти релаксационные). Связанное с циклами-утками явление чрезвычайно быстрого увеличения диаметра предельного цикла при малом изменении управляющего параметра иногда называется "ложной бифуркацией" предельного цикла [61], однако никакой

бифуркации здесь не происходит.

В "квазистационарном" варианте (2.2.11) модели (2.1.7) также есть все условия для существования циклов-уток, которые особенно ярко выражены вблизи правой точки бифуркации Андронова-Хопфа Р#2- На Рис. 2.8 показаны некоторые циклы-утки в системе (2.2.11) при S = 0.7. Циклы построены для значений параметра Рсо от 3.0531672 (самый большой цикл) до 3.0531690 (самый маленький цикл) с шагом 2 х 10~⁷. Пунктирной линией изображен цикл для Рсо = 3.0531681. Отметим, что в своей нижней части траектории разных предельных циклов в масштабах рисунка неразличимы. Изменение цикла от малого до большого происходят на интервале значений Рсо длиной всего 2 х Ю^-6. В исходной трехмерной системе (2.1.7) вблизи

0.175-

0.170-

Z 0.165-

0.160-

0.155-

0.5 0.6 0.7 0.8 0.9

Рис. 2.8: Циклы-утки в системе (2.2.11) при 8 = 0.7.

границ интервала автоколебаний также наблюдаются циклы-утки. Причем чрезвычайно быстрое изменение диаметра цикла с изменением параметра Рсо происходит только вблизи правой точки бифуркации Риг- В следующей главе будет показано, что в паре связанных химических осцилляторов, каждый из которых представлен точечной моделью (2.1.7), в окрестности точки Ря2 наблюдается детерминированный хаос. Странный аттрактор системы связанных осцилляторов представляет собой апериодическую траекторию-утку.

3 Основные результаты главы

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

Показано, что модель принадлежит к классу моделей с быстрыми и медленным переменными. Из трех переменных модели две являются быстрыми и одна - медленной. Показано, что динамика трехмерной модели полностью аналогична динамике двумерной модели, в которой одна из быстрых переменных выражается через остальные на поверхности нулевой скорости ([58]). Описаны циклы-утки в двумерной и трехмерной моделях.

3 Модели колебаний в реакции окисления СО в зерне палладиевого цеолитного катализатора

3.1 Распределенные модели гетерогенного катализа

Рассмотренные в Главе 2 макроскопические модели реакции окисления СО на пал-ладиевом цеолитном катализаторе не дают полного описания экспериментально наблюдаемой динамики реакции. Эти модели более или менее адекватно описывают стационарные режимы и периодические колебания скорости реакции. Однако в экспериментах наблюдаются более сложные "турбулентные" режимы. В конечномерных диссипативных динамических системах "турбулентный" режим представляет собой детерминированный хаос ([10], [34], [20], [36]). В связи с этим, первой задачей данной главы является задача построения и изучения конечномерной модели реакции в зерне цеолита, которая имеет решения типа детерминированного хаоса. По-видимому, наиболее простой способ получения сложных динамических режимов в простых системах состоит в объединении нескольких простых систем в одну составную систему с помощью подходящих связей. В химических системах простые модели, входящие в составную модель, описывают реакции в отдельных (пространственно разнесенных) частях реактора, которые связаны процессами массообмена и энергообмена. Таким образом, составную модель больше нельзя рассматривать как точечную, - она является простейшим примером распределенной модели. В данной работе рассматривается реакция, которая протекает в изотермических условиях, поэтому учитывается только массообмен посредством диффузии. Примеры и анализ условий существования сложной динамики в системах, составленных из пары достаточно простых моделей, описывающих релаксационные колебания в химических, биологических и других системах, приводятся, в частности, в [35].

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

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

X = R(X,P),

где X - вектор динамических переменных, а Р - вектор внешних параметров, R описывает кинетику реакции. В общем виде дискретную распределенную модель можно записать в виде системы уравнений

X_k = R_k(X_k,P) + Y^D_kj(Xj-X_k), k = l,...,N;

Р = а(Р_ехі-Р)-'*Гд_к(Х_к,Р).

k-i

Для определенности можно считать, что X представляет поверхностные концентрации реагентов, а Р - парциальные давления реагентов в газовой фазе. Сумма в

уравнении для Xk представляет диффузионную связь между осцилляторами. Если Dkj Ф 0 только для ближайших соседей, то связь локальна, в противном случае - нелокальная. Уравнение для Р представляет глобальную связь. Изменение Р обусловлено обменом с внешним резервуаром реагентов P_ext и обменом веществом между газовой фазой и поверхностью катализатора посредством адсорбции и десорбции, выраженным суммой членов ди. В непрерывных моделях распределенных систем уравнения для Xk перейдут в систему уравнений в частных производных типа реакция-диффузия, а в уравнении для Р сумма заменится интегралом по поверхности катализатора.

Для описания поверхностных явлений в гетерогенном катализе часто не учитывают связь по газовой фазе и рассматривают модели вида

Щ±$- = R(X,P) + V(DVX(r,t)),

где г - вектор пространственных координат, D - матрица диффузии [26], [19], [62]. С помощью моделей такого типа удается описать множество интересных явлений в гетерогенных и гомогенных системах, начиная от знаменитых структур Тьюринга [24] и до сложных анизотропных пространственно-временных структур [63].

Влияние глобальной связи по газовой фазе на динамику реакции окисления СО на Pt(llO) изучалась в [64]. Показано, что роль глобальной связи по газовой фазе не сводится только к количественным изменениям. Глобальная связь может приводить к качественным изменениям динамического поведения. Численное исследование математических моделей реакций с глобальной связью обнаруживает существование некоторого порогового значения силы (интенсивности) глобальной связи, выше которого локальные осцилляторы полностью синхронизованы и глобальная скорость реакции проявляет регулярные колебания, т.е. в модели наблюдаются пространственно однородные колебания скорости реакции. При уменьшении силы глобальной связи ниже этого порогового значения, в модели наблюдается переход от однородных регулярных колебаний к нерегулярным колебаниям малой амплитуды. Динамика модели с глобальной связью ниже порогового значения зависит от величины локальной связи между осцилляторами. Если локальная связь отсутствует, то

наблюдается кластеризация локальных осцилляторов. В каждом кластере осцилляторы колеблются в фазе, а между колебаниями в разных кластерах устанавливается фиксированная разность фаз. При наличии поверхностной диффузии кластеризация может наблюдаться по крайней мере вблизи бифуркации Андронова-Хопфа. Влияние глобальной связи на динамику реакции изучалась также в абстрактной модели Гинзбурга-Ландау [21, 65].

В случае цеолитного катализатора диффузия в газовой фазе происходит в порах цеолита с весьма малыми коэффициентами диффузии ~ 10^-4-10~⁵ см²/с [66]. Поэтому в данном случае характерные времена диффузии в зерне цеолита сравнимы с характерными временами реакции и связь по газовой фазе не может быть глобальной. Система с локальной связью по газовой фазе может быть представлена в виде

= V(D₃VX(r,t)) + R(X,P), = V(A,VP(r,*))+G(X,P),

dX(r,t)

Ы dP(r,t)

где D_a и D_g матрицы коэффициентов поверхностной диффузии адатомов и диффузии в газовой фазе. Системы такого типа рассматривались в работах [67, 68] в качестве модели абстрактной реакции с переменными по пространству параметрами.

Пусть D_B = 0 и пусть задано граничное условие для Р на поверхности сферы радиуса R, P(r,t)\_r-R = Р. Тогда в стационарном режиме реакции в зерне установится неоднородное распределение концентрации реагентов, которые будут монотонно убывать от поверхности к центру зерна. В колебательном режиме сила связи между локальными осцилляторами будет определяться коэффициентом диффузии D_g. При D_g —> оо давление во всем объеме зерна будет равно Р и осцилляторы будут независимы. При уменьшении D_g влияние соседних осцилляторов друг на друга будет возрастать. При D_g -> 0 реакция будет происходить только у поверхности зерна. Как будет показано ниже, при некоторых средних значениях коэффициента диффузии Dg (точнее, характерной скорости диффузии D_g/R²) динамика реакции будет достаточно сложной.

64 3.2 Строение цеолитного катализатора

Рассматриваемые в данной главе модели предназначены для описания динамики реакции окисления СО на палладиевом цеолитном катализаторе и опираются на соответствующие экспериментальные данные. Поскольку в этой главе представлена модель реакции в зерне такого катализатора, то сначала кратко опишем особенности этой каталитической системы. Технология изготовления палладиевого цеолитного катализатора описана в работе [69] (об атомной структуре цеолитов и их роли в катализе см., например, обзор [57]). В катализаторе цеолит используется как матрица, в которую встраиваются наночастицы (кластеры) палладия. Эти кластеры являются активным элементом катализатора. Молекулы реагентов СО и Ог проникают внутрь цеолита через поры, диффундируют по каналам и адсорбируются на поверхность кластеров палладия. Реакция происходит именно на поверхности кластеров Pd.

Цеолиты являются природными алюмосиликатными минералами и известны уже почти 250 лет. В настоящее время цеолиты представляют большой интерес для катализа, однако природные формы цеолитов имеют существенные недостатки. Они почти всегда содержат нежелательные примеси и их химический состав может меняться от образца к образцу и даже в пределах одного образца. Существенную роль в катализе цеолиты стали играть с начала пятидесятых годов, когда научились изготавливать синтетические цеолиты. Например, важнейшей областью применения синтетического faujasite Y, на основе которого изготавливается палла-диевый цеолитный катализатор, является нефтехимическая промышленность.

Атомарная кристаллическая структура цеолита faujasite показана на Рис. 3.1. Элементарными единицами, из которых построены цеолиты, являются тетраэдральные структуры Si0₄ и АЮ₄ (атом Si или А1 лежит в центре тетраэдра, а атомы О-в его вершинах). Тетраэдры могут соединяться друг с другом вершинами через общий атом кислорода. В результате такого соединения может получиться устойчивая структура, которая служит элементарной ячейкой второго уровня и называется клеткой. В свою очередь клетки образуют кристаллическую решетку

Рис. З.Г. Кристаллическая структура цеолита Faujasite. Вверху слева - элементарная ячейка (Sodalite Unit), - клетка, состоящая из 24-х тетраэдров Al(Si)0₄ с общими атомами кислорода. Вверху слева - кольцо из шести клеток, - основной элемент системы пор цеолита. Внизу - суперклетка, образованная четырьмя кольцами.

- макромолекулу цеолита или матрицу. На Рис. 3.1 показаны основные элементы структура матрицы цеолита faujasite. Клетка этого цеолита состоит из 24-х тетраэдров, соединенных общими вершинами и образующими кубо-октаэдр, который называют ячейкой содалита (sodalite unit). Кубо-октаэдр имеет шесть квадратных граней и восемь шестиугольных. В цеолите faujasite клетки соединяются между собой шестиугольными гранями, образуя тетраэдральную решетку. Основой пористой структуры цеолита является кольцо, состоящее из шести клеток. Это кольцо не плоское, поэтому из таких колец можно составлять пространственные структуры. В частности, на Рис. 3.1 показана суперклетка цеолита, состоящая из четырех колец,

которые молено связать с гранями тетраэдра. Внутри этого тетраэдра находится сферическая полость диаметром 1.3 нм, соединенная четырьмя окнами диаметром 0.74 нм с соседними полостями.

В процессе изготовления катализатора кластеры палладия кристаллизуются внутри цеолита, причем в результате получаются кластеры примерно одного размера, который можно достаточно точно регулировать изменяя условия кристаллизации. Размер кластеров может значительно превышать размер суперклетки цеолита, при этом кластеры локально разрушают матрицу цеолита. Приведем некоторые количественные данные из [70], характеризующие строение катализатора.

Общая масса палладия (загрузка палладия) составляет несколько процентов от массы цеолита. При загрузке палладия, равной 5%, один атом палладия приходится примерно на 10 клеток цеолита. Средний размер кристаллитов (зерен) цеолита, из которых состоит катализатор, равен 5 мкм. Предполагая, что зерно цеолита имеет сферическую форму и учитывая, что плотность цеолита равна 1.9 г/см³, находим массу зерна равной примерно 1.2-10^-10 г. При 5-ти процентной загрузке палладия в зерне цеолита будет содержаться примерно 6-Ю^-12 г палладия. Если предположить, что частицы палладия имеют сферическую форму, то с учетом того, что удельный вес палладия равен ~12 г/см³, находим, что в одном зерне цеолита должно содержаться примерно 1.5-10⁷ частиц палладия диаметром 4 нм, или 9.5-10⁵ частиц диаметром 10 нм. При этом в первом случае среднее расстояние между центрами частиц будет равно ~20 нм и во втором случае, примерно 50 нм. Площадь активной поверхности в зерне цеолита примерно равна 7.5 Ю^-6 см² в случае 4-х нанометровых частиц и 3 Ю^-6 см² в случае 10-ти нанометровых частиц. (Площадь активной поверхности на единицу объема равна, соответственно, 1.15 10⁵ см^-1 и 4.6-10⁴см-¹.)

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

микропористых материалах (диаметр пор < 2 нм) в настоящее время не существует. В литературе обсуждаются различные подходы к моделированию диффузии в пористых средах. В частности, на основе модели многокомпонентного решеточного газа с множественной адсорбцией в клетках цеолита выводятся нелинейные уравнения диффузии [71], которые хорошо описывают процесс при низких концентрациях. Другой подход к аналитическому описанию диффузии газовой смеси в цеолите основан на теории Максвелла-Стефана [38].

3.3 Двухслойная модель реакции в зерне цеолита

Для демонстрации влияния диффузии газа в порах цеолита на динамику реакции была предложена следующая простая модель [51, 72]. Как уже говорилось, для простоты зерно (кристаллит) цеолита можно представить в виде шара радиуса Д.. Будем считать, что все частицы палладия имеют одинаковую форму и размер и равномерно распределены по объему зерна. В ходе реакции газовая смесь, содержащая СО и Ог проникает по порам внутрь зерна, молекулы СО и 0₂ адсорбируются на поверхность частиц Pd, где происходит каталитическая реакция 2С0 + 0₂ —> 2С0₂. В дальнейшем мы не рассматриваем отдельные частицы Pd, а считаем, что активная поверхность катализатора равномерно распределена по внутренней поверхности пор цеолита. При таких предположениях общая модель реакции в зерне должна быть моделью типа реакция-диффузия, которая описывается системой дифференциальных уравнений в частных производных. К этой системе мы обратимся во второй части этой главы, а здесь рассмотрим дискретный вариант модели реакции.

Модель регулярных колебаний в Pd-цеолитном катализаторе

Каждому динамическому режиму соответствует определенный вид спектра показателей Ляпунова. Если система имеет только стационарные решения, то все показатели отрицательны и спектр имеет вид (—, —,..., —). Если есть периодическое решение, то старший показатель равен 0, а остальные отрицательны, спектр - (О,—,...,—). В случае квазипериодического решения число нулевых показателей равно размерности тора, спектр - (0,..., О, —,..., —). Для хаотического решения один или несколько показателей положительны, один показатель нулевой и остальные отрицательны. Например, для системы четырех уравнений возможные хаотические спектры: (+, 0, -, —) - хаос и (+, +, 0, —) - гиперхаос.

В программу решения систем ОДУ встроен модуль оценки показателей Ляпунова по методу Вольфа [46]. Исходя из определения показателей для их вычисления нужно следить за изменением инфинитезимального объема в фазовом пространстве в ходе эволюции системы. Однако на практике очень трудно определить критерий инфинитезимальности, поэтому удобнее работать в касательном пространстве с системой в вариациях (1.4.9).

Для оценки показателей Ляпунова одновременно с решением нелинейной системы (1.4.2) решается линейная система (1.4.9) для п векторов начальных данных, которые представляют собой произвольно ориентированный ортонормированный репер {vi,v2,...,vn} в фазовом пространстве. Для определенности будем считать, что г = е(к\ С течением времени длины векторов, в направлении которых показатели положительны, будут увеличиваться. Кроме того, каждый вектор стремится развернуться в направлении наибольшего роста длины. Поэтому со временем углы между векторами уменьшаются. Чтобы предотвратить "коллапс"репера к нему периодически применяют процедуру ортогонализации Грамма-Шмидта.

Пусть {vi,V2,...,vn} - репер в некоторый момент времени t. Ортогонализация векторов Vk выполняется по формулам: где (,)- скалярное произведение в Еп. При этом направление первого вектора v[ не изменяется, а каждый вектор v k поворачивается оставаясь в линейной оболочке векторов Vi,г»2, , vk так, чтобы стать перпендикулярным векторам v[, v 2,..., v k_1. Теоретически, частота выполнения ортогонализации не имеет значения, однако на практике рекомендуется выполнять ее не менее одного раза за средний период колебаний.

Предположим, что ортогонализация репера проводится в моменты времени tn, п = 1,2,..., N. Тогда оценки показателей Ляпунова вычисляются по формулам

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

Странный аттрактор в диссипативной системе может существовать в некоторой области значений внешних параметров. Предположим, что имеется только один управляющий параметр. При изменении параметра переход к хаосу может происходить только по вполне определенным сценариям. Таких сценариев перехода к хаосу известно три [47]. В литературе наиболее часто встречается переход к хаосу по сценарию Фейгенбаума [48]. В этом случае странный аттрактор образуется в результате бесконечной последовательности бифуркаций удвоения периода предельного цикла. Причем расстояния по параметру между последовательными бифуркациями убывают в геометрической прогресс, поэтому все бифуркации происходят на конечном интервале значений параметра. Вторым сценарием формирования странного аттрактора является сценарий Рюэля-Такенса. В этом сценарии процесс формирование странного аттрактора включает в себя несколько вторичных бифуркаций Хопфа с последующей трансформацией многомерного тора в странный аттрактор. Таким образом, первые два сценария связаны со сверхкритическими бифуркациями предельного цикла. Подкритические бифуркации связаны с третьим сценарием - перемежаемостью. Перемежаемость бывает нескольких типов, из которых наиболее известны типы I, II и III [49]. Эти типы перемежаемости связаны с подкри-тическими бифуркациями цикла, соответственно, седлоузловой, удвоения периода и вторичной бифуркацией Хопфа. При перемежаемости предельный цикл теряет устойчивость через одну из перечисленных бифуркаций и траектория динамической системы переходит на странный аттрактор. При перемежаемости странный аттрактор не формируется; он должен быть уже сформирован посредством одного из первых двух сценариев. В частности, перемежаемость часто встречается на границах интервалов периодичности, которые можно наблюдать в пространстве параметров в области хаоса.

Изменение структуры и типа аттрактора при изменении управляющего параметра можно наблюдать на бифуркационной диаграмме, которая строится следующим образом. При фиксированном значении параметра система (1.4.2) решается на достаточно большом интервале времени, так чтобы траектория приблизилась к аттрактору на достаточно малое расстояние. Затем эта система решается на заданном интервале времени. В процессе решения отмечаются характерные точки, например, локальные экстремумы определенных динамических переменных. Для стационарного состояния получится одна характерная точка, для предельных циклов число экстремумов связано с числом оборотов. В хаотическом режиме экстремумы решения заполняют некоторые интервала. После того как собрано достаточно информации для данного значения параметра, переходят к другому, близкому значению параметра и повторяют описанную процедуру. Результаты представляют в виде графика, на котором по оси абсцисс откладывают значения параметра, а по оси ординат - значения экстремумов. По построенной таким образом диаграмме можно достаточно легко выделить бифуркации периодического решения, области периодичности, квазипериодичности и хаоса, определить сценарии перехода к хаосу. Конечно, анализ бифуркационной диаграммы дает только предварительные результаты, которые должны быть подтверждены вычислением соответствующих количественных характеристик.

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

Для характеризации хаотического режима часто бывают полезными одномерные отображения [1]. В качестве таковых наиболее часто используются отображения следующего экстремума. Эти отображения представляют собой частный случай функции последования Пуанкаре, для которой роль секущей поверхности играет многообразие экстремумов фазовой переменной. В данной работе отображения следующего экстремума используются при анализе сценария перехода к хаосу через перемежаемость. В программу встроен модуль построения отображения га-го максимума (минимума) вместе с диаграммой Ламерея заданной динамической переменной.

Для анализа сложных колебаний в частотной области широко используется анализ Фурье. Для построения амплитудных и энергетических спектров используются стандартные программы быстрого преобразования Фурье временных последовательностей и построения спектров (периодограмм) [42]. По виду спектров Фурье можно определить тип колебаний. Для хаотических колебаний характерен непрерывный спектр с быстрым (экспоненциальным) спадом энергии гармоник при увеличении частоты.

Переход к хаосу через бифуркации удвоения периода

Рассмотренные в Главе 2 макроскопические модели реакции окисления СО на пал-ладиевом цеолитном катализаторе не дают полного описания экспериментально наблюдаемой динамики реакции. Эти модели более или менее адекватно описывают стационарные режимы и периодические колебания скорости реакции. Однако в экспериментах наблюдаются более сложные "турбулентные" режимы. В конечномерных диссипативных динамических системах "турбулентный" режим представляет собой детерминированный хаос ([10], [34], [20], [36]). В связи с этим, первой задачей данной главы является задача построения и изучения конечномерной модели реакции в зерне цеолита, которая имеет решения типа детерминированного хаоса. По-видимому, наиболее простой способ получения сложных динамических режимов в простых системах состоит в объединении нескольких простых систем в одну составную систему с помощью подходящих связей. В химических системах простые модели, входящие в составную модель, описывают реакции в отдельных (пространственно разнесенных) частях реактора, которые связаны процессами массообмена и энергообмена. Таким образом, составную модель больше нельзя рассматривать как точечную, - она является простейшим примером распределенной модели. В данной работе рассматривается реакция, которая протекает в изотермических условиях, поэтому учитывается только массообмен посредством диффузии. Примеры и анализ условий существования сложной динамики в системах, составленных из пары достаточно простых моделей, описывающих релаксационные колебания в химических, биологических и других системах, приводятся, в частности, в [35].

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

В случае автоколебательной реакции распределенную модель можно представит себе как совокупность локальных химических осцилляторов, взаимодействующих непосредственно или через общую среду. Пусть модель химического осциллятора имеет вид где X - вектор динамических переменных, а Р - вектор внешних параметров, R описывает кинетику реакции. В общем виде дискретную распределенную модель можно записать в виде системы уравненийДля определенности можно считать, что X представляет поверхностные концентрации реагентов, а Р - парциальные давления реагентов в газовой фазе. Сумма в уравнении для Xk представляет диффузионную связь между осцилляторами. Если Dkj Ф 0 только для ближайших соседей, то связь локальна, в противном случае - нелокальная. Уравнение для Р представляет глобальную связь. Изменение Р обусловлено обменом с внешним резервуаром реагентов Pext и обменом веществом между газовой фазой и поверхностью катализатора посредством адсорбции и десорбции, выраженным суммой членов ди. В непрерывных моделях распределенных систем уравнения для Xk перейдут в систему уравнений в частных производных типа реакция-диффузия, а в уравнении для Р сумма заменится интегралом по поверхности катализатора.

Для описания поверхностных явлений в гетерогенном катализе часто не учитывают связь по газовой фазе и рассматривают модели вида где г - вектор пространственных координат, D - матрица диффузии [26], [19], [62]. С помощью моделей такого типа удается описать множество интересных явлений в гетерогенных и гомогенных системах, начиная от знаменитых структур Тьюринга [24] и до сложных анизотропных пространственно-временных структур [63].

Влияние глобальной связи по газовой фазе на динамику реакции окисления СО на Pt(llO) изучалась в [64]. Показано, что роль глобальной связи по газовой фазе не сводится только к количественным изменениям. Глобальная связь может приводить к качественным изменениям динамического поведения. Численное исследование математических моделей реакций с глобальной связью обнаруживает существование некоторого порогового значения силы (интенсивности) глобальной связи, выше которого локальные осцилляторы полностью синхронизованы и глобальная скорость реакции проявляет регулярные колебания, т.е. в модели наблюдаются пространственно однородные колебания скорости реакции. При уменьшении силы глобальной связи ниже этого порогового значения, в модели наблюдается переход от однородных регулярных колебаний к нерегулярным колебаниям малой амплитуды. Динамика модели с глобальной связью ниже порогового значения зависит от величины локальной связи между осцилляторами. Если локальная связь отсутствует, то наблюдается кластеризация локальных осцилляторов. В каждом кластере осцилляторы колеблются в фазе, а между колебаниями в разных кластерах устанавливается фиксированная разность фаз. При наличии поверхностной диффузии кластеризация может наблюдаться по крайней мере вблизи бифуркации Андронова-Хопфа. Влияние глобальной связи на динамику реакции изучалась также в абстрактной модели Гинзбурга-Ландау [21, 65].

В случае цеолитного катализатора диффузия в газовой фазе происходит в порах цеолита с весьма малыми коэффициентами диффузии 10-4-10 5 см2/с [66]. Поэтому в данном случае характерные времена диффузии в зерне цеолита сравнимы с характерными временами реакции и связь по газовой фазе не может быть глобальной. Система с локальной связью по газовой фазе может быть представлена в виде где Da и Dg матрицы коэффициентов поверхностной диффузии адатомов и диффузии в газовой фазе. Системы такого типа рассматривались в работах [67, 68] в качестве модели абстрактной реакции с переменными по пространству параметрами.

Пусть DB = 0 и пусть задано граничное условие для Р на поверхности сферы радиуса R, P(r,t)\r-R = Р. Тогда в стационарном режиме реакции в зерне установится неоднородное распределение концентрации реагентов, которые будут монотонно убывать от поверхности к центру зерна. В колебательном режиме сила связи между локальными осцилляторами будет определяться коэффициентом диффузии Dg. При Dg — оо давление во всем объеме зерна будет равно Р и осцилляторы будут независимы. При уменьшении Dg влияние соседних осцилляторов друг на друга будет возрастать. При Dg - 0 реакция будет происходить только у поверхности зерна. Как будет показано ниже, при некоторых средних значениях коэффициента диффузии Dg (точнее, характерной скорости диффузии Dg/R2) динамика реакции будет достаточно сложной.

Формулировка стохастической модели реакции

В данном случае мы имеем дело с парой связанных осцилляторов, поэтому переход в правую полуплоскость первой пары собственных значений может быть связан с возбуждением первого осциллятора, а второй - с возбуждением второго осциллятора. Сравнение фазовых диаграмм отдельных осцилляторов (Рис. 3.3) и связанной пары (Рис. 3.5) показывает, что линии бифуркации Андронова-Хопфа /iieft и bright для пары связанных осцилляторов практически совпадают с границами области колебаний отдельных осцилляторов. Например, верхняя часть левой границы h\eft области колебаний пары осцилляторов совпадает с верхней частью левой границы int,i области колебаний внутреннего осциллятора до ее пересечения с левой границей hext,i области колебаний внешнего осциллятора. Ниже точки пересечения /iieft идет вдоль линии /&ext,i- Поскольку правые границы областей колебаний отдельных осцилляторов не пересекаются, правая граница rjght идет вдоль линии h-mtj2, которая лежит правее линии /iext,2- Линия Л2 перехода второй пары собственных значений через мнимую ось также проходит вблизи соответствующих участков границ области колебаний отдельных осцилляторов.

Линия h2 на фазовой диаграмме двухслойной модели делит область колебаний на четыре подобласти, помеченные римскими цифрами I, II, III и IV. Сразу отметим, что динамические режимы в нижней части фазовой диаграммы, включая область IV, в данной работе не рассматриваются, поскольку при низких скоростях диффузии радиальное распределение СО будет существенно неоднородным и двухслойная модель будет неадекватной. Сопоставляя области I, II, III с соответствующими областями на фазовой диаграмме для отдельных осцилляторов (Рис. 3.3) можно заключить, что в области I активны оба осциллятора, а в областях II и III активен только один внутренний осциллятор. Характер колебаний в области I будет определяться силой связи между осцилляторами, которая увеличивается с уменьшением скорости диффузии. В общем случае можно ожидать, что в области I будут квазипериодические колебания. В областях II и III вблизи внешних границ активные колебания будут только во внутреннем слое, во внешнем слое будут пассивные (вынужденные) колебания. Неустойчивые динамические режимы могут наблюдаться в окрестности границы между областью I и областями II и III, где влияние внутреннего осциллятора на еще не установившиеся колебания внешнего могут приводить к усложнению колебаний.

Общее представление о характере качественных изменений в режиме колебаний при изменении некоторого параметра модели дает бифуркационная диаграмма колебаний. Бифуркационная диаграмма колебаний представляет собой график, по горизонтальной оси которого откладываются значения параметра, а по вертикальной - значения некоторой характерной величины, которые она принимает в процессе колебаний при каждом значении параметра. В качестве характерной величины часто берут максимальное (или минимальное) значение одной из динамических переменных или некоторой функции от динамических переменных. В данной работе в качестве такой величины выбран максимум глобальной скорости реакции которая является одной из величин, доступной экспериментальному измерению.

На практике бифуркационные диаграммы строятся следующим образом. При фиксированном значении параметра изучаемая система уравнений решается на достаточно большом интервале времени. Длительность интервала определяется временем выхода системы на стабильный режим после изменения значения параметра и временем, необходимым для сбора достаточного количества значений наблюдаемой величины. В частности, в данной работе длительность времени установления задавалась числом Ni моментов, когда Rg принимает максимальное значение (для простых периодических колебаний с одним максимумом это означает Ni периодов колебаний), и числом iV2 последовательных максимумов Rg, которые откладываются на бифуркационной диаграмме. По истечении этого времени параметр модели изменялся и все повторялось при новом значении параметра, причем конечная точка траектории системы при старом значении параметра использовалась в качестве начальной точки при новом значении. Числа Ni и N2 подбираются опытным путем и могут изменяться в ходе построения диаграммы в зависимости от сложности динамики системы. Для определения возможных режимов колебаний в двухслойной модели были построены две бифуркационные диаграммы. Одна диаграмма по параметру D при Рсо = 2 Торр (вертикальная прямая на Рис. 3.5), другая диаграмма по параметру Рсо при D = 1052 с-1 (горизонтальная прямая на Рис. 3.5). Первая диаграмма описывает изменения в динамике системы при изменении скорости диффузии внутри области колебаний, а вторая - изменения в динамике при пересечении линий бифуркации Андронова-Хопфа и линии h2 при изменении давления Рсо Рассмотрим бифуркационную диаграмму по скорости диффузии (Рис. 3.6). Диаграмма в верхней части рисунка построена при изменении lg D с шагом 0.005 в интервале от 1 до 10. Диаграмма в нижней части рисунка построена с шагом 0.001. При каждом значении D числа N% и iV2 были равными 10 и 100, соответственно. В каждом цикле колебаний по вертикальной оси откладывается значения локальных максимумов скорости реакции. Для простых периодических колебаний имеется только один максимум (одно-оборотный цикл). Для более сложных много-обротных колебаний число локальных максимумов за период увеличивается. Из вида бифуркационной диаграммы можно заключить, что в системе (3.3.5) при разных значениях D существуют периодические, квазипериодические и, возможно, более сложные хаотические колебания. Рассмотрим более подробно квазипериодические и периодические (резонансные) колебания. Хаотические колебания, в силу их особых свойств, будут описаны в отдельном разделе.

При высоких скоростях диффузии изменения давления Рсо за период колебаний в каждом слое зерна очень малы и колебания реакции в каждом слое можно считать независимыми. Связь между внутренним и внешним осцилляторами слабая, однако частоты колебаний осцилляторов не совпадают, хотя и очень близки. С течением времени разность фаз колебаний осцилляторов медленно изменяется (растет), т.е. в системе двух осцилляторов происходят биения. При уменьшении скорости диффузии сила связи между осцилляторами возрастает. Однако в данном случае увеличение силы связи приводит не к сближению собственных частот, а, наоборот, к их большему различию. При этом период биений уменьшается.

При больших D в правой части верхней диаграммы Рис. 3.6 из-за длинного периода биений сто последовательных максимумов скорости реакции совпадают в масштабе диаграммы и попадают в одну точку по вертикальной оси. Однако положение этой точки на вертикальной оси постепенно изменяется с уменьшением скорости диффузии. Это изменение связано с биениями, которые получаются в результате сложения колебаний двух почти независимых осцилляторов с близкими частотами (частоту биений можно считать примерно равной разности основных частот осцилляторов). На Рис. 3.7а показан пример таких биений при D = 109 с-1. На этом рисунке цветными линиями показаны графики зависимости от времени максимумов и минимумов колебаний скорости реакции.

Оценка величины временной задержки

В этом случае участок кривой S между экстремумами будет неустойчивым, а остальная часть кривой устойчива. В полной системе (4.1.1) участок кривой медленных движений между экстремумами является отталкивающим, внешние ветви кривой являются притягивающими. Если в системе (4.1.1) есть автоколебания, то они представляют собой последовательные быстрые переходы от одной притягивающей ветви кривой S к другой и медленные движения вдоль устойчивых ветвей кривой по направлению к экстремуму.

Красный кружок на рисунках показывает положение стационарного решения системы на кривой медленных движений. На Рис. 4.13а, на котором показана фазовая траектория стохастической модели для 10-ти нанометровой частицы при Рсо = 2.5 Торр, стационарное решение лежит на притягивающей ветви кривой S в области малых значений х и больших значений у (область высокой скорости реакции -высокоскоростная ветвь). Колебаний в точечной модели нет, но в стохастической модели наблюдаются случайные колебания. Эти колебания вызваны статистическими флуктуациями фазовых переменных модели. В ходе эволюции фазовая точка, изображающая состояние модели, стремится к стационарному состоянию. Однако статистические флуктуации не позволяют системе достаточно долго находится в окрестности стационарного состояния. Если случайное отклонение от стационарного состояния достаточно велико, то для того, чтобы вернуться в окрестность стационарного состояния фазовая точка должна проделать длинный путь до противоположной притягивающей ветви кривой S (низкоскоростная ветвь) и вернуться обратно. Эти случайные переходы от высокоскоростной ветви S к низкоскоростной и обратно и являются причиной колебаний, показанных на Рис. 4.10.

На Рис. 4.13Ь показана ситуация, когда стационарное решение лежит на отталкивающем участке кривой медленных движений и в точечной модели существует автоколебательное решение (частица 10 нм, давление Рсо = 3.5 Торр). Предельный цикл в точечной модели показан на рисунке красной линией. Из рисунка видно, что стохастическая траектория следует вдоль детерминистической траектории. Однако, момент срыва стохастической траектории, как правило, опережает момент срыва детерминистической траектории, поэтому средний период стохастических колебаний меньше периода автоколебаний точечной модели. В частности, на Рис. 4.10 для давления Рсо = 3.5 Торр вместе со стохастическими колебаниями скорости реакции красной пунктирной линией показаны детерминистические колебания. Из этого рисунка можно заключить, что за время равное 500 секундам стохастическая траектория совершает 20 оборотов, а детерминистическая только 16.

Наконец, на Рис. 4.13с стационарное решение точечной модели лежит на притягивающей низкоскоростной ветви кривой S (частица 10 нм, давление Рсо = 5.0 Торр), т.е. в точечной модели нет автоколебаний. Однако, в стохастической модели, также как в случае (а), наблюдаются колебания, вызванные статистическими флуктуациями концентрации реагентов. Соответствующие колебания скорости реакции показаны на Рис. 4.11. Здесь колебания имеют большой средний период (по сравнению с колебаниями при Рсо = 3.5 Торр), так как при данных условиях притяжение к низкоскоростной ветви кривой медленных движений сильнее, чем к высокоскоростной.

На Рис. 4.13d показана стохастическая траектория для частицы диаметром 20 нм при давлении Рсо = 3.0 Торр. Здесь кривая медленных движений также имеет форму кубической параболы и стационарное решение лежит на притягивающей низкоскоростной ветви этой кривой. В этом случае стационарное решение лежит достаточно близко к экстремуму кривой S, а статистические флуктуации достаточно велики для того, чтобы фазовая точка могла время от времени отрываться от низкоскоростной ветви и переходить в окрестность высокоскоростной ветви S. Эти переходы, происходящие в случайные моменты времени, проявляются в виде колебаний скорости реакции, показанных на Рис. 4.12 при Рсо — 3.0 Торр.

В моделях реакции в пористом шаре (3.3.5), (3.6.7) не учитывались особенности диффузии в цеолите. Попытка сделать более реалистичную модель с учетом диффузии многокомпонентной газовой смеси в пористой среде предпринята в работе [86]. В этой работе были получены модельные уравнения для описания реакции на одномерной системе наночастиц (цепочке) встроенных в матрицу цеолита. Модель использовалась для численного исследования динамики реакции в зависимости от числа частиц в цепочке и расстояния между частицами.

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

Кристаллическая структура цеолита типа faujasite была кратко описана в начале Главы 3. Основным элементом структуры является так называемая клетка (см. Рис. 3.1). Клетки соединяются между собой шестиугольными гранями и образуют тетраэдральную структуру, состоящую из суперклеток. Диаметр суперклетки равен 1.3 нм, суперклетки соединяются шестиугольными окнами диаметром 0.74 нм. Если клетки соединить между собой квадратными гранями, то получится цеолит типа LTA, который имеет кубическую структуру суперклеток. Для удобства количественных оценок будем считать, что цеолит имеет кубическую структуру с диаметром суперклетки 1.3 нм. Это допустимое упрощение, поскольку нашей целью является только качественное описание реакции. В модели структура цеолита представлена единственным параметром р, который обозначает число суперклеток в единице объема цеолита, т.е. плотность суперклеток. Одна суперклетка цеолита LTA схематически изображена на следующем рисунке.

Молекулы газа, эффективный диаметр которых несколько меньше диаметра окон между суперклетками, могут проникать внутрь цеолита. При этом в одной суперклетке может располагаться несколько молекул, однако для простоты в дальнейшем предполагается, что в суперклетке может находиться не более одной молекулы газа. Через окна молекулы могут перемещаться в соседние суперклетки. Таким образом, диффузия газа в цеолите обладает характерными чертами диффузии решеточного газа. В частности, эффективный коэффициент диффузии зависит от концентрации газа в цеолите [71].

Предполагается, что частицы металла, на поверхности которого происходит реакция, имеют форму октаэдра. Частицы выращиваются внутри цеолита и имеют размеры, намного превосходящие размеры суперклетки. В процессе роста частица "ломает" структуру цеолита вокруг себя. Технология выращивания частиц позволяет достаточно точно контролировать их размер и получать достаточно узкое распределение частиц по размерам. Будем считать, что каждая частица располагается в центре куба, длина ребра которого равна среднему расстоянию между частицами. Такой куб в дальнейшем называется элементарным кубом катализатора. На Рис. 4.15 показан разрез куба цеолита с частицей металла в центре.

Похожие диссертации на Нелинейная динамика молекулярных процессов в гетерогенных системах

Математическое моделирование процессов нелинейной динамики в замкнутых системах автоматического управления с однополярной реверсивной модуляцией Андриянов Алексей Иванович

Методы математического моделирования динамики процессов окисления в системе биологической очистки водыСмирнов Николай Васильевич

Моделирование динамических процессов в конденсированном веществе методом динамики частиц с использованием многопроцессорных вычислительных системЛе-Захаров, Александр Аневич

Математическое моделирование с помощью многопроцессорных вычислительных систем процессов электронного транспорта в вакуумных и твердотельных микро- и наноструктурахПоляков, Сергей Владимирович

Автоматизированное управление регистрациями процессов в системах менеджмента качества предприятий и организацийСеребряков Юрий Евгеньевич

Моделирование вычислительного процесса в системах навигации летательного аппарата, разработка алгоритмов и комплексов программ для его реализации на программируемых логических интегральных схемахЧумакова Екатерина Витальевна

Интеллектуальные модели качественного анализа слабо формализованных временных процессов в системах поддержки принятия решений Муравский Александр Викторович

Модели и алгоритмы информационно-измерительной системы для процесса формования листового стекла Фоминых Дмитрий Сергеевич

Моделирование переходных и колебательных процессов в системах со сложными запаздываниями Музафаров Салих Мухаррамович

Идентификация моделей совместных случайных процессов для систем контроля горной техникиНовосельцева Марина Александровна