Mathematical modelling of a self-oscillating catalytic reaction in a flow reactor

Толық мәтін

1. ВВЕДЕНИЕ

Гетерогенные каталитические реакции являются нелинейными системами, которые могут протекать вдали от термодинамического равновесия. При их изучении возможно наблюдение сложного динамического поведения, которому соответствуют автоколебания, волновые процессы и детерминированный хаос [1–4]. Пространственно-временные структуры, включающие в себя различные виды волновых процессов, были зафиксированы при проведении каталитических реакций как в высоком вакууме на гранях монокристаллов [3], так и при атмосферном давлении на массивных и нанесенных катализаторах [5]. Для обнаружения пространственных структур на поверхности катализатора при высоком вакууме были разработаны методы фотоэлектронной эмиссионной спектроскопии (PEEM) и электронной микроскопии низких энергий (LEEM) [6, 7]. При атмосферном давлении до недавнего времени применялся только метод инфракрасной термографии (IRT), с помощью которого были зафиксированы тепловые структуры на поверхности катализатора в ходе протекания реакций окисления СО [8], Н₂ [9], NH₃ [10] и C₂H₄ [11]. Было показано, что в ходе колебательных режимов на поверхности катализатора возникают горячие “пятна”, и разница температур между горячей и холодной зонами поверхности катализатора может достигать 190°C [12]. В связи с этим в первых работах по математическому моделированию пространственно-временных структур на поверхности катализатора при атмосферном давлении были использованы неизотермические модели, в которых скорость реакции аппроксимировалась либо уравнением первого порядка [12], либо имела вид полинома [13, 14]. М. Шейнтух с сотр. был первым, кто для математического моделирования пространственных неизотермических структур в ходе реакции окисления СО на поверхности Pd-стекловолоконного катализатора применил неизотермическую модель типа STM [15, 16]. Авторам удалось успешно описать наблюдавшиеся экспериментально в данной системе “дышащие” тепловые структуры и “горячие пятна”.

В последнее время для изучения пространственных структур при атмосферном давлении был успешно использован метод визуального наблюдения за состоянием поверхности металлических катализаторов в ходе колебательных реакций [17–22]. Этот метод оказался наиболее информативен и полезен при наблюдении пространственных структур в реакциях окисления метана [17, 18], этана [19–20] и СО [21, 22] на поверхности никелевых и кобальтовых катализаторов, для которых существует значительное различие в цвете восстановленной и окисленной поверхности. Было показано, что во всех данных системах в ходе автоколебаний скорости реакции имеются кинетические пространственно-временные структуры, не обусловленные неизотермическим протеканием процесса. Природа автоколебаний и волновых явлений в данном случае определяется периодическими процессами окисления и восстановления поверхности катализатора, затрагивающими также объем катализатора. При этом наиболее активной в реакциях окисления была светлая (восстановленная) поверхность катализатора, а темная (окисленная) проявляла меньшую активность. Аналогичные закономерности были получены при изучении автоколебаний скорости реакции окисления пропана на никелевом катализаторе методом in situ XPS [23]. Одновременные видеонаблюдения за цветом его поверхности in situ XPS измерения подтвердили кинетическую природу волновых явлений на поверхности никелевой фольги в ходе реакции окисления пропана и их связь с процессами окисления и восстановления катализатора [24].

Несмотря на большое число реакций, в которых обнаружены кинетические волновые явления на поверхности катализатора, существует только одна работа, посвященная математическому моделированию изотермических волн окисления и восстановления, возникающих в ходе реакции окисления СО на никелевом катализаторе [25], в которой показано, что кинетические волновые явления на его поверхности могут возникать в реакторе идеального смешения (CSTR) при наличии влияния внешней диффузии на скорость реакции. Цель настоящей работы – анализ условий возникновения кинетических волновых явлений в проточном реакторе, где возникают градиенты концентраций реагирующих веществ как в объеме реактора, так и на поверхности катализатора. Для получения наиболее общих закономерностей, относящихся сразу ко многим реакциям окисления, была выбрана абстрактная кинетическая STM-модель, в которой колебания обусловлены процессами окисления–восстановления поверхности катализатора [26]. В этой математической модели нетривиальное пространственно-однородное стационарное решение не существует, и традиционные методы теоретического анализа пространственных структур, описанные в работе [27], неприменимы. В связи с этим в настоящем исследовании будут представлены результаты численного изучения распределенной математической модели.

