Численное исследование стационарных решений математических моделей химического реактора с неподвижным слоем катализатора Текст научной статьи по специальности «Математика»

МАТЕМАТИКА

Вестн. Ом. ун-та. 2009. № 4. С. 60-65.

УДК 517.95+541.124

И.Д. Макарова, С.Е. Макаров

Омский государственный университет им. Ф. М. Достоевского

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ СТАЦИОНАРНЫХ РЕШЕНИЙ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ ХИМИЧЕСКОГО РЕАКТОРА С НЕПОДВИЖНЫМ СЛОЕМ КАТАЛИЗАТОРА

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

Ключевые слова: математическая модель, численное исследование стационарных решений.

1. Постановка задачи

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

вв+7в=гев-5(е-ег), а дх

в е

Т-1Х(*')еП (1)

(е-еи = о, в^ =е,

(в,вг )= - заданы

Здесь П - полуполоса (0,1) х (0, да), в, вг - температура в реакторе и в холодильнике, С - концентрация реагирующего вещества, Р,у,д,в0 - константы, из них первые три положительны.

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

© И.Д. Макарова, С.Е. Макаров, 2009

дС +дС=«(1 -о)

д ск

р— +— = г(1 -С)ев-3(0-0), ді дх

(2)

—г--г- =3(0-0), (кі)єП,

в дв д дх

Сх=0 = 0, (в-в)х=0 = 0,

вх=1 =в0,

(С,в,в)(=0 - заданы!

Как и в первом случае, П - полуполо-са (0,1) х (0, да), С - концентрация реагирующего вещества, в,е - температура в реакторе и в холодильнике, а,Р,у,8,в0 -константы, из них первые четыре положительны, начальные функции - гладкие и удовлетворяют условиям согласования нулевого и первого порядков.

Классическая разрешимость начально-краевых задач (1) и (2) следует из общих результатов работ [3, 4, 5]. В [1] указаны условия на параметры, при которых существуют стационарные решения задач

(1) и (2).

Данная работа содержит результаты численных расчетов по нахождению нестационарных решений моделей (1) и (2). При проведении расчетов применялись разностные схемы для решения систем дифференциальных уравнений в частных производных с использованием пакета МаШСаа 2001.

2. Модель реактора при реакции нулевого порядка

Для нахождения численного решения задачи (1) применяем неявную разностную схему [6, 7] первого порядка точности на равномерной прямоугольной сетке

= (х. = гк, г = 0,..., Ы, к = —Л = пт, п = 0,1,...},

используя трехточечный шаблон (а) для второго и шаблон (б) для первого уравнений (см. рис. 1):

п + 1

(а)

п +1

(б)

і + 1

і -1

Рис. 1. Трехточечные шаблоны

„.п+1 „,п л „,п+\ „,п+\

иі - иі + 1 иі - Щ-1

ї

8

т в к і = 1,..., N,

уп+1 - уп у”+! - уп+1

= ±- еи■ - — (ип+1 - Уи+1), в в і

т к

і = 0,.,N -1,

^ = з« - о,

< = у0, К =в0, п = 0,1,.

и0 = V0 =во, г = 0,...,N.

Здесь и, и VП - разностные решения для в и вг соответственно. Разрешая первое уравнение относительно переменной <+1, а второе уравнение относительно переменной v”+1 и обозначая г = т/к , получим разностные уравнения .п+1 = виП + гиП-+ +туви" +т8^+

и =-

і = 1,..., N,

в + г 8 + г

иТх =+',

+1 V" + ГУ^1 +т8ип +1 0

V,- =-, Ум =)

(3)

(4)

1 + г8 + г , = N -1, .,0, при следующих начальных данных

и0 = V,0 =в, г = 0,., N. (5)

Алгоритм расчета достаточно прост: сначала по рекуррентному соотношению (3) насчитываются V,+1, а затем, используя только что найденные значения функции V,+1, по формулам (1) слева направо находим значения и”+1 . Расчеты проводились в области изменения параметров, удовлетворяющих условиям устойчивости задачи (1), а именно, при выполнении следующих условий:

8 + уев° < 1, А = 8 / уев° < 1., (6)

8 + уев0 +у/е8/2 < 1. (7)

Аналитически достаточно сложно построить такую область. Заметим, что из четырех параметров @,8,у,в0 наиболее важным является последний в0 - температура теплоносителя на выходе из реактора (из-за наличия экспоненты уев в правой части первого уравнения (1)). При больших значениях параметра в0 при фиксированном значении у (для нашей модели в0 > 2, 15 при у = 0,1) наступает быстрый рост (при численных расчетах

п

п

наступает переполнение) температуры 0 на выходе из реактора, что на практике приводит к разрушению катализатора и самого реактора. При 00 < 2, 15 расчеты были проведены при следующем наборе параметров: в = 1; 8 = 0,03; у = 0,1; к = 0,01; значения т варьировались от 1/300 до 0,1 до установления стационарного состояния с точностью 10-4. Характерное поведение температуры в реакторе и температуры теплоносителя приведено на рис. 2, а значения температуры

0 на выходе из реактора в зависимости от параметра 00 приведено в таблице.

0гп

1-Ь

Рис. 2. Распределение температуры в в реакторе и температуры теплоносителя вг при

0 = 1,9; т = 0,01

Заметим, что температура теплоносителя 0Г остается практически постоянной на протяжении всего реактора.

Зависимость температуры при выходе из реактора 0(1) от 00

00 -0,5 -0,1 0,1 0,5 0,9 1,3 1,7 1,9

6(1) -0,44 -0,01 0,22 0,68 1,19 1,76 2,51 3,05

На рис. 3 показано поведение температуры в реакторе при значении параметра 00 = 2,1, нарушающего условие экспоненциальной устойчивости, для различных значений времени

і = 50,100, 200, 300 .

Для небольших моментов времени значения температуры в реакторе ограничены, но с увеличением итераций, например, при і = 300 значение 0(1) быст-

ро растет, что приводит к резкому нагреву реактора.

Рис. 3. Распределение температуры 01, 02, 03, 04 в реакторе при і = 50, 100, 200, 300 соответственно

3. Модель реактора при реакции первого порядка

Для нахождения численного решения задачи (2) применяем неявную разностную схему [6, 7] первого порядка точности на равномерной прямоугольной сетке

®т = {к = ік, і = 0,..., N, к = —; і = пт, п = 0,1,...},

для первого и второго уравнений выбираем шаблон (б), для третьего - шаблон (а) (см. рис. 1):

сп+1 _ сп сп+1 _ сп+1 с с + с,- _ С-1 =«(1 -спе ,

т

, .п+1 „п

к

1 ип+1 - ип+1 у

----+-------^=*- (1 - с )еи-(ии+1 - У+1),

т в к в Р ''

і = 1,..., N.

Уп+1 - У Уп+1 - Уп+1 ,

■ ±. + '+1 і =8(ип -У+1),

к

, = 0,.,N-1, С, = сй, иП = V,, =00, п = 0,1,.

и0 = V =в0, , = 0,.,N с,, и,, V, - разностные решения для с,в,вг соответственно. Разрешая первое уравнение относительно переменной СП+1, второе уравнение относительно переменой и"+1, а третье уравнение относительно ^я+1 и обозначая г = Тк , получим разностные уравнения:

с„+1 = <+_<+/ +т«(1 - с,п) еи

1 + г

сп+1 = с

і = 1,..., N,

ви" + гипЦ + ту{1 - сп) еи" + т8у"

ип+1 =

п+1 п+1

и = V

в + г8 + г і = 1,., N,

(8)

(9)

= V, + у;1 +т8ип +1 =в

1 1 + г8 + г ’ N 0, (10)

1 = N -1, .,0,

при следующих начальных данных

и1 = ^ =в0, 1 = 0,.,N, С0 = °. (11) Схема расчета заключается в следующем: сначала по рекуррентному соотношению (10) насчитываются vIn+1, а затем, с использованием только что найденных значений функции vIn+1, по формулам (9) слева направо находятся значения и,п+1, а по формулам (8) также слева направо находятся значения с,я+1.

Число параметров, определяющих устойчивое стационарное решение задачи

(2), равно трем: а,8,у. Это следует из условия существования и экспоненциальной устойчивости в Ь2 - норме стационарного решения начально-краевой задачи (2): 8 < 2е-(^“+1). В отличие от реакции нулевого порядка в данном случае нет зависимости от параметра в0 , хотя в первом и втором уравнениях задачи (2) также присутствует множитель ев . Это объясняется наличием дополнительного множителя 1 - х , отвечающего за скорость реакции от степени превращения (концентрации) вещества х , и при численных расчетах достаточно быстро степень превращения х ^ 1. Характерное поведение концентрации вещества С , температуры реактора в и температуры холодильника вг приведено на рис. 4; а = 2,3 ; в = 0,2 ; 8 = 0,03 ; у = 1,3 ; в0 = 2,1; т = 0,005 ; к = 0,01.

Интересные выводы при численных экспериментах следуют при изменении параметра а, отвечающего за скорость протекания реакции. На рис. 5 приведено поведение концентрации вещества С в зависимости от скорости реакции при а = 0,7;1,5;2,3.

Рис. 4. Распределение концентрации С, температуры в реакторе в и температуры холодильника вг

Рис. 5. Распределение концентрации С при параметре а = 0,7; 1,5; 2,3

При значении параметра а = 2,3 достигается полная степень превращения (концентрация С = 1) при к = 0,2, если а = 1,5, то полная степень превращения достигается при к = 0,4, и наконец, если а = 0,7, то полная степень превращения достигается при к = 1, т. е. в конце реактора. Последний режим, таким образом, является неэффективным на практике. Заметим, что влияние других параметров 8,у,в0 при фиксированных значениях

а = 1,5; 2,3 не существенно.

4. Подсчет числа стационарных режимов

С понятием устойчивости тесно связан вопрос о числе возможных стационарных режимов при одних и тех же значениях всех параметров, причем некоторые из этих режимов могут быть неустойчивыми [8]. Выясним количество стационарных решений для моделей (1) и (2). Пусть вектор-функция у( к) = (у1( к), у2( к))* -

стационарное решение задачи (1), удовлетворяющее краевой задаче:

сV „ ^ / ч

—1 = Ге1 -881-V2),—2 = 81 -v2), хе(0,1), Ох Ох

(12)

Л(0) = V2(0), V2(1) = в0.

Полагая г = v1 - v2 и вычитая из первого уравнения второе, получим:

Ог Ог , .

— = уе 1, —1 =----------8г, х е (0,1),

Ох Ох Ох V ’ (13)

г(0) = 0, ^(1) = в0.

Интегрируя второе уравнение системы (3), имеем:

1

Vl( х) = г( х) - г(1) +81 г (хО =

х

1

= г (х) + в0 +81 г (хО (14)

и

Ог г(x)+8J г( *) Оя г(x)+гJ г( *) О*

— = у е х евв = уе(в е х . (15) Ох

Таким образом, г - монотонная функция, и из (5) получаем:

Ог Ох ’

полагая w =

2 I I -8г—;

Ох ^ Ох) Ох

имеем: w'w = w2 -8wг и

w = Сег +8(г +1), где С = уев - 8(1 + г1)е~21 ,

г = г(1).

Интегрируя уравнение w = — по х

Ох

от 0 до 1 (г от 0 до г1), получим уравнение на г1, которое необходимо для определения стационарных решений задачи (1):

} Ог

= 8.

(16)

0 (е6>0 - (1 + г1)е г)ег +1 + г)

Каждое решение уравнения (16) г1 = v2 (1) - v1 (1) = v2(1) - в0 > 0 позволяет найти 1 (1) и решить задачу Коши для исходной задачи (12) по начальным условиям при х = 1. Таким образом, каждый положительный корень г1 уравнения (16) дает новое стационарное решение задачи (1), т. е. число различных положительных корней уравнения (16) определяет число стационарных решений уравнения (1).

В случае реакции первого порядка пусть (р, v1, v2) - стационарное решение (2):

~~ = а(1 -р) е\

Ох

= Л(1 -Ж18, - V, (17)

Оv2 , ч

— = -8(Vl -v2), х е (0,1),

р(0) = 0, ^(0) = V2(0), V2(1) = в0.

Аналогичная замена переменных г = v1 - v2 приводит систему (17) к виду:

О р Ох Ог

Ох.............. (18)

ОV

—1 = у(1 - р) е 1 - 8г,

Ох

р(0) = 0, г(0) = 0, v2(1) = в0.

Из первого и второго уравнений (18) получаем соотношение на р :

= а(1 -р) е \

= Г(1 -р)е'1, х е (0,1),

а

р = — г .

У

(19)

Второе уравнение системы (18) с учетом (14), (19) преобразуем к виду:

Ог г( х )+8х | г (*) О*

— = (у-аг)ев ■ е х Ох

(20)

Функция г является монотонно возрастающей, Ог]Ох > 0, логарифмируя (20) и дифференцируя полученное уравнение, получим следующее уравнение для г :

(у/а - г)г" + (г')2 = (у/а - г)(г' - 8г)г'.

При г" = w , имеем линейное уравнение относительно w:

(

w

1

Л

-1

w = -8г.

,Ф- г ,

Находя решение этого уравнения и интегрируя г " = w (г от 0 до г1 = г(1)), получаем функциональное уравнение для нахождения г1, определяющего стационарное решение задачи (2):

}-------------------* , - =8. (21)

0 (у/а- г )ег (аев / 8+ I"——О*

** -у/а

Для нахождения числа стационарных решений уравнений (1), (2) будем численно строить графики левых частей уравнений (16), (21) при фиксированных для уравнения (16) значениях ye0a / 8, а для уравнения (21) - значениях у/а, аe0a / 8 и разных значениях 8 . Число параметров, характеризующих режим в реакторе в случае реакции нулевого порядка, равно двум: уев / 8 и 8 . На рис. 6 приведенные расчеты при уев / 8 = 0,05;0,1;0,2;0,5;1 показывают, что данная задача может иметь один или два стационарных режима, а при некоторых значениях параметров даже ни одного.

5

Ц х,0.05)

Г(х,0.13 Г(х,0.23 Дх,0.5)

В1(х)

Рис. 6. Решение уравнения (16), определяющего количество стационарных режимов в реакторе при реакции нулевого порядка

При 8 = 2 из условия уев° / 8 = 0,1, выбрав у = 0,01, а в0 = 2,996 с точностью до 0,001, получаем корни уравнения (16): г11 = 0, 276 и г12 = 2 . Решая дважды задачу Коши (12) с начальными условиями ^(1) = 3,272 , ^(1) = 2,996 и ^(1) = 4,996 , ^(1) = 2,996, получаем два различных стационарных решения задачи (1), причем первое из них, соответствующее меньшему значению г11 , совпадает с решением, полученным по разностной схеме

(3)-(5).

Число параметров, характеризующих режим в реакторе в случае реакции пер-

вого порядка, равно трем: у/а, аев / 8 и 8 . На рис. 7 приведены графики правой части уравнения (21) при фиксированном значении параметра у/а = 10 и значениях выражения

аев / 8 = 0,005 ; 0,008; 0,01; 0,03; 0,1; 0,5. Как и в первом случае, задача (2) имеет одно или два стационарных решения.

Ф(х, 0.005,(3)

Ф(х, 0.008,С»

...... 2

Ф(х,0.01,<3)

Ф(х,0.03,<3)

Ф(х,0.1,<3)

Ф(х,0.5,<Э) 1

5ВД

° 0 2 4 6 8 10

X

Рис. 7. Решение уравнения (21), определяющего количество стационарных режимов в реакторе при реакции первого порядка

ЛИТЕРАТУРА

[1] Зеленяк Т.И. Об устойчивости решений одной

смешанной задачи // Доклады АН СССР. 1966. Т. 171. № 2. С. 266-268.

[2] Зеленяк Т.И. К вопросу об устойчивости реше-

ний смешанных задач для одного квазилинейного уравнения // Дифференциальные уравнения. 1967. Т. 3. № 1. С.19-29.

[3] Аболиня В.Э., Мышкис А.Д. Смешанная задача

для почти линейной гиперболической системы на плоскости // Мат. сб. 1960. Т. 50. № 4. С. 423-442.

[4] Годунов С.К. Уравнения математической физи-

ки. М.: Наука, 1979.

[5] Воробьева Е.В., Романовский Р.К. Об устойчи-

вости решений задачи Коши для гиперболических систем с двумя независимыми переменными // Сиб. матем. журнал. 1998. Т. 39. № 6. С. 1290-1292.

[6] Калиткин Н. Н. Численные методы. М.: Наука,

1978.

[7] Методы моделирования каталитических процессов на аналоговых и цифровых вычислительных машинах / М.Г. Слинько и др. Новосибирск: Наука, 1972.

[8] Число и устойчивость стационарных режимов на непористом зерне катализатора для сложной реакции / М. Г. Слинько и др. // Доклады АН CCCР. 1972. Т. 204. № 6. С.1399-1402.