CC BY
16АПВПМ-2019
ПРИМЕНЕНИЕ МАТЕМАТИЧЕСКИХ МЕТОДОВ ДЛЯ РЕШЕНИЯ СИСТЕМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ, ОПИСЫВАЮЩИХ ПРОЦЕСС ОКИСЛЕНИЯ
ИЗОПРОПИЛБЕНЗОЛА
М, К, Вовденко1, С, А. Габитов1,2, И, М, Губайдуллин1,2
1 Уфимский Федеральный исследовательский центр РАН Институт нефтехимии и катализа, 450075, Уфа 2 ФГБОУ ВО Уфимский государственный нефтяной технический университет, 450062, Уфа
УДК 004.942
Б01: 10.24411/9999-016А-2019-100 13
Окисление изоиропилбеизола (кумола) кислородом воздуха является одной из промежуточных стадий получения фенола и ацетона в кумольном методе, самым распространенном методе получения данных веществ на сегодняшний день. Реакция протекает по радикально-цепному механизму, при этом наряду с целевыми продуктами получаются нежелательные побочные компоненты. Различные авторы указывают разные направления развития цепей реакции, механизмы инициирования и степень влияния тех или иных факторов на ход реакции. Ключевые слова: кумол, жесткая задача, метод Гира, метод имитации отжига
Введение
Процесс окисления изоиропилбеизола (ИПБ) кислородом воздуха является промежуточной стадией технологического процесса получения фенола и ацетона в т.н. кумольным способом. На сегодняшний день более 90% производимого в мире фенола и ацетона получают данным способом [1]. В ходе данной стадии происходит химическое превращение изоиропилбеизола в гидроперекись изоиропилбеизола (ГП ИПБ), которая впоследствии распадается на фенол и ацетон на следующей технологической стадии.
1 Модель
Брутто уравнение целевой реакции окисления имеет вид (1).
ДЯ + 02 ^ Я ООН (1)
где И — кумил-радикал (С6Н5СН(СН3)2) [1].
Однако фактически реакция окисления является радикально-цепной, в ходе которой могут проходить одновременно десятки стадий [2,3]. На основании анализа работ, посвящённых изучению реакционного механизма окисления ИПБ, составим реакционную схемы окисления ИПБ, представленную ниже, в таблице 1 [2-4].
Ключевые вещества в таблице: Н11 — изопропилбензол, ТЮН — диметилфенилкарбинол, Н'О — ацето-фенон.
Реакции (1)-(5) относятся к реакциями инициирования, в ходе которых образуются радикалы, которые в последствии зарождают основную реакционную цепь. Существуют ращличные механизмы инициирования. В данной реакционной схеме представлены как реакции инициирования через распад гидроперекиси на радикалы (3)-(6) [1-3], так и инициирование без участия гидроперекиси [2].
Реакции (6)-(7) относятся к основной реакционной цепи, которая приводит к образованию основного количество гидропероксида. Практически во всех исследованиях реакционного механизма данные стадии не меняются [2,3].
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 18-37-00015.
!ЯВ.\ 978-5-901548-42-4
Таблица 1: Реакционная схема окисления ИПБ
№ Тип реакции Реакция
1 Инициирование ИН + 02 ^ТЬ + Н02•
2 Инициирование ИН + Н02 • ^ТЬ + Н2О2
3 Инициирование тюон^ш> + жь
4 Инициирование + Т1Н ^ ти + тюн
5 Инициирование •ОН + ИН ^ И • + Н2О
6 Основная цепь ть + о2 ^ ш2 •
7 Основная цепь ЕО2 • + Т1Н ^ ти + ШОН
8 Побочная ^ И'О + СНз^
9 Побочная 2110 2 • + 2 СНз^ + О2
10 Побочная СНз• + 02 ^ СНз02^
И Побочная СЩ02^ + ИН ^ СНзООН + ТЬ
12 Побочная СН3ООН ^ Н2СО + ОН2
13 Побочная СЩ02^ ЮН + НСОН + О2
14 Побочная 2НСОН + 02 ^2НСООН
15 Побочная 2ТЬ ^ И-И
16 Побочная ТЬ + Ю2» ^ 110011
17 Побочная 21102• ^ М)ОТ1 + О2
18 Побочная 2110 2 • ^ 2Ш> + О2
19 Побочная + ШОН^ШИ + Ю2»
20 Побочная + ин^ а-мэ + що + ТЬ
21 Побочная 2Е00Н^02^ + Н2О
Теоретически основная реакционная цепь может продолжаться бесконечно. На практике, как правило, основная реакционная цепь обрывается с образованием побочных продуктов, или же приводит к возникновению побочной цепи. В ходе окисления может образовываться большое количество различных типов побочных продуктов, и в различных исследованиях представлены различные типы реакций приводящие к их образованию [1-4]. Для нашей реакционной схемы мы выбрали прежде всего реакции, приводящие к образованию ацетофенона (8), диметилфенилкарбинола (13,19), органических кислот (10-14), перекиси кумола (16,17), альфаметилстирола (20), а также распад гидроперекиси изопропилбензола (21) [1-4].
Для идентификации математической модели, описывающей реакцию, на основе известных данных о поведении многокомпонентной смеси во времени [2,3, 5] необходимо решить обратную задачу химической кинетики [6], математически представляющуюся задачей глобальной условной оптимизации. Необходимо минимизировать критерий отклонения (2) расчетных данных концентраций от экспериментальных:
N М
р = ЕЕ К - хЬ12 ^ тт (2)
= 1 3=1
Для вычисления расчетных концентраций необходимо решить прямую задачу химической кинетики, т.е. на основе известных начальных концентраций и заданных констант скоростей элементарных стадий получить кривые зависимости концентрации от времени. Для этого необходимо решить систему обыкновенных дифференциальных уравнений, которая в простейшем случае строится как раз на основе закона действующих масс. На основе предложенной схемы реакции была составлена следующая система ОДУ, которая оказалась нелинейной (3).
=
<И ~ <1X2 =
<и =
<1X3, = М
=
<И ~
йх^ = М ~
йХб = <И ~
<1X7 -
йх% =
йхд = М ~ <1хго
<и
dxll <И
<1X12
<и
<1X13
<и
<1x14 <И
<1X15
<и
<1X16
<и
йХ17
м
<1X18
<и
<1X19
<и
dX20
<и
<1X21
<и
<1X22
<и
—■¡VI — "Ш2 — 'Ш4 — — — и)ц — 1Л20 —•Ю! — те + тд — т14 + юц- + ^18
IV! + -Ш2 + -Ш4 + -ш5 — -ш6 + т7 + -Шц — 2 -Ш15 — + Ы20 + -ш21 — ш2
■Ш2
+ — т1д — 'Ш21 т3 — 1^4 — -Ш8 — Ы13 + -Ш18 — -Ш1д — -Ш20 + т21 ■шз —
ш4 + Ш13 + Ы1д = Ы5 + -Ш20
= те — -№7 — 2тд — т1е — 2и>17 — 2-Ш18 + "Ш1д + т21 = w8 + 2тд = Ы8 + 2тд — тю
= 1Л10 — Ы11 — 1Л13 = тц — 'Ш12 = ™12 = ™12
= -Ш13 — 2-Ш14 = 2ти
= ть
= ™17
= ы20
(3)
где слагаемые в правой части уравнений выражаются согласно (4).
и)1 = = к1Х1Х2
■Ш2 = = к2Х1Х4
т = = кзхе
и)4 = = к4Х7Х1
= = к^Х8Х1
ые = = кеХ2Хз
и)7 = = к7ХцХ1
™8 = = к8Х2
"Шд = = кдхцхц
™10 = кюХ1зХ2
™11 = кцХ14Х1
'Ш12 = к12Х15
и>13 = к1зХ14Х7
и>14 = кЫХ18Х18Х2
те = к15хз Хз
те = к1ехз хц
т7 = к1бХцХ11
и>18 = к1бХцХ11
тд = ктх7 хе
Ы20 = к1еХ7 Х1
и>21 = к1ехе хе
(4)
2 Численный метод
Как и большинство систем ОДУ, описывающих реальные процессы, приведенная выше система является жесткой, т.е. содержит как быстро, так и медленно изменяющиеся компоненты [7]. Таким образом, метод, решающий систему должен «не потерять» медленно изменяющие составляющие и при этом вести интегрирование с шагом, который не затянет вычисления на очень долгое время, то кроме того чревато накоплением ошибки ЭВМ в решении. С подобной задачей по сравнению с явными методами лучше справляются неявные методы. При этом неявные методы требуют больше вычислительных ресурсов, т.к. на каждом шаге решения ОДУ необходимо вместо вычисления простой формулы решать систему нелинейных уравнений. Примерами неявных методов являются: неявный метод Эйлера, метод Розенброка, метод Гира, метод Адамса-Моултона [8,9]. В настоящей работе был выбран метод Гира (5), который является многошаговым (от 2 до 5 шагов), т.е. требует для вычислений в очередной точке времени не только значение в предыдущей точке, а значения в двух предыдущих точках [10]. Соответственно, несколько начальных шагов приходится решать более простыми методами (например, для 1-го шага нужен одношаговый метод и т.д.).
4 12
Уп+2 - 3Уп+1 + 3Уп = 3Ь/(¿„+2,Уп+2) (5)
Выбор метода объясняется рядом причин. Метод Гира обладает лучшей сходимостью по сравнению с ранее указанными, а также не требует расчёта производной, в отличие от метода Розенброка. В построенной математической модели матрица якобиана плохо обусловлена — содержит много нулей, и поэтому методы, не требующие вычисления производной, не сталкиваются с проблемой расхождения решения.
Для решения систем нелинейных уравнений был выбран метод Вегстейна (6), который показал более быструю сходимость по сравнению с методами итераций, Эйткена-Стефенссона и Ньютона [11]. Он является двухшаговым, поэтому первый шаг при решении системы необходимо вычислить методом простых итераций.
/(хк-\)(хк-1 - Хк-2)
Хк = Хк-1----г-Т,-— (6)
I (Хк-1) - I (Хк-2)
Для решения обратной задачи необходимо многократно решить прямую задачу. Ввиду большого количества локальных минимумов, необходимо применять глобальные методы, среди которых для решения был выбран метод имитации отжига [12]. Этот метод основывается на имитации физического процесса, который происходит при кристаллизации вещества, в том числе при отжиге металлов. Предполагается, что атомы уже выстроились в кристаллическую решётку, но ещё допустимы переходы отдельных атомов из одной ячейки в другую. Предполагается, что процесс протекает при постепенно понижающейся температуре. Переход атома из одной ячейки в другую происходит с некоторой вероятностью, причём вероятность уменьшается с понижением температуры. Устойчивая кристаллическая решётка соответствует минимуму энергии атомов, поэтому атом либо переходит в состояние с меньшим уровнем энергии, либо остаётся на месте [6]. Достоинством метода является отсутствие необходимости вычисления производных, т.к. требуется лишь вычисление значения критерия минимизации. Однако, время и точность получения решения зависят от параметров отжига (функций температуры и перехода к новому состоянию).
3 Результаты
Для расчётов математической модели была реализована программа на языке программирования С++, который хорошо зарекомендовал себя для решения задач численными методами [12], несмотря на наличие готовых математических пакетов (Matlab, MathCad, Maple и др.). Перечисленные среды позволяют управлять многими параметрами моделирования, однако не застрахованы от получения некорректного решения и ошибки в ходе расчётов. Стоит отметить, что реализация неявных методов в этих средах как правило базируется на градиентных методах решения систем нелинейных уравнений, которые требуют дифференцирования, что малоэффективно при плохо обусловленном якобиане системы (что часто встречается в задачах химической кинетики).
Расчётами получены следующие смоделированные графики зависимости концентраций ключевых веществ от времени — на рис. 1 представлено изменение концентрации III IB во времени, а на рис. 2 — ГП III IB во времени.
Рис. 1: Изменение концентрации ИПБ во времени
Рис. 2: Изменение концентрации ГИ ИПБ во времени
Список литературы
[1] Закошанский В.М. Фенол и ацетон: анализ технологий, кинетики и механизма основных реакций. СПб.: ХИМИЗДАТ. 2009. 608 с.
[2] Kazuo Hattori. Yuxi Tariaka. Hiroyuki Suzuki. Tsnneo Ikawa arid Hiroshi Kubota. Kinetics of liquid phase oxidation of curiicrie in bubble column// Journal of chemical Engineering of Japan 1970 P.72 78.
[3] Макалец Б.И.. Кириченко Г.С.. Стрыгин Е.И. и др. Кинетическая модель жидкофазного окисления кумола в гидроперекись// Нефтехимия. 1978 Т. 18. №2. С. 250 255.
[4] Bhattacharya. A. Kinetic modeling of liquid phase autoxidation of cnniene / Bhattacharya. A // Chemical Engineering Journal 2008. P. 308 319
[5] Губайдуллин И.М.. Сайфуллина Л.В.. Еникеев M.I«Информационно-аналитическая система обратных задач химической кинетики». Учебное пособие. Изд-е Башкирск. Ун-та. Уфа. 2003. 89 с.
[6] Барский А.Б. Параллельные процессы в вычислительных системах: планирование и организация. М.: Радио и связь. 1990. 255 с.
[7] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. — М.: Мир, 1999. — 685 с.
[8] Калиткин H.H. Численные методы решения жестких систем // Математическое моделирование. — 1995. - Т. 7. - No 5. - С. 8-11.
[9] Розенброк X., Стори С. Вычислительные методы для инженеров-химиков. — М.: Мир, 1968. — 443 с.
[10] Современные численные методы решения обыкновенных дифференциальных уравнений: Пер. с англ. / под ред. Дж. Холла, Дж. Уатта. — М.: Мир, 1979. — 312 с.
[11] Вержбицкий В.М. Численные методы (математический анализ и обыкновенные дифференциальные уравнения). — М.: Высшая школа, 2001. — 382 с.
[12] Каханер Д., Моулер К., Нэш С. Численные методы и программное обеспечение. — М.: Мир, 1998. — 320 с.
Вовденко Михаил Константинович — аспирант Института нефтехимии и катализа УФИЦ РАН;
e-mail: mikhail_vovdenko@rambler.ru;
Габитов Салих Азатович — магистр Уфимского государственного нефтяного технического университета;
e-mail: salihgabitov@yandex.ru;
Губайдуллин Ирек Марсович — д. ф.-м. «., ст. науч. сотр. Института нефтехимии и катализа УФИЦ РАН,
проф. Уфимского государственного нефтяного технического университета;
e-mail: irekmars@mail.ru. Дата поступления — 1 июня 2019 г.
