Об устойчивости режимов каталитического окисления газовой смеси (двухфазная модель процесса) Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 3, № 4, 1998

ОБ УСТОЙЧИВОСТИ РЕЖИМОВ КАТАЛИТИЧЕСКОГО ОКИСЛЕНИЯ ГАЗОВОЙ СМЕСИ (ДВУХФАЗНАЯ МОДЕЛЬ ПРОЦЕССА)

И. С. Вержбицкая, П. Г. ИцковА, А. Т. Лукьянов Казахский государственный национальный университет им. Аль-Фараби, Алматы, Казахстан

A two-dimensional two-phase model of heat and mass transfer under the oxidation of a gaseous mixture in a cylindrical continuous reactor with a fixed-base layer of the catalyst is suggested, allowing for the radial distribution inhomogeneity of the layer porosity. The dynamic behavior of a system with an exothermic reaction is studied approximately analytically and numerically on the basis of the zero-dimension analogue of this model. The equations of bifurcational diagrams and parametric equations of the non-uniqueness and neutral stability boundaries on the plane of the effective temperature - effective concentration have been obtained. The effect of flow-catalyst mass-transfer coefficient on the domains with different number and stability of stationary regimes has been found and unique and non-unique steady and oscillatory oxidation regimes have been numerically implemented. The results obtained are in good agreement with the earlier data on the catalytic oxidation stability of a pseudohomogeneous model and the published solutions of unidimensional heterogeneous problems.

Математическое моделирование окисления газовой смеси в неподвижном слое катализатора на основе квазигомогенных (двумерной и нульмерной) моделей продемонстрировало существование неединственных стационарных состояний и автоколебательных стационарных режимов, позволило проанализировать устойчивость окисления к внешним периодическим воздействиям [1]. Показано согласие полученной на плоскости определяющих параметров области неединственности с установленной экспериментально и рассчитанных распределений температуры [2] с измеренными [3].

Однако квазигомогенные модели не позволяют проанализировать влияние тепло- и массообмена между потоком и поверхностью зерен катализатора на режимы каталитического окисления. Для выполнения такого анализа естественно рассмотреть гетерогенную модель.

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

© И. С. Вержбицкая, П. Г. Ицкова, А. Т. Лукьянов, 1998.

1. Математическая формулировка задачи

В безразмерных переменных процесс описывается следующей системой уравнений:

F0 > 0: 0 < х < 1, 0 < r < 1, д$к Ак д2§к Ак 1 д ( д§к\

Fl(e)Щ = AG+ nAGГдГ [rnrr) - B(^K - ^g) +0УКГ*(1) £д$с = д^ - Pe ^g + n IА (r^ ^ + B (^ - (2)

дFo дх2 дх r дг \ дг ) '

дУ = LexЙ2 - Pe£ + nLer11 () + ЗД* - У), (3)

дF0 дх2 дх r дг \ дг __

£Ш = -Bi(yK - y) + DayKrv(^к). (4)

Выполняются следующие граничные условия: на входе в реактор (х = 0) —

t = BW - Л),

д^

дх

Pe(tf - tf*),

дУ Pe / ч

дх = Lex(У - У*);

на выходе из него ( х = 1 ) —

д$К д^ ду

дх дх дх

на оси (r = 0) —

д#к = 0 ^^G = дУ = 0

дr ' дr дr

и стенках реактора (r = 1 ) —

0 -те = ^ = 0; (5)

= B;W

дr