2. ТОЧЕЧНАЯ КИНЕТИЧЕСКАЯ МОДЕЛЬ

Модель типа STM включает в себя стадии адсорбции, десорбции, реакции, окисления и восстановления поверхности, которые представлены в табл. 1, где в качестве X используется обозначение для окислителя-кислорода, а Y – для восстановителя в реакционной смеси. V₀, W₀ – концентрации окислителя и восстановителя в газовой фазе, x, y – адсорбированные частицы окислителя и восстановителя, x, y – степени покрытия поверхности окислителем и восстановителем, z – степень покрытия поверхности оксидом металла, k_i – константа соответствующей стадии. Согласно механизму колебаний по типу STM, адсорбция окислителя происходит только на восстановленной поверхности.

 

Таблица 1. Механизм реакции

Стадия	Описание	Скорость стадии
Адсорбция X	X + 2M → 2x	R1 = V0k1(1 – x – y)2(1 – z)2
Адсорбция Y	Y + M → y	R2 = W0k2(1 – x – y)
Десорбция Y	y → Y + M	R3 = k3y
Реакция	x + y → XY + 2M	R4 = k4xy
Окисление катализатора	M + x → Mx	R5 = k5x(1 – z)
Восстановление катализатора	Mx + y → M + XY	R6 = k6yz

 

Система кинетических уравнений, соответствующая механизму реакции, представленному в табл. 1, имеет следующий вид:

$\begin{array}{l}\mathrm{d}x/\mathrm{d}t=V_0{k}_1{(1-x-y)}^2{(1-z)}^2-k_{4}xy-k_{5}x(1-z),\\ \mathrm{d}y/\mathrm{d}t=W_0{k}_2(1-x-y)-k_{3}y-k_{4}xy-k_{6}yz,\\ \mathrm{d}z/\mathrm{d}t=k_{5}x(1-z)-k_{6}yz.\end{array}$ (1)

Систему можно привести в безразмерный вид, введя безразмерное время τ = V₀k₁t, как это было сделано в работе [26]:

$\begin{array}{l}\mathrm{d}x/\mathrm{d}\tau ={(1-x-y)}^2{(1-z)}^2-k_{4}xy-k_{5}x(1-z),\\ \mathrm{d}y/\mathrm{d}\tau =p(1-x-y)-k_{3}y-k_{4}xy-k_{6}yz,\\ \mathrm{d}z/\mathrm{d}\tau =k_{5}x(1-z)-k_{6}yz.\end{array}$ (2)

где $p=\left(W_0{k}_2\right)/\left(V_0{k}_1\right)$ является параметром и ${\kappa }_i=k_i/\left(V_0{k}_1\right)$.

Константы адсорбции k_i (i = 1, 2) находили на основе кинетической теории идеального газа по формуле (3):

$k_i=\frac{S_i}{N_i}\sqrt{\frac{RT}{2\pi M_i}},$ (3)

где R, Дж K^–1моль^–1 – универсальная газовая постоянная, T, K – температура газовой смеси, M_i – масса одного моля газа, кг/моль, N_i – количество активных центров на поверхности катализатора, моль/м², и s_i – коэффициент прилипания. В настоящей работе константы адсорбции были рассчитаны для O₂ (M₁ = 0.032 кг/моль) и CO (M₂ = 0.032 кг/моль). Величина N_c = 2.11 × 10^–5 моль/м² соответствует поверхности металлического Pd. Предполагается, что коэффициенты прилипания s₁ = s₂ = 0.1.

В качестве начальной аппроксимации величины констант других стадий механизма реакции были взяты из работы [26], при которых безразмерные параметры модели (2) принимают следующие значения: κ₃ = 10^–1, κ₄ = 10², κ₅ = 10–6, κ₆ = 2 × 10^–7.

