CC BY
30
Автовслны. Самоорганизация
Изв. вузов «ПНД», т. 14, № 2, 2006 УДК 519.245+517.9+517.39
ФОРМИРОВАНИЕ И РАЗВИТИЕ ПРОСТРАНСТВЕННЫХ СТРУКТУР В СИСТЕМЕ ХИМИЧЕСКИХ РЕАКЦИЙ НА КАТАЛИТИЧЕСКОЙ РЕШЕТКЕ: МОДЕЛИРОВАНИЕ МЕТОДОМ МОНТЕ-КАРЛО
А.В. Ефимов, А.В. Шабунин
Рассматривается формирование пространственных структур (кластеров) в ходе циклических превращений в модели (4+1)-Lattice Lotka-Volterra. Система моделируется с помощью разновидности метода Монте-Карло на поверхности двумерной решетки. Устанавливаются закономерности в распределении кластеров по размерам и в спаде пространственной автокорреляции. Исследуются другие зависимости, характеризующие пространственную динамику модели, а также выясняется влияние кластеров на динамику системы. В работе показывается, что добавление в систему внешнего перемешивания приводит к появлению периодических автоколебаний.
Введение
Синхронизация и образование пространственных структур остается одним из основных направлений исследований в нелинейной динамике. Как правило, объектом исследования являются ансамбли взаимодействующих осцилляторов, то есть динамические системы, демонстрирующие детерминированное поведение. При этом наличие внешних шумов, хотя и играет иногда существенную роль в динамике таких систем, все же не является определяющим фактором. Однако иногда случайные процессы, протекающие на микроуровне, сами являются «строителями» сложных пространственно-временных структур, возникновение которых приводит к появлению «глобальных» колебаний на макроуровне, аналогичных колебаниям детерминированных осцилляторов. Типичным примером здесь могут являться системы, описывающие химические превращения, особое место среди которых занимают модели гетерогенного катализа на поверхности кристаллических решеток различных металлов. Такие реакции, как окисление газа CO на поверхности Pt и восстановление NO в присутствии H2 на Pt, имеют особую важность в различных технологических процессах. Кроме того, гетерогенный катализ - это отрасль, в которой возникновение колебаний, пространственно-временных структур и различные волновые явления часто наблюдаются в натурных экспериментах [1-8].
Пожалуй, первое упоминание о колебаниях в гетерогенных системах было опубликовано еще в 1970 году. После этого регулярные и хаотические колебания были обнаружены более чем в тридцати различных реакциях, практически на всех типах катализаторов. За последние двадцать лет были предложены различные подходы к численному моделированию таких реакций [9-12] и получены новые экспериментальные результаты [3,4,13].
Большинство моделей вышеупомянутых систем можно условно разделить на две группы. В первой группе системы взаимодействующих частиц рассматриваются на макроуровне, например, с помощью феноменологических моделей среднего поля - «Mean Field», которые считаются вполне приемлемыми для гомогенных химических сред с сильным перемешиванием. Уравнения среднего поля дают некоторые представления о ходе химических реакций в таких системах, как брюсселятор и оре-гонатор [14], но не могут адекватно описывать динамику химических превращений на поверхности кристаллической решетки катализатора, так как не учитывают локальный характер взаимодействий. Ко второй группе можно отнести моделирование на более глубоком, микроуровне, позволяющее наблюдать такие интересные явления, как формирование пространственных структур, волновые процессы и сложные колебания концентраций реагентов [9,11,12,15-21].
Один из широко применяемых численных методов моделирования локальных взаимодействий - метод Монте-Карло (Kinetic Monte Carlo simulation), являющийся разновидностью вероятностного клеточного автомата. Возможность моделировать локальные взаимодействия элементов ансамбля с помощью этого метода оказалась полезной для исследования гетерогенных неперемешивающихся сред и других систем, в которых локальные флуктуации имеют особое значение. Именно этот метод применяется в настоящей работе для моделирования динамики системы реакций, которая при использовании метода среднего поля описывается уравнениями, аналогичными системе Лотки - Вольтерры, но с большей размерностью фазового пространства (D = 4). В дальнейшем эту систему будем называть (4+1)-LLV по аналогии с системой (2+1)-Lattice Lotka-Volterra [20] (цифры в скобках означают число реагентов, принимающих участие в реакциях).
В настоящей работе ход процессов моделируется методом Монте-Карло, проводится детальное исследование пространственного поведения системы (4+1)-LLV, исследуются процессы образования и эволюции пространственных кластеров, проводится их статистический анализ, рассматривается зависимость в поведении системы от параметров и начальных условий. В работе также изучается влияние внешнего перемешивания компонентов реакций и показывается, что рост интенсивности перемешивания ведет к рождению периодических колебаний концентраций, сходному с бифуркацией Андронова - Хопфа для динамических систем.
1. Исследуемая система и ее описание методом среднего поля
Рассмотрим реакции четырех веществ X, Y, Z, U, протекающие на поверхности катализатора. Поверхность представляет собой регулярную квадратную решетку, каждая ячейка которой может содержать только одну из частиц реагентов или оставаться «пустой». Вакантное место на решетке, то есть пустая ячейка - «дырка», будет обозначаться как частица S. Взаимодействие между частицами является локальным, то есть происходит между реагентами, находящимися в соседних ячейках решетки,
и описывается следующей кинетической схемой:
X + У кх 2У,
У + Z Н 2Z,
Z + и кз 2и,
и + 5 к4 25,
X + 5
где параметры к представляют собой кинетические константы, пропорциональные вероятностям соответствующих превращений. Так, в соответствии с первой реакцией (1), соседние частицы X и У взаимодействуют друг с другом с вероятностью, пропорциональной к\, при этом частица X превращается в У. Данная реакция является автокаталитической: наличие в точке пространства вещества У вовлекает X в реакцию, в результате которой число частиц У увеличивается. Похожим образом взаимодействуют и остальные частицы У, Z и и. Фактически мы имеем набор циклических превращений, в которых реагенты преобразуются друг в друга. Четвертая строка схемы (1) описывает удаление продукта реакции и с поверхности катализатора, а пятая реакция - процесс адсорбции реагента X на решетку. Важным свойством исходной схемы превращений является то, что в каждом шаге реакции участвуют две частицы и две частицы остаются после взаимодействия. Это приводит к тому, что общее число взаимодействующих частиц в ходе реакции остается неизменным.
В предыдущих работах [20-22] рассматривалась система (2+1)-ЬЬУ, содержащая два типа взаимодействующих частиц (X и У), то есть исследовался более короткий цикл превращений. В рамках модели среднего поля ее уравнения имеют интеграл движения, то есть являются консервативными. Здесь следует отметить, что наличие этого интеграла характерно только для циклических систем типа (2к+1)-ЬЬУ, то есть для систем с четным числом компонентов. Системы с нечетным количеством реагентов не имеют интеграла движения и поведение их моделей среднего поля качественно отличается от поведения систем (2к+1)-ЬЬУ.
Минимальное усложнение системы за счет увеличения числа реагентов, не вызывающее изменения ее класса, дает систему (4+1)-ЬЬУ (1), рассматриваемую в настоящей работе. Ее детальный анализ методом среднего поля был проведен в [23], где были получены и исследованы дифференциальные уравнения для относительных концентраций х, у, г, и
(х
М
йуу (И (г
( (и
(
Уравнения (2) описывают консервативный осциллятор в четырехмерном фазовом пространстве, для которого определен интеграл движения
= -к\ху + к5х(1 — х — у — г — и),
= к\ху — к2уг,
= к2уг — кзги, (2)
= к3ги — к4и(1 — х — у — г — и).
хк2к4 укзкб гк1к4 ик2кб (1 — х — у — и — г)к1кз = К. (3)
Фазовое пространство осциллятора устроено следующим образом. Оно содержит:
• тривиальное состояние равновесия Р\ = (0; 0; 0; 0), соответствующее системе без реагентов;
• четыре линии, образованные из состояний равновесия Р2 = (0;0; z;0), Р3 = (0; у;0;0), Р4 = (0; у; 0; 1 — у), Рб = (х;0;1 — х;0), соответствующие случаю, когда часть реагентов «вымывается» из пространства реакций, что приводит к их прекращению («отравление» решетки);
• нетривиальное состояние равновесие типа центр
= (к3к5 к\кд к2к5)
Р7 = ( А ; А ; А ; А ^ (4)
где А = к\к3 + к2к4 + к3к5 + кк + к2к5.
Интересующая нас область фазового пространства ограничена инвариантными многообразиями: х = 0, у = 0, г = 0, и = 0 и х + у + г + и = 1, образующими замкнутый контур. При выборе начальных условий на одном из ограничивающих многообразий траектория по нему переходит на одно из состояний равновесия Р2 - Рб, что приводит к прекращению реакций. Если начальные условия выбрать внутри контура, то в зависимости от начальных условий и параметров система (2) демонстрирует квазипериодическое или хаотическое поведение [23]. Однако описание методом среднего поля является слишком грубым, и хорошо описывает лишь некоторый средний уровень концентраций веществ, вокруг которого происходят колебания. Непосредственное моделирование реакций в каждой точке пространства методом Монте-Карло дает возможность как детального пространственного анализа процессов на катализаторе, так и получения более реалистичной динамики для усредненных концентраций.
2. Описание используемого алгоритма метода Монте-Карло
Метод Монте-Карло является разновидностью вероятностных клеточных автоматов и моделирует процесс реакции в каждой точке решетки в зависимости от находящейся там частицы и ее ближайших соседей. Используемый в работе алгоритм метода заключается в следующем. Предположим, что мы моделируем реакции (1) на двумерной квадратной решетке со стороной Ь, то есть содержащей М = Ь2 элементов. Заданы значения параметров реакций (кинетические константы кг) и начальные концентрации
Хо Уо Zo ио
Хо = М,уо = М,го = М,ио = М,
где Хо, У0, Zо и ио - количество частиц соответствующих реагентов на решетке в начальный момент времени. До начала моделирования решетка равномерно заполняется этими частицами в соответствии с выбранными начальными концентрациями. При моделировании реакций на каждом микрошаге алгоритма случайным образом выбирается ячейка решетки и один из ее четырех ближайших соседей. В зависимости от того, какие частицы находятся в выбранных ячейках, может осуществиться одна из пяти реакций (1) с вероятностью, пропорциональной соответствующей кинетической константе рг = кг/шах {к^ }5. Пусть, например, в выбранной ячейке находится частица X, а в выбранной случайным образом соседней ячейке - У. Тогда,
в соответствии с первой реакцией схемы (1), частица X с вероятностью pi может превратится в Y. Если же соседние ячейки оказываются заполненными частицами, реакция между которыми невозможна, то состояние решетки на данном микрошаге не изменяется. Единичный временной интервал моделирования, «Monte-Carlo step» (далее - MCS), содержит M описанных выше микрошагов, то есть в среднем каждая из ячеек посещается один раз. В результате от шага к шагу моделируется эволюция пространственного распределения на решетке, по которому можно строить различные усредненные характеристики, то есть определять макросостояния системы.
3. Общие свойства динамики системы
Проведенное моделирование системы (1) при различном выборе начальных концентраций и параметров, а также на решетках разных размеров позволило заключить, что в целом возможны два качественно различных вида эволюции системы:
а) после переходного процесса усредненные по решетке концентрации реагентов совершают установившиеся нерегулярные шумоподобные колебания с малой амплитудой или
б) происходит вытеснение одного или нескольких реагентов, в результате чего все реакции «замораживаются». При этом на поверхности решетки остаются только невзаимодействующие частицы.
Прекращение реакций после переходного процесса происходит, если:
• значения выбранных скоростей реакций ki значительно (в несколько раз) отличаются друг от друга;
• выбран малый размер решетки. В частности, при равных значениях кинетических констант, колебания «замораживаются» на решетках меньших, чем размер M ~ 70 х 70 элементов.
При сопоставимых значениях кинетических констант и при линейных размерах решетки порядка сотен элементов колебания на ней могут наблюдаться неограниченно долго. При этом изменение параметров ki и выбор различных начальных условий не приводит к качественным изменениям в динамике системы, за исключением описанного выше случая «замораживания» и изменения характерного временного масштаба (среднего периода) колебаний. Модель не демонстрирует фазовых переходов при вариации параметров, что вполне согласуется с отсутствием бифуркаций в модели среднего поля подобных систем [20,23]. Интенсивность колебаний средних концентраций спадает с ростом линейного размера решетки L по закону 1/л/М, где M = L2. Это означает, что глобальные колебания системы есть усредненный вклад множества независимых случайных процессов, протекающих локально в каждой точке поверхности решетки. В результате усреднения концентраций реагентов по ансамблю этих осцилляторов «амплитуда» глобальных колебаний стремится к нулю с ростом числа осцилляторов, а значит, и с ростом размера решетки. Фактически мы наблюдаем флуктуации вокруг средних концентраций хср, уср, гср и иср. Значения этих средних концентраций с высокой степенью точности совпадают с координатами нетривиального состояния равновесия (4) уравнений (2). Таким образом, в среднем решетка ведет себя как набор пространственно распределенных стохастических осцилляторов, колебания в каждом из которых происходят вокруг фиксированного состояния равновесия, положение которого определяется кинетическими константами.
4. Пространственно-временная динамика модели Монте-Карло
Рассмотрим локальное распределение частиц на решетке в зависимости от времени. Характерной чертой в поведении модели Монте-Карло является наличие переходного процесса от первоначального равномерного распределения к некоторому неоднородному пространственному распределению. В предыдущей работе [23] было показано, что в пространстве взаимодействия в течение переходного процесса происходит образование кластеров - областей решетки, занятых частицами одного типа. На рис.1 представлено характерное пространственное распределение частиц реагентов на поверхности решетки катализатора в начале реакции и после переходного процесса. Для удобства мы изобразили частицы только одного вещества X. Частицы других веществ и пустые участки решетки отмечены белым цветом. Начальное распределение частиц по поверхности решетки носит случайный характер и обладает свойством пространственной однородности (рис. 1, а). После начала реакций, в соответствии со схемой (1), частицы разных видов взаимодействуют и образуют кластеры (рис. 1, б), которые постоянно «движутся» по решетке, меняя свою форму и размеры. Как будет показано дальше, возникновение кластеров оказывает сильное влияние на временную динамику системы. В силу локальности взаимодействий, частицы соседних кластеров могут реагировать только на их границах. Вследствие этого количество частиц, участвующих в реакции, после образования кластеров уменьшается. Чем больше размеры кластеров, тем меньше частиц провзаимодействуют на каждом шаге алгоритма Монте-Карло.
Рассмотрим процесс образования кластеров на решетке более подробно. Остановимся сначала на исследовании характерных размеров кластеров. На рис. 2 представлена эволюция распределения кластеров по размерам. Построена вероятность р(М) = им/и появления на поверхности решетки X-кластера размером в N частиц в различные моменты времени. Здесь им - количество кластеров, содержащих N частиц X-типа, а и - общее число кластеров, содержащих частицы X. Сначала все кластеры имеют малый размер, обычно не превышающий 10 ячеек (рис. 2, а), что соответствует равномерному начальному распределению частиц. Затем на решетке появляются кластеры все большего размера (рис. 2, б). Переходный процесс, сопро-
Рис. 1. Распределение частиц X на поверхности решетки: в начальный момент времени (а) и после переходного процесса (б)
N ¿) N е N
Рис. 2. Распределение кластеров по размерам. Размер решетки 200 х 200 ячеек
вождающийся формированием кластеров, завершается примерно к моменту времени Ь = 200 МС8, после чего мы имеем распределение (рис. 2, г, д, е) со средним размером кластеров, осциллирующим около значения 25-30 ячеек. Однако время от времени в решетке появляются кластеры, содержащие более тысячи частиц, и очень редко - кластеры размером 40005000 ячеек.
Для решеток малого размера распределение кластеров нестационарно. Оно изменяется соответственно с изменениями концентраций реагентов. Процесс можно рассматривать как стационарный только в предельном случае решетки бесконечного размера, когда «амплитуда» глобальных колебаний стремится к нулю. Рис. 3 демонстрирует временную эволюцию величин, характеризующих динамику кластеров. В частности, на графике 3, б представлена эволюция среднего размера кластеров (К), содержащих частицы X. Для сравнения на графике 3, а показано изме-
400 600
800 1000
Рис. 3. Временные реализации концентрации частиц X (а); зависимости от времени среднего размера X-кластеров (б) и размера наибольшего кластера X-типа (в). Размер решетки Ь2 = = 100 х 100. Параметры к1,2,з,4,б = 1
нение средней концентрации частиц X в решетке, а на графике 3, в - колебания размера наибольшего кластера, состоящего из X-частиц. Начиная с первоначального одно родного распределения, с течением времени в решетке наблюдается постепенный рост среднего размера кластеров. Примерно к моменту времени £ ~ 200 МС8, когда процессы кластерообразования завершаются, размер кластеров перестает расти и начинает осциллировать почти с постоянной интенсивностью около своего среднего значения. Для маленьких решеток средний размер кластеров и усредненная концентрация вещества X хорошо коррелируют с размером наибольшего кластера. То есть рост концентрации веществ сопровождается увеличением размеров соответствующих кластеров, а пики на всех трех графиках примерно совпадают. Для больших решеток (когда Ь ^ 200) это соответствие начинает исчезать. Взаимосвязь между размером самого большого кластера и концентрацией соответствующего вещества численно можно охарактеризовать коэффициентом взаимной корреляции
С =
((Птах(£) - (Птах))(х(£) ~ (х))) \jDnDx
(5)
Рис. 4. Взаимная корреляция С между размером наибольшего X -кластера и средней концентрацией х частиц X как функция линейного размера решетки Ь
Здесь птах - число частиц в самом большом кластере, а Оп^х = ((птах(£) --(Птах))2) и Бх = ((х(г) - (х))2) соответственно дисперсии птах и концентрации частиц X. Зависимость коэффициента С от размера решетки представлена на рис. 4. Его значение монотонно спадает с ростом линейного размера решетки приблизительно от 0.84 для наименьшей из рассматриваемых решеток (Ь2 = 70 х 70) до 0.35 для решетки размером 500 х 500 ячеек. Для небольших решеток корреляция достаточно сильная, что свидетельствует о том, что рост концентрации того или иного компонен-
та связан с образованием на решетке кластеров большого размера, а не с ростом числа средних кластеров. Слабая корреляция для больших решеток говорит о том, что пространственное взаимодействие на больших расстояниях очень слабое.
Происходящее на поверхности решетки образование структур можно рассматривать как смешанный процесс, включающий в себя две составляющие, имеющие различные пространственные масштабы. Первый характерный масштаб соответствует среднему размеру кластеров и определяет пространственно-временную динамику в выделенной малой области решетки. Этот масштаб тесно связан с размером и формой кластеров. Второй масштаб относится к динамике решетки как единого целого. Он определяется перемешиванием, которое происходит из-за несинхронного образования кластеров в различных частях решетки. Вследствие этого статистические свойства распределения кластеров по размерам должны отличаться для малых и больших решеток. Для больших решеток распределение кластеров должно отражать свойства обоих масштабов. Чтобы проверить это предположение, рассмотрим,
Рис. 5. Интегральная вероятность рс(п) появления на поверхности решетки X-кластера, размером не менее N, ячеек для решеток 70х70 ячеек (а); 100х 100 ячеек (б); 150х 150 ячеек (в); 300x300 ячеек (г); 500 х 500 ячеек (д); 1000 х 1000 ячеек (е). Кружки соответствуют максимумам X-концентрации, крестики - минимумам. Для сравнения жирной линией показано начальное равномерное распределение (для Ь = 0 MCS). Штриховой линией на рис. (а) и (б) показана аппроксимирующая функция
как изменяется распределение кластеров с ростом размера решетки. На рис. 5 в двойном логарифмическом масштабе изображены графики зависимости интегральной вероятности рС(п) появления в решетке X-кластера, состоящего не менее чем из N частиц, от размера N для решеток с разными значениями Ь.
<х
РС(К ) = £ р(г). (6)
г=М
Так как процесс на поверхности решетки не является стационарным, мы проверили
поведение вероятности рС(п) для двух характерных моментов времени: когда концентрация х вещества X достигает своего максимального и минимального значения. На рис. 5 соответствующие точки графиков отмечены кружками и крестиками. Минимальный размер решетки, при котором нам удавалось наблюдать колебания - решетка 70 х 70 ячеек. При меньшем размере решетки происходит быстрое «вымывание» одного из реагентов и реакции прекращаются. Распределение вероятности рС(п) для такой решетки представлено на рис. 5, а. На графиках все точки располагаются вблизи прямых (отмечены штриховыми линиями) с наклоном от 0.35 (для максимумов х) до 0.65 (для минимумов х). Из-за малого количества точки графика сильно разбросаны относительно аппроксимирующих прямых. Для решетки размером 100 х 100 ячеек мы имеем хорошую статистику и отклонение от аппроксимирующих прямых становится меньше. Линейная зависимость в двойном логарифмическом масштабе соответствует степенному закону для интегральной функции распределения вероятности
где степень ^ принимает значения от 0.4 (для максимальной концентрации данного вещества) до 0.6 (для минимальной концентрации). В итоге само распределение р(п) описывается степенным законом
Степенной закон в функции распределения кластеров по размерам говорит об их фрактальных свойствах. При увеличении размеров решетки Ь наблюдаются следующие особенности.
• Флуктуации вероятности становятся меньше, так как количество кластеров возрастает и статистика улучшается (рис. 5, в-е).
• Значения вероятности для моментов времени с различными концентрациями реагента становятся ближе друг к другу, что является следствием исчезновения колебаний на макроуровне в больших решетках. Процесс приближается к стационарному.
• Распределение кластеров по размерам для разных решеток ведет себя по-разному. Для малых решеток в распределении наблюдается степенной закон, в то время как для больших размеров над степенным законом начинает доминировать экспоненциальная зависимость.
Последний эффект, по нашему мнению, связан, во-первых, с процессом формирования фрактальных кластеров в малом масштабе и, во-вторых, с процессом перемешивания в макромасштабе из-за несинхронного появления кластеров в различных частях решетки. Перемешивание приводит к экспоненциальному спаду пространственной корреляции и функции распределения кластеров по размерам. Два независимых процесса, одновременно протекающих в микро- и макромасштабе большой решетки, должны приводить к сложному виду распределения кластеров по размерам, составляющими которого будут степенной и экспоненциальный законы. На самом деле, для больших решеток (рис. 5, е) распределение интегральной вероятности хорошо аппроксимируется функцией
рС(п) ~ п
(7)
р(К) - N""
(8)
рС(п) ~ ехр(—Х(п — 1))п
(9)
где X = 0.002 и ^ = 0.6.
5. Модель (4+1)-LLV с внешним перемешиванием
Как отмечалось выше, формирование кластеров определяющим образом влияет на временную динамику системы. Ход реакций сдерживается за счет ограничения пространства взаимодействий: реакции протекают только на границах между однородными кластерами. Глобальная динамика при этом, в пределе большой решетки, представляет собой стационарное состояние с концентрациями, определяемыми состоянием равновесия (4) из уравнений среднего поля (2). Если выбрать начальные концентрации отличными от (4), то по прошествии переходного процесса концентрации придут к данному состоянию равновесия, сделав предварительно несколько витков вокруг него. Глобальное поведение системы напоминает поведение динамической системы в окрестности устойчивого фокуса, то есть системы с диссипацией. Налицо качественное различие в глобальном поведении модели Монте-Карло и в динамике модели среднего поля.
Рассмотрим теперь, как изменится поведение решетки, если мы введем управляющее воздействие, разрушающее однородные кластеры, например, путем внешнего перемешивания. Очевидно, что такое вмешательство должно стимулировать ход процессов на решетке, поскольку взаимодействие между частицами в однородной системе протекает одновременно во всех точках пространства реакций. Чтобы получить подобный эффект, мы добавили в алгоритм метода Монте-Карло на каждой временной единице перемешивание среды, которое заключается в перетасовывании содержимого ячеек решетки случайным образом. Такое перемешивание не похоже на обычную диффузию, так как частица из одной ячейки не проходит весь путь по решетке до второй, попутно взаимодействуя с окружающими частицами, а мгновенно меняется с ней местами. Подобным образом, вероятно, можно моделировать реакции в газовых и жидких фазах с принудительным перемешиванием среды, однако для семейства моделей ЬЬУ такое перемешивание возможно только если частицы всех типов будут постоянно покидать решетку и адсорбироваться из газовой фазы обратно, причем кинетические константы этих процессов должны быть значительно выше соответствующих констант для самих реакций.
Прежде чем переходить непосредственно к результатам моделирования, введем некий параметр, который будет количественно характеризовать степень перемешивания в системе,
п
р = м (10)
где п - число элементарных актов перемешивания (перестановок) на каждом шаге МС8, а М - общее число ячеек решетки. Этот параметр показывает, какая часть ячеек была перемешана на каждом шаге алгоритма. Если р ^ 1, это означает, что в среднем каждая частица изменила свое положение, и можно говорить о нелокальности превращений, входящих в исходную схему (1).
Обратимся теперь к результатам моделирования. На рис. 6 представлена временная реализация для концентрации частиц X при достаточно малом значении р = 0.0056. Несмотря на слабость перемешивания, результаты значительно отличаются от временных реализаций, представленных на графике 3, а. Размах колебаний заметно возрастает, так же как и степень их регулярности. Колебания становятся почти периодическими: в спектре мощности (рис. 7) появляются ярко выраженные пики на кратных частотах на фоне небольшого шумового пьедестала. Эти колебания,
Рис. 6. Временные реализации для концентрации частиц X в системе (4+1)-ЬЬУ с перемешиванием на решетке размером 300 х 300 ячеек: а - малое перемешивание, р = 0.0056; б - большое перемешивание, р = 1. Начальные концентрации: х0 = у0 = г0 = п0 = 0.2; параметры: кх = к2 = к3 = к4 = к5 = 1
по-видимому, можно считать периодическими колебаниями с шумом, который неизбежно возникает в системе за счет вероятностного характера протекающих реакций. Косвенным подтверждением этого является тот факт, что степень «регулярности» колебаний увеличивается с ростом размера решетки. Так, спектр для решетки размером 500 х 500 элементов (изображен темным цветом на рис. 7) имеет более узкие пики и более низкий шумовой пьедестал, чем спектр колебаний решетки размером 300 х 300 элементов (показан более светлым).
Понять причину таких изменений в поведении системы при добавлении в нее даже слабого перемешивания нетрудно, если обратиться к пространственной динамике. Как и в случае исходной модели, мы создаем однородное начальное распределение частиц с выбранными средними концентрациями реагентов (рис. 8, а). Так как перемешивание слабое, оно не может полностью предотвратить образование кластеров, и после переходного процесса образуется характерное распределение кластеров, «снимок» которого представлен на рис. 8, б. Этот рисунок сильно напоминает распределения, полученные для исходной системы (см. рис. 1). Однако есть существенное отличие. При отсутствии перемешивания после образования кластеров все взаимодействия на решетке носят граничный характер. В реакции принимают участие только частицы, находящиеся на границах соседних кластеров, заполненных реагирующими веществами, а большая часть вещества оказывается отделена (блокирована) от частиц, которые необходимы для протекания реакции. Так как формирование кластеров в различных частях решетки происходит несинхронно, то усредненные по пространству реакции концентрации веществ совершают флуктуационные колебания с амплитудой, спадающей при увеличении размера решетки. Когда в систему на каждом шаге добавляется перемешивание, часть блокированных частиц случайным образом
5", дБ|—>—г~г—I—1—I—'—!—'—I—'—I—'—I—'—I—'—I—г
0.0 0.02 0.04 0.06 0.08 /МСБ"1
Рис. 7. Спектры мощности процесса х(Ь). Параметры и начальные условия те же, что и на рис. 6
перебрасывается из одних кластеров в другие. При этом нарушается однородность последних. Кластеры становятся как бы испещрены отдельными частицами разных типов. Как только переброшенные частицы попадают в «благоприятную» среду, начинается быстрый рост соответствующей концентрации в этой области решетки. Естественно, что из-за вероятностного фактора в таких перебросках чаще всего будут принимать участие частицы, усредненная концентрация которых выше в данный момент. Это приводит к дальнейшему росту концентрации доминирующего вида, пока количество частиц второго вещества, принимающего участие в бимолекулярном взаимодействии, не станет слишком низким. Таким образом, даже слабое перемешивание способствует появлению пространственного синхронизма и, как следствие, росту амплитуды и регулярности колебаний. При этом, как видно из сопоставления рис. 8, б и 8, в, решетка в целом периодически «окрашивается» в разные цвета, то есть ведет себя как единое целое.
Рассмотрим теперь зависимость динамики системы от параметра перемешивания. Для этого выберем значения параметров к = 1 и будем плавно увеличивать р от нуля до единицы. При значениях р близких к нулю, система ведет себя так же, как и без перемешивания. Глобальная динамика сходна с «зашумленным» устойчи-
Рис. 8. Пространственные «снимки» системы с перемешиванием для разных моментов времени г: 0 МСЗ (а); 900 МСЗ (б); 1800 МСЗ (в). Оттенками серого показаны частицы различных видов. Начальные условия и параметры системы те же, что и на рис. 6
г
1.0
0.8
3
0.2
0.0
0.4
0.6
0.06
0.08
0.04
0.02
-0.2
0.0 0.2 0.4 0.6
0.8
х 0.0
0.0 0.002 0.004 0.006
Рис. 9. Проекции «фазовых портретов» глобальных Рис. 10. Зависимость дисперсии для процесса
колебаний на решетке на плоскость X — У в зави- х(Ь) от параметра перемешивания р для решеток
симости от степени перемешивания: 1 - р = 0, 2 - размером 100 х 100 ячеек (линия 1), 200 х 200
р = 0.0028, 3 - р = 0.0056. Прямыми изображена ячеек (2) и 300 х 300 ячеек (3). Начальные кон-
проекция контура из инвариантных многообразий центрации и относительные скорости реакций
для системы среднего поля (2) выбраны такими же, как и на рис. 6
вым фокусом (рис. 9, линия 1). Начиная с р ~ 0.001, наблюдается возникновение и постепенный рост колебаний в системе. «Фазовый портрет» становится похожим на предельный цикл в динамических системах (рис. 9, линия 2). Этот процесс проиллюстрирован на рис. 10, на котором построена зависимость дисперсии процесса от параметра р для нескольких размеров решеток. Видно, что, во-первых, «амплитуда» колебаний практически не зависит от размера решетки, во-вторых, существует пороговое «бифуркационное» значение рождения колебаний и, в-третьих, дисперсия зависит от «бифуркационного» параметра линейно, что означает рост амплитуды колебаний пропорционально квадратному кореню от надкритичного значения параметра, характерный для суперкритической бифуркации Андронова - Хопфа в динамических системах. Таким образом, по своим особенностям динамика рассматриваемой стохастической системы напоминает рождение предельного цикла из состояния равновесия при потере им устойчивости через бифуркацию Андронова - Хопфа. При дальнейшем увеличении управляющего параметра происходит монотонный рост амплитуды колебаний, пока, наконец, «предельный цикл» не «влипает» в замкнутый контур, образованный многообразиями. На рис. 9 линия 3 соответствует подобной ситуации, когда «предельный цикл» располагается в непосредственной близости от контура. После этого «цикл» исчезает, и мы наблюдаем динамику, характерную для неустойчивого фокуса внутри замкнутого контура из многообразий: из любых начальных условий колебания сначала резко возрастают, а затем, достигнув некоторой максимальной амплитуды, прекращаются. Данную временную динамику можно наблюдать на рис. 6, б, построенном для системы с полным перемешиванием.
В работе исследуется пространственная динамика модели Монте-Карло для системы (4+1)-ЬЬУ и рассматривается влияние внешнего перемешивания. В процессе реакций одноименные частицы образуют на поверхности двумерной решетки катализатора фрактальные структуры - кластеры. Наблюдаемые в системе колеба-
Заключение
ния концентраций реагентов являются результатом граничного взаимодействия этих кластеров. Характеристики колебаний определяются формой, структурой взаимодействующих кластеров и их количеством. Законы распределения кластеров по размерам различны для решеток малого и большого размера. Для решеток малого размера характерно распределение по степенному закону, в то время как в больших решетках к степенному закону распределения кластеров по размерам добавляется экспоненциально спадающая зависимость. Такой вид распределения в больших решетках свидетельствует об одновременном протекании двух процессов на поверхности решетки: формирование и взаимодействие кластеров в малом пространственном масштабе и перемешивание и суперпозиция несинхронных колебаний различных частей решетки в больших масштабах. Добавление внешнего перемешивания в систему приводит к синхронизации колебаний различных областей решетки, что в свою очередь ведет к росту амплитуды и регулярности глобальных колебаний в системе. В зависимости от степени перемешивания можно наблюдать аналог бифуркации Андронова -Хопфа, реализующейся в динамических системах, в результате которой рождается зашумленный устойчивый «предельный цикл». Амплитуда «предельного цикла» растет с увеличением параметра перемешивания, пока значения средних концентраций не выходят за рамки пороговых значений, соответствующих в модели среднего поля инвариантным многообразиям, образующим в фазовом пространстве замкнутый контур. После этого колебания в системе прекращаются, а поверхность решетки оказывается «отравленной» одним или двумя компонентами реакций.
При добавлении в систему сильного перемешивания (p ^ 1) реакции протекают по всей поверхности катализатора, а взаимодействия на решетке становятся нелокальными. Несмотря на то, что при составлении уравнений среднего поля основным приближением является нелокальность взаимодействий, добавление перемешивания не приводит к динамике, описанной в рамках макроскопического подхода модели среднего поля: динамика системы как при наличии внешнего перемешивания, так и при его отсутствии соответствует поведению диссипативных моделей. Вопрос о причине такого расхождения остается пока открытым.
Авторы выражают свою признательность профессору Астеро Провате (Институт физической химии, научный центр «Demokritos», Афины, Греция) за проведенные плодотворные обсуждения и консультации по материалам данной работы.
Библиографический список
1. Ertl G. Oscillatory kinetics and spatiotemporal self-organization in reactions at solid surfaces // Science. 1991. Vol. 254. P. 1750.
2. Wintterlin J.Scanning tunneling microscopy studies of catalytic reactions // Adv. Catal. 2000. Vol. 45. P. 131.
3. Ertl G., Norton P.R. and Rustig J.Kinetic oscillations in the platinum-catalyzed oxidation of CO // Phys. Rev. Lett. 1982. Vol. 49. P. 177.
4. Imbihl R. and Ertl G. Oscillatory kinetics in heterogeneous catalysis // Chem. Rev. 1995. Vol. 95. P. 697.
5. Shvartsman S.Y., Schutz E., Imbihl R. and Kevrekidis I.G. Dynamics on microcomposite catalytic surfaces: The effect of active boundaries // Phys. Rev. Lett. 1999. Vol. 83. P. 2857.
6. Voss C. and Kruse N. Chemical wave propogation and rate oscillations during the NO2 /H2 reaction over Pt // Ultramicroscopy. 1998. Vol. 73. P. 211.
7. Slinko M., Fink T., Loher T., Madden H.H., Lombardo S.J., Imbihl R. and Ertl G. The NO+H2 reaction on Pt(100) - steady-state and oscillatory kinetics // Surface Science. 1992. Vol. 264. P. 157.
8. Hartmann N. and Madix R.J. Dynamical rearrangements of the (2 x 1) O adlayer during CO oxidation on Cu(110) // Surface Science. 2002. Vol. 516. P. 230.
9. Ziff R.M., Gulari E., Barshad Y. Kinetic phase transitions in irreversible surface-reaction model // Phys. Rev. Lett. 1986. Vol. 56. P. 2553.
10. Brosilow B.J., Gulari E., ZiffR.M. Boundary effects in a surface reaction model for CO oxidation // J. Chem. Phys. 1993. Vol. 98. P. 674.
11. Zhdanov V.P. Surface restructuring and kinetic oscillations in heterogeneous catalytic reactions // Phys. Rev. E. 1999. Vol. 60. P. 7554.
12. Zhdanov V.P. Surface restructuring, kinetic oscillations, and chaos in heterogeneous catalytic reactions // Phys. Rev. E. 1999. Vol. 59. P. 6292.
13. Voss C., Kruse N.Field ion microscopy during an ongoing surface reaction: NO/H2 on Pt // Applied Surface Science. 1994. Vol. 87/88. P. 127.
14. Nicolis G. and Prigogine I. Self-organization in Nonequilibrium Systems. New York.: Wiley, 1977.
15. Albano E.V. Monte Carlo simulations of surface chemical reactions: Irreversible phase transitions and oscillatory behaviour // Computer Physics Communications. 1999. Vol. 121-122. P. 388.
16. Albano E.V. and Marro J.Monte Carlo study of the CO- poisoning dynamics in a model for the catalytic oxidation of CO // J. Chem. Phys. 2000. Vol. 113. P. 10279.
17. Tammaro M. and Evans J.W. Chemical diffusivity and wave propagation in surface reactions: lattice-gas model mimicking CO-oxidation with high CO-mobility // J. Chem. Phys. 1998. Vol. 108. P. 762.
18. Liu D.J. and Evans J.W. Symmetry-breaking and percolation transitions in a surface reaction model with superlattice ordering // Phys. Rev. Lett. 2000. Vol. 84. P. 955.
19. De Decker Y., Baras F., Kruse N. and Nicolis G. Modeling the NO + H2 reaction on a Pt field emitter tip: Mean-field analysis and Monte-Carlo simulations // J. Chem. Phys. 2002. Vol. 117. P. 10244.
20. Provata A., Nicolis G. and Baras F. Oscillatory dynamics in low dimensional lattices: A lattice Lotka-Volterra model // J. Chem. Phys. 1999. Vol. 110. P. 8361.
21. Tsekouras G.A. and Provata A. Fractal properties of the lattice Lotka-Volterra model // Phys. Rev. E. 2002. Vol. 65. art. no 016204.
22. Frachebourg L., Krapivsky P.L. and Ben-Naim E. Spatial organization in cyclic Lotka-Volterra systems // Phys. Rev. E. 1996. Vol. 54. P. 6186.
23. Efimov A., Shabunin A., Astakhov V. and Provata A. Chaotic dynamics of chemical reactions in low-dimensional substrates: Mean-Field and Monte-Carlo approaches // Изв. вузов. Прикладная нелинейная динамика. 2003. Т. 11, № 2. С. 72.
Саратовский государственный Поступила в редакцию 12.01.2005
университет
FORMATION AND EVOLUTION OF THE SPATIAL STRUCTURES IN THE SYSTEM OF CHEMICAL REACTIONS ON THE CATALITYC SURFACE: MONTE CARLO SIMULATION
A.V. Efimov, A.V. Shabunin
The cluster formation in the cyclic (4+1)-Lattice - Lotka-Volterra model is studied by Kinetic Monte Carlo simulations on a square lattice support. The features of cluster size distribution, spatial autocorrelation function and other dependences of the spatial dynamics of the system are under consideration. The role of cluster formation process and it effect on the systems dynamics is studied in this work. We show that the external mixing added to the initial scheme leads to the periodic self-oscillations appearance.
Ефимов Антон Викторович - окончил физический факультет СГУ по специальности радиофизика и электроника (2004). Работает инженером кафедры радиофизики и нелинейной динамики, с осени 2005 года - аспирант кафедры. Автор нескольких научных публикаций. Область научных интересов: нелинейная динамика, хаос, образование структур, хаотическая синхронизация, ее количественный анализ, клеточные автоматы. Е-таП^тоу@сЬао8.88и.гиппе1;.ги
Шабунин Алексей Владимирович - окончил Саратовский государственный университет (1990). Доцент кафедры радиофизики и нелинейной динамики СГУ, кандидат физико-математических наук (1998). Научные интересы - нелинейная динамика, теория колебаний, синхронизация и управление хаосом. Автор более 40 научных публикаций. Е-таП:а1ехеу@сЬао8.88и.гиппе1;.ги
