CC BY
76
УДК 544.43
ИДЕНТИФИКАЦИЯ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ РЕДУЦИРОВАННЫХ СХЕМ РЕАКЦИЙ
В.А. Вайтиев, Е.В. Степашина, С.А. Мустафина
Стерлитамакский филиал Башкирского государственного университета E-mail: stepashinaev@ya.ru
Разработан алгоритм решения обратной задачи для редуцированных схем реакций, позволяющий определить кинетические параметры схем реакций меньшей размерности с сохранением динамики концентраций целевых веществ. Проведена идентификация математической модели редуцированной схемы реакции димеризации a-метилстирола, для которой найдены константы скоростей стадий и значения энергий активации.
Ключевые слова:
Обратная задача, редукция кинетических схем, механизм реакции.
Введение
В химической промышленности широко распространены нестационарные технологические процессы, связанные с проведением химических реакций, содержащие большое число реагентов. Такие процессы описываются в общем случае системами нелинейных дифференциальных уравнений, в которых число уравнений определяется количеством участвующих в реакции веществ. Однако экспериментальным способом измерению доступны концентрации не всех веществ, так как время их образования и расходования может быть малым по сравнению со временем, за которое происходит существенное изменение состава реагирующей смеси. При этом для анализа механизма реакции требуется, как правило, исследовать динамику изменения концентраций лишь некоторых веществ. Помимо этого, высокая размерность систем дифференциальных уравнений, описывающих математические модели реакций, влечет за собой необходимость выполнения большого объема вычислений. Поэтому поиск рациональных подходов к замене исходной системы на эквивалентную ей систему с сохранением динамики концентраций целевых веществ представляет научный и практический интересы.
В этом случае возникает необходимость создания комбинированного численного алгоритма «сокращения» (редукции) кинетической схемы сложной химической реакции на основе метода анализа графа прямых связей с распространением ошибки для заданного момента времени протекания реакции как с постоянным, так и с переменным реакционным объемом. Построение математических моделей редуцированных схем реакций влечет за собой необходимость решения задачи идентификации математической модели реакции, то есть обратной задачи химической кинетики.
Постановка задачи
Кинетический анализ химических систем на основе математической обработки экспериментальных данных, полученных в ходе протекания химической реакции, условно включает следующие основные этапы:
1) разработка математического описания химического процесса - построение кинетической модели исследуемой химической реакции;
2) решение прямой задачи химической кинетики выбранным методом - расчет состава многокомпонентной реагирующей смеси и скорости реакции на основе заданной кинетической модели с известными параметрами (константами скоростей, энергиями активации каждой стадии химической реакции);
3) разработка алгоритма решения обратной кинетической задачи и ее решение с целью восстановления (уточнения) на основе экспериментального материала вида кинетической модели и кинетических параметров исследуемого механизма;
4) сравнение расчетных и экспериментальных данных (как правило, концентраций одного или нескольких исходных веществ и продуктов реакции), полученных в ходе кинетического эксперимента с помощью различных методов анализа;
5) вывод об адекватности модели и кинетических параметров исследуемой реакции. При необходимости уточнение кинетической модели, поиск новых кинетических параметров.
Решение обратной кинетической задачи тесно
связано с формулировкой прямой кинетической задачи, то есть с разработкой математического описания для расчета состава реакционной смеси и скоростей стадий реакции на основе кинетической модели.
Обратная кинетическая задача представляет собой задачу минимизации функционала отклонения между расчетными и экспериментальными данными:
Q=Е Е|хр> - х#| ^ т*п (1)
,=1 ] =1
где хр - расчетные значения концентраций веществ; хЭ - значения концентраций веществ, полученные экспериментальным путем; I - количество точек эксперимента; п - количество веществ.
В результате сокращения схемы реакции система уравнений, описывающих ее математическую
модель, имеет меньшую размерность по сравнению с исходной. Динамика концентраций веществ во времени редуцированной схемы реакции описывается системой обыкновенных дифференциальных уравнений:
Сх.
— = Г (к, Х).
Ж
где х=(х1,х2,^,хг) - вектор концентраций веществ; к=(к1,к2,.,кт) - вектор констант скоростей стадий, вид функции / определяется схемой реакции.
В качестве неизвестных параметров выступают константы скоростей стадий к, определяемые на основе уравнения Аррениуса:
( E. \
kj (T) = k0J exp —j , (j = 1,..., m)
(2)
где k0 - вектор значений предэкспоненциальных множителей (их размерность совпадает с размерностью к); Ej - значение энергии активации j-й стадии (кДж/моль); T - температура протекания реакции (К); R - универсальная газовая постоянная (8,31 Дж/(моль-К)).
Для решения задачи идентификации математической модели редуцированной схемы реакции необходимо рассчитать значения кинетических констант k0j и энергии активации E. Решение обратной задачи для редуцированной схемы реакции дает возможность оценить корректность замены исходной схемы реакции ее сокращенной схемой на основе сравнения численного решения прямой задачи для исходной и редуцированной схем реакций.
Процедура решения обратной задачи состоит в поиске констант скоростей стадий, минимизирующих функционал (1). Для минимизации функционала (1) будем использовать метод Хука-Дживса, который представляет собой комбинацию исследующего поиска с циклическим изменением переменных и ускоряющего поиска по образцу.
Сформулируем алгоритм поиска кинетических констант.
Шаг 1. Ввод исходных данных: t - время протекания реакции; xdwn - вектор начальных концентраций веществ; х8ксп - матрица экспериментальных концентраций в фиксированные моменты времени; s - минимальное значение функционала (критерий остановки поиска); ht - шаг вариации по каждой координате (константе скорости) в ходе исследующего поиска (i=1,...,m); et - минимальное значение шага по каждой координате (константе скорости) в ходе исследующего поиска; maxi -максимальное количество итераций поиска; [kdlim,kup] - область поиска по каждой координате (константе скорости).
Шаг 2. Задаются стартовая точка k1=(k11,k21,km1), где m - количество кинетических параметров, номер текущей константы s:=1, номер итерации i:=0. Вычисляется значение функционала (1) Q(k1) для набора констант k1.
Шаг 3. Решается прямая задача химической кинетики при двух наборах кинетических пара-
метров к^к/ЛЧ^,...,^1) и ^(к/А1-^,..,^1). Вычисляются значения функционала (1) <2(к/+) и <2(к/_) при этих наборах параметров; г':=г'+1.
Шаг 4. Если Q(k1)<Q(k1s+) и Q(k1)<Q(k,!-), то шаг Н, уменьшается. В противном случае к1 присваивается один из векторов к1, к1+, соответствующий минимальному значению из функционалов Q(ks1+), Q(ks1-).
Шаг 5. Если э<т, то я:=8+1 и переход на Шаг 3, иначе путем исследующего поиска по всем координатам (Шаг 2-Шаг 4) получен новый набор констант к2=(к12,к22,кт2) с наименьшим значением функционала (1). Если условия окончания работы алгоритма не выполняются (значение функционала (1) меньше е или достигнуто максимально допустимое значение количества итераций алгоритма поиска 1=тах1), то переход на Шаг 6.
Шаг 6. Поиск по образцу. Рассчитывается вектор констант к3=(к13,к23,кт3) по формуле к3=к‘+Я(к2-к1), где Я - параметр алгоритма, обычно выбирается равным 2.
Шаг 7. Исследующий поиск (Шаг 2-Шаг 4) для нового набора констант к3 за исключением того, что шаги вариации по константам Нв на этой фазе метода не уменьшаются. В результате получается набор констант к4=(к14,к24,к,4).
Шаг 8. Если вектор к4 отличен от к3, то к‘:=к2, к2:=к*и переход на Шаг 6. В противном случае к‘:=к2и переход на Шаг 3.
Вектор значений энергии активации Е=(Е1,Е2,.,Ет) рассчитывается исходя из уравнения Аррениуса (2):
Е
1п к = 1п к0-------.
0 КТ
если известны значения констант скоростей при разных температурах. Пусть найдены константы скорости реакции при пТ различных значениях температуры. Они удовлетворяют системе уравнений:
Е
1пкт = 1пко(; = 1,...,п ). (3)
КТ ,
ln k0 = a, - r = b, lnkT = yj, T = xj (j = 1,...,nT).
Введем обозначения:
е
К
J
Тогда система уравнений (3) примет вид:
yj = а + Ьх} (] = 1,...,Пт ).
Для определения значений энергий активации по стадиям реакции воспользуемся методом наименьших квадратов, согласно которому необходимо минимизировать функционал:
Q (а, Ь) = £((а + Ьх}) - у) )2.
1=1
Исходя из найденных коэффициентов а и Ь определяются значения энергии активации Е и константы к0:
к0 = еа, Е = -ЬК.
Вычислительный эксперимент
Проведем идентификацию кинетических параметров промышленно значимой редуцированной схемы реакции димеризации а-метилстирола в присутствии цеолитного катализатора МаНУ. Продукты реакции (линейные и циклические димеры) находят практическое применение в качестве пластификаторов, модификаторов полимеров, каучу-ков в производстве синтетических масел и др. Совокупность химических превращений, описывающих данную реакцию, и соответствующие кинетические уравнения имеют вид [1]:
2 X1 о X2, ш1 (C, Т) = k1C12 - k10C2
2 ^ о X 3, ш2(С, Т) = к2С12 - к11С3
2 XI ^ X4, ш3(С, Т) = к3С12,
X2 оXз, шА(С,Т) = к4С2 -к^Сз
X2 ^X4, ш5(С, Т) = к5С3,
Xз ^ X4, ®б(С, Т) = кбС2Сі,
X1 + X2 ^ X5! ш7(С, Т) = к7С1С3,
X1 + X, ^^, ш8(С, Т) = к8С1С„
XI + X4 ^ X5, ш9(С, Т) = к9С2С4,
где введены обозначения Х1 - а-метилстирол; Х2 -а-димер; Х3 - ^-димер; Х4 - циклический димер; Х5 - тримеры, где а^х) - скорость £-й стадии (кмоль/(м3-ч)) (£=1,...,9); С=(С!,.,С5) - вектор концентраций компонентов (кмоль/м3); к=(к1,.,к12) -вектор кинетических констант скоростей ]-й реакции (м3/(кмоль-ч)) (/=1,...,12).
Значения кинетических констант и энергии активации (табл. 1) были рассчитаны в лаборатории приготовления катализаторов Института нефтехимии и катализа РАН (г. Уфа) с учетом наличия цеолитного катализатора. Константа скорости ;-й реакции рассчитывается через выбранную опорную температуру Топ=373 К по формуле:
к1 (Т) = к1 (Топ)ехр (Е- ^ .
Таблица 1. Кинетические параметры процесса димеризации а-метилстирола в присутствии катализатора ЫаНУ при температуре 373 К
№ к (373 К), м3/(кгкат-ч) Е, кДж/моль № к (373 К), м3/(кгкат-ч) Е, кДж/моль
1 61,357 196 7 0,019308 247
2 8,9534 263 8 41,556 194
3 7,7916 259 9 0,03662 115
4 1,1693 238 10 0,04547 279
5 0,11922 275 11 0,0995 204
6 0,12041 127 12 0,05132 138
Кинетическая модель димеризации а-метил-стирола с учетом изменения числа молей в ходе протекания химической реакции представляется системой [2]:
Шхі = ^ (х,Т) - х^т (х,Т)
Ж N ’
где ^ =£уЖ , і =
к=1
^ (х, Т), где ^ = £ Wj £ у ,
Ш j=l і=і
с начальными условиями:
х. (0) = х°, і = 1,...,5, N(0) = 1,
где х1 - концентрация і-го компонента (мольная доля); Ы=С/С0 - относительное изменение числа молей реакционной среды, где С0 - начальная суммарная концентрация реагентов (кмоль/м3); V) -матрица стехиометрических коэффициентов (і=1,...,5; 7=1,...,9); Wj=aj/C0 - приведенные скорости химических реакций (7=1,..,9) (1/ч).
Сокращенная схема данной реакции получена в работе [3] на основе комбинированного алгоритма редукции схемы реакции в диапазоне времени и температуры протекания реакции. Данный алгоритм основывается на методе анализа графа связей с распространением ошибки [4, 5]. Редуцированная схема реакции имеет вид:
2 Х1 о X 2,
2 X1 о X 3,
2 XI ^ X 4!
X2 оXз,
X2 ^ X4,
Xз ^ X4. (5)
Скорости стадий реакции димеризации а-ме-тилстирола выражаются уравнениями:
ш 1 (С, T) = кС2 - k7C2,
ш 2 (С, T) = к2С1 - к8С3,
ш3(С, T) = kзC12!
®4(С,T) = k4C2 -k9Cз,
ш5(С, T) = к5С2, ffl6(C> T) =
где С=(С!,С2,С3,С4) - вектор концентраций компонентов; к=(к!,...,к9) - вектор кинетических констант скоростей стадий реакции (5).
Таблица 2. Кинетические параметры сокращенной схемыы реакции димеризации а-метилстирола в присутствии катализатора ЫаНУ при температуре 373 К
№ к, (373 К), м3/(кгкат-ч) Е, кДж/моль № к, (373 К), м3/(кгкат-ч) Е, кДж/моль
1 62,788 197,6 6 0,70168 320,9
2 6,037 231,4 7 0,00121 301,3
3 9,055 263,3 8 0,00847 242,4
4 1,092 311,9 9 0,00468 184,0
5 0,0012 573,3
С помощью построенного алгоритма решения обратной задачи химической кинетики рассчитаны кинетические параметры сокращенной схемы реакции а-метилстирола (табл. 2).
Результаты и обсуждение
вычислительного эксперимента
В результате решения обратной кинетической задачи рассчитаны значения энергий активации Ej и констант ^ (/=1,...,9) для редуцированной схемы реакции димеризации а-метилстирола (табл. 2), на основе которых решена прямая кинетическая задача. Относительная разница между расчетными и экспериментальными значениями концентраций веществ составила не более 11 %, что укладывается в погрешность измерений при проведении эксперимента. На рисунке представлена динамика концентраций целевых веществ сокращенной схемы и концентраций этих же веществ в исходной схеме при температуре Т=353 К. Как видно из рисунка, сокращение схемы реакции (4) не изменило общую динамику изменения концентраций веществ во времени. Относительные погрешности векторов концентраций веществ Х1, Х2, Х3, Х4 для сокращенной схемы реакции димеризации а-ме-тилстирола составили: 5(х1')=1,35 %,
ад=1,68 %, <5(х3)=10,24 %, ад=7,93 %.
Отсюда видно, что точность описания динамики концентраций целевых веществ сокращенной схемой реакции димеризации а-метилстирола находится в пределах погрешности количественного анализа. Поэтому редуцированная схема реакции димеризации а-метилстирола является эквива-
лентной схеме реакции (4) меньшей размерности, что позволяет использовать ее при решении других задач, основанных на анализе кинетической модели схемы реакции.
t. ч
Рисунок. Динамика изменения концентраций целевых веществ при Т=353 К (X, - вещества исходной схемы реакции, X’ - вещества сокращенной схемы реакции, i=1,2,3)
Для автоматизации процедуры поиска параметров кинетической модели химической реакции в среде визуального программирования Borland Delphi на языке Object Pascal разработан программный комплекс, который позволяет осуществить построение кинетической модели реакции, решить прямую задачу химической кинетики, рассчитать значения констант скоростей стадий реакции на основе метода Хука-Дживса, вычислить значения энергий активации стадий.
СПИСОК ЛИТЕРАТУРЫ
1. Степашина Е.В., Мустафина С.А. Формирование математической модели каталитических процессов с переменным реакционным объемом на основе теоретико-графового подхода // Известия Томского политехнического университета. - 2012. -Т. 320. - №3. - С. 31-36.
2. Байтимерова А.И. Математическое моделирование и численное исследование каталитических процессов в каскаде реакторов: дис.... канд. физ.-мат. наук. - Уфа, 2009. - 127 с.
3. Степашина Е.В. Алгоритмы редукции кинетических схем сложных химических процессов: дис. ... канд. физ.-мат, наук. - Уфа, 2013. - 160 с.
4. Степашина Е.В., Мустафина С.А. Графовый подход при моделировании каталитических процессов с переменным реакционным объемом // Системы управления и информационные технологии. - 2012. - Т. 47. - № 1. - С. 14-19.
5. Xia G., Michelangeli D.V., Makar P.A. Mechanism reduction for the formation of secondary organic aerosol for integration into a 3-dimensional regional air quality model: a-pinene oxidation system // Atmospheric Chemistry and Physics. - 2009. - V. 9. -P. 4341-4362.
Поступила 19.06.2013 г.
UDC 544.43
IDENTIFICATION OF MATHEMATICAL MODELS OF REDUCED REACTION SCHEME
V.A. Vaytiev, E.V. Stepashina, S.A. Mustafina Sterlitamak Branch of Bashkir State University
The authors have developed the algorithm to solve an inverse problem for the reduced reaction schemes which allows determining kinetic parameters of reaction schemes with smaller dimension preserving dynamics of concentration target substances. The mathematical model of the reduced scheme of a-methylstyrene dimerization was identified. Rate constants of stages and activation energy values were found for this model.
Key words:
Inverse problem, reduction of the kinetic schemes, reaction mechanism.
REFERENCES
1. Stepashina E.V., Mustafina S.A. Bulletin of the Tomsk Polytechnic University, 2012. 320, 3, pp. 31-36.
2. Baytimerova A.I. Matematicheskoe modelirovanie i chislennoe is-sledovanie kataliticheskikh processov v kaskade reaktorov. Diss. ... cand. fiz.-mat. nauk. (Mathematical modeling and numerical investigation of catalytic processes in reactor column): Diss. ... Cand. fys-math. Sci. Ufa, 2009. 127 p.
3. Stepashina E.V. Algoritmy reduktsii kineticheskikh skhem slozhnykh khimicheskikhprotsessov. Diss. ... cand. fiz. and mat. nauk. (Algorithms of reduction of kinetic diagrams in complex chemical processes) Diss. ... Cand. fys-math. Sci. Ufa, 2013. 160 p.
4. Stepashina E.V., Mustafina S.A. Sistemy upravleniya i informat-sionnye tekhnologii, 2012. 47, 1, pp. 14-19.
5. Xia G., Michelangeli D.V., Makar P.A. Atmospheric Chemistry and Physics, 2009. 9, pp. 4341-4362.
УДК 544.452.2
ЗАВИСИМОСТЬ ВЫХОДА НИТРИДА АЛЮМИНИЯ ОТ МАССЫ НАВЕСКИ И ДАВЛЕНИЯ ВОЗДУХА ПРИ ГОРЕНИИ НАНОПОРОШКА АЛЮМИНИЯ
Роот Л.О., Звягинцева Е.С., Ильин А.П.
Томский политехнический университет E-mail: tolbanova@mail.ru
Экспериментально обоснованы методы повышения выхода нитрида алюминия в продуктах сгорания нанопорошка алюминия в воздухе: повышение массы навески нанопорошка алюминия и повышение давления воздуха. Установлено, что с помощью увеличения массы навески сжигаемого нанопорошка алюминия до 15 г или при повышении давления воздуха до 120 кПа увеличи-вется выход нитрида алюминия примерно на 30 %. Экспериментально подтвержден ранее предложенный механизм формирования нитевидных кристаллов из газовой фазы.
Ключевые слова:
Нитрид алюминия, оксид алюминия, нанопорошок алюминия, горение, продукты сгорания, газофазная реакция, высокотемпературная реакция, синтез сжиганием, нитридсодержащие керамические материалы.
Введение
Изучение горения в воздухе нанопорошков (НП) алюминия, бора, кремния, титана, циркония, гафния, ниобия, тантала и хрома, а также смесей их оксидов с НП алюминия показало, что в составе продуктов их сгорания в воздухе содержались в виде самостоятельных кристаллических фаз соответствующие нитриды [1-3]. Среди различного вида материалов интерес представляет нитрид алюминия и нитридсодержащие композиты на его основе, что связано с его уникальными свойствами: нитрид алюминия является хорошим изолятором, устойчив в среде азота до высоких температур, в то же время по теплопроводности сравним с серебром. В настоящее время алюмони-тридная керамика широко используется в качестве
подложек в электронике. Нитрид алюминия не смачивается жидким алюминием, поэтому изделия из него широко применяются в литейном производстве.
Одним из наиболее перспективных методов получения нитридсодержащих композитов является синтез сжиганием порошкообразного алюминия или нанопорошка алюминия в воздухе [4]. Это направление в настоящее время интенсивно развивается, так как синтез протекает в самоподдерживающемся экзотермическом режиме за относительно короткое время. Для масштабного производства нитридсодержащих материалов и технической керамики с помощью синтеза сжиганием в воздухе необходимо разработать методы повышения выхода нитрида алюминия в продуктах сгорания.