Для проведения расчетов давление окислителя в исходной реакционной смеси было принято постоянным и равным 2 кПа (V₀ = 2.83 × 10^–1 моль/м³) и температура газа T принималась 850 K, что соответствует условиям работы [26]. В дальнейшем параметр p будет использоваться в качестве управляющего параметра, определяемого внешними условиями. Продолжение нетривиального стационарного состояния по параметру p позволило определить интервал существования колебательного решения (предельного цикла), величина которого составляет 0.165 < p < 0.525. При p < 0.165 поверхность катализатора находится в окисленном состоянии, а при p > 0.525 в восстановленном состоянии.

На рис. 1 показано колебательное решение системы (2), полученное при p = 0.5, при использовании начальных условий x(0) = y(0) = z(0) = 0. Период колебаний равен 1.2 × 10⁷ безразмерных единиц, что соответствует 50 с. Из рисунка видно, что колебания имеют релаксационный характер. Во время колебаний происходит как окисление поверхности $\frac{\mathrm{d}z}{\mathrm{d}t}>0$, так и восстановление поверхности $\frac{\mathrm{d}z}{\mathrm{d}t}<0$, причем во время окисления скорость реакции выше, чем во время восстановления. Период колебаний определяется величинами констант к₅ , к₆, тогда как отношение данных констант к₅ / к₆ обуславливает длительность стадий окисления и восстановления во время предельного цикла.

 

Рис. 1. Колебательное решение системы (2) при p = 0.5 (а); колебания скорости реакции R при p = 0.5 и остальных параметрах, соответствующих величинам (3) (б).

 

3. РАСПРЕДЕЛЕННАЯ МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Моделью проточного реактора является прямоугольный канал шириной L (м) и высотой L/2. Плоский квадратный катализатор размером L × L (м²) расположен в нижней части канала. Газовая смесь, содержащая восстановитель Y и окислитель X, пропускается через канал с линейной скоростью u (м/с).

Пусть ξ обозначает координату вдоль реактора по направлению потока газа, а катализатор расположен на участке 0 < ξ < L. Предположим, что распределения концентраций реагентов V(ξ, t) и W(ξ, t) (моль/м³) по длине реактора газа описываются с помощью системы уравнений типа реакция–диффузия–конвекция. В безразмерных переменных эта система имеет следующий вид:

$\begin{array}{l}\partial v/\partial \tau +\frac{1}{{\tau }_L}\partial v/\partial \zeta =\\ =d_v{\partial }^{2}v/\partial {\zeta }^2-a_{v}v{(1-x-y)}^2{(1-z)}^2,\\ \partial w/\partial \tau +\frac{1}{{\tau }_L}\partial w/\partial \zeta =\\ =d_w{\partial }^{2}w/\partial {\zeta }^2-a_w\left[pw(1-x-y)-{\kappa }_{3}y\right],\\ \partial x/\partial \tau =j(x,s)+v{(1-x-y)}^2\times \\ \times {(1-z)}^2-{\kappa }_{4}xy-{\kappa }_{5}x(1-z),\\ \partial y/\partial \tau =j(y,s)+pw(1-x-y)-\\ -{\kappa }_{3}y-{\kappa }_{4}xy-{\kappa }_{6}yz,\\ \partial z/\partial \tau ={\kappa }_{5}x(1-z)-{\kappa }_{6}yz,\end{array}$ (4)

где ζ = ξ/L – безразмерная длина, а v = V/V₀ и w = W/W₀ – безразмерные концентрации реагирующих веществ в газовой фазе. Параметры системы (4) в безразмерной форме имеют следующий вид: ${\tau }_L=V_0{k}_1{t}_L,d_v=\frac{D_X}{V_0{k}_1{L}^2},d_w=\frac{D_Y}{V_0{k}_1{L}^2},$ $a_v=\frac{2N_c}{V_{0}L},a_w=\frac{2N_c}{W_{0}L}$, где $t_L=L/v$ – время контакта; D_X и D_Y – коэффициенты диффузии реагентов X и Y в газовой фазе.

