Формирование математической модели каталитических процессов с переменным реакционным объемом на основе теоретико-графового подхода Текст научной статьи по специальности «Математика»

СПИСОК ЛИТЕРАТУРЫ

1. Branzoi I.V., Iordoc M., Codescu M. Electrochemical studies on the stability and corrosion resistance of new zirconium-based alloys for biomedical applications // Surf. Interface Anal. - 2008. - V. 40. -№ 3-4. - P. 167-173.

2. Rosalbino F., Maccio D., Giannoni P., et al. Study of the in vitro corrosion behavior and biocompatibility of Zr-2.5 Nb and Zr-1.5 Nb-1 Ta (at %) crystalline alloys // J. Mater. Sci.: Mater. Med. - 2011. - V. 22. - № 7. - P. 1293-1302.

3. Дятлова В.Н. Коррозионная стойкость металлов и сплавов / Справочник. - М.: Машиностроение, 1964. - 352 с.

4. Балашова Н.А., Кабанов Б.Н. Электрохимическое поведение циркония в серной кислоте // Доклады АН СССР. - 1957. -Т. 121. - №1. - С. 126-128.

5. Гильман В.А., Колотыркин Я.М., Малкина Р.И. Растворение циркония в концентрированной соляной кислоте // Защита металлов. - 1966. - Т. 2. - № 4. - С. 490-492.

6. Колотыркин Я.М., Гильман В.А. Влияние ионов хлора на электрохимическое и коррозионное поведение циркония // Доклады АН СССР. - 1961. - Т. 137. - № 3. - С. 642-645.

7. Поляков С.Г., Григоренко Г.М., Оноприенко Л.М. и др. Коррозионно-электрохимическое поведение сплава Zr-2,5 Nb и его сварных соединений в серной кислоте // Защита металлов. -1989. - Т 25. - № 3. - С. 419-426.

8. Шавкунов С.П., Толкачев А.Б. Электрохимическое поведение моно- и поликристаллических электродов циркония в растворах серной кислоты // Защита металлов. - 2003. - Т. 39. -№ 3. - С. 253-258.

9. Атанасьянц А.Г, Никитин К.Н. Активация и пассивация при электрохимической обработке циркония в водных и неводных растворах хлоридов // Электрохимия. - 2001. - Т. 37. - № 9. -С. 1049-1056.

10. Глухова А.И., Андреева В.В. Влияние добавок окислителей в растворы соляной кислоты на коррозионное и электрохимическое поведение циркония и сплавов системы цирконий -титан // Журнал прикладной химии. - 1962. - Т. 35. - № 3. -С. 567-572.

11. Антропов Л.И. Теоретическая электрохимия. - М.: Высшая школа, 1984. - 519 с.

12. Felmy A., Girvin D., Jenne E. MINTEQ: A computer program for calculating aqueous geochemical equilibria. - Washington: US Environmental Protection Agency, 1984. - 98 p.

13. Божко П.В., Коршунов А.В., Ильин А.П. и др. Электрохимическое поведение пластически деформированного титана в растворах серной кислоты // Известия Томского политехнического университета. - 2011. - Т 319. - № 3. - С. 17-24.

14. Лурье Ю.Ю. Справочник по аналитической химии. - М.: Химия, 1989. - 448 с.

15. Самсонов Г.В., Борисова А.Л., Жидкова Т.Г. и др. Физико-химические свойства окислов. - М.: Металлургия, 1978. - 472 с.

16. Черняева Т.П., Стукалов А.И., Грицина В.М. Поведение кислорода в цирконии // Вопросы атомной науки и техники. Сер. Физика радиационных повреждений и радиационное материаловедение. - 2000. - № 2. - C. 71-85.

17. Головин А.М., Потапов В.С. К расчету коэффициента диффузии кислорода в циркониевых сплавах методом растворения оксидной пленки // Вестник Моск. ун-та. Сер. 2. Химия. -2006. - Т. 47. - №6. - С. 421-424.

Поступила 06.03.2012 г.

УДК 544.43

