РАСЧЕТ ПРОЦЕССА ГИДРИРОВАНИЯ α-ПИНЕНА В ТРУБЧАТОМ РЕАКТОРЕ С ПРЯМОТОКОМ ТЕПЛОНОСИТЕЛЯ СРЕДСТВАМИ ПАКЕТА MATHCAD PRIME Текст научной статьи по специальности «Математика»

УДК 66.011, 519.61

А. Н. Иванов, С. А. Мустафина

РАСЧЕТ ПРОЦЕССА ГИДРИРОВАНИЯ а-ПИНЕНА В ТРУБЧАТОМ РЕАКТОРЕ С ПРЯМОТОКОМ ТЕПЛОНОСИТЕЛЯ СРЕДСТВАМИ ПАКЕТА MATHCAD PRIME

Ключевые слова: математическое моделирование, гидрирование, пинен, жесткие системы, численные методы, Mathcad.

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

Key words: mathematical modeling, hydrogenation, pinene, stiff equations, numerical methods, Mathcad.

The article describes the development of a mathematical model for the hydrogenation of а-pinene in a tubular reactor with a direct feed of the refrigerant, as well as its subsequent implementation by the Mathcad Prime package. The mathematical formulas describing in a differential form the change in the system parameters by the relative length of the reactor are given.

Введение

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

осуществлением лабораторных исследований и опытно-промышленных испытаний,

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

Цель и объект исследования

В описываемой работе предлагается составление и расчет математической модели процесса гидрирования а-пинена (2,6,6-

триметилбицикло[3. 1. 1]гепт-2-ена). Данное

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

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

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

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

Процесс гидрирования а-пинена происходит при предварительном нагреве реакционной смеси, а реакция сопровождается выделением тепла. Химизм процесса, а именно целевая реакция гидрирования (1) и побочные реакции изомеризации (2) а-пинена, представлены ниже:

(1)

(2)

При этом следует отметить, что при перегреве реакционной массы наблюдается интенсивное протекание побочных реакций, в связи с чем данный процесс следует проводить в режиме, близком к изотермическому. С этой целью гидрирование ведут в реакторах теплообменного типа -кожухотрубчатых или типа «труба в трубе» [5]. В межтрубное пространство подается вода в качестве теплоносителя для отвода тепла реакции, а внутритрубное пространство заполнено

соответствующим пористым катализатором гидрирования.

Математическая модель процесса гидрирования пинена

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

1) кинетику химических реакций, протекающих как в жидкой, так и в газовой фазе;

2) термодинамику процесса в жидкой и газовой фазах;

3) процессы межфазного массообмена -испарения и конденсации компонентов смеси;

4) процесс теплопереноса между жидкой и газовой фазой, а также теплопередачи через стенку теплоносителю.

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

Математическое описание материального баланса процесса основано на уравнениях неразрывности потоков жидкости и газа, а также на уравнениях химической кинетики [6]. Значения энергий активации и предэкспоненциальных множителей в уравнении Аррениуса для вышеприведенных реакций были получены путём решения обратной кинетической задачи на основе имеющихся справочных данных. Направление подачи воды для охлаждения совпадает с направлением подачи сырья, то есть аппарат представляет собой теплообменник прямоточного типа.

Материальный баланс включает следующие зависимости параметров процесса в форме дифференциальных уравнений:

- уравнение неразрывности потока жидкости:

dL :_ Fl • Ур , где dL - скорость изменения мольного расхода жидкости, кмоль/ч; Fl - объемная скорость испарения, кмоль/ч-м3; Ур - реакционный объем (межтрубное пространство за вычетом объема, занимаемого катализатором, м3.

- скорости изменения концентраций компонентов в жидкой фазе:

dx(/):_ F(/)-р 1 •Я • Ур, Р1

где dx(/) - скорость изменения мольной доли /-го компонента в жидкой фазе; F(/) - скорость изменения количества /-го компонента в жидкой фазе в единице объема аппарата, кмоль/ч-м3; р -

вектор параметров процесса, а именно: р;.+1 -мольная доля /-го компонента в жидкой фазе, р1 -мольный расход жидкости, кмоль/ч.

Компоненты жидкой фазы: 1 - а-пинен, 2 - цис-пинан, 3 - изомеры пинена. Компоненты газовой фазы: 1 - а-пинен, 2 - цис-пинан, 3 - изомеры пинена, 4 - водород.

- уравнение неразрывности потока газов:

dG := ФО Ур, где dG - скорость изменения мольного расхода газа, кмоль/ч; ФG - объемная скорость изменения количества газовой фазы, кмоль/ч-м3.

- скорости изменения концентраций компонентов (за исключением водорода) в газовой фазе:

. Ф(/) - р+5 • ФО

dy(/ ):=^——--Ур ,

р5

где dy(/) - скорость изменения мольной доли /-го компонента в газовой фазе; Ф(/) - скорость изменения количества /-го компонента в газовой фазе в единице объема аппарата, кмоль/ч-м3; р/+5 -мольная доля /-го компонента в газовой фазе, р5 -мольный расход газа, кмоль/ч.

- скорости изменения концентрации водорода в газовой фазе:

^ _-(1 -ф)• м0)-ф• ш(1)-р9 • ФО Ур