Модель (4) учитывает диффузию адсорбированных частиц x и y по поверхности катализатора, а поверхностный оксид Мх считается неподвижным. В настоящей работе рассматривается миграция прыжкового типа, т.е. частицы x или y могут перепрыгивать из занятого активного центра поверхности в другой свободный активный центр. Таким образом, интенсивность миграции определяется не только концентрацией поверхностных частиц, но и числом свободных центров на поверхности катализатора [28]. Уравнение такого типа миграции можно представить в следующем виде:

$\partial x/\partial \tau =H_x(s{\partial }^{2}x/\partial {\varsigma }^2-x{\partial }^{2}s/\partial {\varsigma }^2)(\equiv J(x,s\left)\right),$

где H_x, м²/с – коэффициент, аналогичный коэффициенту диффузии, а s = 1 – x – y обозначает долю свободных активных центров на поверхности катализатора. Аналогичное уравнение записывается также для миграции частиц y. Безразмерное выражение для J записывается как:

$j(x,s)=\frac{H_x}{V_0{k}_1{L}^2}\left(s{\partial }^{2}x/\partial {\varsigma }^2-x{\partial }^{2}s/\partial {\varsigma }^2\right),$

Для расчета изменения концентраций реагентов в газовой фазе модель (4) была дополнена следующими граничными условиями Данквертса на входе в реактор (ξ = 0):

$\begin{array}{l}u\left(V-V_0\right)\left|{}_{\xi =0}\right.=D_X\partial \xi V\left|{}_{\xi =0}\right.,\\ u\left(W-W_0\right)\left|{}_{\xi =0}\right.=D_Y\partial \xi W\left|{}_{\xi =0}\right.;\end{array}$ (5)

и граничными условиями Неймана на выходе из реактора (ξ = L):

${\partial }_{\xi }V\left|{}_{\xi =L}\right.=\mathrm{0,}{\partial }_{x}W\left|{}_{\xi =L}\right.=0.$ (6)

Для степеней покрытий поверхности x, y, z на входе и выходе из реактора использовали граничные условия Неймана.

Коэффициенты диффузии реагирующих веществ в газовой фазе принимали равными D_X = D_Y = D = 5 × 10^–5 м²/с. Коэффициенты диффузии реагирующих веществ по поверхности принимали равными H_x = D/10, H_y = 2H_x. Время контакта t_L принимали равным 0.05 с, что соответствовало линейной скорости газа 0.2 м/с.

4. РЕЗУЛЬТАТЫ РАСЧЕТОВ

4.1. Численный метод

Для численного решения уравнений (4) распределенной модели был применен метод конечных разностей. Производные по ζ в правой части системы были заменены конечными разностями на равномерной сетке из N узлов, ζ_k = h/2 + (k − 1)h, k = 1, 2, ... , N, где h = 1/N является шагом сетки. Производную первого порядка ∂_ζf заменяли разностью (f_k − f_k−1)/h (1-й порядок точности). Производную второго порядка аппроксимировали центральной разностью ∂_ζζf ≈ (f_k−1 − 2f_k + f_k+1)/h² (2-й порядок точности), где f_k = f(ζ_k).

Члены f₀ и f_N+1 определяли с использованием граничных условий (5), (6). Из условий Данквертса при ζ = 0 получаем

$u(f_0-1)=\frac{D}{Lh}(f_1-f_0).$

Из условий Неймана при ζ = 0, 1 следует, что

$f_0=f_1,f_{N+1}=f_N.$

Таким образом, система (4) преобразовывалась в систему 5N обыкновенных дифференциальных уравнений для 5N-мерного вектора неизвестных q = {v_k, w_k, x_k, y_k, z_k, k = 1, 2, ... , N}.

$\mathrm{d}q/\mathrm{d}\tau =F\left(q\right(\tau \left)\right).$ (7)