ФОРМИРОВАНИЕ МАТЕМАТИЧЕСКОЙ МОДЕЛИ КАТАЛИТИЧЕСКИХ ПРОЦЕССОВ С ПЕРЕМЕННЫМ РЕАКЦИОННЫМ ОБЪЕМОМ НА ОСНОВЕ ТЕОРЕТИКО-ГРАФОВОГО ПОДХОДА

Е.В. Степашина, С.А. Мустафина

Стерлитамакская государственная педагогическая академия им. Зайнаб Биишевой E-mail: [email protected]

На основе теории графов доказаны свойства решения кинетической модели химического процесса. Разработан алгоритм уточнения механизма химической реакции на основе метода анализа графа прямых связей с распространением ошибки. Численный алгоритм апробирован на каталитическом процессе димеризации а-метилстирола.

Ключевые слова:

Граф реакции, дифференциальные уравнения на графах, механизм реакции.

Key words:

Graph of the reaction, the differential equations on graphs, reaction mechanism.

Введение

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

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

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

В настоящее время широко применяется геометрическая трактовка механизма реакции. На основе теории двудольных графов возможен анализ существования и свойств решения системы обыкновенных дифференциальных уравнений кинетической модели схемы реакции [3]. Одним из методов, реализующий графовый подход к механизму реакции, является метод анализа графа прямых связей [4]. Для того чтобы применять этот алгоритм необходимо знать скорости всех реакций в анализируемые моменты времени, то есть иметь только численное решение прямой кинетической задачи. Однако данный метод рассматривает только прямые связи между веществами. Но вещества могут быть связаны через промежуточные вещества, то есть косвенно. В связи с этим требуется использовать метод анализа графа прямых связей с распространением ошибки [5], который учитывает косвенное влияние веществ.

В ходе протекания химических реакций возможно изменение реакционного объема или числа молей реакционной среды. Например, при проведении реакции в цилиндре с подвижной стенкой (поршнем) объем системы можно произвольно менять вне всякой связи с происходящими в системе химическими превращениями. На практике с протеканием газовых реакций при переменном объеме приходится сталкиваться в тех случаях, когда в реакции изменяется число молекул, а давление в системе поддерживается постоянным. Поэтому при разработке математического описания сложного процесса необходимо учитывать изменение числа молей реакционной смеси (или реакционного объема).

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

Геометрическая трактовка механизма реакции

Пусть химическая реакция протекает в т стадий, в которых участвует п веществ А^А^.Д. Схему реакций можно представить в виде:

Еа«А ’ ' = 1>->т’ (1)

к=1 к=1

где ак, - стехиометрические коэффициенты.

Если Ск - концентрация вещества Ак (к=1,...,п), то закон изменения во времени функций С(£),..,СЛ(?) записывают в виде дифференциальных уравнений

—С т

-Сг = ^.к®1 , =Рк -ак, к = 1,-> п (2)

-г ,=1

Здесь оI - скорость г-й реакции. По закону действующих масс

п

о ' = кП с“'к, ' = 1,-, т, (3)

к=1

где к - константа скорости реакции.

Зависимость константы скорости реакции от температуры определяется по уравнению Аррениуса:

_ _Е_

к = к0 ■ в*Т, (4)

где Е- энергия активации; Т- температура протекания реакции; Я - универсальная газовая постоянная.

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

Механизму (1) соответствует двудольный граф Г: А={А1,.,Ап} - множество вершин, соответствующих веществам, участвующим в реакциях; В={Ьь...,Ьт} - множество вершин, соответствующих реакциям. Тот факт, что ак единиц вещества Ак вступает в реакцию Ь, условимся обозначать стрелкой, идущей от вершины Ак к вершине Ь с весовым коэффициентом а|к. Аналогично с помощью стрелок, идущих от вершины Ь к вершине Ак, будем отображать, что рл единиц вещества Ак является продуктом реакции Ь . Каждой вершине Ак соответствует функция концентрации С(), каждой вершине Ь; - функция скорости г-й реакции ю;(/,С). Тогда кинетическая модель схемы реакции (1) представляется системой дифференциальных уравнений на графе Г: т

-Ст =ССу'к°' (г, С), к = п (5)

-г =1 к