BiW(^G -

д^К = BiW(^K - *»),

дr

W

BiD (У - У~)-

дУG „.W/

дr

Начальные условия взяты в виде

Fo = 0: 0 < r < 1, 0 < х < 1, ^к = ^g = У = Ун, Ук = Ун- (6)

Здесь

§ = ЯТ/Е, х = х/1, г = г/п, Г0 = аг/!2,

Г\(е) = 7(1 - е), 7 = Скрк/рсро, У = С/Со*, П = 12/г2,

Ба = 12к0/а, 8 = Daq, д = QЯ/рКcКE, Ре = и!/а,

Ьех = Ох/а, Ьег = Бг/а, В = аБуд12/Ао, В1 = вБуд1/а,

Б1Кх = а^/Ак, = а2п/Хо, = а3п/Ак, Б1С = ав1/Б, (7)

где — температура; С — концентрация реагента; и — скорость потока; , к0, Q — энергия активации, предэкспоненциальный множитель, тепловой эффект реакции; Я — универсальная газовая постоянная; р, с — плотность и теплоемкость; О, А, а = А/рсро — коэффициенты диффузии, тепло- и температуропроводности; г, X, г — время и пространственные координаты; ук — концентрация газовой смеси на входе в реактор; Со* — концентрация реагента у поверхности катализатора; а, а1, а2, аэ — коэффициенты теплоотдачи между газом и катализатором, на входе в реактор, между газом и стенкой, между стенкой и катализатором соответственно; ап, в — коэффициенты массообмена с окружающей средой и на поверхности катализатора соответственно; Буд = Б(1 — е)/(2Угр) — удельная поверхность зерен в единице объема слоя; гр — радиус частиц катализатора; г, I, Б, У — радиус, длина, площадь боковой поверхности и объем реактора; у*, §* — концентрация и температура газа до входа в реактор; Ре, Ье, Ба, Б1 — числа Пекле, Льюиса, Дамкеллера, Био; д — тепло реакции. Индексы О и К относятся к газовой и твердой фазам, г, х — радиальной и продольной составляющим, вх, /ш, н и то — к входу и стенкам реактора, началу процесса и окружающей среде; Т, Б — к тепловому и диффузионному процессам. Коэффициенты переноса и тепло- и массообмена определяются эмпирическими зависимостями [4]. Пористость катализатора изменяется в радиальном направлении по закону

е(г) = е {1 + ехр[3г4/гр(г/г4 — 1)]} , (8)

где е — среднее значение пористости, г и гр — соответственно радиусы реактора и частицы катализатора.

2. Нульмерная модель

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

Р1(е) ЩО = т12 (§эф! — — В (¿к — ¿о) + Ук8^(4) = Р (4Л,У,Ук) , (9)

о

¿Го

0

тэ4 (§эф2 — ¿о) + В (ёк — ¿о) = Р2 (ёк, ёо,У, Ук) , (10)

¿У ( ё ё \

Т = т56 (Уэф — У) + В1 (Ук — У) = Рэ \ ,^о,У,Ук^ , (11)

е Л = —В (Ук - У) + УкВаг^ (¿к) = Р ¿с, У, Ук) , (12)

где ¿к, ¿с, У> Ук — значения температуры, концентрации реагента в центре полуслоя, где предположительно находится фронт реакции. В дальнейшем индекс """ опускается. Здесь

Ш12 = т1 + т2; тз4 = т3 + т4; т,5б = т5 + тб;

т1 = Ак/Ас2Б1кх/(1 + Б1кх/2), т2 = Ак/АспБ^/(1 + Б^/2), тз = Ре(2 + Ре/2)/(1 + Ре/2), т4 = пБ^/(1 + Б^/2)[2 + 1/(2г»)], т5 = Ре[2 + Ре/(2Ьел)]/[1 + Ре/(2Ье)л)], т6 = ЬегБЬгс[2 + 1/(2г)]/(1 + Б10/2) — (13) эффективные коэффициенты тепло- и массообмена,

¿эф! = (т^* + т^те)/т,12, ¿эф2 = (т^* + m4¿те)/mз4,

Уэф = (т5У* + теУ^)/т5б —

эффективные значения температуры и концентрации, г = 0,1/2,1.

Течение процесса зависит от 10 безразмерных эффективных физико-химических параметров процесса. Расчеты выполнялись при фиксированных значениях параметров: т12 = 321.1; тз4 = 152.9; т5б = 127.85; В = 0.325; Уэф = 0.00375; ¿эфх = 0.0515; ¿эф2 = 0.05; Ба = 2.18 • 108; 6 = 2.04 • 1010. Исследовалось влияние массообмена на поверхности зерна катализатора, эффективных температуры и концентрации смеси на возможные режимы процесса.

3. Координаты и число стационарных состояний

Координаты стационарного состояния $ks, Vks, Vs в системе (9)-(12) являются решениями четырех трансцендентных алгебраических уравнений:

Pi(^Ks,^Gs,VKs,Vs) = 0, i =1, 4. (14)

Для построения параметрических зависимостей в качестве бифуркационного параметра выбираем эффективную температуру. Тогда, согласно (14), уравнение бифуркационной диаграммы на плоскости , $ks запишется в виде

о о , Втз4 /0 q \ BiVs^r v (^Ks) /. - \

** = ЙК' + mi2(m34 + B) (ЙК' - - (Bi + B2rv№<,))mi2' (15)

где Vs = т5бУэф/ {^55 + B^rv(#Ks)/[Bi + B2rv(#Ks)]}.

На рис. 1 представлены возможные виды бифуркационных диаграмм при различных значениях концентрации уэф (B1 = 490.99). Температура в стационарном состоянии определяется графически, абсциссой точки пересечения прямой

= const (16)

с бифуркационной кривой (15). Как видно из рисунка, в системе (9)-(12) возможны от одного (кривая 1) до трех (кривые 2-4) стационарных состояний; высокотемпературное, низкотемпературное и неустойчивое среднетемпературное, которое реализуется только

Рис. 1. Влияние концентрации уэф на стационарные состояния: уэф = 0.000375 (1), 0.001 (2), 0.002 (3), 0.0025 (4).

принудительно. Критические условия "воспламенения" и "потухания", соответствующие экстремумам кривых 2-4, зависят от параметра уэф. Концентрация уэф слабо влияет на условия воспламенения $кв и весьма значительно на условия потухания $кп.

Границы области неединственности стационарных состояний на плоскостях параметров определяются уравнением (15) совместно с условием касания прямых (16) и бифуркационных диаграмм в экстремальных точках д$эф1 /д§Кз = 0, т. е.

Уэф = {1 + Втз4/[шх2 (тз4 + В)]} х

х {Ш56 - В1В2Ту(^)/[В1 - В2Ту(^)]2}/ (в25т'у(&квт§6)) . (17)

Уравнения (15), (17) являются параметрическими уравнениями границы неединственности стационарных состояний на плоскости д^эф1, уэф.

4. Параметрические уравнения границ устойчивости

В предположении квазистационарности массообмена между потоком и катализатором из уравнения (12) получим

УК = В1 у/В + БаТу (^к)]. Тогда порядок динамической системы (9)-(12) понизится до третьего:

d#

К

dFо

т12 (^эф1 - #к) - В (^к - #о) +

£

+В1у6Ту (^к)/[В1 + БаТу (0к)] = Р1 (#к, ^о, у), -г^ = тз4 (^эф2 - #о) - В (0К - ^о) = Р2 (#к, У),

dFг

о

dFо

= т5б (уэф - у) - В1 уБаТу($к)/В + БаТу ($о)] = Рз ($к, , у)

(18)

Устойчивость стационарных состояний исследуется на основе линейного приближения системы (18), т.е.

АС,

(19)

где

дк - дкз дс - дс, У — Уз

А = (аз)

/ дР1 дР1 дР1

ддк ддс дУ

дР2 дР2 дР2

ддк ддс дУ

дРз дРз дРз

\ ддк ддс дУ )

Локальная устойчивость зависит от знака действительной части собственных чисел матрицы , которые задаются корнями характеристического уравнения

А3 + а А2 + ДА + 0 = 0,

где

а = —эрА,

0 = — аеА

33

Д = ^ — аи аи)-

¿=1 3=2

Границы а = 0, Д = 0, 0 = 0 на плоскости дэфг, уэф определяются параметрическими уравнениями (15) и

где

Уэф = {1 + ВВтъ(дк)/ Кб (В1 + В2Г.0(дк))]}Уз,

Уяа = {т12 + тз4 + Ш56 + 2В+ +ВхВ2Тъ (дк, )/ [В1 + В2Т.и (дк, )] [В1 + В2Т.и (дк, )]2}/ [в216т'ь (дк) , УяА = {^12^34 + Ш12 В + 1П34В + (Ш12 + тз4 + 2В )(Ш56 + В1В2 т^ (дк, )/[В1 + +В2Т*(дк,)]} /[В^(дк,)] (тз4 + В + т5б)/ [В1 + В2Т*(дк)]2 ,

У^ = [В1 — В2Т„(дк,)] /[(тм + В)5В1ТЬ(дк)] [В(т12 + тз4) + +Ш12Ш34] [В1 — В2ТУ (дк,)] / [В1 — В2ТУ (дк )(т12Шз4 + В (Ш12 + тз4))/т5б].

(20)

(21) (22)

(23)

Расчеты показали, что граница 0 = 0 (клин) совпадает с границей неединственности (15), (17). Граница нейтральной устойчивости совпадает с границей Д = 0 (линия с петлей). Граница а = 0 находится полностью внутри границы неединственности 0 = 0 и не вносит ничего нового в деление плоскости на области возможных режимов окисления (рис. 2, кривая 3).

Для получения наиболее полной информации о возможных состояниях системы (9) -(12) необходимо проанализировать взаимное расположение кривых Д = 0, 0 = 0. На

г/эфхю4

.В !. 1 \ 1 V 4 \ / \\

— III \Ш\\ \ А \ * \ \

•А V ^ ^ Л

\ I \ > , \ " >4. V -<2 \ \ \

\ \ •V. II \ N У

0.052 0.060 иэф1

Рис. 2. Области возможных режимов окисления при 1 = 490.99.

рис. 2 приведены области возможных режимов окисления на плоскости параметров , уэф. Видно, что граница однозначности В = 0 (кривая 1), нейтральной устойчивости Д = 0 (кривая 2) делят параметрическую плоскость на области единственных, множественных, автоколебательных режимов окисления. В областях I, II стационарный режим единственный: без колебаний (кривая 1), автоколебательный (II). Области Ш-1У отвечают трем положениям равновесия, среднее из которых является неустойчивым и не реализуется. Два других положения равновесия в области III устойчивы, а в области У — неустойчивы. В области IV неустойчив высокотемпературный и устойчив низкотемпературный, в области VI, наоборот, устойчив высокотемпературный и неустойчив низкотемпературный.

На рис. 3 области возможных режимов окисления при В1 = 22.05 сравниваются с полученными ранее на основе псевдогомогенной модели [1]. Как видно, области неединственных стационарных режимов (III, У) и области единственных автоколебательных режимов (II), полученные на основе гетерогенной модели, значительно больше, и, кроме того, учет массообмена между потоком и катализатором приводит к появлению областей (IV, VI) с неустойчивым высокотемпературным (IV) или низкотемпературным (VI) стационарными состояниями в области неединственности. Численная реализация этих режимов показала их автоколебательный характер. Из расчетов следует, что существует интервал значений коэффициента массобмена на поверхности катализатора, при которых наблюдается наилучшее согласование результатов для квазигомогенной и двухфазной моделей: В1 = 14-24.

Рис. 3. Границы неединственности и нейтральной устойчивости для двухфазной (линии 1, 2) и псевдогомогенной (12') моделей.

УэфХЮ4

" N \ > I; Г / \

' \\ 4 \\ 4 \ \ х

2 1 ^ V "Л \ \ \ \ \ \ *

\\ ч ^ 1 V ч \ _ у

0.05 0.058 иэф]

Рис. 4. Влияние коэффициента массообмена 1 на границы неединственности (линии 1, 1') и нейтральной устойчивости (2, 2'); 1 = 490.7 (1), 41.4 (2).

Рис. 4 иллюстрирует влияние коэффициента массообмена В1 на положение границы

УК

0.08

0.04

а

£ •

0.08

0.04

.5 1 1 б

Б» 2 - -б1

0.001

0.003 УК

0.001

0.003 УК

Рис. 5. Фазовые портреты для устойчивых единственного (а) и неединственных (б) высокотемпературного (линия 1) и низкотемпературного (линия 2) режимов.

неединственности © = 0 и нейтральной устойчивости А = 0. Видно, что с увеличением интенсивности массообмена Б\ область неединственности расширяется (см. кривые 1, 1'), при этом потухание происходит при более низких значениях концентрации уэф, а область единственных автоколебательных режимов внутри петли уменьшается (кривые 2, 2). Последнее связано с ослаблением условия возникновения автоколебательных режимов, когда эффективная скорость массопереноса меньше эффективной скорости теплопереноса [2], [5]. Расчеты показали, что теплообмен на поверхности катализатора незначительно влияет на положение границ неединственности и нейтральной устойчивости.

5. Результаты численной реализации

На рис. 5, 6 приводятся результаты численного решения задачи (9)-(12) методом Рунге— Кутта. Стационарные состояния, полученные из бифуркационной диаграммы, показаны точками на фазовых портретах (см. рис. 5 и 6, б) и пунктиром при изменении температуры во времени (рис. 6, а).

Фазовый портрет на рис. 5, а иллюстрирует единственный устойчивый стационарный режим окисления, соответствующий варианту на рис. 2. Видно, что со временем система приходит в стационарное состояние (точка Б).

На рис. 5, 6 показан выход системы на неединственные стационарные состояния, соответствующие вариантам Б и на рис. 2. Как и ожидалось, в обоих случаях высокотемпературный режим устойчив (линии 1 на рис. 5, б, 6), низкотемпературный режим устойчив для варианта и автоколебательный для варианта (линии 2 на рис. 5, б, 6). Изменение температуры во времени варианта показывает, что стационарная температура в устойчивом высокотемпературном режиме (линия 1, рис. 6, а) близка к полученной из бифуркационной диаграммы, а автоколебания происходят вблизи стационарного состояния низкотемпературного режима (линия 2, рис. 6). Расчеты показали, что температура катализатора ниже, чем в потоке газа (см. линии 1', 2! на рис. 6, а). Это объясняется тем, что зажигание происходит вблизи входа в реактор [2] и движущимся потоком тепло переносится быстрее, чем по слою катализатора. Этот результат, полученный на нульмерной

0.04

Ус а

= —Л=Г — . -г.

Г vк г

да* 2' Ш /г- 1у|П

1 "К 12 1 !

ук 0.12

0.04

б

\

Б VI V«

0.04

£

г/кх

Рис. 6. Изменение со временем температуры (а), катализатора (1, 2), газа (12') в автоколебательном низкотемпературном и устойчивом высокотемпературном режимах и фазовый портрет типа предельного цикла для низкотемпературного режима (б).

модели, описывающей поведение процесса во фронте реакции, согласуется с рассчитанными значениями температуры катализатора и потока во фронте и за фронтом химической реакции [6, 7].

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

Список литературы

[1] Вержбицкая И. С., ИцковА П. Г., Лукьянов А. Т. ТОХТ, 24, 1990, 412-416.

[2] Вержбицкая И. С., ИцковА П. Г. Математическое моделирование нестационарных процессов. Алма-Ата, 1988, 36-42.

[3] Каьтнорр О., Уоктмеуек Б. ОНвш. Engng вег. 35,1979, 49.

[4] Чудновский А. Ф. Теплофизические характеристики дисперсных материалов. Физматгиз, М., 1962.

[5] Лукьянов А. Т., Артюх Л.Ю., ИцковА П. Г. Математическое моделирование задач теории горения. Наука, Алма-Ата, 1981.

[6] Фатеев А. Г. Тепломассоперенос при фазовых и химических превращениях. Минск, 1990.

[7] Стегасов А.Н., ШигАРов А. Б., Кириллов В. А. ТОХТ, 23, №3, 1989, 351-356.

[8] werzhbitzkaya i.s., itskova p.g., lukyanov a.t. Mathematical simulation combustion of gas mixture in chemical reactor. Model calculation including deactivation and regeneration of catalyst. In: Int. Symp. "Chemistry of Flame Front". Almaty: Казак Университет^ 1997, 50.

Поступила в редакцию 3 марта 1998 г.