Матрица Якоби системы (7) является 11-диагональной, поэтому разумно использовать метод разреженных матриц для численного решения системы. Для выбора числа N систему (7) решали с разными значениями N и сравнивали полученные решения. Оказалось, что до величины N = 1000 решения системы заметно различаются, хотя эти различия постепенное уменьшались. Начиная со значения N = 1000 и до значения N = 10000 решения системы были практически неразличимы. Следовательно, все численные решения системы (7), представленные ниже, были получены при N = 1000 (h = 0.001).

Для решения системы (7) использовали пакет программ Matlab “ode15s”. Расчеты проводили с относительной точностью 10^–6, а в областях качественного изменения поведения системы при изменении параметра p точность решения составляла 10^–9.

4.2. Результаты численного моделирования

Динамическое поведение модели (7) изучали при варьировании параметра p с шагом 0.001, начиная со значения p₁ = 0.162 (точка зарождения колебаний в точечной модели). Колебания с заметной амплитудой возникали в распределенной системе при достижении параметром p величины равной 0.235. Увеличение параметра p приводило к возрастанию периода и амплитуды колебаний. Более того, колебания в распределенной системе (4) существуют и при значении параметра p > p₂, т.е. за пределами области колебаний в точечной системе. На рис. 2 и 3 представлены результаты численного решения системы (4) при p = 0.8. Из рис. 2 видно, что решение для z имеет вид “дышащей” пространственно-временной структуры. Сначала фронт окисления распространяется против потока газа от конца пластины катализатора (ζ = 1) в ее начало (ζ = 0). Однако он не доходит до начала пластины, постепенно замедляется и, наконец, останавливается. Далее фронт окисления начинает двигаться в направлении потока газа и возвращается к концу пластины (ζ = 1). Линии уровня функции z(ζ, t), приведенные на рис. 2, показывают участки “наступления” фронта, где точка z = c движется справа налево, а также участки “отступления” фронта, где эта точка перемещается в противоположном направлении. Из рисунка видно, что чем меньше величина z, тем дальше по пластине происходит распространение волны окисления.

 

Рис. 2. Пространственно-временная диаграмма степени окисления катализатора z(ζ, t) с представленными линиями z = const для z = 0.3, 0.5, 0.7, 0.9 при p = 0.8.

 

Рис. 3. Пространственно-временные диаграммы для концентраций x(ζ, t) (а) и y(ζ, t) (б). Черная линия – линия уровня z(ζ, t) = 0.7.

 

Соответствующие пространственно-временные диаграммы для концентраций x(ζ, t) и y(ζ, t) вместе с линией уровня для z = 0.7 представлены на рис. 3. Линию уровня z можно разделить на две части, связанные с особенностями распространения фронта окисления. В нижней части линии уровня z фронт окисления распространяется с конца пластины катализатора в ее начало против потока газа, а в верхней части линии уровня z фронт окисления движется в противоположном направлении. На рис.3 также видно, что точка z = 0.7 перемещается против потока газа по поверхности, покрытой преимущественно окислителем (верхняя панель рис. 3), при этом степень покрытия поверхности окислителем x снижается с уменьшением ζ, что вызывает замедление скорости стадии окисления и продвижения волны окисления. С другой стороны, уменьшение ζ приводит к росту концентрации восстановителя (нижняя панель рис. 3), и обратное движение точки z = 0.7 происходит уже по поверхности, покрытой восстановителем.

Временная зависимость безразмерных концентраций окислителя и восстановителя на выходе из реактора показана на рис. 4. Из рисунка следует, что в газовой фазе при p = 0.8 происходят регулярные периодические колебания, которые сохраняются при возрастании параметра p до 0.94. Последующее увеличение параметра p приводит к усложнению колебательного режима, и, как видно из рис. 5, возникают сложные периодические колебания, когда за один период колебаний наблюдаются три последовательных пика. Дальнейший рост параметра p вызывает нарушение регулярности колебаний. Динамическое поведение при p = 1.01, показанное на рис. 6, представляет собой нерегулярные колебания, которые обуславливаются тем, что фронт окисления в этом случае произвольным образом останавливается на разных участках пластины катализатора.

 