При этом функции (О ¡(¡, С) непрерывны по / и С (>) и непрерывно-дифференцируемы по С, причем

о(г, С) > 0 при г > 0, Ск > 0, к = 1,..., п (6)

Анализ кинетической модели химической реакции

На основе свойств теории двудольных графов, сформулируем и обоснуем некоторые утверждения для полученной системы (5).

А-вершину Ак графа Г будем называть непосредственно предшествующей В-вершине Ь, если а1к>0. Функция о ¡(¡,С подчинена вершине Ак, если

ю;(/,С)=0при Ск=0. Таким образом, все функции о ¡(¡,С) подчинены всем А-вершинам, непосредственно предшествующим вершине Ь .

Утверждение. Если Ск(0)=Ск°>0 (к=1,...,п), то «гладкое» решение задачи (5), (6) положительно:

Ск(г) > 0 (г> 0; к = 1,..., п)

на интервале, где это решение существует.

Следствие. Если Ск°>0, к=1,...п, то «гладкое» решение задачи (5), (6) неотрицательно:

Ск(г) > 0 (г > 0; к = 1,...,п)

на интервале, где это решение существует.

Для дальнейших исследований свойств решения введем в рассмотрение линейные формы

о,

Ж =^

1 С

С0

связанные с графом Г.

Лемма. Пусть С(/) - решение задачи (5), (6). Тогда если Я=(ЯЬ...,Я„) является решением системы неравенств

Ц(Я) < 0 (' = 1,...,т), (7)

то функция

(8)

0 < С,(ї) <£ЯкСк0.

(9)

Если суммарная концентрация С = ^ С,.

изме-

няется со временем, то система (2) преобразуется к виду

й (Сх,)

йї

= Хуу® у,(і = 1>->п)

(10)

у=1

с начальными условиями:

X' (0) = X'0, С (0) = С0,

где х=С/С - концентрация г-го компонента в мольных долях.

Систему уравнений (10) замыкает условие нормировки по компонентам реакционной среды:

(11)

Разделим (10) на начальную суммарную кон центрацию С0:

й (Ых,)

йї

'=ХУУЖ7’(І = 1>-> п)

(12)

=1

где М=С/С0 - относительное изменение числа молей реакционной среды; Щ=а /С0 - приведенные скорости химических реакций (=1,...,т). Обозначим

является невозрастающей функцией и

Теорема 1. (Априорная оценка). Пусть существует неотрицательное решение системы неравенств (7) Як>0 (к=1,.,п), причем при некотором /: Я;>1. Пусть С(7) - решение системы (5), (6). Тогда имеет место оценка

Р =£ УуЖу,(, = 1,..., п).

(13)

=1

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

й (Ых,) ..Ох йЫ

у = N—L + х.---.

Теорема 2. (Теорема существования.) Если существует положительное решение Як>0 (к=1,.,п) системы неравенств (7) (или уравнений (8)), то решение задачи (5), (6) с произвольными неотрицательными начальными данными существует на полуоси />0.

Доказательство. Умножив положительное решение системы (7) на достаточно большую константу, всегда можно получить решение, удовлетворяющее условиям Як>1 (к=1,...,п). Тогда оценка

(9) имеет место при каждом / (/=1,...,п). Таким образом, решение существует во всей области 00.

Чтобы определить, происходит ли изменение реакционного объема в ходе реакции, строят матрицу, элементами которой являются стехиометрические коэффициенты веществ Уд ( г=1,...,т,у=1,...,п). Затем находят компоненты вектора 8д по следующему правилу

т

81=Суа (1 = ^п).

'=1

Если хотя бы одна из компонент вектора 8д отлична от нуля, то реакция протекает с изменением числа молей (или реакционного объема).

йї йї йї С учетом (13) и (14) система (12) примет вид:

(14)

ЛГОх,. йЫ „

Ы—L + х,-------= Р.

йї йї

(15)

Отсюда, с учетом условия нормировки (11), получим:

-№ _ ^ Т7 Я и/

~-Г = т, т = § 11.

Выразим — из выражения (15). Тогда система -г

уравнений (12) примет вид обыкновенных нелинейных дифференциальных уравнений:

^ ,(' = п);

йї Ы

йЫ = Р (х, т) йї

(16)

с начальными условиями:

X,.(0) = х“ (' = 1,...,п); N(0) = 1. (17)

Полученная система уравнений (16) с начальными условиями (17) является кинетической мо-

,= 1

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

По теореме 2 решение кинетической модели

(10) существует. В силу утверждения переменная Ы=С/С0>0 на полуоси о0. Учитывая, что кинетическая модель (16), (17) получена из (10) равносильными преобразованиями, можно сделать вывод, что существует решение кинетической модели (16), (17).

Алгоритм уточнения механизма химической реакции

Приведем численный алгоритм уточнения химической реакции, основанный на анализе графа прямых связей с распространением ошибки.

Шаг 1. Выбор целевых веществ.

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

Шаг 2. Построение графа связей между веществами.

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

Шаг 3. Расчет коэффициентов зависимости между веществами.

Чтобы определить степень зависимости между всеми веществами, рассчитывается коэффициент гав [4, 5]:

т

X КОЛ/1

r = -

' AD

где о)1 - скорость г-й реакции; уА1 - стехиометрический коэффициент вещества А в г-й реакции (стехиометрический коэффициент положительный, если А - продукт, и отрицательный, если А - реагент); 8В=1, если вещество В участвует в г-й реакции, 8В=0, в противном случае; т - количество реакций в системе.

Коэффициент зависимости между каждой парой веществ является элементом матрицы. Каждый ее элемент удовлетворяет соотношению 0<ГАВ<1.

Если вещество А связано с веществом С косвенно (реакция вида А^В, В^С), то коэффициент зависимости между веществами А и С рассчитывается по формуле:

где гАВ - коэффициент зависимости между веществами А и В; гВС - коэффициент зависимости между веществами В и С.

Шаг 4. Расчет обобщенньх коэффициентов связи между веществами.

Если в схеме реакций вещество А связано с веществом С как прямо, так и косвенно, тогда рассчитывается обобщенный коэффициент зависимости между веществами А и С:

rac = max {fAC }.

all path i

где максимум определяется по всем путям графа связей, соединяющим каждое целевое вещество с остальными. Любое вещество Х будет выбрано в качестве вещества, связанного с целевым веществом А, если

rax

где е является определенным пороговым значением (0<е<1).

Шаг 5. Формирование сокращенного механизма реакции.

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

Вычислительный эксперимент

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

2X о X2,

2X о X3.

2X ^ Х4,

X2 о X3,

X2 ^ X4.

X3 ^ X4.

Xj + X2 ^X5.

Xj + X3 ^ X5.

Xj + X4 ^ X5. (18)

где введены обозначения веществ X1 - а-метилсти-рол; Х2 - а-димер; Х3 - в-димер; Х4 - циклический димер; Х5 - тримеры.

Кинетические уравнения, соответствующие схеме химических превращений, выписаны согласно (3) и имеют вид:

©ДС, Т) = ^(Т )С12,

©2 (С, Т) = к2(Т )С2,

®э(С, Т) = кз(Т )С12,

©4 (С, Т) = к4(Т )Сз,

©5(С, Т) = к5(Т )С2,

®в(С, Т) = к6(Т )С2,

©7 (С, Т) = к7(Т )Сз,

©8(С, Т) = к8(Т )С2,

©9(С, Т) = к9(Т )С3,

©10 (С, Т) = кю(Т )С1С2,

«11 (С, Т) = кц(Т )С1Сз,

©12 (С, Т) = к12(Т)С1С4.

Кинетическая модель димеризации а-метил-стирола с учетом изменения числа молей в ходе протекания химической реакции представляется системой:

0х1 = р (x, Т) - хЛ(x, Т) йї Ы ,

где

= F (x, T), где N = ,

dt

(19)

(20)

Результаты и обсуждение

Исходя из полученных значений обобщенных коэффициентов связей и точности расчета из схемы реакции были исключены вещества Х4 и Х5. В соответствии с этим сокращенная схема реакции имеет вид:

2X1 о X2,

2Х о X3,

X2 о Xз.

с начальными условиями:

X (0) = х0, ' = 1,...,5, N(0) =1,

где XI - концентрация г-го компонента в мольных долях (номер г компонента соответствует веществу X, г=1,...,5); N - переменный реакционный объем.

Уточним механизм димеризации а-метилстиро-ла с учетом изменения реакционного объема. Решим систему дифференциальных уравнений, представляющую кинетическую модель схемы реакции (19), (20) с учетом изменения реакционного объема, явным методом Рунге-Кутта четвертого порядка. Поскольку основной целью процесса является получение линейных димеров, в качестве целевых веществ были выбраны Х2 и Х3. Расчет проводился с точностью £=0,001 в момент времени /=0,02 ч при температуре Т=373 К. Рассчитаны коэффициенты зависимости глВ для всех веществ. С помощью алгоритма поиска графа в глубину определены все пути, соединяющие каждое целевое вещество с остальными. По найденным путям рассчитаны обобщенные коэффициенты связи веществ:

1 0,000630 0,000177 0 0

0,976226 1 0,000173 0 0

0,661918 0,095594 1 0 0

0,980872 0,016596 0,002157 1 0

1 0,095164 0,995495 0,000911 1

(21)

На рисунке представлено изменение концентраций компонентов в мольных долях сокращенного механизма и концентраций этих же компонентов в исходном механизме.

Как видно из рисунка, уточнение механизма реакции (18) не изменило общую динамику изменения концентраций целевых веществ во времени. Поэтому для анализа кинетических и термодинамических моделей реакции (18) можно использовать схему реакции (21), которая является эквивалентной исходной.

Для автоматизации процедуры уточнения механизма химической реакции в среде визуального программирования Delphi на языке Object Pascal разработан программный комплекс, который позволяет осуществить построение кинетической модели как с изменением, так и без изменения реакционного объема, решить прямую задачу химической кинетики и на основе анализа графа прямых связей с распространением ошибки исключить из механизма вещества, мало влияющие на образование целевых веществ.

Программное обеспечение зарегистрировано в Объединенном фонде электронных ресурсов «Наука и образование» (ОФЭРНиО), Всероссийском научно-техническом информационном центре (ВНТИЦ).

Выводы

Разработана математическая модель каталитического процесса. Принципиальным отличием модели является учет изменения числа молей реакционной среды в ходе химического превращения. На основе теоретико-графового подхода доказана математическая корректность постановки задачи.

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

Разработано программное обеспечение для расчета каталитических процессов. При замене блока реакций программный продукт может быть адаптирован к другим процессам в жидкофазной реакционной системе.

0.00 0.01 0.02

г , ч

Рисунок. Динамика изменения концентраций ¡-го компонента в мольных долях (х, соответствует исходному механизму, х' - сокращенному механизму. Индексы 1 - а-метилстирол, 2 - а-димер, 3 - р-димер)

СПИСОК ЛИТЕРАТУРЫ

1. Brown N.J., Li G. Mechanism Reduction via Principal Component Analysis // International Journal of Chemical Kinetics. - 1997. -V. 29. - P. 393-414.

2. Dickinson R.P, Gelinas R.J. Sensitivity Analysis of Ordinary Differential Equation Systems - A Direct Method // Journal of Computational Physics. - 1976. - V. 21. - P. 123-143.

3. Вольперт А.И., Худяев С.И. Анализ в классах разрывных функций и уравнения математической физики. - М.: Наука, 1975. -395с.

4. Lu Т., Law C.K. A Directed Relation Graph Method for Mechanism Reduction // Proceedings of the Combustion Institute. - 2005. -V.30. - P. 1333-1341.

5. Xia G., Michelangeli D.V., Makar PA. Mechanism reduction for the formation of secondary organic aerosol for integration into a 3-dimensional regional air quality model: а-pinene oxidation system // Atmospheric Chemistry and Physics. - 2009. - V. 9. -P. 4341-4362.

6. Байтимерова А.И., Степашина Е.В., Мустафина С.А. Математическая модель процесса в РИС на двудольном графе // Обозрение прикладной и промышленной математики. - 2010. -Т 17. - Вып. 3. - С. 462-463.

Поступила 27.10.2011 г.