' р5 '

где dyh - скорость изменения мольной доли водорода в газовой фазе; ф - мольная доля газовой фазы; М (/) - скорость /-ой реакции в жидкой фазе; ш (/) - скорость /-ой реакции в газовой фазе; р9 -

мольная доля водорода в газовой фазе,

Приложения к материальному балансу:

- расчет мольной доли газовой фазы -фактической и для равновесной системы жидкость-газ:

Рб

Ф :=

Pi + Р5

Ф2 := 0.5 - начальное приближение для равновесной системы.

a/fa (/) := pi • pi+/ + p5 • p5+/ , Pi + P5

где a/fa (i) - мольная доля /-го компонента в жидкой и газовой фазах.

a/fa (/ )•( k/sp (/)-1) pg

РФ2 (ф2) :=£■

1 + Ф2 •(k/sp (/)-1) (Р1 + Р5 )• Ф2

где

k/sp (/) - константа равновесия процесса

испарения для /-го компонента.

Значение равновесной мольной доли газовой фазы вычисляется численным способом c помощью функции root. Предварительно производится проверка на наличие действительного корня в интервале [0, 1]:

ф2 := if (Fф2 (TOL) • Fф2 (1) < 0)

|| root (Fф2 (ф2), ф2,TOL,1) else

if Fф2 (TOL)> 0

if Fф2(TOL)> Fф2 (1)

I11

else

II 0

else

if Fф2(TOL)> Fф2 (1)

II 0

else

III

где TOL - точность расчета в среде Mathcad. - скорость испарения Vnap (кмоль/ч):

Vnap := ф2 •(Р1 + p5)-p5,

Fl :=

-Vnap Vp ;

Y (/) :=if Vnap > 0

|| Р/+1 else

II Р/+5

где Y (/) - мольная доля /-го компонента в жидкой фазе при испарении или в газовой фазе при

конденсации.

и :=

-1 1 0 -1 0 1

матрица

стехиометрических коэффициентов для основной и побочной реакций.

F(/") := (1 -ФMfc •w(j))

j=1

Vnap

"Vp"

Y (/),

0G := -(1 - Ф) • W(1) - ф • w(1) +

Vnap

"VT

Ф(/) := ф• £(Uj,, • w(j)) + VnpP• Y(/).

j=1

Vp

Тепловой баланс: зменение темпе

1000 • Vp • ¿(-Q (j) • ((1 - Ф) • W (j) + Ф • w( j)))

- изменение температуры реакционной смеси:

2

" .........) +

dT :=-

j=1

С^ • р + CpG • р5

-\/лар• £^(у)•У(у)) + ах• Зх, • (р„ -Р10)• \/р

+-^-,

С^ • Р, + CpG • р5

где dT - скорость изменения температуры реакционной смеси, К; Q (у) - тепловой эффект для у-ой реакции, кал/моль; CpL - молярная теплоемкость жидкой фазы, кал/моль-К; СрС -молярная теплоемкость газовой фазы, кал/моль-К; ах - коэффициент теплопередачи, ккал/(м2-ч-К); р10 - температура реакционной смеси, К; р,, -температура теплоносителя, К; Зх, - внутренняя поверхность теплообмена на единицу объема, м-1; dHV (у) - теплота испарения у-го компонента, кал/моль.

- изменение температуры теплоносителя:

,= ах • Зх2 •(Р10 - Р11 )^Р ' Сх • Сх '

где dTx - скорость изменения температуры теплоносителя, К; Зх2 - внешняя поверхность теплообмена на единицу объема, м-1; Сх -теплоемкость хладагента (воды), ккал/(кг-К); Сх -массовый расход хладагента, кг/ч. Приложения к тепловому балансу:

- площадь поверхности теплообмена на единицу объема аппарата:

Зх1 := 3, 1 ^р'

Зх2:= ^, 2 ^р'

где З1 и З2 - внутренняя и внешняя поверхности

теплообмена, м2.

- тепловой эффект химических реакций:

Q(j) :=-ЕК • dH0,):

