Содержание к диссертации
Введение
1 Обзор численных методов анализа автономных систем 15
1.1 Численное интегрирование жестких систем 15
1.2 Построение диаграмм стационарных решений методом продолжения по параметру 19
1.3 Исследование устойчивости стационарных решений 29
1.4 Одно- и двупараметрический анализ стационарных решений 33
2 Базовые алгоритмы пакета программ STEP 42
2.1 Алгоритм метода продолжения по параметру 43
2.2 Исследование устойчивости стационарных решений по методу Годунова - Булгакова в зависимости от параметра . 50
2.3 Численные схемы интегрирования Розенброка и Гира . 54
2.4 Построение бифуркационных кривых «овражным методом» . 57
3 Численное исследование математических моделей из приложений 68
3,1 Примеры модельных задач учебного курса 69
3.2 Математическая модель поверхностных процессов в реакции окисления СО на иридии 82
3.3 Математическая модель регуляции внутриклеточного биосинтеза холестерина 99
4 Организация численного эксперимента в пакете STEP 108
4.1 Ввод и вывод информации 108
4.2 Работа с моделью в диалоговом режиме 111
4.3 Раздел «Задача Коши» 116
4.4 Раздел «Нелинейные системы» 120
Заключение 125
Приложение 128
Список литературы 136
- Построение диаграмм стационарных решений методом продолжения по параметру
- Исследование устойчивости стационарных решений по методу Годунова - Булгакова в зависимости от параметра
- Математическая модель поверхностных процессов в реакции окисления СО на иридии
- Работа с моделью в диалоговом режиме
Введение к работе
При изучении многих прикладных задач важную роль играет вычислительный эксперимент, предоставляя возможность более глубокого понимания существа изучаемых процессов через исследование соответствующей математической модели. Кроме того, в ряде случаев моделирование процессов в лабораторных или натурных условиях сложно и дорого.
Для проведения численного исследования физического процесса необходимо:
- построить физическую модель процесса, т. е. найти закономерности, позволяющие судить о характере протекания процесса (это, например, схемы кинетических реакций, уравнения теплового баланса, равновесия, законы сохранения и т. п.);
- построить математическую модель, в нашем случае в терминах обыкновенных дифференциальных уравнений (ОДУ);
- найти решения системы дифференциальных уравнений, применяя методы численного интегрирования;
- провести численное исследование поведения решений в зависимости от параметров системы;
- сопоставить полученные результаты с известными качественными представлениями о характеристиках изучаемого физического процесса.
Как отмечается, например, в [73], достаточно полное аналитическое исследование автономных систем (в том числе и бифуркационный анализ) возможно, как правило, лишь для систем специального вида. Поэтому важна роль численных методов, составляющих основу организации численного эксперимента при изучении качественных характеристик поведения решения, и разработок программного обеспечения, использующего эти методы.
Востребованность разработок программного обеспечения такого рода состоит в том, что многие математические модели из приложений формулируются в терминах автономных систем дифференциальных уравнений. Из литературы известен ряд программных комплексов типа пакета STEP, реализующих численные алгоритмы качественного анализа автономных систем. В частности, много работ было посвящено изучению бифуркаций в гидродинамике (см., например [5]). Для изучения бифуркаций в биологических системах были созданы (НИИВЦ АН СССР, Пущино) такие программные комплексы, как LINLBF [98] и более поздние его модификации, L0CBIF [99] и CONTENT [101]. Можно также назвать такие широко известные пакеты программ, как DERPAR [83] и AUTO [85]. Для интегрирования задачи Коши применяются, например, комплекс программ MADONNA, см. http://www.berkeleymadonna.com/, пакеты GEAR [82] и DIFFSUB [87].
Опыт работы с приложениями показывает, что решение конкретных задач на практике оказывается затруднительным в случае применения готовых комплексов программ, являющихся, как правило, коммерческим продуктом, закрытым для доступа к текстам программ. В этом случае приме нение так называемого «черного ящика» может вызвать затруднения при организации численного эксперимента. Кроме того, в рамках используемого комплекса программ может оказаться проблематичным и привлечение новых алгоритмов.
На основе вычислительных алгоритмов, разработанных в Институте математики им. С. Л. Соболева СО РАН, было создано программное обеспечение численного исследования автономных систем, включающее в себя, в частности, оригинальный вариант метода продолжения для решения систем нелинейных уравнений, зависящих от параметра, определение устойчивости стационарных решений по методу Годунова - Булгакова в зависимости от параметра, алгоритмы проведения одно- и двупараметрического анализа стационарных решений автономной системы.
Постановка задачи
При изменении какого-либо параметра системы может произойти качественное изменение - бифуркация - установившегося режима при некотором значении параметра, которое будем называть критическим. Очевидно, что фазовый портрет изучаемой системы при этом меняется. Так, в случае, когда стационарное решение теряет устойчивость при критическом значении параметра, может возникнуть замкнутая траектория (предельный цикл) - бифуркация Андронова - Хопфа или, иначе, комплексная бифуркация. Другим возможным проявлением нелинейности модели является существование критического значения параметра, при котором возникает множественность решений.
Говорят, что разбиение пространства параметров Rm на области с качественно различными типами динамического поведения определяет бифуркационную диаграмму.
Качественные, топологические перестройки фазовых портретов динамических систем уравнений при непрерывном, плавном изменении параметров системы описываются теорией бифуркаций. Изложение основных методов теории бифуркаций и результатов, полученных в этой области, начиная с основополагающих работ А. М. Ляпунова, А. Пуанкаре и А. А. Андронова, можно найти, например, в работах [2], [3], [б], [7], [73].
При проведении численного эксперимента наибольший интерес представляет качественная картина, сохраняющаяся при изменении фазового портрета системы. Это прежде всего вопросы о характере режимов, которые устанавливаются в системе по завершении переходного процесса: число состояний равновесия, отвечающих стационарным режимам данной системы, и их устойчивость; возможность существования предельных циклов, отвечающих режимам периодических колебаний и т. п. Очевидно, что построе ниє бифуркационной диаграммы в совокупности с соответствующими фазовыми портретами позволит ответить на вопросы о возможных в системе динамических режимах и их качественных перестройках.
Предметом диссертационной работы являются вопросы организации численного эксперимента с помощью вычислительных средств, которые позволяют проводить качественное исследование автономных систем вида (0.0.1) в зависимости от параметров и систем нелинейных уравнений, не обязательно являющихся правыми частями автономных систем. Применение пакета программ STEP при проведении численного эксперимента предоставляет следующие возможности.
• Численное интегрирование задачи Коши для автономной системы с целью выхода на стационарный режим, либо с целью получения колебательного режима. Определение периода колебаний в режиме диалога и нахождение решения задачи Коши при изменении значения параметра в процессе счета,
• Построение связных ветвей стационарных решений в зависимости от выбранного параметра модели методом продолжения по параметру, с учетом возможного существования точек ветвления типа поворот.
• Численное исследование устойчивости найденных стационарных решений на основе теоремы Ляпунова об устойчивости по первому приближению.
• Обнаружение на диаграмме стационарных решений точек, в которых возможна бифуркация Андронова - Хопфа.
• Проведение двупараметрического анализа стационарных решений автономной системы. Цели и задачи работы
• Разработка пакета программ, предназначенного для численного исследования автономных систем ОДУ общего вида и систем нелинейных трансцендентных уравнений.
• Применение пакета STEP в качестве математического обеспечения учебного процесса в курсе «Инженерная химия каталитических процессов».
• Организация численного эксперимента с помощью пакета программ STEP для выявления качественных свойств конкретных математических моделей из приложений.
Научная новизна
В работе представлено программное обеспечение численного исследования автономных систем ОДУ общего вида, оформленное в виде пакета программ STEP. В пакете STEP использованы вычислительные алгоритмы, ориентированные на изучение систем произвольной размерности. Применение этих методов дает возможность эффективно организовать численное исследование математических моделей в зависимости от параметров. Кроме того, указанные вычислительные алгоритмы учитывают возможность проявления нелинейных эффектов, которые, как правило, присутствуют в математических моделях, описывающих «нелинейные» процессы (гистерезис, сильная параметрическая чувствительность, возникновение автоколебаний и т.д.).
Предложен и программно реализован новый алгоритм, расширяющий возможности численного построения бифуркационных кривых для автономных систем из п уравнений. Представленный в работе пакет программ STEP применялся для решения многих прикладных задач (см. [9], [26], [28], [33], [50], [57] - [68], [79], [88], [89], [91] - [95], [104]). Как показала практика, пакет легко адаптируется к моделям из различных приложений благодаря использованию в пакете алгоритмов универсального типа.
Применение пакета STEP в учебном процессе показано на примерах мо дельных задач вычислительного практикума учебного курса «Инженерная химия каталитических процессов». Вычислительный практикум, в качестве ф1 программного обеспечения которого используется пакет STEP, разработан на кафедре катализа и адсорбции НГУ.
Проведено численное исследование математической модели, описывающей поверхностные процессы, протекающие в реакции окисления СО на иридии. Полученные для этой модели результаты являются новыми.
Проведено численное исследование математической модели, описывающей регуляцию внутриклеточного биосинтеза холестерина, в зависимости от параметров модели и построена матрица параметрической чувствительности.
Структура диссертации
Диссертация состоит из введения, четырех глав, заключения, приложения и списка литературы. Общий объем диссертации составляет 150 страниц. Библиография содержит 108 наименований.
Во введении обоснована актуальность темы, сформулированы цели и задачи численного исследования автономных систем ОДУ, описывающих различные физические процессы, показана научная новиза. В первой главе дан обзор известных численных методов, которые используются для анализа нелинейных динамических моделей. В первом параграфе рассматривается проблема численного интегрирования задачи Ко-ши для автономной системы ОДУ, в частности, вопросы, связанные со свойствами устойчивости численной схемы. Во втором параграфе рассмотрены численные методы построения диаграммы стационарных решений автономных систем в зависимости от параметра, в третьем - методы исследования устойчивости стационарных решений. В четвертом параграфе даны схемы проведения одно- и двупараметрического анализа стационарных решений автономных систем.
Вторая глава диссертации посвящена изложению алгоритмов, которые являются базовыми в пакете программ STEP. К ним относятся вариант метода продолжения по параметру, предложенный С. И. Фадеевым [23], [56], метод Годунова - Булгакова для определения устойчивости стационарных решений [17], алгоритмы нахождения точек локальной бифуркации положений равновесия автономной системы (точек ветвления типа поворот и точек бифуркации Андронова - Хопфа) и построения бифуркационных кривых в плоскости двух параметров. Для численного интегрирования систем ОДУ предлагаются многошаговый метод Гира переменного порядка точности [82] и полунеявный метод Розенброка 2-го порядка, вариант которого изложен в [22].
В третьей главе рассматриваются задачи из приложений и из учебного курса «Инженерная химия каталитических процессов». В первом параграфе этой главы приводятся примеры модельных задач, которые подбирались по следующим работам, освещающим проблемы гетерогенного катализа: [1], [10], [11], [30], [31], [37], [42], [52], [53], [70]. Во втором параграфе рассмотрена математическая модель реакции окисления окиси углерода на иридии. В третьем - математическая модель регуляции биосинтеза холестерина в клетке.
В четвертой главе дается структура пакета STEP по разделам. В первом параграфе описан ввод и вывод данных и раздел «Создание модели». Во втором параграфе рассмотрены схемы проведения численного исследования моделей с применением пакета. Во третьем параграфе приводятся блок-схемы алгоритмов численного интегрирования раздела «Задача Ко-ши», дано описание основных подпрограмм и параметров метода расчета для каждого алгоритма. В четвертом параграфе дано описание структуры раздела «Нелинейные системы», включая описание основных подпрограмм и параметров метода продолжения.
В приложении к третьей главе, § 2 приведена физико-химическая модель поверхностных процессов в реакции СО + Ог на иридии и результаты экспериментальных исследований каталитических процессов в ходе реакции. В приложение к § 3 вынесена формулировка математической модели регуляции биосинтеза холестерина в клетке.
Практическая ценность
Приведенный выше список публикаций, свидетельствует о большой практике использования пакета STEP в приложениях. Пакет программ передан в ряд научно-исследовательских учреждений РАН, таких как Институт катализа им. Г. К. Борескова, Институт цитологии и генетики, (Новосибирск); Институт вычислительного моделирования (Красноярск); Государственная Академия тонкой химической технологии, Научно-исследовательский физико-химический институт им. Л. Я. Карпова, Институт химической физики им. Н. Н. Семенова, (Москва), где применяется для проведения чис ленного исследования прикладных задач.
Пакет программ STEP является составной частью учебного курса «Инженерная химия каталитических процессов», который ведется на факультете естественных наук НГУ. В качестве руководства для проведения численного исследования этого класса математических моделей и документации по пакету программ было выпущено учебное пособие [67]. Многооконный интерфейс дает возможность проводить исследование в режиме диалога, что делает пакет удобным инструментом проведения численного эксперимента как для научных работников, занимающихся численным анализом нелинейных проблем, так и для студентов математического, физического факультетов и факультета естественных наук.
Апробация работы
Материалы работы докладывались на Школе-семинаре «Динамика процессов и аппаратов непрерывной технологии», 1991, Яремча; Международной конференции «Современные проблемы прикладной и вычислительной математики», 1995, Новосибирск; Международной научно-методической конференции «Новые информационные технологии в университетском образовании», 1996, 1997 г.г., Новосибирск; XIII Международной конференции по химическим реакторам, 1996, Новосибирск; 2-м и 3-м Сибирском конгрессе по прикладной и индустриальной математике, 1996, 1998 г.г., Новосибирск; Международной конференции «Всесибирские чтения по математике и механике», 1997, Томск; XIV Международной конференции по химическим реакторам, 1998, Томск; Международной конференции посвященной академику С. К. Годунову, 1999, Новосибирск; XV Международной конференции по химическим реакторам, 2001, Хельсинки, Финляндия; Российско-голландском симпозиуме "Catalysis for Sustainable Development", 2002, Новосибирск; 2-й и 3-й Международных конференциях по биоинформатике "Bioinformatics of Genome Regulation and Structure", 2000, 2002 r.r., Новосибирск.
Автор выражает глубокую благодарность научному руководителю С.И.Фадееву, а также А.Ю.Березину, Ю.С.Волкову, Е.П.Волокитину, А.И.Воронину, В.И.Елохину, В.А.Лихошваю и А.В.Ратушному, принявшим участие в обсуждении работы.
Построение диаграмм стационарных решений методом продолжения по параметру
Рассмотрим систему из п нелинейных уравнений где х - n-мерный вектор переменных, а - скалярный параметр, fi(x,a),.. .,fn(x, а) - достаточно гладкие функции по совокупности аргументов. Обозначим через х{а) решение системы (1.2.1), зависящее от параметра а. В дальнейшем используются обозначения: J — fx(xta) — (п х п)-матрица Якоби системы, J = [/х(ж,а),/а(х, а)], - (п х (п + 1))-матрица, составленная из матрицы Якоби и вектор-столбца производных по параметру.
Согласно [73], диаграммой стационарных решений системы (1.2.1) называется график кривой х(а). Пусть (х ,а ) -точка на диаграмме стационарных решений, т. е. Точка (х ,а ) множества решений системы (1.2.1), для которой выполнено условие называется регулярной точкой.Рассмотрим случай, когда det(J) = 0 и rank(J) = п. Тогда, заменив один из столбцов матрицы J (допустим, fc-й столбец) последним столбцом матрицы J, можно добиться того, чтобы полученная квадратная матрица имела ранг п. В этом случае систему (1.2.1) можно представить в виде: где z = (х\}..., avj-i), а = xn+i, хк не принадлежит z. То есть, все компоненты вектора z, в том числе и а, будут функциями Хк. Роль параметра а теперь будет играть переменная Xk- Как показано в [67], в точке (а; , а ) справедливо равенство:
Если d2a(dx\ 0 - это означает, что при а а пара стационарных решений изчезает. Если же fta/dx1 0, то при а. а возникает пара стационарных решений. Такую особую точку (х ,а ) называют точкой ветвления типа поворот. В том случае, когда cPa/dxl 0, точка (х ,о ) -особая точка типа перегиб.
Если rank (/(а; , а )) п, то точку (ж , а ) называют сингулярной или существенно особой точкой. В дальнейшем мы будем предполагать, что диаграмма стационарных решений не содержит существенно особые точки.
Пусть при некотором значении параметра a — ao известно решение х системы (1.2.1), т. е. выполняется уравнение и определитель матрицы Якоби системы не равен нулю в некоторой окрестности точки (x,Qfo). Тогда, как следует из теоремы о неявных функциях, существует решение системы (1.2.1) в окрестности а = ао- Требуется построить зависимость х(а) для рассматриваемой системы в указанной окрестности по а.
Заметим, что задаваясь некоторым интервалом по а, допустим [а аг], мы, вообще говоря, заранее не знаем, в какой окрестности ао существует решение. Кроме того, зависимость х(а) может быть построена регулярным образом только там, где т. е. между двумя существенно особыми точками, если они существуют.
Рассматриваемая постановка задачи имеет следующий геометрический смысл. Введем в качестве универсального параметра длину дуги 5 пространственной кривой в (п + 1)-мерном пространстве {х, а). Элемент дуги ds определяется по формуле
По теореме о неявных функциях, графиком решения системы (1.2.1) в пространстве Rn+1, проходящим через точку (х,о о), является единственная гладкая пространственная кривая S, которая имеет параметрическое представление
Требуется построить гладкую пространственную кривую Е, выходящую из точки (xQ,ао), для значений а из некоторой окрестности ао Є [а\, а ]- Для построения связных ветвей стационарных решений системы (1.2.1) предлагается использовать метод продолжения по параметру, который дает возможность изучать зависимость поведения решения от параметров модели, и, в том числе, находить бифуркационные значения параметров.
Идея продолжения решения известна и используется в математике и механике довольно давно. В литературе метод продолжения фигурирует под самыми разными названиями: метод последовательных нагружений, метод гомотопии и др. (см., например, [44]). Метод продолжения реализует построение диаграммы стационарных решений как шаговый процесс по параметру с использованием информации о решениях, полученных на предыдущих шагах, для задания начального приближения на текущем шаге. Решение затем уточняется с помощью итерационной процедуры, например, методом Ньютона - Рафсона (см. [23], [56], [73], [74], [84]). Остановимся подробнее на различных вариантах реализации метода продолжения по параметру.
В методе Д. Ф. Давиденко [19] продолжение решения системы (1.2.1) формулируется как задача Коши для системы ОДУ с начальными условиями X(CXQ) х, где J = [fx{x,a)] - матрица Якоби, det(J) ф 0 в интервале изменения параметра а, а Є [аь г], х(ао) - стартовое решение системы (1.2.1) при a = CXQ. ДЛЯ построения решения х(а) используется какая-либо численная схема интегрирования задачи Коши.
В задаче (1.2.3) предполагается, что в рассматриваемом интервале значений параметра а, а Є [а аг], определитель det(J) матрицы Якоби системы (1.2.1) отличен от нуля. Пусть при некотором а — а нарушается условие теоремы о неявных функциях, а именно, матрица Якоби становится вырожденной, т.е. det(/x(j;(a+),Qr ) — 0. Это означает, что х(а) многозначная функция, а точка (а:(а ),а: ), принадлежащая , разделяет ветви однозначности. Для того, чтобы продолжать решение и в окрестности точек «поворота» гладкой пространственной кривой S, реализуется идея, высказанная Кубичеком [73], которая позволяет построить интегральную кривую задачи (1.2.3).
Обозначим xn+i = a, z = (жі,. ..,хтхп+\). Будем рассматривать неиз вестные х\ хп и параметр а как равноправные переменные, т. е., в качестве параметра продолжения можно выбирать либо а, либо одну из компонент вектора х. Это возможно сделать, поскольку, по предположению, ранг расширенной матрицы J = [fx{xya)}fa(x,a)] равен п и, следовательно, всегда можно указать такую компоненту х составного вектора z = (х, а), для которой матрица
Исследование устойчивости стационарных решений по методу Годунова - Булгакова в зависимости от параметра
Рассмотрим общую схему проведения одно- и двупараметрического анализа стационарных решений автономной системы уравнений вида (0.0.1). Пусть сначала изменяется лишь один параметр а — р изучаемой системы, а остальные параметры фиксированы, т. е. рассматривается диаграмма зависимости стационарных решений от параметра а. Однопараметрический анализ состоит в исследовании диаграммы стационарных решений с целью выявления на ней бифуркационных точек. Допустим, что исследуемая диаграмма содержит критическое значение параметра а, отвечающее локальной бифуркации стационарных решений системы: точку ветвления типа поворот или точку бифуркации Андронова - Хопфа.
Сделаем теперь «активным» другой параметр /3 = рк, к ф j, и на плоскости двух параметров (а,/3), построим бифуркационную кривую (линию кратности или линию нейтральности), проходящую через найденную бифуркационную точку. По определению, кривая в параметрической плоскости (а,Р), состоящая из координат точек поворота, называется линией кратности, В этом случае бифуркационная кривая разбивает (а,/3)-плоскость на области с определенным числом стационарных решений. Соответственно, кривая в параметрической плоскости (a,j3), состоящая из точек бифуркации Андронова - Хопфа, называется линией нейтральности. Построение бифуркационных кривых является задачей двупараметрического анализа.
Однолараметрический анализ системы Рассмотрим нелинейную систему из п уравнений вида где х є IVх, а Є R активный параметр модели. Кратко остановимся на численных методах нахождения критических значений параметра а системы (1.4.1) на некотором интервале изменения параметра, предложенных в работах [71], [73], [84], [99]. Общая схема этих методов состоит в присоединении к исходной системе так называемых условий вырождения [71] и решении полученной расширенной системы уравнений. Способ выбора условий вырождения обуславливает выбор метода для определения критического значения параметра. Рассмотрим различные варианты.
Присоединив к исходной системе уравнений (1.4.1) условие равенства нулю якобиана системы, получим систему из (п + 1) уравнений относительно неизвестных (х\}... ,хп,а). Равенство det(J) = 0 нарушает условие разрешимости системы (1.4.1) относительно х, которое содержится в известной теореме о неявных функциях. Таким образом, из системы (1.4.2) определяются координаты точки ветвления типа поворот. Решение этой системы может быть найдено с помощью итерационного метода Ньютона - Рафсона. Отметим, что для систем больших размеров возникает проблема нахождения начального приближения для стартового решения системы (1.4.2) в итерационной процедуре. 2. Условие det(J) = 0 можно записать иначе: где v - не равный нулю вектор. Вектор v определен с точностью до последнего множителя, поэтому к этому соотношению добавляется условие «нормировки» вектора v: при заданном к. Эти уравнения образуют систему из (2д + 1) нелинейных уравнений относительно (2п + 1) неизвестных [х, v, а). Решение системы также можно искать с помощью метода Ньютона. Построение матрицы Яко-би для такой системы требует вычисления вторых производных функций fi, г = 1,...,п, или использования соответствующих разностных формул.
Другой метод определения координат точек поворота основан на использовании условия где А; выбирается таким образом, чтобы матрица J , образованная из матрицы J = [fx(x,o ),fa(xfa)] вычеркиванием к-го столбца, была невырождена. Производную [da/dxk) можно найти, решая систему линейных алгебраических уравнений
В пакете программ STEP (см. [67]) на каждом шаге продолжения, начиная со 2-го, решение на элементарном отрезке по текущему параметру \ь — rrjt представляется в виде эрмитовой кубической параболы. Это позволяет, продифференцировав кубическую параболу по Xk, на каждом интер-вале проверять условие (1.4.3) при построении диаграммы стационарных решений системы (1.4.1).
Приведем теперь краткое описание алгоритмов определения координат точки бифуркации Андронова - Хопфа. 1. В работе [73] предлагается следующая схема нахождения двух комплексно-сопряженных чисто мнимых собственных чисел матрицы Яко-би: Составим собственный многочлен Р{Х) матрицы Якоби и вычислим его коэффициенты, используя, например, метод Леверье [35]. После нахождения значений коэффициентов ai,a2,... ,оп многочлен Р(Х) можно записать: где РП_2(А) - многочлен степени (п-2). Для того, чтобы найти в, перепишем Р{Х) в форме где коэффициенты pi, ..., рп-2, А, В могут быть вычислены по рекуррентным формулам
Математическая модель поверхностных процессов в реакции окисления СО на иридии
В данном параграфе рассматривается математическая модель реакции окисления окиси углерода на иридии. Как известно (см., например, [76]), реакция окисления СО на металлах платиновой группы является модельной реакцией, которая позволяет более глубоко понять сущность физико-химических процессов, протекающих в системе «реакционная среда - катализатор». Особенно это относится к изучению реакций окисления на металлах в области высокого вакуума, где возможно эффективное применение широкого спектра физических методов исследования.
В работах [12], [13], [25], [80] изложены результаты экспериментальных исследований и теоретического анализа изучаемого процесса, на основе которых была построена детальная схема каталитического цикла поверхностных процессов в ходе реакций окисления на платиновых металлах (см. рис. 4.1, Приложение) и предложен механизм реакции. Экспериментальные результаты были получены с помощью современных методов фотоэлектронной спектроскопии: РФЭС, УФФЭС и масс- спектроскопии. В качестве образцов использовались монокристаллы и поликристаллическая фольга иридия.
Следуя работам [12], [13], [25], [80], приведем формулировку проблемы. Согласно концепции Борескова [11] в условиях воздействия реакционной среды происходит изменение химического состава и структуры поверхности катализаторов. В свете этой концепции авторы приходят к выводу о необходимости двуслойного («двухуровневого») рассмотрения адсорбированного слоя. Показано (см., например, [25]), что в ходе реакции на поверхности 1г могут быть реализованы пять типов промежуточных соединений - кислород (cvi - О = ZO) и СО («і - СО = ZCO) на исходной поверхности (ИП), кислород (а2 - О = (ZyOZ) О) и оксид углерода (а2 - СО = (ZvOZ)CO) на поверхностном окисле и собственно кислород в виде поверхностного окисла (Р - СО = (ZyOZ)), здесь Z, Zy - активные центры поверхностного и приповерхностного слоя иридия.
Предложенный механизм отражает образование и взаимодействие промежуточных соединений в исследуемой каталитической системе с учетом влияния изменения свойств поверхностного слоя (ИСПС) на адсорбционные и кинетические характеристики реагентов. Для каждой реакции приводят- ся константы скорости, которые были экспериментально определены [13] и затем уточнялись в процессе моделирования. Предэкспоненциальные множители являются стандартными для молекулярных взаимодействий [27].
Здесь обозначено: Яд - газовая постоянная, Т - температура. В таблице значения Еа приводятся в [ккал/моль], а Р(СО) и Р(02) - [Торр]. Пред-экспоненты для реакций как первого, так и второго порядка приведены к размерности [сек-1].
Первые четыре реакции представляют собой традиционный адсорбционный механизм, на основе которого можно качественно интерпретировать множественность стационарных состояний и критическое замедление скорости реакции в кинетической области (см., например, [24], [75], [76]). Однако для описания обнаруженной модификации поверхностных центров рассматривается более детальная модель, для которой принят ряд предположений (см. [25]). В силу этих предположений и с учетом двух уравнений баланса по активным центрам поверхностного и приповерхностного слоя иридия: была сформулирована математическая модель данной каталитической реакции в виде системы из 5-ти дифференциальных уравнений: где Р(02), Р(СО) - парциальные давления реагентов 02 и СО; х;, і = 1,..., 5, - безразмерные концентрации промежуточных веществ ZO, ZCO, (ZvOZ), (ZvOZ)0, (ZvOZ)CO соответственно; kj, j = 1,...,12, - константы скорости элементарных реакций изучаемого механизма.
Цель исследования математической модели (3.2.1) состояла в том, чтобы найти области значений параметров модели, где существуют автоколебания, построить области множественности стационарных решений системы, а также получить информацию как о скоростях отдельных стадий, так и о суммарной скорости образования СОг:
С помощью пакета STEP была проведена серия численных экспериментов, направленных как на воспроизведение данных численного исследования модельной реакции, приведенных в [25], так и получения новых результатов [95]. К ним относится, в частности, обнаружение в узкой области варьирования параметров существования автоколебаний, а также областей множественности стационарных режимов при реальных физических условиях.
Как показано в работе, численное моделирование стационарной кинетики реакции дало удовлетворительное совпадение с экспериментальными данными во всем исследуемом диапазоне изменения параметров реакции:
В работе [25] приведены результаты расчетов при условиях Р(02)/Р{СО) да 9/1, Р(СО) = 1.1 Ю-7 (рис. 4.2 в Приложении). Расчеты проводились методом установления, а именно, для каждого значения температуры решалась задача Коши до выхода на стационар. Это стационарное решение бралось затем в качестве точки диаграммы стационарных решений. Построение диаграммы, таким образом, потребовало проведения большого числа расчетов и больших затрат машинного времени.
Численное исследование, проведенное по пакету STEP, для представленных на рис. 4.2 условий, было организовано следующим образом. Решив задачу Коши с нулевых начальных данных при тех же соотношениях давлений, Т = 700 К, мы выходим на стационарное решение, которое затем берется в качестве стартового в методе продолжения по параметру. Методом продолжения по параметру строится зависимость концентраций ХІ, і — 1,..., 5, от температуры на интервале [300 — 700 К], двигаясь по параметру в сторону уменьшения температуры (шаг по параметру задается с минусом). Затем, во время дополнительного исследования вычисляются скорости реакций. (Выражения для их вычисления записываются в программу FK0M при создании или коррекции изучаемой модели.)
На рис. 3.1, 3.2 представлены полученные результаты численных расчетов. Рис. 3.1 представляет график зависимости суммарной скорости реакции Rsum (кр. 1) от температуры и разложение зависимости суммарной скорости на составляющие скорости на ИП иридия (кр.2 - Лип) и ПО (кр.З - Япо)- Как следует из рис. 4.2 и рис. 3.1, математическая модель (3.2.1) хорошо описывает практически все особенности эксперимента, а именно: 1) величину максимальной скорости R$Um(T) , 2) резкий рост скорости после снятия торможения (десорбции СО); 3) небольшое падение скорости после выхода ее на максимум; 4) гистерезис скорости реакции при попеременном увеличении и уменьшении температуры; 5) s-образную форму «активной» (верхней) ветви скорости реакции.
Работа с моделью в диалоговом режиме
Для работы с моделью пользователь переходит в соответствующий раздел: «Задача Коши» или «Нелинейные системы», где выбирает модель для исследования из банка моделей. После этого происходит автоматический переход в режим трансляции и генерирования exe-файла, т. е. исполняемой программы, дающей пользователю возможность проведения численного исследования интересующей его модели. При трансляции программ информация, которая может помочь проверить правильность сформированной модели, записывается в файл, который доступен в окне просмотра. Далее представлены блок-схемы организации численного исследования моделей в рамках пакета. Схема численного исследования
Предположим, что в разделе «Задача Коши» мы вышли на стационарное решение при некотором значении параметра ао- Это решение можно взять в качестве стартового в методе продолжения. Оно записывается в файл исходных данных для раздела «Нелинейные системы»).
В этом разделе методом продолжения по параметру строится диаграмма стационарных решений системы т.е. графики компонент вектор-функции х(а), являющихся решением рассматриваемой системы.
Точки ветвления решения типа поворот фиксируются при этом следующим образом. На каждом шаге продолжения, начиная со 2-го, решение на элементарном отрезке по текущему параметру р представляется в виде эрмитовой кубической параболы. Это позволяет, продифференцировав кубическую параболу по ц, на каждом интервале проверять условие da/dp. = 0 при построении диаграммы стационарных решений системы (4.2.1).
После завершения процесса продолжения формируется таблица, отражающая изменение характера устойчивости матрицы J = }х{х{а) а) на решениях системы при изменении параметра а. Иными словами, для каждого значения а вычисляется характеристика устойчивости k(J). В таблицу выдается сообщение STABIL, если k(J) К, (где & - достаточно большое заданное число, определяемое машинной точностью вычислений); UNSTABIL, в противном случае.
Кроме того, таблица содержит значения определителя матрицы J. Эта информация позволяет высказать ряд суждений относительно полученных стационарных решений.
Рассмотрим ситуацию, когда в процессе продолжения решения по параметру обнаружилось, что в соседних точках а cv+ и a = at определитель матрицы J меняет знак. Следовательно, в некоторой точке а = а, а а. ort+, имеет место ветвление решений типа поворот. В результате характер устойчивости стационарных решений либо меняется: UNSTABIL при а = а & STABIL при а = а+ . Это означает, что некоторые из собственных чисел матрицы J, гурвицевой при а = а или а а , выходят на мнимую ось комплексной плоскости при некотором а = а, ос Є [а , а ]. Как отмечалось, при этом возможно зарождение периодических колебаний, описываемых автономной системой уравнений (автоколебания), в которые переходит сколь угодно мало возмущенное неустойчивое стационарное решение в достаточно малой окрестности точки а. = а (бифуркация Андронова - Хопфа). Значение а на интервале [а , а ] уточняется с помощью метода деления отрезка пополам.
Во многих случаях это явление можно обнаружить непосредственно, если решать задачу Коши с начальными данными в виде «возмущенного» неустойчивого стационара. Таким же способом, через использование диаграммы стационарных решений для задания начальных условий задачи Коши, может быть установлена граница интервала по а, за которой колебания отсутствуют. При этом проверяется, что стационар, полученный из решения задачи Коши, принадлежит построенной ранее диаграмме.
Раздел включает в себя методы численного интегрирования систем обыкновенных дифференциальных уравнений в форме Коши: полунеявный метод Розенброка и метод Гира, численные схемы которых приведены во второй главе, 3.
С помощью меню можно выбрать для решения задачи Коши один из методов и затем задать вектор начальных данных х, значения параметров модели и алгоритмические константы соответствующего метода. Вся входная информация сохраняется в файле исходных данных, который заполняется с экрана. Эти файлы составляют «банк исходных данных», с помощью которого происходит накопление информации для каждой новой модели.
