CC BY
12
УДК 530.182
ТЕОРЕТИЧЕСКОЕ ИССЛЕДОВАНИЕ ОБРАЗОВАНИЯ НЕОДНОРОДНЫХ СТРУКТУР НА ФРОНТЕ ГОРЕНИЯ НА ОСНОВЕ РЕДУЦИРОВАННОЙ МОДЕЛИ
Э.О. Якупов, В. В. Губернов, А. А. Полежаев
В химических системах, в частности, в системах горения, наблюдается множество пространственно-временных режимов. Так, в зависимости от числа Льюиса, определяющего отношение эффективности переноса тепла к диффузии химических компонентов, на фронте горения могут возникать либо автоволновые, либо ячеистые структуры. На основе простой математической модели, полученной редукцией полной модели, описывающей кинетику горения водорода, были исследованы особенностии формирования такого типа структур. В данной работе нас интересовала возможность описания в рамках модели ячеистых (тьюринговских) структур на фронте горения. Был проведён параметрический анализ, определяющий критерии возникновения тьюринговской неустойчивости. Оказалось, что необходимым условием является число Льюиса меньше единицы, что согласуется с экспериментальными наблюдениями. В найденной области параметров были проведены численные расчеты, а в качестве результатов были продемонстрированы различные варианты сформированных структур.
Ключевые слова: неустойчивости на фронте реакции, структуры Тьюринга, нелинейные системы, математическое моделирование.
Введение. С точки зрения нелинейной динамики химические системы являются удобными объектами для изучения разнообразных пространственно-временных структур, наблюдаемых в системах различной природы как неживых, так и живых. В химических системах удалось наблюдать множество пространственно-временных режимов ФИАН, 119991 Россия, Москва, Ленинский пр-т, 53; e-mail: edward.thrst@gmail.com.
[1-5], в частности, распространяющиеся волны, пульсирующие, хаотичные, кольцевые, спиральные, квазипериодические, стационарные неоднородные структуры. В некоторых случаях химические реакции протекают не равномерно в пространстве, а в виде распространяющегося фронта. Среди таких - реакция между изначально разделенными реагентами, которые в результате реакции образуют неравномерно распределенный осадок (так называемые структуры Лизеганга) [6], и самораспространяющийся высокотемпературный синтез (СВС) [7]. Другим ярким примером является распространение фронта пламени. В системах горения возникновение нелинейных волновых структур может быть обусловлено диффузионно-тепловыми [8] или гидродинамическими неустойчивостями [4]. Как известно, возникновение диффузионно-тепловых неустой-чивостей связано с нарушением равновесия между потоками тепла и реагентов, которое выражается числом Льюиса [9]. В случае, когда число Льюиса меньше единицы, возникает ячеистая структура или, в терминах нелинейной динамики, диссипативная (тьюринговская) структура, что проявляется в формировании неоднородного распределения температуры в плоскости фронта пламени. Для числа Льюиса больше единицы на фронте пламени возникают неустойчивости распространения и пульсации, которые приводят к формированию различных нестационарных режимов горения.
Известно, что пульсирующие, концентрические, спиральные, ячеистые нелинейно-волновые структуры могут возникать при определенных условиях в пламенах предварительно перемешанных смесей углеводородных топлив, стабилизированных на пористой горелке [10-15]. Динамика этих нелинейных структур также может усложняться и приводить к формированию квазипериодических и хаотических режимов горения. Количественно и качественно с помощью асимптотического [12, 13] и численного анализа [14, 15], соответственно, было показано, что формирование таких режимов обусловлено возникновением диффузионно-тепловой неустойчивости.
Ранее мы продемонстрировали, что как автоволновые, так и ячеистые структуры на распространяющемся фронте горения могут быть описаны в рамках блочной модели, в которой один блок описывает волну горения, а другой - возникающие на ней структуры [16-18]. При этом мы для каждого из блоков использовали простые феноменологические модели, никак не отражающие конкретные детали процесса. Кроме того, используя известную кинетическую схему элементарных реакций окисления водорода, мы показали, что в результате редукции соответствующей детальной математической модели вначале до четырёх, а затем и до двух уравнений, соответствующих критическим про-
цессам на фронте волны, можно получить модель, способную объяснить возникновение на фронте спиральных волн, наблюдаемых в экспериментах [19].
В настоящей работе мы продолжаем исследования полученной ранее редуцированной модели, а именно демонстрируем, что при изменении экспериментальных условий, приводящих к изменению параметров модели, она способна объяснить возникновение на фронте горения не только волновых, но и ячеистых (тьюринговских) структур.
Модель. Построение модели в данном исследовании аналогично тому, что было проделано в работе [19]. Повторим здесь основные моменты. Редукция детальной кинетической схемы горения водорода в случае богатых смесей приводит к модели с четырьмя уравнениями - для безразмерной концентрации радикалов Н, Н02, кислорода 02 и температуры в. Численное исследование этой модели продемонстрировало, что фронт горения не однороден по распределению радикалов, а формирует два слоя - условно высокотемпературный, в котором доминируют радикалы Н, и низкотемпературный, где в основном присутствуют радикалы Н02. Предполагая, что структуры на фронте волны горения формируются именно в низкотемпературном слое, мы провели дальнейшую редукцию модели и получили два уравнения - для радикалов Н02 и температуры в:
ди
— = р — т ехр(—Е/в)и + БиАи, (1а)
дв
— = 7 (т ехр(—Е/в)и — в(в — во)) + £еДв. (1Ь)
дЬ
Следует отметить, что в качестве характерных масштабов длины и времени брались тепловая толщина фронта горения и тепловая толщина, деленная на ламинарную скорость горения, соответственно. Температура масштабировалась на адиабатическую температуру горения, а концентрация реагентов - на характерное значение во фронте горения. Подобный подход является общепринятым [9]. Более детально этот вопрос будет раскрыт в нашей будущей работе, в которой будет проведена строгая редукция полного кинетического механизма реакции к задаче вида (1).
Качественно смысл членов этих уравнений следующий. В первом уравнении для радикала Н02, концентрация которого обозначена переменной и, первый член соответствует их производству с постоянной скоростью р за счёт подпитки радикалами Н из высокотемпературного слоя. Данный путь реакции главным образом определяется элементарной стадией Н + 02 + М ^ Н02 + М, которая является доминирующим каналом разветвления цепи при высоких давлениях и обладает нулевой энергией активации. Второй член в уравнении (1а) соответствует рекомбинации радикала Н02 и определяется в основном реакциями Н02 и Н, имеющими относительно невысокую энергию активации
[20]. Скорость реакции рекомбинации подчиняется закону Аррениуса mexp(-Е/в) с энергией активации Е. Эта реакция экзотермическая, поэтому аналогичный член входит и во второе уравнение для температуры в. Поскольку мы будем исследовать эти уравнения в двумерном случае (в плоскости низкотемпературного слоя), то в уравнение (1b) для температуры мы добавили линейный член з(в — во), который описывает релаксацию температуры к равновесному для рассматриваемого слоя значению и по сути является следствием диффузии тепла поперёк слоя, которой мы пренебрегли. Множитель y соответствует отношению характерных времён, возникающих в результате редукции исходной кинетической схемы и обезразмеривания. Отметим, что полученная нами модель, по сути, аналогична феноменологической модели Сальникова [21], которая применялась для объяснения автоволн в пламёнах [22].
Как отмечалось выше, для того чтобы в системе горения возникали тьюринговские (ячеистые) структуры при наличии диффузии, число Льюиса должно быть меньше единицы. Исходя из этого, потребуем, чтобы меньше единицы было отношение коэффициентов диффузии D = D@/Du < 1. Как известно [23], для развития тьюрингов-ской неустойчивости нужно, чтобы переменная, играющая роль активатора, диффундировала медленнее, чем другая - ингибитор. Таким образом, роль активатора играет температура в, а ингибитора - радикалы HO2 (переменная u). Для того чтобы модель типа "реакция-диффузия" демонстрировала нетривиальное поведение, в частности, возникновение диссипативных структур, необходимо, чтобы нуль-изоклина активаторной переменной - в нашем случае в - была немонотонной. Это имеет место, когда E > 4в0. На рис. 1 изображены нуль-изоклины первого уравнения ui(e) = p ехр(Е/в) и второго уравнения и2(в) = p(в — в0)ехр(Е/в). Их точка пересечения, стационар системы, имеет координаты (в st,u ¡i) = (во + p, p ехр(Е/(во + f))).
Проведём линейный анализ уравнений (1) и найдем область параметров, отвечающую достаточным и необходимым условиям для диффузионной неустойчивости и возникновения тьюринговских структур. Для этого, линеаризовав систему этих уравнений вблизи стационара (вst,ust), найдём матрицу A:
A = | fu fe\ = f—mexp(—Е/в^) —pE/в2st
ge) u e l YmexP(—E/вst) y(РЕ/в^ — s)
Соответствующее дисперсионное уравнение имеет вид:
IA — A +I1 0 ]k2 0 D
Рис. 1: Фазовый портрет системы уравнений (1). Нуль-изоклина первого уравнения: п\(0) = р ехр(Е/0), второго уравнения - и2(0) = р (© — @0)ехр(Е/О). Стационарная точка: = (О0 + р, р ехр(Е/(О0 + р))). Немонотонность одной из изоклин и
расположение стационара на средней неустойчивой ветви являются необходимыми, но недостаточными условиями для образования структур в нашем случае.
или
Л2 + Л[к2(1 + Б) — (/ + д©)] + Бк4 — (Б ¡и + д@)к2 + /ид© — /©ди = 0, (2)
где Л (к2) - собственные числа задачи, и для возникновения диссипативной неустойчивости необходимо, чтобы было И,е(Л(к2)) > 0 в некотором ненулевом диапазоне волновых чисел к. Соответствующие условия образования диссипативных структур в двухкомпо-нентной системе хорошо известны (см., напр., [23]). Применительно к нашей задаче они приобретают вид:
/«д© — /©д« > 0, (За)
/и + д© < 0, (ЗЬ)
Б/и + д© > 0, (Зс)
(Б/и + д©)2 — 4Б(/ид© — /©ди) > 0. (Зё)
Условие (3а) выполняется в нашем случае всегда, поскольку единственная стационарная точка модели не может быть седлом.
Модель, представленная уравнениями (1), допускает большой набор варьируемых параметров. Зафиксируем часть из них: Е = 1.0, во = 0.2, 7 = 0.25, т = 1.2, О = О©/Ои = 0.02. Таким образом, в нашем случае мы оставляем свободными только два параметра системы: р - скорость рекомбинации радикала из первого уравнения (1а), и в - поперечный коэффициент диффузии из второго уравнения (1Ь). Соотношения (3) определяют следующие области на плоскости параметров (р, в), представленные на рис. 2, а именно:
2.00 1.75 1.50 1.25 1.00 0.75 0.50 0.25 0.00
0.0 0.2 0.4 0.6 0.8 1.0
Р
Рис. 2: Области параметров, отвечающие набору условий (3). Внутри синей области стационарная точка находится на средней ветви изоклины. Внутри зеленой области выполняются условия тьюринговской неустойчивости. Внутри оранжевой области состояние системы отвечает бифуркации Пуанкаре-Андронова-Хопфа. В области между зеленой и оранжевой границами выбран ряд точек, и для соответствующих наборов параметров был проведен численный счет модели (1). Результаты представлены на рис. 4.
• Из критериев (3Ь) и (3с) следует, что члены на главной диагонали матрицы А (/и и д©) должны иметь противоположные знаки, что в свою очередь определяет положение стационарной точки на средней, неустойчивой ветви изоклины и2(в) на рис. 1. Это соответствует области внутри синей границы на рис. 2.
• Соотношение (3ё) гарантирует наличие неустойчивых пространственных мод. При этом свободный член квадратного уравнения (2) отрицателен в некотором диапазоне волновых чисел, что гарантирует существование в этом диапазоне положительного собственного числа Л (к2). Соответствующая область находится внутри зеленой границы на рис. 2. Заметим, что между зеленой и синей границами, несмотря на то, что стационар находится на неустойчивой ветви, система ещё сохраняет свою устойчивость.
• Требование (3Ь) на след матрицы А гарантирует наряду с условием (3а) (выполняющимся тождественно), что в отсутствие диффузии стационарная точка устойчива. Соответствующая область лежит вне оранжевой границы на рис. 2. При нарушении этого условия будет иметь место бифуркация Пуанкаре-Андронова-Хопфа, приводящая к автоколебаниям.
Таким образом, объединив все вышеперечисленные условия, мы определяем область параметров (р, в) внутри зеленой границы, за исключением оранжевой области. Множество значений р и в, полученные так вместе с прочими зафиксированными параметрами системы (1), полностью отвечает требованию тьюринговской неустойчивости, что должно обеспечить образование диссипативных структур. Чтобы исследовать особенности возникающих структур, были выбраны несколько точек, отмеченные на рис. 2, для которых было проведено численное моделирование.
Результаты численного моделирования. Численный счет проведен с помощью метода разделения на физические процессы, при котором диффузионная часть вычислялась методом Писмена-Рачфорда (Реасешап-И,асЬГогё), а реакционная - методом Кранка-Николсона.
Уравнения (1) рассматривались на двумерной квадратной области с размерами сторон 100x100. В качестве граничных использовались условия нулевых потоков. Начальное состояние системы задавалось вблизи стационара случайным шумом в промежутке [0, 0.25], для которого использовался генератор псевдослучайных чисел с посевом, основанным на текущем времени.
Параметры счета: шаг по пространству = = 0.1 и шаг по времени ^ = 0.01.
Точки на параметрической плоскости (р, в), для которых были проведены расчеты, выбраны таким образом, чтобы равномерно заполнить широкую область параметров, что позволило продемонстрировать различные виды получаемых структур. На рис. 3(а) показаны изоклины для каждого случая, обозначенные одним цветом, где видно, как по
Рис. 3: Изоклины (й) и дисперсионные кривые Яе(А) (Ъ) для выбранных точек.
мере роста параметров р и в стационар (точка пересечения) смещается вдоль средней ветви одной из изоклин. Справа показаны дисперсионные кривые И,е(А(к2)) (рис. 3(Ь)). Видно, что для каждого набора параметров р и в существует область значений волновых чисел к, при которых собственное значение задачи на линейную устойчивость является действительным положительным числом и стационарные состояния теряют устойчивость по отношению к тьюринговского типа неустойчивости.
Результаты численного счета модели (1) представлены на рис. 4 для тех же наборов параметров.
Можно заметить, что при движении в параметрической области (в,р) слева направо, что соответствует перемещению стационарной точки вдоль неустойчивого участка изоклины д0/д£ = 0 от её максимума до минимума (см. рис. 1 и 3(а)), меняется вид получаемых структур. Сначала структуры организуются в виде "горячих" ячеек на "холодном" фоне (рис. 4(а), 4(Ь)). Это наблюдается, когда точка пересечения изоклин и^в) и и2(0) находится в верхней части неустойчивой ветви. В то время, как точка пересечения нуль-изоклин смещается в нижнюю часть ветви, а значит с ростом температуры 0, характер структур, по сути, становится инвертированным - а именно, имеет вид "холодных" ячеек на "горячем" фоне (рис. 4(£)). В области параметров, при которых стационар располагается в средней части неустойчивой ветви, наблюдаются лабиринтные или полосные структуры (рис. 4(ё)), которые, по сути, являются переход-
О О О О 0-0о 5 о о О о °о ос
> о оо0°ооооо о О о00°0000 о ° о о ° )00о0°о0о оо0о
> О о о О ° ° ° О О С
°о° о
> п о
о ° о О О
О ° о о О о
ООО
о о с
. о О
О и „ _ О и
"оооооооо о 0 0 000°00 00 С °0 ° 0000000 о°с
0оо 0 0 0 о00 о °о° с
Й|1Ш|И!Й|1|5
00°00Л°0 о°о°0Оооо 0 00о00000о0°0о°000 О О п о° о О О °о О °о° о о о
50°0°00°0 О О ОО О о о 5 о О О г->° 00 ° 2 ° о О
> 2 0 о ° о о о 0°о° о о° ° о С
° о° ° ООО ° °
^ Л ^ л г\ ^ ^ ^ ^ ^ г\ г
ш
зхо^и I
Инн
°о о°о° оЪ о°о°Ь О о°о°пс 00л0п0000л° ° о П° 0 ° О о0 С 0°ОоО°оОо0°Оп°Оо°35о
)00оо°о00^°0л00л0п00 зооооОо0п°о0Оп0о^о2оо ООооО°о00°Оо0Ол0Оп00 )00000°00000°00°00°0( 0°000°00°0000000000< 3 О О О Оо Ог1 Ол^П ^о °
0000Оп0Оо0Ол0о0°о9Г)0л0с о О О О о0 о о°ОпО 00 о о°оо О 3 о О О о о О о° о° О О О о°о0°
3 О О о о о о о °о °п 0 °п О °о о ° о О О О О О О ОоО °п о °п О °о О 0
ГоГо®о0°о° оОоОс
Рис. 4: Структуры Тьюринга, полученные в результате численного счета модели (1) с параметрами Е = 1.0, 7 = 0.25, т = 1.2, в0 = 0.2, О = О©/Ои = 0.02, проведённого методом деления по физическим процессам в области 100x100 с шагом по пространству вх = 0.1 и шагом по времени вЬ = 0.01. Нумерация рисунков соответствует, точкам, отмеченным на рис. 2 и 3.
ными от одного предельного случая к другому. В промежуточных случаях получаются такие состояния системы, когда пятна (ячейки) и полосы сосуществуют (рис. 4(с), 4(е)), что приводит к долгому, по сравнению с остальными результатами, дрейфу структур в плоскости и их превращению из одного вида в другой.
Структуры, представленные на рис. 4(а), 4(Ь), 4(ё), 4(£), ячеистые и лабиринтные, приходят к устойчивому стационарному состоянию за время, не превышающее 2000 без-
размерных единиц. В промежуточных же режимах, представленных на рис. 4(c), 4(e), эволюция структур - дробление полос на ячейки и обратный процесс - продолжается даже после 4000 единиц безразмерного времени, что может говорить о том, что в этих случаях в принципе не возникает стационарного состояния.
Таким образом, в результате численного моделирования было показано, что в полученной области параметров можно наблюдать возникновение тьюринговских структур. Структуры имеют вид как пятен, так и полос, в зависимости от того, в каком месте изоклины находится стационар.
Заключение. В данной работе было продемонстрировано, что последовательной редукцией полной модели, основанной на детальной схеме окисления водорода, может быть получена простая модель, которая, как показали численные эксперименты, способна описать формирование двумерных ячеистых структур на фронте волны горения. Модель по своей структуре аналогична модели Сальникова [21], предложенной им в 40-х годах прошлого века для объяснения осцилляций в процессах горения, однако, в отличие от последней, она учитывает реальные химические процессы в пламени. Мы показали, что нелинейно-волновые структуры могут формироваться в лидирующей -низкотемпературной - зоне фронта горения, где доминирующим процессом являются формирование и рекомбинация гидропероксильных радикалов. В рамках полученной модели, описывающей динамику радикалов и поля температуры в низкотемпературном слое фронта, проведен параметрический анализ условий возникновения неустойчивости Тьюринга. Численное интегрирование модельных уравнений для параметров, соответствующих возникновению этого типа неустойчивости, показало, что возможно формирование структур различной морфологии: горячие ячейки, полосатые структуры и холодные ячейки. В дальнейшем мы планируем показать существование подобных решений для более полной модели, включающей как низко-, так и высокотемпературные зоны фронта горения.
Работа выполнена при финансовой поддержке РНФ (грант № 22-21-00546).
ЛИТЕРАТУРА
[1] D. Walgraef, Spatio-temporal pattern formation: with examples from physics, chemistry, and materials science (Springer Science & Business Media, 2012).
[2] M. Cross, P. Hohenberg, Reviews of Modern Physics 65, 851 (1993). DOI: 10.1103/RevModPhys.65.851.
[3] S. Camazine, J. Deneubourg, N. Franks, et al., Self-organization in biological systems (Princeton University Press, 2003).
[4] G. Sivashinsky, Annual Review of Fluid Mechanics 15, 179 (1983). DOI: 10.1146/annurev.fl.15.010183.001143.
[5] Y. Ju, K. Maruta, Progress in Energy and Combustion Science 37, 669 (2011). DOI: 10.1016/j.pecs.2011.03.001.
[6] R. Liesegang, Naturwissenschaftliche Wochenschrift 11, 353 (1896).
[7] A. Merzhanov, I. Borovinskaya, Doklady Akademii Nauk SSSR, Seriya Khimiya 204, 366 (1972).
[8] B. Matkowsky, D. Olagunju, SIAM J. Appl. Math. 39, 290 (1980). DOI: 10.1137/0139025.
[9] Y. Zeldovich, G. Barenblatt, V. Librovich, G. Makhviladze, The mathematical theory of combustion and explosions (Consultants Bureau, 1985).
10] M. Gorman, C. Hamill, M. El-Hamdi, K. A. Robbins, Combustion Science and Technology 98, 25 (1994). DOI: 10.1080/00102209408935394.
11] M. Gorman, M. El-Hamdi, K. A. Robbins, Combustion Science and Technology 98, 47 (1994). DOI: 10.1080/00102209408935396.
12] J. Buckmaster, SIAM Journal on Applied Mathematics 43, 1335 (1983). DOI: 10.1137/0143089.
13] V. Kurdyumov, M. Sanchez—Sanz, Proceedings of the Combustion Institute 34, 989 (2013). DOI: 10.1016/j.proci.2012.06.039.
14] V. Gubernov, V. Bykov, U. Maas, Combustion and Flame 185, 44 (2017). DOI: 10.1016/j.combustflame.2017.07.001.
15] B. Stefen, Combustion Science and Technology 22, 143 (1980). DOI: 10.1080/00102208008952379.
16] E. Yakupov, A. Polezhaev, Bulletin of the Lebedev Physics Institute 45, 165 (2018). DOI: 10.3103/S1068335618060015.
17] E. Yakupov, A. Polezhaev, V. Gubernov, T. Miroshnichenko, Physical Review E 99, 042215 (2019). DOI: 10.1103/PhysRevE.99.042215.
18] E. Yakupov, V. Gubernov, A. Polezhaev, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 083147 (2020). DOI: 10.1063/5.0012435.
19] E. Yakupov, V. Gubernov, A. Polezhaev, Izvestiya VUZ. Applied Nonlinear Dynamics 29, 538 (2021). DOI: 10.18500/0869-6632-2021-29-4-538-548.
[20] G. Stahl, J. Warnatz, Combustion and Flame 85, 285 (1991). DOI: 10.1016/0010-2180(91)90134-W.
[21] И. Сальников, Журнал Физической Химии 23(3), 258 (1949).
[22] S. Scott, J. Wang, K. Showalter, J. Chem. Soc., Faraday Trans. 93, 1733 (1997). DOI: 10.1039/A608474E.
[23] J. II. Murray, Spatial models and biomedical applications (Springer, 2003).
Поступила в редакцию 25 июня 2022 г. После доработки 29 августа 2022 г. Принята к публикации 30 августа 2022 г.
