Введение
Актуальность работы. Работа посвящена математическому моделированию критических явлений в химических системах на основе новых методов геометрической теории сингулярных возмущений.
Сингулярно возмущенные системы обыкновенных дифференциальных уравнений используются для моделирования процессов различной природы. Так, в моделях химической кинетики наличие малого параметра связано с тем, что в химической системе одновременно происходят резко отличающиеся по скорости процессы.
Основы теории сингулярных возмущений были заложены в работах Тихонова А. Н. Наиболее широкое распространение получил метод пограничных функций Васильевой-Тихонова. Дальнейшее развитие теория получила в работах Аносова Д. В., Андронова А. А., Боголюбова Н. Н., Бутузова В. Ф., Васильевой А. Б., Вишика М. И., Дмитриева М. Г., Дородницына В. А., Емельянова С. В., Ильина А. М., Кащенко С. А., Кобрина А. И., Коровина С. К., Красовского Н. Н., Крей-на С. Г., Крылова Н. М., Ломова С. А., Люстерника Л. А., Мартынен-ко Ю. Г., Маслова В. П., Митропольского Ю. А., Мищенко Е. Ф., Моисеева Н. Н., Найфэ А. X., Нейштадта А. И., Новожилова И. В., Понт-рягина Л. С, Розова Н. X., Хапаева М. М., Чернышова К. И., Bonet С, Chang Н. С, Cole J., Howes F. A., Kelley A., Miller J., O'Malley R. E., Schneider K. R., Smith D. R., Van Dyke M. и многих других авторов (см. [2, 3, 11, 12, 15, 16, 17, 18, 19, 20, 22, 32, 33, 43, 48, 49, 54, 56, 58, 64, 66, 67, 68, 69, 70, 71, 74, 92, 96, 124, 154, 169, 197]).
Следует отметить, что методы теории сингулярных возмущений применялись для исследования моделей химической кинетики в работах Бутузова В. Ф., Васильевой А. Б., Вольперта А. И., Калачева Л. В., Худя-ева С И. и др. [13, 18, 25].
Обычное предположение теории сингулярных возмущений состоит в
том, что основной функциональный определитель быстрой подсистемы отличен от нуля. Однако во многих прикладных задачах, в частности в моделях химических систем, это условие нарушается, и возникают различные критические ситуации. Различные критические случаи рассматривались в работах Бутузова В. Ф., Васильевой А. В., Волосова В. М., Нефедова Н. Н., Соболева В. A., Gu Z., O'Malley R. E., Schneider К. R., Williams F. [18, 21, 23, 24, 52, 82, 146].
Нарушение этого условия может привести к возникновению траекторий-уток. Термин "утка" возник в математической литературе в связи с применением нестандартного анализа к исследованию дифференциальных уравнений. Первое упоминание об утках принадлежит, по-видимому, J. L. Callot, M. Diener, F. Diener (1978) [129]. Исследование траекторий—уток для различных классов систем проводилось в работах многих авторов. Следует отметить работы Арнольда В. И., Горелова Г. Н., Звонкина А. К., Ильяшенко Ю. С, Колесова А. Ю., Колесова Ю. С, Мищенко Е. Ф., Покровского А. Н., Розова Н. X., Сам-борского С. Н., Соболева В. А., Шубина М. A., Benoit E., Eckhaus W. [5, 10, 30, 37, 51, 65, 73, 77, 78, 120, 121, 124, 130, 131, 133, 142, 143]. Заметим, что если первоначально термин "утка" употреблялся применительно к предельным циклам уравнений типа уравнения Ван-дер-Поля (так называемые циклы-утки), то позднее речь идет об объектах более общей природы — о траекториях-утках как одномерных устойчиво-неустойчивых медленных интегральных многообразиях.
Приложение траекторий-уток к моделированию критических явлений в химических системах позволило решить ряд интересных задач. Заметим, что исследованиям критических явлений в химических системах посвящено множество публикаций, среди которых следует отметить работы Абрамова В. Г., Бабкина В. С, Бабушка В. И., Баренблатта Г. И., Барзыкина В. В., Быкова В. И., Гольдштейна В. М., Горелова Г. Н., Дубо-вицкого Ф. И., Зельдовича Я. Б., Кащенко С. А., Мержанова А. Г., Семенова Н. И., Соболева В. А., Тодеса О. М., Франк Каменецкого Д. А., Ху-дяева С И., Gray В. F., Griffiths J. F., Kassoy D., Linan A., Mclntosh A. C, Sivashinsky G. I. [7, 14, 28, 29, 38, 45, 59, 60, 61, 62, 63, 80, 81, 84, 85, 86, 89, 90, 95, 116, 117, 128, 140, 141, 145, 152, 153, 158, 162, 163, 170, 195].
Следует отметить, что при моделировании процессов горения, сопровождающихся резким ростом температуры, при помощи дифференциальных уравнений с частными производными в работах Галактионо-ва В. А., Курдюмова С. П., Малинецкого Г. Г., Михайлова А. П., Пове-
щенко Ю. А., Попова Ю. П., Самарского А. А. использовались режимы с обострением [6, 27, 39, 55, 75, 76]. Явление потери устойчивости в сложных системах химической кинетики рассматривалось в работах Ивановой А. Н., Тарнопольского Б. Л. и других авторов [40, 41].
В последние годы появилось значительное число публикаций, посвященных применению траекторий-уток в различных задачах биологии, механики, химии, экономики и электроники. Не претендуя на полноту, среди них можно выделить работы Колесова А. Ю., Колесова Ю. С, Мищенко Е. Ф., Покровского А. Н., Розова Н. X., Соболева В. А., Ваег S. M., Bar-Eli К., Braaksma В., Br0ns M., Dumortier F., Erneux Т., Freire E., Gamero E., Gaspar V., Guckenheimer J., Hoffman К., Krupa M., Mazzotti M., Milik A., Moehlis J., Morbidelli M., Peng В., Rodriguez-Luis A. J., Roussarie R., Serravalle G., Showalter K., Szmolyan P., Weckesser W. [51, 65, 83, 87, 118, 119, 125, 126, 132, 134, 137, 139, 147, 148, 157, 161, 164, 165, 170, 185, 186].
Одним из основных методов исследования сингулярно возмущенных систем является метод интегральных многообразий Боголюбова-Митропольского. Под интегральным многообразием здесь понимается гладкая инвариантная поверхность дифференциальной системы. Интегральное многообразие называется медленным, если движение по нему осуществляется со скоростями порядка единицы (в полной системе есть движения со скоростями порядка отрицательной степени малого параметра). Использование медленных интегральных многообразий позволяет понижать размерность изучаемых моделей и избавляться от вычислительной жесткости. Теория интегральных многообразий применялась для исследования сингулярно возмущенных систем в работах Ба-риса Я. С, Задираки К. В., Лыковой О. В., Митропольского Ю. А., Самойленко А. М., Сидоренко В. В., Соболева В. А., Стрыгина В. В., Фодчука В. И., Fenichel N., Hale J., Henry D., Knobloch H., Kokotovic" P. и др. [8, 9, 36, 64, 79, 88, 93, 94, 135, 149, 150, 155, 156].
Данная работа посвящена исследованию интегральных многообразий со сменой устойчивости систем обыкновенных нелинейных сингулярно возмущенных дифференциальных уравнений и применению полученных результатов для моделирования критических явлений в химических системах.
Как уже отмечалось, для широкого круга химических процессов характерно резкое различие скоростей превращения веществ, участвующих в этих процессах и резкое различие скоростей тепловых и концентраци-
онных изменений.
В каталитических системах скорости реакций, происходящих на поверхности катализатора, на порядок или даже на несколько порядков выше, чем скорости реакций в газовой фазе. Как правило, реакции на поверхности катализатора идут при сравнительно небольших температурах и скорость выделения тепла существенно ниже, чем скорости изменения концентраций реагирующих веществ.
Для задач теории горения является характерной высокая скорость тепловыделения при сравнительно низкой скорости расходования горящего вещества. Это различие носит настолько резкий характер для газофазных систем, что явление самовоспламенения приобрело название теплового взрыва.
Более того, в химической кинетике устоявшимся принципом является различение участвующих в процессе веществ по скоростям их превращения. В газофазной кинетике по этому принципу различают активные центры и основные вещества.
В последнее время существенно вырос интерес к изучению влияния на основной процесс так называемых буферных явлений. Для каталитических систем это может быть учет старения катализатора, адсорбции нереакционного продукта реакций на поверхности катализатора. Естественно, скорости буферных явлений на порядок ниже скорости реакции.
Исследования химических систем, разделенных на медленную и быструю подсистемы в силу различия скоростей, производится, как правило, методом квазистационарных концентраций Боденштейна-Семенова. Идея метода проста: предполагается, что система "подстраивается" под медленную подсистему за счет быстро приходящей к равновесию (квазистационарному режиму) быстрой подсистемы. Это позволяет учитывать при анализе значительно меньшее число параметров, что приводит к существенным упрощениям при исследовании модели. Формализм метода квазистационарных концентраций сводится к замене части дифференциальных уравнений в модели процесса алгебраическими, что позволяет понизить размерность модели.
При исследовании моделей химических процессов особенно интересен качественный анализ. Если размерность дифференциальной системы, описывающей поведение химического процесса, больше двух, то ее исследование методами качественной теории дифференциальных уравнений, как правило, связано с существенными трудностями. Особенно
8
это проявляется в неизотермических моделях из-за сильной нелинейности системы. Возможности численного анализа таких моделей при наличии существенно разномасштабных переменных также ограничены. Причиной этого является вычислительная жесткость системы, то есть высокая чувствительность к малым погрешностям вычислений. Метод квазистационарных концентраций пригоден только для анализа грубых ситуаций, когда дифференциальная система имеет притягивающее медленное интегральное многообразие.
Наличие в системе быстрых и медленных переменных позволяет использовать асимптотические методы, которые, как правило, предназначены для решения краевых и начальных задач на конечных промежутках и не приспособлены для качественного исследования. В данной работе используется метод качественного исследования систем с быстрыми и медленными переменными, являющихся типичными для моделей химических систем, так, в частности, для рассматриваемых моделей лазеров, каталитических реакторов и теплового взрыва.
Основная задача математической теории теплового взрыва заключается в исследовании динамики процесса горения при заданных размерах реакционного сосуда, теплофизических и кинетических характеристиках реагирующего вещества, коэффициенте теплоотдачи. Для классической модели теплового взрыва эти характеристики отражает некоторый параметр, значение которого определяется начальным состоянием химической системы. В зависимости от значения этого параметра происходит либо переход реакции к медленному режиму, что ведет к затуханию реакции, либо реакция переходит в режим самоускорения, что приводит к взрыву. Численные расчеты для сосредоточенной двумерной модели показывают, что переход от медленного режима к взрывному происходит в чрезвычайно узком промежутке изменения параметра, характеризующего начальное состояние системы. При некотором значении этого параметра, которое называется критическим, реакция идет максимально долго, не срываясь в режим взрыва и не переходя в медленный режим выгорания. Соответствующий режим будем называть критическим.
Задачи определения критических значений параметров и моделирования критических режимов для многофазных и многотемпературных систем и являются главными в рамках исследования моделей. Формализм решения этих задач сводится к обоснованию существования и построения асимптотических разложений медленных интегральных поверхностей со сменой устойчивости.
Ряд важных прикладных задач биологии, биофизики, механики, лазерной оптики, экономики и теории управления также приводит к необходимости изучения медленных процессов со сменой устойчивости.
Хорошо известно, что сингулярно возмущенные системы обыкновенных дифференциальных уравнений используются для моделирования процессов различной природы. В работе рассматриваются системы вида
х = f(x,y,z,e),
У = g(x,y,z,e), (1)
ez = p(x,y,z,a,e),
где е — малый положительный параметр, а — скалярный параметр, х — скалярная переменная, у и z — векторные переменные размерности п и т+1, соответственно, /, д яр — достаточно гладкие функции. Уравнение для скалярной переменной х выделяется по техническим соображениям. Обычно в качестве переменной х используется одна из медленных переменных, строго монотонная по времени.
Под медленной поверхностью системы (1) понимается поверхность, описываемая уравнением
p(x,y,z,a,0)=0. (2)
Основное предположение обычно состоит в том, что
det — (х, у,ф{х, у, а),а,0) ф 0,
где z = ф(х,у,а) — изолированное решение уравнения (2). Нарушение этого условия может привести к возникновению так называемых траекторий-уток.
Лист медленной поверхности устойчив, если собственные числа матрицы
dp
— (х,у,ф(х,у,а),а,0) (3)
имеют отрицательные вещественные части. Если хотя бы у одного из собственных чисел этой матрицы вещественная часть становится положительной, то лист теряет устойчивость. Листы медленной поверхности
10
разделяются так называемыми поверхностями срыва, имеющими размерность вектора у, на которых
др
det dz^Xt У: ^У) а^ а> °)= 0-
В е-окрестности устойчивого и неустойчивого листов медленной поверхности лежат притягивающее (устойчивое) и отталкивающее (неустойчивое) медленные интегральные многообразия. Медленное интегральное многообразие представляет собой инвариантную поверхность, движение по которой осуществляется со скоростью порядка единицы.
Наличие дополнительного скалярного параметра а обеспечивает условия для того, чтобы устойчивое и неустойчивое интегральные многообразия можно было "склеить" в одной точке поверхности срыва. Именно через эту точку проходит траектория, которая является уткой. Под траекторией-уткой обычно понимается траектория сингулярно возмущенной системы, которая проходит вначале вдоль устойчивого интегрального многообразия, а затем вдоль неустойчивого, причем оба раза расстояния порядка единицы.
Следует отметить, что для случая dim у = О, dim z = 1 уже в первых работах, посвященных траекториям-уткам, было установлено существование единственной траектории-утки, отвечающей единственному значению параметра а = а* (точнее, значение параметра а* существует на интервале порядка О(е~ 1/?)).
В случае же dim у = 1 картина принципиально иная: в работе установлено существование однопараметрического семейства траекторий-уток, так как наличие параметра а позволяет выбрать точку склейки устойчивого и неустойчивого интегральных многообразий.
Из вышесказанного следует, что траектория утка представляет собой одномерное медленное интегральное многообразие, "склеенное" из неустойчивой и устойчивой частей.
Рассмотрим несколько простых примеров.
П р и м е р 1. (dimу = 0).
В качестве простейшей системы с траекторией-уткой можно предложить систему вида
х = 1,
ez — 2xz + a. Ясно, что при а — 0 система имеет траекторию-утку z = 0.
11
П р и м е р 2. (dim у = 0). Для системы на плоскости
х ~ f(x,z,e),
ez = p(x, z, а, е)
устойчивые и неустойчивые части медленной кривой разделены точками срыва, в которых выполняется равенство
°
Обычно исследование таких систем в окрестности нерегулярных точек проводится в предположении, что в этих точках выполняется неравенство
Z
\t \/
ч ч 1 \
\ \ /\
ч ч ч /
о/ | У ч %ч\ 1 х
\/ ч [
ч ч
л
Рис. 1: Траектория-утка (сплошная линия) и "ложная" утка (пунктирная линия) системы из Примера 2.
Однако существует класс задач, для которых это условие нарушается при некотором значении а. Например, в системе
12
dx _
е— = z -х +а, at
при а = ±е линии z = ±? проходят вдоль медленной кривой бесконечно долго (см. рис. 1). Заметим, что траекторией-уткой является только z = х. В этом примере [142] точка х = 0, z = 0 является точкой самопересечения медленной кривой при а = 0. Подобные ситуации рассматривались в работах [5, 30, 130, 142]. Аналогичные системы возникают в моделях теплового взрыва в случае автокаталитической реакции. При этом траектории-утки являются естественным математическим объектом, позволяющим моделировать критические явления и находить критические значения параметров в виде асимптотических разложений по целым степеням малого параметра е.
Рис. 2: Медленная кривая (пунктирная линия) и траектория-утка (сплошная линия) системы из Примера 3.
Как было отмечено ранее, работы, посвященные исследованию траекторий-уток обычно содержат утверждения типа "Жизнь уток коротка". Под этим подразумевается, что значения параметра а, при
13
которых существуют траектории-утки, отличаются друг от друга на величину порядка О(е~1^К?), где к — некоторое положительное число. При этом только в работе [5] приведены достаточные условия, при которых этот факт имеет место, в то же время, нетрудно предложить примеры
"вечно живых уток".
П р и м е р 3. Рассмотрим систему [52]
О О О
ez = х + z - а .
Окружность
представляет собой траекторию-утку, движение по которой осуществляется по часовой стрелке. Верхняя полуокружность неустойчива, нижняя — устойчива (см. рис. 2) . Эта траектория существует при всех а2 > е2 /А. Рассмотрим обобщение Примера 3 на случай dim у > 0. Система
х = z,
ez = bz2 + f(x, e)
имеет особые точки при z = 0, f(xa,e) = 0. В соответствии с матрицей Якоби
' 0
1\
о)
особая точка линеаризованной системы является седлом (центром), если f'(xs,e)>0(f'(xs,e)<0).
Для исходной системы тип особой точки сохраняется. В случае седла это хорошо известный факт. В случае центра мы можем применить теоремы Ляпунова об особой точке типа центр благодаря существованию аналитического первого интеграла исходной системы. Этот первый интеграл может быть получен из следующего уравнения
?-(z2)' = bz* + f(x,e). (4)
Решение (4) описывается выражением
z2 = Се2Ъ*1е + - Г e2b^-T^?f{r,e)dr. ? Jo
14
Таким образом, аналитический по ж и г первый интеграл находится в виде
Г
е Если функция f{x,e) является многочленом
N
t(r F) _ Y^n rn
71=0
то частное решение уравнения (4) удовлетворяет соотношению
N
z2 = У qnxn.
71=0
Подставляя это соотношение и выражение
N ЛГ-1
(г2)' = ф(х,е) = ^nqnxn'1 = ^{п + l)qn+lxn
П=1 71=0
в (4), получим
ЛГ-1 N N
2 Z-j ' n+1 zl^ n Z^
тг=О n=0 n=0
Отсюда имеем
bqN + Рлг = 0,
т{п + l)qn+i = bqn + pn, n = N- l,...,0; 2
или
&v = —Pn/Ь,
Qn = (~Pn + |(n + l)tfn+i) A w = N - 1,..., 0. Инвариантное многообразие уравнения (4), описываемое уравнением
z2 = ф(х,е), устойчиво при bz < 0 (неустойчиво при bz > 0).
15
П р и м е р 4. При N = 1 имеем
дг = -a/b, qo = - (a+ ^a) /b, и медленное интегральное многообразие описывается уравнением
bz2 + ах + а -\—-а = 0. 26
Соответствующая кривая является параболой с устойчивой (неустойчивой) ветвью при z > 0 (z < 0), см. рис. 3.
Рис. 3: Медленная кривая (пунктирная линия) и траектория-утка (сплошная линия) системы из Примера 4.
П р и м е р 5. При N = 2 имеем f() 2+p1x+p0:
д2 = -а/Ь, qi = -ea/b2, qQ = ~ ( а + ~а) J /Ъ.
16
Рис. 4.
Рис. 5.
Медленное интегральное многообразие описывается уравнением
еа
е2а
bzl + ах2 + —х + а + —- = О
или
bz2
17 |
