Теоретические основы химической технологии, 2023, T. 57, № 5, стр. 606-611
Параметрический анализ математической модели каталитического осциллятора
А. Я. Наимов a, *, С. Л. Назанский a, В. И. Быков b
a МИРЭА – Российский технологический университет
Москва, Россия
b Институт биохимической физики им. Н.М. Эмануэля РАН
Москва, Россия
* E-mail: alexandermitht@gmail.com
Поступила в редакцию 16.05.2023
После доработки 11.06.2023
Принята к публикации 23.07.2023
- EDN: MGGZKB
- DOI: 10.31857/S0040357123050160
Аннотация
В ряде случаев колебательный режим позволяет провести реакцию с большей селективностью по целевому продукту. Для организации производства в данном режиме необходимо определить условия, при которых возникают колебания, а также рассмотреть сам характер колебаний. В работе проведен параметрический анализ базовой кинетической модели колебательной реакции без автокатализа. Определены границы параметров, при которых система имеет колебательный характер. Построены фазовые портреты системы, бифуркационные кривые. Проведен анализ стационарных состояний системы. Определены тип и количество стационарных состояний. Показано, что система при определенных параметрах имеет три стационарных состояния: два неустойчивых узла и седло. Параметрический анализ базовых моделей позволит подобрать начальные приближения для расчетов более сложных моделей реальных реакций.
ВВЕДЕНИЕ
Исследования колебательных процессов в химии показали, что для наличия колебаний в системе необходимо соблюдение некоторых условий:
1. реакция имеет хотя бы два различных пути и может переключаться между ними;
2. существует некий промежуточный продукт, количество которого направляет реакции по одному из путей;
3. в ходе реакции количество этого продукта будет расходоваться до определенной поры, после чего характер реакции меняется и интермедиат уже начинает накапливаться в системе, и так происходит переключение с одного пути на другой.
Сложные механизмы, описывающие колебательные реакции, содержат большое количество параметров. Математическая модель для данного механизма будет выглядеть как система дифференциальных уравнений высоких порядков. Решение подобных систем напрямую невозможно из-за их сложности. Однако существует метод параметрического анализа нелинейных систем, который позволяет от простого к сложному решать подобные задачи и предполагать механизм колебательных реакций [1].
Все дискуссии, которые велись до сих пор, не касались экспериментальной практики, за исключением упоминаний об открытых и закрытых системах. Колебания концентрации и цвета могут возникать в обычных реакторах идеального смешения или периодического действия, и эти конфигурации обычно сочетаются с механизмом смешивания, чтобы избежать неоднородностей. Введение диффузионных явлений переноса вещества в сочетании с петлей обратной связи реакции и возмущением приводит к морфогенезу в форме волн и пространственных структур [2].
Химические осцилляторы предлагают множество вариантов с точки зрения состава и поведения, но можно ли сконструировать такую систему в соответствии с определенными потребностями? В 80-х годах по этому поводу были проведены некоторые исследования. Разработка генераторов, использующих большое количество соединений со схожими свойствами, помогла бы достичь более конкретных целей [3].
В настоящее время активно изучаются колебательные реакции окислительного карбонилирования [4, 5]. Селективность образования продукта в реакции окислительного карбонилирования фенилацетилена, катализируемой палладием, при прочих равных условиях в колебательном режиме выше. Было установлено, что возникновение колебаний pH совпадает с увеличением скорости потребления фенилацетилена и образования сопутствующих продуктов.
Таким образом, изучение механизмов колебательных реакций позволяет выявить режимы, при которых обеспечивается максимальная селективность. В связи с этим является актуальной задача определения областей рабочих параметров реактора, при которых возникают колебания. Знание данных областей позволит провести оптимизацию реактора по какому-либо технологическому или экономическому критерию.
Экспериментатор в ходе работы неизбежно сталкивается с тем, что все измерения производятся с некоторой погрешностью. Это может быть следствием случайных ошибок или же того, что сам прибор не может быть идеально точным. В результате измерений получаются кривые, характеризующиеся случайными колебаниями измеряемых концентраций. Наличие данных колебаний (шумов) затрудняет создание математической модели системы и ее анализ.
Помимо этого для системы, которая находится в критическом состоянии, малые флуктуации могут приводить к резкому изменению ее характера. Поэтому необходимо рассмотреть влияние колебаний самих параметров на характер поведения всей системы.
Цель настоящей работы: проведение параметрического анализа базовой кинетической модели химической реакции, определение значений параметров, при которых система проявляет колебательный характер, исследование поведения системы при изменении данных параметров.
ТЕОРЕТИЧЕСКИЙ АНАЛИЗ
Настоящая статья посвящена параметрическому анализу каталитического осциллятора.
Рассмотрим механизм реакции, содержащий четыре стадии:
(1)
$\begin{gathered} Z\mathop {\overset {{{k}_{1}}} \longleftrightarrow }\limits_{k - 1} X,\,\,\,\,2Z\mathop {\overset {{{k}_{2}}} \longleftrightarrow }\limits_{{{k}_{{ - 2}}}} {{Z}_{2}}, \\ X + {{Z}_{2}}\xrightarrow{{{{k}_{3}}}}3Z,\,\,\,Z\mathop {\overset {{{k}_{4}}} \longleftrightarrow }\limits_{{{k}_{{ - 4}}}} Y. \\ \end{gathered} $Следуя закону действующих масс, схеме отвечает кинетическая модель:
(2)
$\begin{gathered} \frac{{dx}}{{dt}} = {{k}_{1}}z - {{k}_{{ - 1}}}x - {{k}_{3}}x{{z}_{2}},\,\,\,\,\frac{{dz}}{{dt}} = - {{k}_{1}}z + {{k}_{{ - 1}}}x - \\ - \,\,2{{k}_{2}}{{z}^{2}} + 2{{k}_{{ - 2}}}{{z}_{2}} + 3{{k}_{3}}x{{z}_{2}} - {{k}_{4}}z + {{k}_{{ - 4}}}y, \\ \frac{{dy}}{{dt}} = {{k}_{4}}z - {{k}_{{ - 4}}}y,\,\,\,\,{{z}_{2}} = \frac{{1 - x - y - z}}{2}. \\ \end{gathered} $Константы скоростей химических реакций k1, k–1, k2, k–2, k3, k4, k–4 выступают в качестве параметров модели (2). После приравнивая к нулю правых частей системы (2) получим уравнение стационарности относительно одной переменной х:
(3)
$\begin{gathered} ({{k}_{{ - 2}}} + {{k}_{3}}x)(2{{k}_{1}} + {{k}_{3}}\alpha x)((1 - x)(2{{k}_{1}} + {{k}_{3}}\alpha x) - \\ - \,\,\alpha (2{{k}_{{ - 1}}}x + {{k}_{3}}x(1 - x)) - \\ - \,\,2{{k}_{2}}{{\left( {2{{k}_{{ - 1}}}x + {{k}_{3}}(1 - x)} \right)}^{2}} = 0, \\ \alpha = \left( {1 + \frac{{{{k}_{4}}}}{{{{k}_{{ - 4}}}}}} \right). \\ \end{gathered} $Также при приравнивании к нулю правых частей системы уравнений (2) были найдены зависимости z(x), z2(x), y(x):
(4)
$\begin{gathered} z(x) = \frac{{2{{k}_{{ - 1}}}x + {{k}_{3}}x(1 - x)}}{{2{{k}_{1}} + {{k}_{3}}\alpha x}},\,\,\,\,{{z}_{2}}(x) = \frac{{1 - x - \alpha z(x)}}{2}, \\ y = \frac{{{{k}_{4}}}}{{{{k}_{{ - 4}}}}}z(x). \\ \end{gathered} $Из уравнения стационарности (3) выделим два параметра k2, k–2.
Данное уравнение F(x, k1, k–1, k2, k–2, k3, k4, k–4) = 0 линейно относительно указанных параметров. Поэтому из него можно получить параметрические зависимости k2(x) и k–2(x).
(5)
${{k}_{2}} = \frac{{{{k}_{3}}x{{z}_{2}}(x) + {{k}_{{ - 2}}}{{z}_{2}}(x)}}{{{{z}^{2}}(x)}} = \frac{{{{k}_{3}}x + {{k}_{{ - 2}}}}}{{{{z}^{2}}(x)}}{{z}_{2}}(x),$(6)
${{k}_{{ - 2}}} = \frac{{{{k}_{2}}{{z}^{2}}(x) - {{k}_{3}}x{{z}_{2}}(x)}}{{{{z}_{2}}(x)}} = {{k}_{2}}\frac{{{{z}^{2}}(x)}}{{{{z}_{2}}(x)}} - {{k}_{3}}x.$Пример построения параметрических зависимостей k2(x) приведен на рис. 1.
Рис. 1.
Параметрические зависимости x(k2) при k1 = 1, k–1 = 0.1, k–2 = 0.05, k4 = 0.0032, k–4 = 0.002: 1 – k3 = 1; 2 – k3 = 10; 3 – k3 = 100.

Начальным приближением для подбора параметров k1, k–1, k2, k–2, k3, k4, k–4 служили параметры автокаталитической кинетической модели [9 ] . Также условились, что константа скорости нелинейной стадии k3 превышает значения остальных параметров; константы скоростей прямой и обратной реакции для буферной стадии k4, k–4 имеют наименьшие значения; константы скорости k1, k–1 взяты в отношении 10 : 1, как и в автокаталитической модели; константы скорости k2 : k–2 = = 100 : 1, что предполагает практическую необратимость реакции образования Z2.
Так как фазовое пространство системы имеет размерность три, то процедура анализа устойчивости стационарных состояний (ст. с.) и их бифуркаций становится более трудоемкой по сравнению с предыдущей моделью [9 ] .
Для проверки устойчивости ст. с. воспользуемся первым методом Ляпунова. Необходимо составить матрицу Якоби для системы дифференциальных уравнений и определить корни характеристического уравнения [10 ] .
Элементы матрицы Якоби aij для (2) имеют вид:
Характеристическое уравнение при анализе устойчивости ст. с. является кубическим:
гдеНеобходимое и достаточное условие того, что все характеристические корни имеют отрицательные действительные части, состоит в выполнение условий Рауса–Гурвица σ, ξ, ∆ [11 ] . Так как всегда σ > 0, то бифуркации числа и устойчивости ст. с. определяются из равенства нулю величин ∆ и ξ. Зная параметрические зависимости, можно численно построить соответствующие бифуркационные кривые L∆, Lσ (рис. 2 и рис. 3 соответственно).
Рис. 2.
Кривая кратности L∆(k2, k–2) при k1 = 1.1, k–1 = = 0.12, k4 = 0.05, k–4 = 0.0145: 1 – k3 = 50; 2 – k3 = 78; 3 – k3 = 100.

Рис. 3.
Кривая нейтральности Lσ(k2, k–2) при k1 = 1.1, k–1 = 0.12, k4 = 0.05, k–4 = 0.0145: 1 – k3 = 50; 2 – k3 = = 78; 3 – k3 = 100.

Стоит заметить, что при увеличении значения параметра k3, отвечающего за нелинейную стадию, область множественности стационарных состояний смещается в сторону увеличения соотношения параметров k2 : k–2, что говорит о преобладании скорости прямой реакции над обратной.
Параметрический портрет представлен на рис. 4. В данном случае кривая нейтральности Lσ и кривая кратности LΔ расположены таким образом, что делят пространство на характерные области. При этом, варьируя значения параметров k2, k–2 и оставляя значения остальных параметров неизменными, можно переводить ход реакции в осциллирующий режим. Либо, обладая этими знаниями, можно кратчайшим путем вывести систему из колебательного режима.
Рис. 4.
Параметрический портрет на плоскости (k2, k–2) при k1 = 1.1, k–1 = 0.12, k3 = 78, k4 = 0.05, k–4 = 0.0145: 1 – область с единственным устойчивым ст. с.; 2 – область с единственным неустойчивым ст. с.; 3, 4 – области с тремя ст. с.

Область с тремя ст. с. можно наблюдать и на параметрической зависимости концентрации компонента Х от варьируемого параметра k2. При этом данная область совпадает с областью значений k2 на плоскости (k2; k–2), что еще раз подтверждает адекватность построения областей множественности ст. с. Для значения k–2 = 0.4 построена кривая х(k2) на рис. 5.
Рис. 5.
Параметрическая зависимость х(k2) при k1 = 1.1, k–1 = 0.12, k–2 = 0.4, k3 = 78, k4 = 0.05, k–4 = 0.0145.

Фазовые траектории системы для случая единственного неустойчивого ст. с. (существует предельный цикл) представлены на рис. 6. С любых начальных данных фазовые траектории стремятся к устойчивому предельному циклу. Траектории наматываются на цикл по часовой стрелке.
Рис. 6.
Проекция фазового портрета для системы (1) при k1 = 1.1, k–1 = 0.12, k2 = 6, k–2 = 1, k3 = 78, k4 = 0.05, k–4 = 0.0145: А – стационарное состояние системы.

Из уравнения ст. с. (3) при заданных значениях параметров была рассчитана концентрация компонента Х. Значения Z, Z2, Y были найдены из уравнений их зависимости (4) от переменной Х. Ст. с. отмечено на проекции фазового портрета и имеет следующие значения по каждой из переменных: Х = 0.436106; Z = 0.125632; Z2 = 0.002524; Y = 0.433213.
Для области параметров с тремя ст. с. фазовый портрет представляет собой предельный цикл с тремя особыми точками, две из которых являются неустойчивыми узлами и одна точка седловидного типа (рис. 7).
Координаты особых точек в фазовом пространстве были рассчитаны так же, как и для случая с единственным неустойчивым ст. с. При подборе решения выбирались различные начальные значения переменной Х. Результаты подсчетов приведены в таблице.
Таблица 1.
Стационарные состояния системы.
Стационарное состояние | X | Z(x) | Z2(x) | Y |
---|---|---|---|---|
N2 | 0.116096 | 0.189072 | 0.021429 | 0.651974 |
C1 | 0.206717 | 0.173699 | 0.010312 | 0.598961 |
N1 | 0.392562 | 0.135066 | 0.003314 | 0.465745 |
На рис. 8 видно, как фазовые траектории выходят из окрестности ст. с. и стремятся к предельному циклу. Траектория движения проходит таким образом, что наматывание на предельный цикл проходит по часовой стрелке (как и в случае, когда фазовые траектории выходят из точек концентрационного пространства вне предельного цикла).
Временные зависимости лучше всего демонстрируют колебательный характер системы. Общая площадь под кривой характеризует количество образовавшегося продукта, поэтому от характера колебаний зависит селективность реакции по данному компоненту. Для случая единственного неустойчивого ст. с. представлены на рис. 9. Колебания устойчивые и незатухающие, амплитуда колебаний заметна.
ЗАКЛЮЧЕНИЕ
Таким образом, в ходе параметрического анализа были установлены области изменения констант скорости, при которых исследуемая система реализует одно, два и три стационарных состояния. Полученные состояния отличаются составом реакционной смеси, и, как следствие технологическими показателями. Определены тип и количество стационарных состояний.
Для систем, характер поведения которых сильно зависит от внешних факторов, необходимо исследование их влияния на нее. Колебания устойчивы лишь при некотором наборе и соотношении констант скоростей. В данной работе показано, что их изменение может вывести систему из колебательного режима. Работа в данном направлении позволит в будущем рассмотреть влияние шума на характер колебаний. А что более важно, отделить шумы и составить адекватную кинетическую модель.
Работа выполнена при финансовой поддержке гранта РНФ № 19-19-00620-П.
ОБОЗНАЧЕНИЯ
Список литературы
Bykov V., Tsybenova S., Yablonsky G. Chemical Complexity via Simple Models. De Gruyter Graduate, 2018
Gray C., An analysis of the Belousov-Zhabotinskii reaction // Rose-Hulman Undergrad. Math. 2002. V. 3. P. 1.
Kurin-Csörgei K., Epstein I.R., Orbán M. et al. Systematic design of chemical oscillators using complexation and precipitation equilibria // Nature. 2005. V. 433. P. 139.
Novakovic K., Bruk L., Temkin O. History, versatility and future prospects of oscillatory carbonylation reactions of alkynes // RSC Advances. 2021. V. 11(39). P. 24336.
Parker J., Novakovic K. The Effect of Temperature on Selectivity in the Oscillatory Mode of the Phenylacetylene Oxidative Carbonylation Reaction // ChemPhysChem. 2017. V. 18(15). 1981–1986 (2017)
Наимов А.Я., Быков В.И. Параметрический анализ базовой кинетической модели, содержащей автокатализ // Эл. сб. тр. Первая научно-техническая конференция Московского технологического университета. М., 2016. С. 682
Вольтер Б.В., Сальников И.Е. Устойчивость режимов работы химических реакторов. 2-е изд., перераб. и доп. М.: Химия, 1981.
Слинько М.Г. Некоторые пути развития методов моделирования химических процессов и реакторов // Теорет. основы хим. технологии. 1976. Т. 10. № 2. C. 171.
Дополнительные материалы отсутствуют.
Инструменты
Теоретические основы химической технологии