CC BY
16
УДК 53.01
ИССЛЕДОВАНИЕ МЕХАНИЗМА ВОЗНИКНОВЕНИЯ АВТОВОЛНОВЫХ СТРУКТУР НА ФРОНТЕ РЕАКЦИИ
Э. О. Якупов, А. А. Полежаев
Предложен качественный механизм образования волновых структур на фронте реакции. Предполагается, что структуры образуются в результате взаимодействия двух подсистем, одна из которых ответственна за формирование фронта, а вторая - за возникновение самих структур. Рассмотрены три модели, в каждой из которых численно продемонстрированы двумерные аналоги концентрических и спиральных волн. В качестве подсистем были выбраны модели Фицхью-Нагумо, Фишера-КПП и Орегонатор.
Ключевые слова: автоволны, неустойчивости на фронте реакции, математическое моделирование.
Введение. Есть много примеров, когда химическая реакция протекает не однородно в пространстве, а в виде распространяющегося фронта. К таковым можно отнести реакцию между исходно пространственно разделёнными реагентами, образующими в результате реакции нерастворимое вещество, выпадающее в осадок (так называемые структуры Лизеганга) [1], а также СВС - самораспространяющийся высокотемпературный синтез [2].
Другой яркий пример - это распространение фронта горения. Этим явлениям посвящено множество экспериментальных и теоретических исследований, в которых продемонстрировано, что на таких фронтах могут возникать различные неустойчивости, приводящие как к осцилляциям, так и к формированию различных пространственных структур [3, 4]. Относительно недавно экспериментально было продемонстрировано, что на фронте пламени могут возникать спиральные и концентрические автоволновые структуры [5]. Такие структуры наблюдались как на расширяющемся сферическом фронте [6, 7], так и в микрореакторах [8]. Механизм их формирования до сих пор не понят.
ФИАН, 119991 Россия, Москва, Ленинский пр-т, 53; e-mail: edward.thrst@gmail.com, apol@lpi.ru.
В настоящем исследовании мы предлагаем качественное объяснение механизма формирования этих структур на основе построения и численного исследования соответствующих математических моделей. Идея предлагаемого механизма возникновения структур на фронте волны реакции заключается в следующем: исходную систему мы рассматриваем как комбинацию двух взаимодействующих подсистем, первая из которых формирует фронт волны реакции, а вторая отвечает за возникновение структур на этом фронте. Каждая из подсистем описывается моделью типа "реакция-диффузия".
Вообще говоря, интересующие нас структуры существенно трёхмерны, однако предлагаемые механизмы можно продемонстрировать на относительно простых пространственно двухмерных моделях.
Математические модели. Предположим, что источник концентрических волн расположен в исходном трехмерном пространстве на некоторой линии, причём фронт волны перпендикулярен этой линии. В качестве двумерной области, в которой мы будем моделировать волны, аналогичные концентрическим, возьмём сечение, проходящее через эту линию, причём она расположена на краю области и может являться источником автоволн (рис. 1). В качестве граничных мы используем условия нулевых потоков переменных на границах области.
Двумерным аналогом спиральной волны на фронте реакции является спиральная винтовая линия на поверхности цилиндра, причём она формируется при движении фронта несущей волны вдоль оси цилиндра. Соответственно, в численных экспериментах в прямоугольной области мы выбираем в качестве граничных периодические условия на боковых сторонах области, по которым происходит склейка цилиндра (рис. 2). На двух других сторонах в качестве граничных условий мы используем нулевые потоки.
Начальные условия для обоих случаев, а также качественный механизм формирования структур будут описаны в разделе, касающемся конкретных моделей.
Реализация обоих видов структур будет продемонстрирована на примере трёх блочных моделей. В качестве элементарных блоков мы выбрали хорошо известные модели Фицхью-Нагумо [9], Фишера-Колмогорова-Петровского-Пискунова (Фишера-КПП) [10, 11] и Орегонатор [12].
Первый блок модели - автономная система уравнений, описывающая распространяющийся плоский фронт реакции. Второй блок параметрически зависит от переменных первого блока, при этом вне фронта реакции он находится в устойчивом невозбудимом состоянии, а на фронте переходит в возбудимое состояние (для спиральной волны) ли-
бо на основной части фронта также в возбудимое состояние, а в небольшой области у края - в колебательное состояние (для аналога концентрической волны).
Для получения искомых структур весьма важно правильно задать начальные условия. Если для второго случая достаточно задать произвольное отклонение от стационарного состояния, то в первом, чтобы возникла спиральная волна, надо задать её зародыш, в качестве которого можно взять одномерное решение в виде движущейся волны в соответствующей модели.
Для численных экспериментов использовалась собственная программа. Вычисления проводились на квадратной сетке 20 х 20 с пространственным шагом 0.1 и с шагом по времени 0.0025. Для проверки сходимости варьировали шаги по пространству и времени. Полученные решения сходятся, поскольку при уменьшении шагов по пространству и времени в два раза результат существенно не меняется. При численном счете применялся метод расщепления по физическим процессам. Кинетическая часть уравнения решалась методом Рунге-Кутты четвертого порядка и затем включалась в решение уравнения диффузии методом переменных направлений.
Фицхью-Нагумо и Фицхью-Нагумо. Рассмотрим первый случай, когда в качестве обоих блоков выступают пары уравнений Фицхью-Нагумо
^ = -[(2 - 0.1)(1 - д) - К] + Бд Ад,
д1 = 1 дь т
дК , ^ .,
— = д - К + £>ЛАК,
— = - [и(и - а)(1 - и) - ь] + БиАи, дЬ £
дь
дЬ
— = и - ь + Б„Аь.
Здесь д, К, и, V - переменные, описывающие систему, которые в случае горения или химической реакции имеют смысл концентраций реагентов и температуры. Связь между подсистемами задается в виде а = а - Ьд, где параметры а и Ь выбираются исходя из требований к параметру а и наличия источника. Неоднородность величины а по пространству определяет, в каком состоянии находятся элементы среды (в невозбудимом, возбудимом или колебательном). Кроме того, для возникновения автоволн принципиально важно, чтобы коэффициент диффузии первой переменной и(д), "активатора", был существенно больше коэффициента диффузии второй переменной ь(К) - "ингибитора": Ди ^ ^ (Дд ^ Оставшиеся параметры в уравнениях, т и £ (много меньшие 1), определяют скорости и размеры получаемых волн и подбираются исходя из размеров области, на которой ведутся вычисления.
Рис. 1: Результаты моделирования системы Фицхью-Нагумо и Фицхью-Нагумо. Здесь и далее показана переменная, и степень затемнения обозначает величину этой переменной. Стрелки с пунктиром указывают на траектории движения отдельных волн: (а) аналог концентрических волн с источником у левой границы, (б) аналог спиральных волн на поверхности цилиндра.
В численном эксперименте были выбраны следующие параметры модели: т = 0.005, Од = 0.5, О = 0, £ = 0.002, Ои = 1, О = 0. Для получения структур, имитирующих концентрические волны, как на рис. 1(а), у левой границы задавалась узкая полоса (уравнение связи а = 0.5 — 0.8$), служащая источником автоволн, которые затем двигались по возбудимой области, возникающей на фронте реакции (уравнение связи а = 0.5 — 0.5$). В качестве граничных использовались условия нулевых потоков, в качестве начального условия выбиралось произвольное отклонение от стационара.
Для получения структур, аналогичных спиральным волнам, как на рис. 1(б), в начальный момент необходимо было задать зародыш волны на всей области. Эта волна выживала только на фронте, задаваемом первым блоком, на котором возникала возбудимая полоса (уравнение связи а = 0.5 — 0.5$), и двигалась по ней при наличии периодических граничных условий.
Модели Фишера-КПП и Фицхью-Нагумо. В следующей рассматриваемой системе первая пара уравнений заменяется уравнением Фишера-КПП, решение которого представляет собой движущуюся волну переключения. Такой выбор замены исходит из идеи
отойти от четырехкомпонентной системы к более простой трехкомпонентной. К тому же, волна переключения соответствует более реалистичной модели, например, простейшей модели волны горения
дд дЬ
— = - [и(и - а)(1 - и) - ь] + БиАи, дЬ £
дь дь
-£ = кд(1 - д) + БдАд,
— = и - ь + Б„Аь.
Параметр к определяет скорость роста переменной д в уравнении Фишера-КПП. Связь между подсистемами вводится схожим с предыдущим, но нелинейным образом: а = а - Ьд(1 - д). На систему налагается то же требование, что и в предыдущем случае -Би » .
В численных экспериментах для данной модели использовались следующие параметры: к = 6, Бд = 0.2, £ = 0.002, Би = 1, Бь = 0. Качественно схема возникновения структур аналогична описанной в предыдущей модели. На рис. 2(а) колебательная область у левой границы задавалась уравнением связи а = 0.6 - 3д(1 - д) (вертикальный пунктир), фронт реакции определялся как а = 0.6 - 2д(1 - д). Начальное условие -небольшое отклонение от стационара, граничные условия - нулевой поток. На рис. 2(б) фронт реакции задавался уравнением связи а = 0.6 - 2д(1 - д).
Модели Фишера-КПП и Орегонатор. Уравнения Фицхью-Нагумо являются простейшей моделью, описывающей автоволновые структуры, однако они мало подходят для описания химической системы. Чтобы ещё более приблизить модель к реальной системе, заменим блок уравнений Фицхью-Нагумо на Орегонатор [12] - модель, предложенную для реакции Белоусова-Жаботинского. Фронт реакции по-прежнему определяется решением уравнения Фишера-КПП
дд
дЬ ди 1
дЬ = £
-£ = кд(1 - д) + Бд Ад,
(Л \ ^ t д- и
и(1 - и) + /ь-
д + и
+ БиАи,
дь
— = и - ь + БьАь. дЬ
Параметры £, /, д определяют поведение Орегонатора, и требования к ним подробно описаны в соответствующей литературе [13]. Уравнение связи между подсистемами мы определим как / = а - Ьд(1 - д).
Рис. 2: Результаты моделирования системы Фишера-КПП и Фицхью-Нагумо. Стрелки с пунктиром указывают на траектории движения отдельных волн: (а) аналог концентрических волн с источником у левой границы, (б) аналог спиральных волн. 20.01
Рис. 3: Результаты моделирования системы Фишера-КПП и Орегонатор. Стрелки с пунктиром указывают на траектории движения отдельных волн: (а) аналог концентрических волн с источником у левой границы, (б) аналог спиральных волн.
В численных экспериментах для данной системы использовались следующие параметры: к =1, Бд = 0.1, £ = 0.02, д = 0.02, Би =1, Бь = 0. Кроме некоторых особенностей при расчете, таких как контроль за скоростью и шириной волны или чувствительность к параметрам системы, результат моделирования аналогичен предыдущим случаям. Решение автономного уравнения Фишера-КПП формирует фронт, на котором Орегонатор переходит в возбудимое состояние. Возникающая в последнем волна движется вдоль этого фронта.
Так для получения аналога кольцевых волн, представленных на рис. 3(а), мы определили источник на левом краю области с помощью уравнения связи между подсистемами: / = 5 - 13д(1 - д), а возбудимую область - / = 5 - 10д (1 - д); как и в предыдущих случаях в качестве граничных использовали условия нулевого потока.
Для получения аналога спиральной структуры, представленной на рис. 3(б), фронт волны реакции мы определили как / = 4 - 7д(1 - д); задали начальную волну на всей области и ввели периодические граничные условия.
Обсуждение. Цель данного исследования - разработать механизм возникновения автоволновых структур на распространяющемся фронте реакции. Мы предполагаем, что такие структуры могут возникать как результат взаимодействия двух подсистем: одна соответствует распространяющемуся фронту, а другая ответственна собственно за автоволны на этом фронте. Соответствующая математическая модель состоит из двух блоков для каждой из подсистем. В настоящей работе реалистичность данного механизма продемонстрирована в численных экспериментах в двумерном пространстве для простых базовых моделей. Следующими этапами будут исследование реальных моделей химических реакций, демонстрирующих такой тип поведения, в частности, волн горения, а также переход к численному моделированию в трёхмерном пространстве.
Работа частично поддержана грантами РФФИ № 17-01-00070 и 16-03-00758.
ЛИТЕРАТУРА
[1] R. E. Liesegang, Naturwiss. Wochenschr. 11, 353 (1896).
[2] А. Г. Мержанов, в: Физическая химия. Современные проблемы. Ежегодник, под ред. Акад. Я. М. Колотыркина (М., Химия, 1983), с. 6.
[3] G. I. Sivashinsky, Annual Review of Fluid Mechanics 15(1), 179 (1983).
[4] G. Wang, Y. Li, W. Yuan, et al., Combustion and Flame 191, 368 (2018).
[5] S. K. Scott, J. Wang, and K. Showalter, J. Chem. Soc., Faraday Trans. 93, 1733 (1997).
[6] G. Jonaas, J. K. Bechtold, and C. K. Law, Proceedings of the Combustion Institute 31(1), 1039 (2007).
[7] G. Jomaas and C. K. Law, Physics of Fluids 22(12), 124102 (2010).
[8] S. Kumar, K. Maruta, S. Minaev, and R. Fursenko, Physics of Fluids 20(2), 024101 (2008).
[9] R. FitzHugh, Bull. Math. Biophysics 17, 257 (1955).
[10] R. A. Fisher, Ann. Eugenics 7, 353 (1937).
[11] А. Н. Колмогоров, И. Г. Петровский и Н. С. Пискунов, Бюл. МГУ. Математика и механика 1(6), 1 (1937).
[12] R. J. Field and Richard M. Noyes, J. Chem. Phys. 60, 1877 (1974).
[13] John J. Tyson and Paul C. Fife, J. Chem. Phys. 73, 2224 (1980).
Поступила в редакцию 6 марта 2018 г.
