CC BY
46
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ОПТИМАЛЬНОЕ УПРАВЛЕНИЕ
УДК 534; 517.925.41/42
О БИФУРКАЦИЯХ В ДИНАМИЧЕСКОЙ СИСТЕМЕ ЧУМАКОВА - СЛИНЬКО © 2011 г. В.И. Потапов
Норильский индустриальный институт
PotapovVI@norvuz. т
Поступила в редакцию 07.02.2011
Проведено качественное исследование системы Чумакова - Слинько, адекватно описывающей химические процессы в реакции окисления молекулярного водорода на поверхностях металлических катализаторов. С помощью программного средства Ма&САС выявлены всевозможные перестройки фазовых портретов динамической системы при изменении значений ее параметров. Показано, что в такой нелинейной системе реализуется бифуркация Андронова - Хопфа рождения предельного цикла из седло-фокуса (1, 2), а также бифуркация появления сложного предельного множества фазовых траекторий (аттрактора Чумакова - Слинько). Наблюдалась бифуркация удвоения периода цикла в окрестности седло-фокуса (1, 2) с отрицательными седловыми величинами.
Ключевые слова: система Чумакова-Слинько, состояния равновесия, циклы, бифуркации, хаос, осцилляции.
Введение
Одной из простейших схем реакции окисления молекулярного водорода на поверхности никелевого или платинового катализатора [1] является:
—
(I) Н7 + 2М 2НМ,
—К2>
(II) О + 2М 20М,
^12-
(III) 2НМ + ОМ Кз >3М + НО, (1)
(IV) Н2 + ОМ К4 > М + НО,
—К5>
(V) М п + ОМ (М п О) + М.
< К_5
Первые две стадии обратимые и представляют собой адсорбцию молекул водорода и кислорода, третья основная необратимая стадия
- окисление водорода по рекомбинационному механизму Лэнгмюра - Хиншельвуда, четвертая необратимая стадия - окисление водорода по ударному механизму Или - Ридила. Пятая обра-
тимая стадия - растворение кислорода в приповерхностном слое катализатора, её называют буферной [2-4], то есть промежуточный комплекс (М п О) не участвует в основной стадии
(III). Приповерхностный слой (подложка) выступает в качестве «депо» для сброса излишка адсорбированного кислорода. Иначе идет обмен поверхностными и приповерхностными центрами адсорбции через окисление.
В проточном реакторе, где происходят эти стадии, поддерживаются постоянными концентрации молекул водорода и кислорода, отводится Н2О, соблюдается постоянство температуры в газовой фазе. В этих условиях имеем открытую гетерогенную систему, которой на основании закона действующих масс соответствует математическая модель:
= к • (1 — х — х2) — * х| — 2к3 • х| • х2,
■ х2 = к2 * (1 — Х1 — х2)2 — к—2 * х2 — к3 * х12 * х2 — к4 * х2> (2)
хз = к5 * х2 * (1 — х3) — к_5 * х3 * (1 — х — х2),
где х1; х2, х3 - безразмерные концентрации (степени покрытия поверхности) промежуточных
веществ НМ, ОМ, (МО) соответственно, кг-, к_г- - константы скоростей стадийных реакций (I)-(V), х4 = 1-х1 -х2, х5 = 1 -х3 - степени покрытия поверхности М и подложки Мп катализатора свободными центрами адсорбции соответственно.
Математическая модель (2) представляет собой нелинейную динамическую систему с 8 параметрами, компьютерное исследование которой показало существование одного состояния равновесия типа топологический узел с тремя отрицательными характеристическими корнями. Например, положим
к = 1.5, к—1 = 0.008, к2 = 20,к-2 = 0.02,
к3 = 100, к4 = 20, к5 = 0.0024; к_5 = 0.019, х1 = 0.1, х2 = 0.5, х3 = 0.
Используя MathCAD, находим координаты состояния равновесия О1 (0.0261245; 0.3669598; 0.0709552) и корни характеристического уравнения
=-0.0124121 Х2 =-44.5543457 ^ =-3.8053669
Замечаем, что точка О] является топологическим устойчивым узлом в Я3 с ведущим одномерным многообразием (соответствующим ^2) и неведущим двумерным многообразием 1¥^ (соответствующим Х2 и Х3) [5]. В этом случае динамика достаточно простая.
Если допустить, что константы скоростей основных стадий (III), (IV) схемы (1) экспоненциально зависят от концентраций поверхностных компонент [4, 6]:
к3 = к30е^3'х2, к4 = к40е^4'х2-Ц5'Х3, (3)
то решения динамической системы демонстрируют богатое поведение [1, 7]:
х = к (1 — х — х2) — к_2 х2 — 2к3 ^е ^3 2 х2 х2,
х2 = к2(1 - х1 - х2)2 - к40е ^х2 ^х3 х2 -
<
— кое ^3 2х2 х2 — к_2х2, хз = е(х(1 -х3)-ох3(1 -х -х2)),
где е = к5 - малый параметр, а а = к—5 / к5 -отношение констант скоростей 5-й стадии схемы (1). Систему (4) мы называем системой Чумакова - Слинько.
В широкой области параметров к , к у системы (4) имеется единственный устойчивый
предельный цикл, окружающий неустойчивое состояние равновесия в 3-мерном фазовом пространстве. Убеждаемся, что в симплексе реакции (1):
£ = {(х1, х2’ хз)1 х1 ^ 0 х2 ^ 0 1 — х — х2 ^ о, о ^ х ^ 1} дивергенция векторного поля системы (4)
,. . . „ дх, дх1 дх-,
х, х2, х3) = т~ + т— + — = дх дх2 дх3
= —2к^ (1 — х — х2) — 2к_^ х — 4к3 де 2 х ^ —
— 2к2(1 — х1 — х2) — 2к_2х2 + к40Ц4Є ^4 2 3 х2 — (6)
— к Є_Ц4х2 —^5х3 — к Є~^3х2 х2 +
к40е к30е х1 ^
+ ц3к30ех^ х2 — єх2 — ає(1 — х1 — х2)
всегда отрицательна. Это означает, что система (4) диссипативная, фазовые траектории не покидают симплекс (5), ее фазовый поток сжимается. Это было справедливо и для гетерогенной каталитической реакции Яблонского окисления угарного газа на платиновой поверхности [3, 8].
Устойчивость состояний равновесия системы Чумакова - Слинько и предельные циклы
Исследуем устойчивость состояния равновесия (стационарного решения) системы Чумакова - Слинько (4). Замечаем, что в общем случае явные значения компонент стационарного решения найти невозможно, поэтому решаем нелинейную трансцендентную систему для определения состояний равновесия с помощью пакета МаШСАБ при каждом наборе параметров [8].
Положим к = 15, к_2 = 0.008,
к2 = 20, к_2 = 0.02,к30 = 100,к40 = 20, є = 0.0024 а = 7.88, х = 0.1,х2 = 0.5,х = 0,
ц3 = 2, ц4 = 7, ц5 = —12.
Находим координаты состояния равновесия (4) О(00631010.6287730.205696. Далее строим матрицу Якоби линеаризованной системы (4) и решаем характеристическое уравнение. В результате получаем
^ = —7.801106, Х2 = —0.060694 0.275525,
Я3 = —0.0606940.275525, и предельный цикл не наблюдается.
Уменьшаем ц5 до -12.320462, находим стационарную точку О (0.064321 0.621117,
Рис. 1. Предельный цикл в системе Чумакова-Слинько при ^5 = -12.320463
„<1>
Рис. 2. Автоколебания адсорбированных атомов водорода
0.20036!), и так как \= -7.919414Х2 =
= -1.42678х 10-7 + 0.289214, Х3 = -1.42678Х
х10-7 -0.289214, то состояние равновесия О становится слабо устойчивым фокусом. Замечаем, что Яе Х2 3 « 0, то есть в интерпретации
А. Тондла [9] состояние равновесия «близко» к устойчивому пространственному центру.
Далее при уменьшении Ц5 до -12.320463 состояние равновесия
о (0.0643210.6211170.20036^) становится седло-фокусом, так как
^ = -7.919414 Х2 = 4.956545<10-8 +
+ 0.2892144,Х3 = 4.956545<10-8 -0.289214.
Таким образом, произошел переход характеристических корней через мнимую ось: реализуется бифуркация Андронова-Хопфа рождения устойчивого предельного цикла из седло-фокуса. Конфигурация цикла показана на рис. 1, а динамика адсорбированных компонент - на рис. 2 (колебания с периодом Т = 75).
При дальнейшем уменьшении параметра ц5 до -13 состояние равновесия
О (0.0668560.6052210.189769
остается седло-фокусом типа (1, 2), для которого Хх = -8.167587^2 = 0.1342891- 0.272772,
Х3 = 0.134289-0.272772.
Динамика остается прежней, но с периодом Т = = 65.
Далее, при Ц5 = -13.911121952 мнимая
часть характеристических корней седло-фокуса становится малой величиной
\ = -8.494498Х2 = 0.323481 7.798317 10-6 г, Х3 = 0.32348-7.798317<10-6г.
Вновь с помощью программного средства МаШСАБ решаем систему (4) и замечаем, что предельный цикл сохраняется, как на рис. 1. Динамика поведения адсорбированных компонент НМ аналогична (рис. 2), но период колебаний становится еще меньше Т = 55.
При ц5 = -13.911121951 состояние равновесия 01 становится седлом, для которого \ = -8.494498, Х2 = 0.323489, Х3 = 0.323473. Однако конфигурация цикла не меняется -остается такой же как, на рис. 1.
Рис. 3. Седло-узловой цикл в системе Чумакова-Слинько
Рис. 4. Релаксационные колебания свободных центров адсорбции на поверхности катализатора
При дальнейшем уменьшении параметра ц5 вплоть до значения -28 состояние равновесия О остается седлом. В симплексе реакции (5) по-прежнему сохраняется предельный цикл, как на рис. 1. Все промежуточные концентрации испытывают колебания, как на рис. 2, но с периодом Т = 30 при и5 = -28.
Возьмем следующие значения параметров из
[1]:
к = 0.005, £_1 = 0.001,к2 = 0.109, к_2 = 0.002, к30 = 500, к40 = 1.5, е = 0.00065 а = 10,
х1 = 0.1,х2 = 0.5,х3 = 0, и3 = 80,и4 = 10,и5 = 0, при этом стационарное состояние такое: (0.1733023, 0.6919165, 0.3392199), и характеристические корни равны
X = -0.0013258 К = -0.0250664 Х3 = 2.9504152<10-11.
Замечаем, что один из корней практически равен нулю, когда критическое значение параметра и5 = 0, а два других отрицательные, то есть имеем негрубое седло с вырождением. При ц5 < 0 X] и Х2 остаются отрицательными, Х3
становится положительным и седловая величина ст = шах(Х1,\2) + ^3 принимает отрицательное значение. Тогда в симплексе реакции (1) существует единственная периодическая орбита (рис. 3). Мультипликаторы цикла удовлетворяют соотношению |и • ц21 = 1.
Все адсорбированные компоненты испытывают колебания (рис. 4).
Покажем, что при критических значениях параметров к = 0.13582426к40 = 2 [10] в
системе (4) рождаются два цикла. Пусть параметры принимают значения:
к = 0.13582426 к_! = 0.008, к2 = 20, к_2 = 0, к30 = 10С, к^ = 2, е = 0.0024 а = 7.88, X = 0.1, х2 = 0.5, X = 0, Из = 30, ц.4 = 12, ^5 = -10.
При этом стационарное решение имеет вид (0.3252505, 0.595804, 0.4892084) и X =
= -1.9004282 Х2 = -0.0001976ъ0.0050586/, Х3 = -0.00019760.0050586/ (устойчивый узел-фокус). Далее интегрируем систему (4) методом Рунге - Кутта с различными начальными условиями и получаем два цикла (рис. 5). Удвоение цикла соответствует мульти-
пликатору ц = +1. Зеленая траектория скручивается по спирали к состоянию равновесия. Красная траектория движется по циклу, а черная удаляется и вновь приближается к предельному циклу (рис. 5).
Фазовая точка системы Чумакова - Слинько совершает несколько осцилляций вокруг красного цикла (рис. 5) и затем совершает колебания по черному циклу длительное время (точнее, по черной замкнутой орбите). Таким образом, красный цикл является неустойчивым, а черный - устойчивым. В этом случае существуют два цикла, как в системе Ван-дер-Поля с нелинейностью пятого порядка [8].
Сразу появляются автоколебания свободных центров адсорбции с красным начальным усло-
вием (рис. 6а) с периодом Т = 1000 и медленная в течение t = 1200 нерегулярная осцилляция при движении по черной (рис. 6б) траектории вокруг красного цикла. Затем в течение t = 6000 идут осцилляции черной траектории вокруг предельного цикла, и происходит выход на новую черную орбиту с временным удалением от неустойчивого красного цикла (рис. 6б).
С момента t = 7300 начинается движение по черной траектории вдали от красного цикла с возвратом в окрестность красного, то есть идут релаксационные колебания (рис. 6б). Выход на стационар по зеленой траектории очень медленный.
Замечаем, что идет медленный выход в течение t = 1200 динамической компоненты X
Рис. 5. Неустойчивый (красный) и устойчивый (черный) предельные циклы в системе Чумакова-Слинько (а)
Х4
У4
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 1 -10
X1)
(б)
74
Т 4
Рис. 6. Колебания свободных центров адсорбции на поверхности катализатора
(степени покрытия поверхности атомами водорода) на черную орбиту (рис. 7б), спирально раскручивающуюся вокруг красного цикла в течение ? = 6000. Затем с момента ? = 7300 фазовая траектория выходит на устойчивый черный цикл, удаляющийся и возвращающийся в окрестность красного неустойчивого цикла (рис. 7б). Идет очень медленный выход на стационар с зеленого начального условия (рис. 7б).
Структура аттрактора Чумакова - Слинько
С помощью программного средства МаШСАБ [8] подтвердим наличие в системе Чумакова - Слинько (4) бифуркации удвоения
периода цикла в окрестности седло-фокуса с отрицательными седловыми величинами. Для этого берем из [10] критические значения параметров
к = 0.145, к_! = 0.008,к2 = 20,к_2 = 0, к30 = 10С, к40 = 1.9 3 3 4 е = 0.0024 а = 7.88,
X = 0.1, х = 0.5, х = 0,
И3 = 30, И4 = 12, И5 = -10, находим стационарную точку
(0.3414096, 0.5783812, 0.4778317) и X = -1.9192539,
Х2 = 0.0004210.0051248/,
X = 0.000421 0.0051248.
(а)
(б)
X ^
Рис. 7. Колебания адсорбированных атомов водорода
Рис. 8. Многообходный предельный цикл в системе Чумакова-Слинько
Далее решаем систему (4) методом Рунге-Кутта при различных начальных условиях, и в симплексе реакции (1) в окрестности состояния равновесия рождается цикл (рис. 8).
Ему соответствуют колебания с удвоением периода (рис. 9) и мультипликатором и = -1.
Кинетика окисленных центров подложки испытывает сложные колебания с удвоением периода Т = 2000 при выходе из красных, зеленых и черных начальных условий, а из синих - медленно удаляется от состояния равновесия (рис. 10) с нарастанием амплитуды колебаний, т.е. соответствующее состояние равновесия неустойчивое, окруженное многообходным предельным циклом (рис. 8).
Покажем, что сценарий Фейгенбаума [8, 11]
- последовательного удвоения периода цикла -
приводит к бифуркации сложного предельного множества в системе (4).
В этом случае берем критические значения параметров:
к = 0.144, = 0.008,к2 = 20, к_2 = 0,к30 = 100,
к40 = 1.92, є = 0.0024 а = 7.88, хг = 0.1, х2 = 0.5, х3 = 0,
Ц = З0, Ц4 = 12, Ц = -10,
находим стационарную точку (0.3378532,
0.5824988, 0.4813541) и
= -1.9130324,
Х2 = 0.000281- 0.0051073/,
Х3 = 0.0002810.0051073/.
Результат численного решения системы (4) представлен на рис. 11.
Рис. 9. Колебания свободных центров адсорбции на поверхности катализатора с элементами релаксации
Рис. 10. Изменения концентрации окисленных центров подложки
Рис. 11. Аттрактор Чумакова - Слинько
(а)
(б)
Рис. 12. Колебания адсорбированных атомов водорода
Детальное компьютерное исследование динамики адсорбированных компонент показало, что компонента НМ с красным начальным условием длительное время (0 < t < 4000) совершает не] б улярные колебания, но затем выходит на незатухающие устойчивые колебания (рис. 12а) с постоянной амплитудой и периодом Т= 1000.
Динамическая переменная х с синим начальным условием медленно выходит на колебания (рис. 12а) с нарастающей амплитудой (стационарное состояние системы Чумакова -Слинько в этом случае неустойчиво).
Движение по зеленой траектории до момента t = 4000 идет с нарастающей амплитудой (рис. 12б) и мягко переходит на орбиту постоянного радиуса и периода Т = 2000. Двигаясь по черной траектории длительное время (рис. 12б) 0 < t < 630(0 система (4) совершала нерегулярные колебания и затем мягко перешла на автоколебательный режим (рис. 12б) с постоянным периодом Т = 1000.
Заметим, что при кх = 0.144; к40 = 1.92 система Чумакова - Слинько длительное время находится в хаотическом состоянии и затем выходит на регулярный режим (рис. 12а, б).
Заключение
В настоящей работе рассматривались актуальные в области водородной энергетики вопросы математического моделирования процес-
сов гетерогенной каталитической реакции окисления водорода (1). Изучалась система Чумакова - Слинько с алгебраической (квадратичной и кубической) нелинейностью. С помощью компьютера проверено, что при всевозможных значениях её параметров стационарное состояние является пространственным топологическим узлом с тремя отрицательными характеристическими корнями. Ему соответствует устойчивая динамика поведения адсорбированных компонент НМ, ОМ и (М„О) с быстрым выходом на стационар без колебаний.
Однако если принять допущение, что энергия активации стадий (III) и (IV) реакции (1) линейно зависит от степеней покрытия поверхности и подложки катализатора атомами кислорода, то константы скоростей этих стадий будут изменяться по экспоненциальному закону (3). В этом случае приходим к системе (4) с трансцендентной нелинейностью. У такой системы в широкой области значений параметров имеется единственный устойчивый предельный цикл, окружающий неустойчивое стационарное состояние. Установлено, что в симплексе (5) реакции (1) дивергенция векторного поля системы (4) всегда отрицательна. Это означает, что система (4) диссипативная, фазовые траектории не покидают симплекса (5), её фазовый поток сжимается. При этом в системе (4) Чумакова - Слинько выделена широкая область параметров, в которой реакция (1) длительное время находится в автоколебательном режиме. Это позволяет решать проблему
а
аккомодации тепловой энергии, выделяемой в рекомбинационных стадиях.
При уменьшении параметра и5 (т.е. при изменении константы скорости основной стадии
(IV) схемы (1)) происходит переход характеристических корней через мнимую ось, появляется седло-фокус типа (1, 2) с одномерным устойчивым многообразием Ж5 (соответствующим X < 0) и двумерным неустойчивым многообразием Ши (соответствующим комплексным корням Яе Х2 з > 0). В системе Чумакова - Слинько (4) реализуется бифуркация Андронова -Хопфа рождения цикла из состояния равновесия. Поскольку седловая величина ст = X + Яе Х23 отрицательна, то
рождается единственная устойчивая периодическая орбита [5, 12].
При этом химическая реакция (1) выходит на автоколебательный режим, все адсорбированные компоненты НМ, ОМ и (М„О) которой испытывают незатухающие колебания с постоянным периодом.
С помощью программного средства MathCad подбором параметров в системе Чумакова -Слинько (4) в окрестности устойчивого фокуса удалось наблюдать два предельных цикла (меньший неустойчивый и больший устойчивый), как в системе Ван-дер-Поля с нелинейностью 5-го порядка [8].
Также мы визуализировали сложное притягивающее инвариантное множество фазовых траекторий - аттрактор Чумакова-Слинько. При этом интегрировали систему методом Рун-ге - Кутта на больших временных интервалах и заметили, что все динамические переменные длительное время испытывают нерегулярные осцилляции, но затем мягко выходят на колебания с постоянными амплитудами и периодами.
Подобные исследования ранее проводились Г.А. Чумаковым и М.Г. Слинько [1, 4, 6, 7]. Здесь мы с помощью прикладного программного пакета MathCad детализировали и визуализировали колебательные процессы в 3Б-графике. Также подтвердили наличие сценария Фейген-баума - удвоения периода цикла в системе Чумакова - Слинько (используя критические значения параметров из работы [10]), приводящего к рождению сложного притягивающего множества фазовых траекторий (аттрактора Чумакова - Слинько).
Динамика аттрактора в течение длительного времени носит нерегулярный характер, а затем выходит на автоколебательный режим. В этом случае говорят о переходном хаосе.
Выражаю благодарность Г.А. Чумакову, любезно предоставившему свои новые работы по слабо устойчивой динамике в трехмерной кинетической модели каталитического окисления водорода, а также студенту Норильского индустриального института Е.Ю. Радченко за активное участие в компьютерном исследовании системы Чумакова - Слинько.
Огромная благодарность профессору А.Д. Морозову за полезные замечания и исправления рукописи.
Список литературы
1. Чумаков Г.А., Слинько М.Г., Беляев В.Д. // ДАН СССР. 1980. Т. 253.
2. Eigenberger G. // 4th Intern. Symp. Chem. React. Engng. Heidelbery, 1976.
3. Быков В.И., Яблонский Г.С., Ким В.Ф. Об одной простой модели кинетических автоколебаний в каталитической реакции окисления СО // ДАН СССР. 1978. Т. 242. В. 3.
4. Slinko M.M., Slinko M.G. // Catalysis Rev. Sci. Engng. 1978. V. 17(1).
5. Methods of qualitative theory in nonlinear dynamics. Part 1 and 2 / L. P. Shilnikov, A.L. Shilnikov, D.V. Turaev, L.O. Chua // World Sci., Singapore. V. 5. 1998 and 2001 [Перевод: Методы качественной теории в нелинейной динамике. Часть 1. М.-Ижевск, НИЦ РХД, ИКИ, 2004. Часть 2. 2009].
6. Слинько М.М., Слинько М.Г. Автоколебания скорости гетерогенных каталитических реакций // Успехи химии. 1980. Т. 49. № 4.
7. Чумаков Г.А., Слинько М.Г. Кинетическая турбулентность скорости реакции взаимодействия водорода с кислородом на металлических катализаторах // ДАН СССР. 1982. Т. 266. В. 5.
8. Потапов В.И. Математические модели нелинейных динамических явлений. Их численный и качественный анализ. Норильск, 2005.
9. Тондл А. Нелинейные колебания механических систем. М.: Мир, 1973.
10. Chumakov G.A., Chumakova N.A. Weakly Stable Dynamics in a Three-Dimensional Kinetic Model of Catalytic Hydrogen Oxidation // Chemistry for Sustainable Development. 2003. 11. P. 63-66.
11. Guckenheimer J., Holmes Ph. Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Field. Springer, 1983 [Перевод: Нелинейные колебания, динамические системы и бифуркации векторных полей. М.-Ижевск: ИКИ, 2002].
12. Ilyashenko Yu., Li W. Nonlocal Bifurcations. American Mathematical Society, 1999.
ON BIFURCATIONS IN DYNAMICAL SLINKO - CHUMAKOV SYSTEM
V.I. Potapov
Qualitative research has been carried out of Slinko - Chumakov system which adequately describes the chemical processes in the reaction of oxidation of molecular hydrogen on the metal catalyst surfaces. Using MathCAD software, all possible phase portraits of the dynamical system have been revealed while changing its parameters. It has been shown that in such a nonlinear system the Andronov-Hopf bifurcation is realized: birth of a limit cycle from saddle-focus [1, 2], as well as bifurcation of the complex limit set of phase trajectories (Slinko - Chumakov attractor). A period-doubling bifurcation has been observed in the vicinity of the saddle-focus [1, 2] with negative saddle values.
Keywords: Slinko - Chumakov system, equilibrium states, cycles, bifurcations, chaos, oscillations.
