CC BY
154
УДК 51.74
Профессор В.К. Битюков, профессор С.Г. Тихомиров,
(Воронеж. гос. ун-т. инж. технол.) кафедра информационных и управляющих систем. тел. (473) 255-38-75
доцент Д.В. Арапов,
(Воронеж. гос. ун-т. инж. технол.)
кафедра информационных технологий, моделирования и управления. тел. (473) 255-25-50 E-mail: SavvinSS@yandex.ru
аспирант С.С. Саввин
(Воронеж. гос. ун-т. инж. технол.) кафедра информационных и управляющих систем. тел. (473) 255-38-75
Professor V.K. Bityukov, professor S.G. Tikhomirov,
(Voronezh state university of engineering technologies)
Department of information and control systems department. phone (473) 255-38-75
associate professor D.V. Arapov,
(Voronezh state university of engineering technologies)
Department of information technologies, modeling and control department, phone (473) 255-25-50 E-mail: SavvinSS@yandex.ru
graduate S.S. Savvin
(Voronezh state university of engineering technologies)
Department of information and control systems department. phone (473) 255-38-75
Моделирование процесса пиролиза прямогонного бензина с использованием генетического алгоритма
Modeling of naphtha pyrolysis with using genetic algoritm
Реферат. При эксплуатации промышленных печей пиролиза основной задачей является выбор оптимального режима термического разложения исходного сырья в зависимости от выхода целевых продуктов в условиях технологических ограничений на процесс. Для решения данной проблемы для действующего реактора, в качестве которого в работе рассматривается промышленная крупнотоннажная печь SRT-VI, разработана математическая модель процесса пиролиза, использующая кинетическую схему, состоящую из первичной реакции разложения сырья и вторичных элементарных реакций взаимодействия рассматриваемых компонентов смеси, уравнения теплового баланса и гидродинамики потока в змеевике. Сырьем для установки выбранного типа служит прямогонный бензин. Выходными параметрами модели являются мольные расходы товарных углеводородов. Реактор описывается уравнением идеального вытеснения в статическом режиме функционирования. Принято, что все реакции имеют Аррениусовскую зависимость от температуры. Энергии активации химических процессов оценивались с помощью уравнения Поляни-Семёнова, а идентификация предэкспоненциальных множителей проводилась с помощью генетического алгоритма (ГА). Данная задача требует одновременного решения системы дифференциальных уравнений, описывающих процесс пиролиза и поиска большого количества неизвестных параметров, в связи с чем предложено модифицировать ГА. Оптимальная схема включает операторы кодировки по Грею, турнирной селекции, с рангом турнира больше 4, кроссовера с частичным случайным выбором аллей, мутации с высокой вероятностью ее свершения и эллитизма с конкурентным глобальным состязанием. С использованием предложенного подхода осуществляется параметрическая идентификация модели процесса. Анализ результатов моделирования с данными действующего реактора показал его пригодность для использования в целях управления процессом пиролиза.
Summary. In operation of industrial pyrolysis furnaces, the main task is the selection of the optimal mode of thermal decomposition of the feedstock, depending on the yield of the desired products under conditions of technological limitations on the process. To solve this problem for an operating reactor, this paper considers the SRT-VI Large-Capacity industrial Furnace , the mathematical model of the pyrolysis process was constructed, using a kinetic scheme which consists of primary reaction of decomposition of raw materials and secondary elementary reactions of interaction of the considered mixture components, the heat balance equation and hydrodynamics of flow in the coil. The raw material for the selected installation type is naphtha (straight-run petrol). Output parameters of the model are the molar costs of marketable hydrocarbons. The reactor is described by the equation of ideal displacement in the static mode of operation. It is assumed that all reactions have a temperature dependence that follows the Arrhenius law. The activation energies of chemical processes were estimated using the Polanyi-Semenov equation and identification of pre-exponential factors was carried out using a genetic algorithm (GA). This task requires solving simultaneous system of differential equations describing the pyrolysis process and a search for a large number of unknown parameters, and therefore it is proposed to modify the GA. Optimal scheme includes Gray encoding arithmetic operators, tournament selection, with tournament ranking more than 4, crossover with partial random choice of alleys, mutations with a high probability of occurring and elitism with competitive global competition. Using the proposed approach, the parametric identification of model process is accomplished. The analysis of the simulation results with the data of operating reactor showed its suitability for use in order to control the pyrolysis process.
Ключевые слова: математическая модель, пиролиз, генетический алгоритм, тепловой баланс.
Keywords: mathematical model, genetic algorithm, pyrolysis, heat balance.
© Битюков В.К., Тихомиров С.Г., Арапов Д.В., Саввин С.С., 2015
ВестникВГУИТ, №3, 2015_
Пиролиз углеводородного сырья - важный источник для производства олефинов и ароматических соединений, являющихся основой нефтехимической промышленности. Задачи оптимизации процесса высокотемпературного пиролиза углеводородов, представленные в работе, направлены на достижение максимального выхода целевых продуктов пиролиза и связаны с определением оптимальной совокупности режимных параметров печи.
Разработанная математическая модель узла пиролиза учитывает взаимосвязь протекающих в реакторе физических процессов - гидродинамических, тепловых, массопередачи, с уравнениями химической кинетики. Принимается, что в качестве сырья для пиролизного реактора используется прямогонный бензин, все компоненты, участвующие в реакциях, находятся в газовой фазе, а змеевик пиролизной печи является реактором идеального вытеснения, в связи с высокой турбулентностью потока реагентов.
Для описания химизма процесса предложена кинетическая схема, состоящая из первичной реакции разложения углеводородного сырья, рассмотренной в [1], и вторичных элементарных химических реакций:
1. Naphtha = 0.5Н2 + 0.76СН4 + 1.16С2Н4 +
0.13С2Н6 + 0.38С3Н6 + 0.09С3Н8 + 0.008С4Н1О
0.245С4Н8 + 0.113С4Н6 + 0.08С+
2. С2Н4 о С2Н2 + Н2
3. С2Н6 о С2Н4 + Н2
4. С3Н8 о С3Н6 + Н2
5. С4Н8 о С4Н6 + Н2
6. С4Н10 о С4Н8 + Н2
7. С^Ню ^ 2С2Н4 + Н2
8. С3Н6 о С2Н2 + СН4
9. 2С2Н6 ^ С3Н8 + СН4
10. С3Н8 ^ С2Н4 + СН4
11. С4Н8 ^ С2Н4 + С2Н4
12. С4Н10 ^ СзН6 + СН4
13. C4HW ^ С2^4 + С2^6
14. 2С3Н6 ^ 3С2Н4
15. С3Н6 + Н2 о С2Н4 + СН4
16. С4Н8 + Н2 о С3Н6 + СН4
17. С2Н2 + С2Н4 ^ С4НЬ
18. С2Нв + С2Н4 о С3Н6 + СН4
19. С3Н8 + С2Н4 ^ С3Н6 + С2Н6
20. С3Н6 + С2Н6 о С4Н8 + СН4
21. С2Н4 + С4Н6 ^В + 2Н2
22. С3Н6 + С4Н6 ^Т + 2Н2
23. С4Н8 + С4Н6 ^ЕВ + 2Н2
24. С4Н6 + С4Н6 ^ ST + 2Н2
25. Т + Н2^В + СН4
26. ЕВ + Н2^Т + СН4
27. С4Н8 + СН4 ^ С2Н4 + С3Н8
28. С4Н6 + СН4 ^ С3Н6 + С2Н4
29. СН4 + СН4 ^ С2Н4 + 2Н2
где В - бензол, Т - толуол, ЕВ - этилбензол, ST - стирол, С+- углеводороды в составе которых больше 4 атомов углерода.
Система включает в себя 29 прямых элементарных молекулярных реакций и 10 обратных, протекающих между 15 продуктами. В работе принято, что константы скоростей химических реакций имеют Аррениусовскую зависимость от температуры:
Ei(T)
K(T)i = U • exp
(1)
где ¿ = 1. .39, К(Т)1 - константа скорости ьой реакции схемы (с-1, м3/(моль*с)), - пред-экспоненциальный множитель >ой реакции схемы (с-1, м3/(моль*с)), R - газовая постоянная (Дж/(моль*К), Е1 - энергия активации >ой реакции схемы (Дж/моль), Ь - текущая точка по длине змеевика, Т - температура пирогаза (К).
С целью уменьшения размерности задачи в работе проводится оценка энергий активации по правилу Поляни-Семёнова [2], описывающему соотношения между энтальпией элементарной реакции и энергией активации: Для экзотермических реакций:
Е1(Т)=А-0.25^И1(Т),
Для эндотермических реакций:
Е1(Т)=А + 0.75^И1(Т),
где И - изменение энтальпии >ой элементарной химической реакции (Дж/моль), А = 48 (кДж/моль) для реакций отрыва атома и замещения, А = 42 (Дж/моль) для реакций присоединения.
Для расчета тепловых эффектов реакций применяется уравнение Кирхгофа:
т
И1(Т),= АИ1(ТН) + ^АСринд.^йТ,
где СрИнд. - индивидуальные теплоемкости компонентов, входящих в >ую реакцию (Дж/(моль*К)), знак А обозначает разность между продуктами и реагентами >ой реакции с учетом стехиометрических коэффициентов, Тк - температура смеси в текущей точке разбиения змеевика (К), Тн - температура смеси в предшествующей точке (К).
СрИНДк(Т) = Ак + ВкТ(Ь) + СкТ2(Ь) + ОкТ3(Ь),
где Ак, Вк, Ск, - константы в уравнении идеально-газовой теплоемкости для к-ого компонента схемы.
На основе кинетической схемы составлен материальный баланс, описываемый системой обыкновенных дифференциальных уравнений первого порядка:
ль
<и
р(0
Я • Т(Ь)
лй„
где к = 1. .15; i = 1. .39, - мольный расход компонентов пирогаза (моль/с), - стехио-метрический коэффициент молекулярной реакции, Р - давление в текущей точке змеевика (Па), йвн - внутренний диаметр змеевика (м).
В систему входят 39 неизвестных параметра - по числу предэкспоненциальных множителей в уравнении Аррениуса (1).
Тепловой баланс процесса описывается дифференциальным уравнением следующего вида:
йТ ф • а(Т) -л • й
<И (С + Ср) • Ср(Т)
15
С
(Тст(Т) — Т(Ь))
(С + Ср)-Ср(Т) {-
^(НкСо • хл>тк),
(2)
'Р) -ГЧ-, к = 1
где ф - коэффициент неравномерности обогрева, а(Т) - коэффициент теплоотдачи от стенки змеевика к движущемуся потоку (Дж/(с*К*м2), d - наружный диаметр змеевика (м), Ср(Т) - теплоемкость реакционной смеси (Дж/(кг*К)), G и Gp - расход сырья и пара соответственно (кг/с), х_у1 - мольные доли влажных компонентов пирогаза, Тст(Т) - температура стенки змеевика (К).
X 1(1) = т°'С' т<(Ь)к
~ к вр • тй(Ь)к +т0 • С
где т0 и тк - молекулярные массы воды и кого компонента соответственно, тй - мольные доли компонентов пирогазовой смеси.
тй(1)к = у 15^, )
Для расчета уравнения (2) оценка термодинамических параметров осуществлена с использованием следующих зависимостей:
1
а(Т) = —-т-Т,
ии(тух(т) лст хк где Ыи(Т) - критерий Нуссельта, А - теплопроводность пирогазовой смеси (Дж/(м*с*К)), Зст - толщина стенки (м), <5к - толщина отложений кокса (м), Аст - теплопроводность материала трубы (Дж/(м*с*К)), Ак - теплопроводность кокса (Дж/(м*с*К)).
Ыи(Т) = 0.0214 • Яе(Т)0 8 • Рг(Т)04,
где Яе (Т) - критерий Рейнольдса, Рг( Т) - критерий Прандтля.
р(Т) • йвн • р
Яе(Т) =
КТ) • д
где р(Т) - скорость потока пирогаза (м/с), р - плотность пирогаза (кг/м3), //(Т) - динамическая вязкость пирогазовой смеси (кг*с/м2), д - ускорение свободного падения (м/с2).
= д • КТ • ср(Т
к ) А(Т) '
Скорость потока пирогаза рассчитывается по формуле [3]:
р(Т) =
УаРоТ(Ь)(т0С + Ср ^тт^к) 5Т0Р(Ь~)т0 Ук5=1 ткт<1(Ь)к ' где Уа - молярный объем газов при нормальных условиях (л/моль), Р0 и То - нормальные условия (Па и К соответственно), 5 - площадь сечения трубы (м2).
Теплоемкость пирогазовой смеси описывается зависимостью:
1
Ср(Т) =
Ср 1,1=1 ткт^;)к + тоС
15
mkmd(L)k (ср11Н^(Т))
0 к=1
+
15
+т,
>с(!
к=1
md(L)k тк
(СРиндк(Т))
где СринДо, - индивидуальная теплоемкость водяного пара (Дж/(моль*К)). СРинд0(Т) = Ао + ВоТ(Ь) + СоТ2(Ь) + ОоТ3(Ь),
где .0, В0, С0, Б0 - константы в уравнении идеально-газовой теплоемкости для водяного пара.
Вязкость индивидуального компонента пирогазовой смеси рассчитывается по уравнению Чэмпена-Энскога с использованием потенциала Леннарда-Джонса [4]:
/индк(Т) = 26.69 •
Vтк • Т(Ц)
п2 а к
а,
(3)
где а - радиус твердой сферы (м), А - интеграл
столкновений.
а
2,3551 — 0,0874 • ш
Ф13 ,
1с
где Рс, Тс - критические давление и температура (Па, К), ш - фактор ацентричности.
А =
1.16145 0.52487
+ ■
+
•7^0.14874 0.7732ТрО1
2.16178
(4)
е2Л3787^Тро1-где Тр01 - потенциальная температура.
2
4
ВестникВГУИТ, №3, 2015
г П1)
где £ - параметр потенциала Леннарда-Джонса £ = (0,7915 + 0,1693 • ш) • Тс,
Для расчёта вязкости полярной молекулы воды (3) интеграл столкновений (4) пересчиты-вается по уравнению Брокау:
_ 0.2 *52
^Н20 = ^ + ^ ,
1Ро*Н20
где 8 - параметр полярности молекулы.
Теплопроводность индивидуального компонента пирогазовой смеси рассчитывается по методу Мисика и Тодоса:
АИНя.,(П = 4,45 • 10-6 • Тг
Ср,
Р
2/3
индк с к
'■индк V
" т1/6к-ш1/2к
(5)
где Тг - приведенная температура.
т Т(1) 1т"к= т ' 1ск
Зависимость (5) справедлива для метана, нафтенов и ароматических углеводородов, во всех остальных случаях применяется выражение:
ЛИНДк = 10-6 • (14,52 • ТГк - 5,14)2/3
СРишдк 'Рс/ к
т1/6к-т1/2к'
Вязкость всей смеси с учетом пара определяется по уравнению Васильевой [4]: 16
' И-ИНД1
где у - мольные доли компонентов смеси и пара, А- параметр уравнения Вильке:
т.;
4,1 =
^ + • (т1)1/4)1/2
№инд; тЬ
(8-(1+т)1/2)
К К тг
Для определения теплопроводности смеси пирогаза с учетом пара проводятся аналогичные вычисления.
Для вычисления температура стенки в зависимости (2) используется уравнение:
ТСТ(Т)4+--
' Тст(Т) Тт
ТО) = о,
где (;ст - параметр, учитывающий зачернен-ность поверхности трубы, д - коэффициент, учитывающий интенсивность теплообмена между источником нагрева и участком трубы, Тт - температура в топке печи (К).
д =
а
ъ
1 горелки
2
*<6=1ргоРелки 2 2=1 а2 + (—--
\ 1зм 6.
где ^гОрелки - расход топлива в горелках (кг/с), 1зм - длина змеевика (м), а - настроечный параметр.
Дифференциальное уравнение описывающее изменение давления по длине змеевика
имеет вид: & - — ----
аь
Iв(т) ^ (с + сраг)т(р
^вн ' Тстанд ^ ^вн 1зм' Мсмеси
где Рстанд и Тстанд - нормальные давление и температура соответственно (Па, К), 9(Т) - коэффициент гидравлического трения потока о стенки змеевика, 2 - коэффициент местного сопротивления калача змеевика, Мсмеси - средняя мольная масса смеси.
0.25
А
68
Дифференциальные уравнения кинетики, записанные с учетом теплового баланса и гидродинамики в совокупности, представляют полное математическое описание процесса пиролиза в змеевике печи.
Оценку предэкспоненциальных множителей и из уравнения (1) предложено проводить с помощью поискового алгоритма, так как данная задача имеет множество локальных экстремумов. На данный момент одним из предпочтительных способов многоцелевой оптимизации являются генетические алгоритмы, основанные на механизмах, аналогичных естественному отбору в природе. Поиск оптимального решения опирается на гипотезу селекции.
Генетический алгоритм [5], оперирует не искомыми переменными, а закодированными отображениями - хромосомами. В качестве функции цели используется модульный критерий отклонения расчетных значений от экспериментальных:
15
ъ
,к=1
(хк Хехрк)
• 100 ) ^ тт
Для корректной работы алгоритма с большим числом неизвестных параметров предложено его модифицировать. В работе применяется кодировка по Грею, что гарантирует соответствие соседних порядков хромосом с ближайшими декодируемыми точками пространства. Используется турнирная селекция, при которой все особи текущей популяции случайным
образом разбиваются на подгруппы с последующим выбором в каждой из них особи с наилучшей приспособленностью, при этом повышение ранга турнира приводит к уменьшению числа итераций всего алгоритма, но снижению скорости вычислений. Применяется кроссовер с частичным случайным выбором аллелей. Данная схема, основанная на следующем правиле: совпадающие аллели родительских хромосом, сохраняются у потомков, а не совпадающие выбираются случайным образом. Этот прием позволяет увеличить поисковую способность генетического алгоритма. Задается высокая вероятность мутации хромосом, что обеспечивает широкие возможности перебора решений. Новая популяция создается на основе принципа эли-тизма. В данной работе выбрана его разновидность - конкурентный подход с глобальным состязанием, при которой все родительские особи "состязаются" со всеми потомками и победители, чья приспособленность выше, независимо от возраста, переходят в следующее поколение.
В качестве экспериментальных данных используются результаты промышленной эксплуатации пиролизной печи SRT-VI. Идентификация модели проводилась по выборке из 14 наборов выходов процесса пиролиза, полученных при различных технологических режимах. Часть из них представлена в таблице 1.
Т а б л и ц а 1 Данные промышленной эксплуатации печи SRT-VI
Продукт (% мас.) № эксперимента
1 2 3 4 5
Н2 1.28 1.39 1.41 1.37 1.45
СН4 13.65 13.44 13.07 14.14 14.62
С2Н4 29.4 28.97 27.44 29.55 30.46
С2Н6 4.87 4.87 4.95 5.01 4.82
СэНб 18.2 18.42 20.14 17.68 17.37
С3Н8 0.47 0.58 0.59 0.59 0.48
С4Н10 0.23 0.23 0.24 0.2 0.19
С4Н8 4.43 4.4 5.3 3.96 3.81
С4Н6 7 7.53 78.01 7.26 6.75
С2Н2 0.58 0.58 0.71 0.49 0.48
В 13.3 12.75 12.01 13.83 13.57
т 3.73 3.71 3.42 3.33 3.47
ЕВ 1.4 1.51 1.41 1.2 1.25
S 1.52 1.62 1.3 1.4 1.54
Т (К) 1113 1108 1098 1118 1118
Р S 0.6 0.6 0.6 0.55 0.6
F S 1540 1540 1540 1540 1540
где P_S - отношение пар/сырье, Т - температура на выходе из змеевика (К), F_S - расход прямогонного бензина (кг/ч).
Разработанная математическая модель процесса пиролиза позволяет получить профили концентраций продуктов пиролиза по длине змеевика (рисунок 1).
10
20
L (т)
Рисунок 1. Профиль концентраций продуктов пиролиза по длине змеевика.
Профиль изменения температуры пиро-газовой смеси по длине реактора представлен на рисунке 2.
1200
£ Н
1100-
100С
900
0
10
20
Ь (т)
Рисунок 2. Профиль изменения температуры пирогаза.
0
Вестпик<ВТУИТ, №3, 2015
Профиль изменения давления реакционной смеси по длине змеевика отображен на рисунке 3.
220
215
рц 210 ^ 205
200 195
0 10 20 L (m)
Рисунок 3. Профиль изменения давления пирогаза.
Среднее расхождение полученных значений от экспериментальных данных не превышает 10 %, что соизмеримо с погрешностью хромотографии. В таблице 2 приведен анализ погрешностей целевых продуктов процесса:
ЛИТЕРАТУРА
1 Haghighi S.S., Rahimpour M.R., Raeissi S. et al. Investigation of ethylene production in naphtha thermal cracking plant in presence of steam and carbondioxide // Chemical Engineering Journal. 2013. V. 228. P. 1158-1167.
2 Долганова И.О., Долганов И.М., Ивашкина Е.Н. и др. Развитие подхода к моделированию процесса получения этилбензола // Вестник науки Сибири. 2012. Т. 2. №1. С.35-44.
3 Битюков В.К., Арапов Д.В., Саввин С.С. Кинетическая модель пиролиза бензина в крупнотоннажной печи // Сборник трудов ММТТ-27. 2014. Т.8. С. 182-185.
4 Разносчиков В.В., Демская И.А. Математическая модель расчета теплофизических свойств синтетического жидкого топлива [Электронный ресурс] // Электронный журнал «Труды МАИ». 2012. № 50. Режим доступа: http://www.mai.ru/science/ trudy/published.php (04 мая 2015 г.)
5 Keyvanloo K., Sedighi M., Towfighi. J. Genetic algorithm model development for prediction of main products in thermal cracking of naphtha: Comparison with kinetic modeling // Chemical Engineering Journal. 2012. V. 209. P. 255-262.
Т а б л и ц а 2 Сравнение результатов моделирования с экспериментальными данными по целевым продуктам пиролиза
Продукт Экспериментальные данные (% масс) Рассчетные данные (% масс.) Погрешность
T= 1118 (K), P S=0.6
С2Н4 25.629 26.14 2.03
C3H6 15.337 14.13 7.84
C4H6 6.301 6.64 5.51
T=1103 (K), P S=0.55
С2Н4 24.9 26.31 5.66
C3H6 16.2 15.67 3.27
C4H6 6.1 6.53 7.04
В работе проведено моделирование процесса пиролиза прямогонного бензина с учетом теплового баланса и гидродинамики. Реализована идентификация модели с помощью модифицированного генетического алгоритма на основе данных с действующего производства. Анализ результатов моделирования показал ее пригодность для дальнейшего использования в целях управления технологическим процессом и научных исследованиях.
REFERENCES
1 Haghighi S. S., Rahimpour M. R., Raeissi S. et al. Investigation of ethylene production in naphtha thermal cracking plant in presence of steam and carbondioxide. Chemical Engineering Journal, 2013, vol. 228, pp. 1158-1167.
2 Dolganova I.O., Dolganov I.M., Ivashkina E.N. et al. The development approaches to modeling of production ethylbenzene. Vestnik nauki Sibiri. [Siberian Journal of Science], 2012, vol. 2, no. 1, pp. 35-44. (In Russ.).
3 Bityukov V.K., Arapov D.V., Savvin S.S. Kinetic model of naphtha pyrolysis in large-capacity furnaces. Sbornik trudov ММТТ-27 [Proceedings of MMTT-27]. 2014, vol. 8, pp. 182-185. (In Russ.).
4 Raznoschikov V.V., Demskaya I.A. Mathematical model for calculating the thermal properties of synthetic liquid fuels. Trudyi MAI. [Proceedings of MAI], 2012, no. 50. Available at: http://www.mai.ru/science/ trudy/ (Accessed 04 May 2015). (In Russ.).
5 Keyvanloo K., Sedighi M., Towfighi. J. Genetic algorithm model development for prediction of main products in thermal cracking of naphtha: Comparison with kinetic modeling. Chemical Engineering Journal, 2012, vol. 209, рр. 255-262.