где dH0( - стандартная теплота образования /-го компонента, ккал/моль.

Вычисление кинетических констант

Обозначим константы скорости реакции в жидкой фазе:

к? - константа скорости реакции гидрирования в жидкой фазе;

к2 - константа скорости реакции изомеризации в жидкой фазе;

/=1

i=1

к3 - константа скорости реакции, обратной реакции изомеризации в жидкой фазе;

и аналогично для процессов, протекающих в газовой фазе:

к4 - константа скорости реакции гидрирования в газовой фазе;

к5 - константа скорости реакции изомеризации в газовой фазе;

к6 - константа скорости реакции, обратной реакции изомеризации в газовой фазе.

В таблице 1 приведены значения кинетических параметров на основании известных экспериментальных данных [7].

Таблица 1 - Данные кинетики процесса гидрирования а-пинена

Константа скорости реакции Температура, °С Энергия активации E, ккал/моль

70 100 130

в жидкой фазе

ki, ч-1 •атм-1 0,0792 0,1792 0,3588 6,92

k2, ч-1 0,0350 0,1088 1,2454 9,63

кз, ч-1 0,5488 1,2454 2,5058 6,96

в газовой фазе

к4, ч-1 -1 3 моль •м 436,83 1580,6 4723,1 10,9

к5, ч-1 32,102 207,07 1012,1 15,8

кб, ч-1 64,966 449,81 2335,1 16,4

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

к, := А • е ,

где А, - предэкспоненциальный множитель для ,-ой реакции; Е, - энергия активации ,-ой реакции,

Дж/моль; Т - абсолютная температура, К; R = 8,314 - универсальная газовая постоянная, Дж/моль-К.

Переведя значения температуры из шкалы Цельсия в шкалу Фаренгейта, и значения энергии активации в Дж/моль, стало возможным вычислить значение предэкспоненциального множителя для каждой из реакций с целью получения непрерывной зависимости между скоростью реакции и температурой смеси (табл. 2).

Таблица 2 - Численные значения кинетических параметров процесса гидрирования а-пинена

Константа скорости реакции Предэкспоненциальный множитель Энергия активации E, Дж/моль

ki 2,037-103 28972,656

k2 4,799-104 40318,884

кз 1,496-104 29140,128

k4 3,867-109 45636,120

к5 3,77240й 66151,440

кб 1,8411012 68663,520

Выбор метода решения системы дифференциальных уравнений

Вышеприведённая математическая модель процесса гидрирования а-пинена в конечном счете представляет собой систему из 11 дифференциальных уравнений - 9 уравнений материального и 2 уравнений теплового балансов. Для вычисления профилей изменения параметров модели по длине аппарата был начат поиск эффективного численного метода её решения.

Как видно из таблиц 1-2, кинетические параметры реакций, протекающих в жидкой и газовой фазах, отличаются на порядки, что является признаком жесткой системы дифференциальных уравнений [8]. Это доказывает и тот факт, что при использовании стандартных решателей - rkfixed (метода Рунге-Кутты 4-го порядка с фиксированным шагом), rkadapt (метода Рунге-Кутты 4-го порядка с адаптационным размером шага), а также решателей Adams (метода Адамса) и Bulstoer (метода Булишера-Штера) не принесло результатов -происходит резкое возрастание погрешности уже в начале вычислений.

Для решения жестких систем дифференциальных уравнений в пакете Mathcad Prime 4.0 предусмотрены несколько решателей:

1) Radau - метод Radaus для жестких систем;

2) BDF (backward differentiation formulas) - метод дифференцирования назад;

3) Stiffb - метод Булишера-Штера для жестких систем;

4) Stiffr - метод Розенброка.

Однако для 2-х последних методов якобиан системы является обязательным аргументом [9], тогда как для методов Radau и BDF он необязателен, хотя и повышает точность численных вычислений. Вычисление якобиана системы дифференциальных уравнений по данной математической модели в символьном виде не представляется возможным по двум причинам. Во-первых, из-за присутствия вычислений равновесной мольной доли газовой фазы, решение которой в символьном виде не найдено. Во-вторых, из-за использования фазовых ограничений на процесс и, как следствие, присутствия условных конструкций в математической модели.

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

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

Core i-3-2330M 2,20 ГГц, 8 ГБ ОЗУ, Windows 8.1, PTC MathCAD Prime 4.0 64-bit выпуск F000. Результаты анализа представлены на рисунке 1. Для наглядности шкала числа интервалов дана в логарифмическом масштабе.

Рис. 1 - Время численного решения системы с помощью методов Radau и BDF

Как видно из графика, наиболее эффективным по времени решателем при одинаковой точности является метод Radau, который впервые появился в Mathcad 2001i. При это наблюдается высокая устойчивость алгоритма даже при большом числе разбиений.

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

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