Рис. 4. Временная зависимость безразмерных концентраций окислителя (v) и восстановителя (w) на выходе из реактора.

 

Рис. 5. Сложные периодические колебания безразмерных концентраций окислителя (v) и восстановителя (w) при p = 0.97.

 

Рис. 6. Нерегулярные колебания безразмерных концентраций окислителя (v) и восстановителя (w) при p = 1.01.

 

При p = 1.069 колебания исчезают, и в системе устанавливается стационарное состояние, которое является неоднородным по длине пластины катализатора и в котором левая часть пластины катализатора восстановлена (z ≈ 0), а правая часть окислена (z ≈ 1). Увеличение p приводит к стационарному состоянию, в котором вся пластина восстановлена. Детали этого перехода и распределение степени окисления по длине катализатора в стационарном состоянии при разных значениях p показаны на рис. 7. Из рисунка видно, что восстановление пластины катализатора начинается с ее заднего края. При возрастании параметра p функция z(ζ) становится унимодальной с максимумом в центральной части пластины. Рис. 7 также показывает, что при определенных значениях параметра p существует пространственно-неоднородное стационарное состояние, в котором центр пластины окислен, тогда как концы пластины полностью восстановлены.

 

Рис. 7. Распределение степени окисления по длине катализатора в стационарном состоянии при разных значениях величины p.

 

5. ЗАКЛЮЧЕНИЕ

В настоящей работе впервые проведен анализ возникновения пространственно-временных структур в проточном реакторе в изотермических условиях с учетом диффузии реагирующих веществ как в газовой фазе, так и по поверхности катализатора. В качестве кинетической модели была выбрана STM-модель, которая является основой для описания колебательных режимов с учетом окисления и восстановления поверхности катализатора.

Представленные выше результаты показывают, что в данном случае могут возникать три типа пространственных кинетических структур. Это неоднородные по пространству: 1) стационарное состояние, 2) периодическая по времени волновая структура, 3) апериодическая по времени “дышащая структура”. Определяющими факторами в возникновении пространственных структур в рассматриваемой модели являются наличие автоколебаний в кинетической модели и неоднородное распределение реагентов в газовой фазе.

Показано, что природа возникновения волн окисления и восстановления поверхности катализатора связана с наличием градиента концентраций реагирующих веществ в проточном реакторе и на поверхности катализатора. Тип наблюдаемой пространственной структуры зависит от параметра p, который равен отношению скоростей стадий адсорбции реагирующих веществ. Наиболее сложная апериодическая “дышащая структура” возникает в области границы района колебаний при больших p. Показано, что нерегулярность колебаний скорости реакции обусловлена случайной позицией места остановки распространения волн окисления.

Одним из важных результатов исследования является установление большего района колебаний в распределенной системе по сравнению с кинетической моделью. Это связано с неоднородностью распределения концентрации реагирующих веществ по объему реактора и по поверхности катализатора.

В настоящей работе приведена зависимость решения модели только от параметра p. Более подробная информация о влиянии других параметров на динамическое поведение системы представлена в препринте [29].

ФИНАНСИРОВАНИЕ

Работа выполнена в рамках государственного задания № 122040500058-1.

КОНФЛИКТ ИНТЕРЕСОВ

Авторы заявляют об отсутствии конфликта интересов, требующего раскрытия в данной статье.

 

Сокращения и обозначения: STM – математическая модель Sales–Turner–Maple; PEEM – фотоэлектронная эмиссионная спектроскопия, LEEM – электронная микроскопия низких энергий, IRT – инфракрасная видео термография, XPS – рентгеновская фотоэлектронная спектроскопия, CSTR – реактор идеального смешения.

Авторлар туралы

N. Peskov

Moscow State University

Хат алмасуға жауапты Автор.
Email: peskovnick@gmail.com

Faculty of Computational Mathematics and Cybernetics

Ресей, Leninskie Gory, Moscow, 119991

M. Slinko

Semenov Institute of Chemical Physics

Email: peskov@cs.msu.ru
Ресей, Kosygina Str., 4, Moscow, 119991

Әдебиет тізімі

Қосымша файлдар