представлены на рисунках 2-4.

Как видно из рисунка 2, на протяжении порядка 45% длины аппарата происходит снижение мольного расхода газовой фазы без изменения мольного расхода жидкости (за исключения небольшой просадки в самом начале):

Рис. 2 - Профили мольного расхода жидкой и газовой фаз по относительной длине реактора

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

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

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

изомеры водород

Рис. 3 - Профили молярных концентраций компонентов по относительной длине реактора

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

Температурные профили процесса представлены на рисунке 4.

реакционная смесь

Рис. 4 - Профили температур по относительной длине реактора

Благодаря оптимальным параметрам теплосъема и подачи теплоносителя [10, 11] удается поддерживать процесс в состоянии, близком к изотермическому, что наиболее благоприятно для каталитических процессов.

Полученные результаты показывают удовлетворительное согласование с аналитическим решением.

Заключение

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

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

Работа выполнена при финансовой поддержке РФФИ в рамках научного проекта № 17-47-020068 и проекта №13.5143.2017/БЧ, выполняемого вузом в рамках государственного задания Минобрнауки РФ.

Литература

1. А.Н. Иванов, С.А. Мустафина, Е.А. Шулаева, Н.С. Шулаев, Объектно-ориентированное моделирование химико-технологических систем с целью повышения безопасности производства, Математическое моделирование процессов и систем, V Всероссийская научно-практическая конференция, приуроченная к 110-летию со дня рождения академика А.Н. Тихонова (Стерлитамак, Россия, 17-19 ноября, 2016), Стерлитамак, 2016, часть 1, С. 137-142.

2. С.А. Мустафина, Р.С. Давлетшин, Некоторые алгоритмы поиска оптимального управления каталитического процесса, Журнал Средневолжского математического общества, 8, 2, 152-158 (2006).

3. В.И. Быков, В.М. Журавлев, Моделирование и оптимизация химико-технологических процессов. ИПЦ КГТУ, Красноярск, 2002. 298 с.

4. Р.С. Давлетшин, С.А. Мустафина, А.В. Бадаев, С.И. Спивак, О моделировании процесса гидрирования а-пинена, Катализ в промышленности, 6, 34-40 (2005).

5. С.А. Мустафина, Р.С. Давлетшин, А.В. Бадаев, С.И. Спивак, Выбор типа реактора для проведения каталитического процесса гидрирования а-пинена, Обозрение прикладной и промышленной математики, 12, 2, 446-447 (2005).

6. Р.С. Давлетшин, С.А. Мустафина, С.И. Спивак, Исследование процесса гидрирования а-пинена в реакторе с неподвижным слоем катализатора, Вестник Башкирского университета, 10, 2, 15-18 (2005).

7. С.А. Мустафина, А.В. Балаев, Р.С. Давлетшин, С.И. Спивак, У.М. Джемилев, Моделирование процесса газожидкостного гидрирования а-пинена в трубчатых реакторах, Доклады Академии наук, 406, 5, 647-650 (2006).

8. Д. Каханер, К. Моулер, С. Нэш, Численные методы и программное обеспечение. Мир, Москва, 1998. 575 с.

9. М.В. Тихонова, И.М. Губайдуллин, С.И. Спивак, Численное решение прямой кинетической задачи методами Розенброка и Мишельсена для жестких систем дифференциальных уравнений, Журнал Средневолжского математического общества, 12, 2, 26-33 (2010).

10. S.A. Mustafina, Yu.A. Valieva, R.S. Davletshin, A.V. Balaev, S.I. Spivak, Optimization of catalytic processes and reactors, Kinetics and Catalysis, 46, 5, 705-711 (2005).

11. С.А. Мустафина, Р.С. Давлетшин, С.И. Спивак, Математическое моделирование и оптимизация процесса гидрирования а-пинена, Обозрение прикладной и промышленной математики, 11, 2, 376 (2004).

© А. Н. Иванов - магистрант кафедры математического моделирования Стерлитамакского филиала ФГБОУ ВО «Башкирский государственный университет», sanekclubstr@mail.ru; С. А. Мустафина - д.ф.-м.н., профессор, зав. кафедрой математического моделирования Стерлитамакского филиала ФГБОУ ВО «Башкирский государственный университет», mustafina_sa@mail.ru.

© A. N. Ivanov - graduate student of the department of mathematical modeling, Sterlitamak Branch of Bashkir State University, sanekclubstr@mail.ru; S. A. Mustafina - doctor of Science, Professor, head of chair of mathematical modeling, Sterlitamak Branch of Bashkir State University, mustafina_sa@mail.ru.