CC BY
31
УДК 517.977.5
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА КАТАЛИТИЧЕСКОЙ ИЗОМЕРИЗАЦИИ ПЕНТАН-ГЕКСАНОВОЙ ФРАКЦИИ
© А. Г. Фасхутдинов1*, И. В. Ахметов2, И. М. Губайдуллин1'2, А. Е. Мусина2
1 Институт нефтехимии и катализа УФИЦ РАН Россия, Республика Башкортостан, 450075 г. Уфа, пр. Октября, 141.
Тел./факс: +7 (347) 284 27 50.
*Етай: azamat_d1989@mail.ru
2 Уфимский государственный нефтяной технический университет Россия, Республика Башкортостан, 450064 г. Уфа, ул. Космонавтов, 1.
Рассмотрены основные схемы превращений каталитической изомеризации пентан-гексановой фракции. С помощью расчетов показано, что протекание всех реакций, входящих в схему превращений, вероятны с точки зрения термодинамики при условиях проведения процесса. Вероятность протекания реакций оценена по значению энергии Гиббса. Обоснован выбор кинетической модели для наиболее полного описания основных химических процессов исходя из механизма протекания реакций. На основе кинетической модели построена математическая модель, заключающаяся в рассмотрении всех компонентов сырьевой фракции. Заложены основы для последующего моделирования всей химико-технологической системы процесса.
Ключевые слова: изомеризация, пентан-гексановая фракция, бифункциональный катализатор, кинетическая модель, математическая модель.
В связи с переходом отечественной нефтепе-
рерабатывающей промышленности на выпуск автобензинов, соответствующих стандарту «Евро-5» и «Евро-6», актуальной задачей становится снижение содержания ароматических углеводородов и, в особенности, бензола в моторном топливе при одновременном сохранении значения октанового числа. Известно, что основное количество арома-тики в автомобильных бензинах образуется в процессе каталитического риформинга. Каталитическая изомеризация легких парафинов в свою очередь позволяет получать высокооктановый компонент автомобильного бензина с минимальным содержанием ароматики.
Целевым назначением процесса каталитической изомеризации пентан-гексановой фракции является получение высокооктановых компонентов автомобильных бензинов. Данный процесс обладает высокой эффективностью, поскольку в качестве сырья используются такие низкооктановые компоненты нефти, как фракции 62-70 °С, а также рафинаты каталитического риформинга. Процесс осуществляется в среде водорода в присутствии бифункциональных катализаторов [1-2].
Следует отметить, в послевоенные годы в ходе проектирования и строительства большинства установок, в том числе реакторов и аппаратов, глубокому изучению механизма процессов уделялось недостаточное внимание. Интенсивное развитие таких современных информационных технологий, как технология параллельных вычислений, искусственные нейронные сети, 3Б-моделиро-вание, имитационное моделирование и т.п., позволяют разработать детализированные кинетические модели сложных процессов. Затем на основе данных моделей существенно увеличить выход целе-
вого продукта, минимально затрачивая при этом энергетические и материальные ресурсы [3].
При рассмотрении процессов нефтепереработки часто возникают проблемы, связанные с описанием кинетики. Это, в первую очередь, можно связать с большим компонентным составом сырья, а также протеканием огромного количества химических реакций, протекающих на кислотных и металлических центрах бифункциональных катализаторов [4-5]. В этой связи, для упрощения описания кинетики, стадии образования промежуточных веществ, отсутствующих в составе целевых продуктов, обычно не рассматривают.
Объектом исследования является реакторный блок установки каталитической изомеризации пентан-гексановой фракции, состоящий из трех последовательных реакторов. В качестве сырья взяли гидроочищенную бензиновую фракцию с пределами выкипания 62-70 °С, содержащую в основном пентан и гексан. Загрузка сырья составляла 47.6 м3/ч. Расход водородсодержащего газа (ВСГ) - 3400 нм3/ч. Давление в зоне реакции составляло 3.2 МПа. Температуры на входе первого, второго и третьего реакторов составляли 147, 160 и 145 °С соответственно. На выходе из реакторов -160, 170 и 146 °С соответственно. В качестве катализатора процесса применяли бифункциональный катализатор отечественного производства СИ-2. Как было отмечено, реакторный блок объекта исследования состоит из трех реакторов, в которых катализатор распределяли по 9000 кг.
Характеристика сырья и продукта (стабильный изомеризат) приведены в табл. 1 (данные эксплуатации промышленной установки каталитической изомеризации пентан-гексановой фракции).
Таблица 1
Характеристика сырья и продукта
Анализируемый показатель Методика Значение
Сырье Стабильный изомеризат
Плотность при 20 0С, кг/м3 ГОСТ Р 51069 653.6 643.2
Октановое число (исследова- ГОСТ 8226 67 81
тельский метод)
Фракционный состав, 0С
начало кипения 42 36
90% выкипает ГОСТ 2177 62 59
97% выкипает 64 63
конец кипения 65 64
Выход, % 98 98
Таблица 2
Состав сырья и продукта процесса
Компонент
Сырье, % мас.
Изомеризат, % мас.
2.2-ДМС4
2.3-ДМС4
2-МС5
3-МС5 Бензол Изопентан Метилциклопентан Н-гексан Н-пентан Изобутан н-Бутан
Пропан
Этан
Метан
0.36 2.15 19.32 14.73 1.25 10.75 2.7 18.19 30.81
21.13 6.57 19.58 11.23
23.56 0.12 7.27 7.80 1.04 0.87 0.61 0.13 0.09
В табл. 2 представлен состав сырья и продукта процесса каталитической изомеризации пен-тан-гексановой фракции (данные промышленной установки).
Очевидно, что для разработки детализированной кинетической модели достаточно моделирование только каскада реакторов. Состав газопродуктовой смеси после реакторного блока рассчитывается из данных по продуктам после блока стабилизации, имеющейся в наличии.
Начальным этапом математического описания объекта является составление схемы реакций углеводородов в ходе процесса. Следует отметить, от степени детализации химических превращений напрямую зависит точность расчетов и адекватность математической модели реальному процессу. Поэтому схема превращений должна в достаточной степени отражать физико-химическую сущность процесса [6-7].
В данной работе была предпринята попытка более детализированного анализа схемы превращений исходя из механизма.
Известно, что реакции в процессе изомеризации в присутствии бифункциональных катализаторов протекают по карбоний-ионному механизму [8-10]. При этом основу составляет термодинамика процесса. Результатом термодинамического анализа является список реакций, которые вероятно протекают в заданных условиях [11-12]. Вероятность протекания реакций при технологических условиях процесса каталитической изомеризации пентан-гексановой фракции (температура - 150 °С (423 К) и давление -3.2 МПа) была оценена по значению энергии Гиббса Дв (табл. 3).
Из табл. 3 видно, что реакции изомеризации, дегидроциклизации, превращения нафтенов, дегидрирования являются обратимыми. Тогда как реакции гидрокрекинга - необратимыми [13].
Согласно уравнениям реакций, протекающих в процессе каталитической изомеризации пентан-гексановой фракции, запишем кинетические уравнения следующим образом (табл. 4).
Таблица 3
Уравнения реакций процесса каталитической изомеризации пентан-гексановой фракции и значения термодинамических характеристик
Реакция | АН, кДж/моль | AG, кДж/моль
1. н-С5Н12 ^ изо-С5Н12 -8.22 -5.76
2. Н-С6Н14 2-МП -6.98 -3.77
3. н-СбН14 3-МП -4.44 -7.60
4. Н-С6Н14 2,2-ДМБ -18.53 -5.60
5. Н-С6Н14 2,3-ДМБ -10.93 -0.99
6. 2-МП 3-МП -2.54 -3.02
7. 2,2-ДМБ 2,3-ДМБ -7.61 -4.61
8. БЗ + Н2 ЦГ -213.94 -52.65
9. БЗ + Н2 МЦП -197.03 -53.83
10. ЦГ МЦП -112.34 49.10
11. Н-С6Н14 ^ ЦГ+ Н2 8.57 7.82
12. 2-МП МЦП+ Н2 -66.96 -29.59
13. 3-МП МЦП+ Н2 -64.42 -26.58
14. 2,2-ДМБ МЦП+ Н2 -78.51 -31.43
15. 2,3-ДМБ МЦП+ Н2 -70.91 -26.82
16. МЦП+ Н2—3-МП -64.42 -26.58
17. Н-С5Н12 + Н2 — С3Н8+С2Н6 -44.49 -49.96
18. Н-С5Н12+Н2 — Н-С4Н10+СН4 -56.57 -61.66
19. Н-С6Н14 + Н2—2С3Н -37.94 -48.88
20. Н-С6Н14 + Н2—-Н-С4Н10+С2Н6 -45.91 -51.80
21. Н-С6Н14 + Н2—Н-С5Н12+СН4 -56.14 -47.34
22. Н-С6Н14 + Н2 — НЗ0-С5Н12+СН4 -64.35 -66.71
23. 2-МП + Н2 — НЗ0-С4Н10+С2Н6 -47.38 -49.82
24. 3-МП + Н2 — НЗ0-С5Н12+СН4 -59.91 -65.95
25. 2,3-ДМБ +Н2— НЗ0-С4Н10+С2Н6 -43.44 -52.63
26. 2,3-ДМБ +Н2— МЗ0-С5Н12+СН4 -53.43 -65.71
27. изо-С5Н12 + Н2 — изо-С4Н10+ СН -56.82 -57.71
28. 2,2-ДМБ + Н2 — U30-C4H10+C2H6 -35.83 -48.05
29. 2,2-ДМБ + Н2 — U30-C5H12+CH4 -45.82 -61.11
Обозначения: 2-МП - 2-метилпентан, 3-МП - 3-метилпентан, 2,2-ДМБ - 2,2-диметилбутан, 2,3-ДМБ - 2,3-диметилбутан, ЦГ - циклогексан, МЦП - метилциклопентан, БЗ - бензол.
Таблица 4
Кинетические уравнения процесса каталитической изомеризации пентан-гексановой фракции
Реакция_|_Кинетические уравнения
1. н-С5Н12 « изо-С5Н12 ^1=i1x1-i2x2
2. Н-С6Н14 ~ 2-МП W2=k3x3-k4x4
3. Н-С6Н14 ~ 3-МП W3=k5x3-k6x5
4. Н-С6Н14 ~ 2,2-ДМБ W4=k7x3-k8x6
5. Н-С6Н14 ~ 2,3-ДМБ W5=k9x3-k10x7
6. 2-МП ~ 3-МП W6=knX4-k12X5
7. 2,2-ДМБ ~ 2,3-ДМБ W7=k13x6-k14x7
8. БЗ+Н2 ~ ЦГ W8=k15X8X93-k16X10
9. БЗ+Н2 ~ МЦП W9=k17x8x93-k18xn
10. ЦГ ~ МЦП W10=k19x10-k20xn
11. Н-С6Н14 ~ ЦГ+ Н2 Wn=k21x3-k22x10x9
12. 2-МП ~ МЦП+ Н2 W12=k23x4-k23xnx9
13. 3-МП ~ МЦП+ Н2 W13=k25x5-k26xnx9
14. 2,2-ДМБ ~ МЦП+ Н2 W14=k27x6-k28xnx9
15. 2,3-ДМБ ~ МЦП+ Н2 W15=k29x7-k30xnx9
16. МЦП+ Н2 — 3-МП W16=k31xnx9
17. Н-С5Н12+ Н2 — С3Н8+С2Н6 W17=k32x1x9
18. Н-С5Н12+ Н2 — Н-С4Н10+СН4 W18=k33x1x9
19. Н-С6Н14+ Н2 — 2С3Н8 W19=k34x9
20. Н-С^Н.14+ Н2 —* Н-С4Н10+С2Н6 W20=k35x3x9
21. Н-С6Н14+ Н2 — Н-С5Н12+СН4 W21=k3fx3x9
22. Н-С6Н14+ Н2 — изо-С5Н1г+СНл W22=k37xx
23. 2-МП+ Н2 — изо-С4Ню+С2Н6 W23=k38xx
24. 3-МП+ Н2 — изо-СзН12+СН4 W24=k39x5x9
25. 2,3-ДМБ+Н2 — изо-С4Ню+С2Н6 W25=k40x7x9
26. 2,3-ДМБ+Н2 — изо-СзН12+СН4 W26=k41x7x9
27. изо-С5Н12+Н2 — изо-С4Н10+ СН4 W27=k42x2x9
28. 2,2-ДМБ+Н2 — изо-С4Ню+С2Н6 W28=k43x6x9
29. 2,2-ДМБ+Н2 — изо-СзН12+СН4 W29=k44x6x9
Обозначения: х1 - концентрации компонентов в мольных долях, х\ - н-пентан, х2 - изопентан, х3 - н-гексан, х4 - 2-МП, х5 -3-МП, х6 - 2,2-ДМБ, х7 - 2,3-ДМБ, х8 - бензол, х9 - водород, х10 - ЦГ, хп - МЦП, х12 - пропан, х13 - этан, х14 - н-бутан, х15 -метан, х16 - изобутан, к/ - кинетическая константа скорости /-ой реакции, л-моль^-ч"1 (/ = 8, 9, 11-29), ч-1 (/ = 1-7, 10), Ш/ -скорость /'-ой реакции, моль/(л-ч).
Скорости реакций, входящие в кинетическую модель (табл. 4), были записаны согласно закону действующих масс [14].
Математическая модель процесса представляет собой систему нелинейных дифференциальных уравнений [15-16]:
dx1
Л dx2 Л
= -Щ7 -Щ8 + W21
= щ + ж22 + ж23 + ж24 + ж26 - ж27 + ж
= —ж —ж—ж—ж—ж —ж —ж —ж —ж
ж2 ж3 ж4 ж5 ж11 '19 ж 20 ' 21 ж 22
dt
= ж2 —жб —ж12 —ж2з
dt ЛХ5 dt ЛХ6 dt ЛХ1 dt
dt
= —3ж8 —зж, + жп + ж!2 + Жз + ^4 + ^5 —^6 —^7 ^8 ^9 — ж2о —
dt
ж ж ж ж ж ж ж ж ж
ж21 ж22 ж23 ж24 ж25 ж26 ж27 ж28 ж29
—10 = ж —ж + ж
ж 8 ж10 ^ ж11
= ж, + жб —^з + ^6 —ж24
=ж ж ж ж ж
= ж4 ж7 ж14 ж28 ж29
= ж5 + ж7 —^5 —ж25 —ж26
= —ж — ж9
dt
dt
dt
= ж, + жт + Ж2 + ^з + ж,4 + ж,5 —ж16
= ж, + 2жс
dt ЛХ, ,
. = ^7 + ж25 + ж27 + ж2
dt
dt
=ж
= ж18
= ^8 + ж2! + ж22 + ж2з + ж24 + ж26 + ж
dt
16 = ж25 + ж27 + ж28
с начальными условиями при t = 0, хг(0) = хг0, где г' - соответствующий углеводород.
Таким образом, размерность системы уравнений математической модели совпадает с количеством веществ, концентрация которых определяются в расчетах. Для применения полученной модели сначала необходимо провести сравнительную оценку кинетических параметров.
При решении обратной кинетической задачи планируется использовать метод параболического спуска. Для решения системы обыкновенных нелинейных дифференциальных уравнений выбран модифицированный метод Кутты-Мерсона пятого порядка точности [17].
Алгоритм решения системы дифференциальных уравнений = Лх) требует пятикратного вычисления правых частей Лх) [18]:
Р = Ахк)Н, Р=Лхк+ВДк, Рз = А(хк+(¥1+¥2)/6)К Р = ¡(хк+(¥1+з¥2)/8)Н
Хк+1 = Хк+(^!-з^2+4^4)/2
Р = /(Хк+\)Ь
Х*к+1 = Хк+(р1+зр4+р5)
(1) (2)
(3)
(4)
(5)
(6) (7)
где к - шаг интегрирования по пространственной или временной координате, хк - значения переменных в начале интервала интегрирования. В конце интервала интегрирования получаются два решения: хк+х - методом Рунге-Кутты с точностью О(к4) и х к+1 - методом Кутты-Мерсона с точностью О(к5).
Выводы
1. С помощью расчетов показали, что протекание всех реакций, входящих в схему превращений, вероятны с точки зрения термодинамики при условиях проведения процесса.
2. Предложена детализированная математическая модель, заключающаяся в рассмотрении всех компонентов сырьевой фракции с пределами выкипания 62-70 °С (пентан-гексановая фракция).
Заключение
Разработка на основе механизма реакций адекватной кинетической модели является важным этапом в создании математической модели. Поскольку путем решения обратных кинетических задач определяются параметры, которые будут служить основой для математической модели. А также учитывающие физико-химические аспекты процесса. В конечном итоге математическая модель позволит проводить прогнозируемые расчеты, подобрать оптимальные технологические условия с целью повышения ресурсоэффективности процесса каталитической изомеризации пентан-гексановой фракции.
ЛИТЕРАТУРА
1. Ахметов С. А. Технология глубокой переработки нефти и газа. Учебное пособие для вузов. Уфа: Гилем, 2002. С. 672.
2. Koncsag C. I. Study of C5/C6 isomerization on Pt/H-zeolite catalyst in industrial conditions // Ovidius Univ. Annal. Chem. 2011. Vol. 22. №2. P. 102-106.
3. Слинько М. Г. Теоретические основы химической технологии. 2007. Т. 41. №1. С. 16.
4. Stijepovic M. Z., Ostojic A. V., Milenkovic I., Linke P. Development of a Kinetic Model for Catalytic Reforming of Naphtha and Parameter Estimation Using Industrial Plant Data // Energy Fuels. 2009. Vol. 23. №2. P. 979-983.
5. Исмагилов И. З., Матус Е. В., Нефедова Д. В., Кузнецов В. В., Яшник С. А. Керженцев М. А., Исмагилов З. Р. // Кинетика и катализ. 2015. Т. 56. №3. С. 397.
6. Белинская Н. С., Силко Г. Ю., Францина Е. В., Ивашкина Е. Н., Иванчина Э. Д. Разработка формализованной схемы превращений углеводородов и кинетической модели процесса гидродепарафинизации дизельных топлив // Известия Томского политехнического университета. 2013. Т. 322. №3. С. 129-133.
7. Chuzlov V. A., Ivanchina V. A., Dolganov I. M., Molo-tov K. V. Simulation of Light Naphtha Isomerization Process // Procedia Chemistry. 2015. №15. P. 282-287.
8. Barron Y., Maire G., Muller J. M., Gault F. G. The mechanisms of hydrogenolysis and isomerization of hydrocarbons on metals: II. Mechanisms of isomerization of hexanes on platinum catalysts // Journal of Catalysis. 1966. Vol. 5. №3. P. 428-445.
9. Ono Y. A survey of the mechanism in catalytic isomerization of alkanes // Catal. Today. 2003. №81. P. 3-16.
10. Liu H., Meng X., Zhang R., Liu Z., Meng J., Xu C. Reaction performance and disappearance kinetics of n-pentane isomerization catalyzed by chloroaluminate ionic liquid // Catal. Commun. 2010. №12. P. 180-183.
11. Dymond, J. H., Young, K. J., Isdale J. D. // Chem. Ther-modyn. 1979. №11, P. 887.
12. Maginn E. J., Alexis T. B., Doros N. T. Sorption Thermodynamics, Siting, and Conformation of Long n-Alkanes in Silicalite As Predicted by Configurational-Bias Monte Carlo Integration // J. Phys. Chem. 1995. Vol. 99. №7. P. 2057-2079.
13. Triwahyono S., Jalil A. A., Mukti R. R., Musthofa M., Razali N.A.M., Aziz M.A.A. Hydrogen spillover behavior of Zn/HZSM-5 showing catalytically active protonic acid sites in the isomerization of n-pentane // Appl. Catal. A: Gen. 2011. №407. P. 91-99.
14. Logan S. R. Fundamentals of Chemical Kinetics // Addison Wesley Longman. 1996. P. 280.
15. Ахметов И. В., Губайдуллин И. М., Коледина К. Ф., Са-фин Р. Р. Математическое моделирование и оптимизация реакций синтеза ароматических соединений // Электротехнические и информационные комплексы и системы. 2015. Т. 11. №2. С. 53-58.
16. Ахметов И. В., Губайдуллин И. М. Кинетическая модель реакции получения метилового эфира 5-ацетил-2-пирролкарбоновой кислоты // Обозрение прикладной и промышленной математики. 2009. Т. 16. №5. С. 80.
17. Полак Л. С., Гольденберг М. Я., Левицкий А. А. Вычислительные методы в химической кинетике. М.: Наука. 1984. С. 280.
18. Арушанян О. Б., Залеткин С. В. Численное решение обыкновенных дифференциальных уравнений на Фортране. М.: изд-во МГУ. 1990. С. 336.
Поступила в редакцию 30.08.2018 г.
MATHEMATICAL MODELING OF THE PROCESS OF CATALYTIC ISOMERIZATION OF PENTAN-HEXANE FRACTION
© A. G. Faskhutdinov1*, I. V. Akhmetov2, I. M. Gubaydullin1'2, A. E. Musina2
1Institute of Petrochemistry and Catalysis, RAS 141 Oktyabrya Avenue, 450075 Ufa, Republic of Bashkortostan, Russia.
2Ufa State Petroleum Technological University 1 Kosmonavtov Street, 450064 Ufa, Republic of Bashkortostan, Russia.
Phone: +7 (347) 284 27 50.
*Email: azamat_d1989@mail.ru
The authors of the article studied the main schemes of reaction networks of the catalytic isomerization of the pentane-hexane fraction. With the help of calculations, it was shown that all reactions occurring in the transformation scheme under the conditions of the process are possible from the point of view of thermodynamics. The probability of reactions was estimated from the value of Gibbs energy AG. The selection of the kinetic model for the most complete description of the basic chemical processes based on the mechanism of the reactions proceeding was substantiated. The mathematical model was constructed on the basis of kinetic model, taking into consideration all components of the raw material fraction. The foundations for the subsequent modeling of the entire chemical-technological process system were laid.
Keywords: isomerization, pentane-hexane fraction, bifunctional catalyst, kinetic model, mathematical model.
Published in Russian. Do not hesitate to contact us at bulletin_bsu@mail.ru if you need translation of the article.
REFERENCES
1. Akhmetov S. A. Tekhnologiya glubokoi pererabotki nefti i gaza. Uchebnoe posobie dlya vuzov [Technology of deep processing of oil and gas. Textbook for universities]. Ufa: Gilem, 2002. Pp. 672.
2. Koncsag C. I. Ovidius Univ. Annal. Chem. 2011. Vol. 22. No. 2. Pp. 102-106.
3. Slin'ko M. G. Teoreticheskie osnovy khimicheskoi tekhnologii. 2007. Vol. 41. No. 1. Pp. 16.
4. Stijepovic M. Z., Ostojic A. V., Milenkovic I., Linke P. Energy Fuels. 2009. Vol. 23. No. 2. Pp. 979-983.
5. Ismagilov I. Z., Matus E. V., Nefedova D. V., Kuznetsov V. V., Yashnik S. A. Kerzhentsev M. A., Ismagilov Z. R. Kinetika i kataliz. 2015. Vol. 56. No. 3. Pp. 397.
6. Belinskaya N. S., Silko G. Yu., Frantsina E. V., Ivashkina E. N., Ivanchina E. D. Izvestiya Tomskogo politekhnicheskogo universiteta. 2013. Vol. 322. No. 3. Pp. 129-133.
7. Chuzlov V. A., Ivanchina V. A., Dolganov I. M., Molotov K. V. Procedia Chemistry. 2015. No. 15. Pp. 282-287.
8. Barron Y., Maire G., Muller J. M., Gault F. G. Journal of Catalysis. 1966. Vol. 5. No. 3. Pp. 428-445.
9. Ono Y. Catal. Today. 2003. No. 81. Pp. 3-16.
10. Liu H., Meng X., Zhang R., Liu Z., Meng J., Xu C. Catal. Commun. 2010. No. 12. Pp. 180-183.
11. Dymond, J. H., Young, K. J., Isdale J. D. Chem. Thermodyn. 1979. No. 11, Pp. 887.
12. Maginn E. J., Alexis T. B., Doros N. T. J. Phys. Chem. 1995. Vol. 99. No. 7. Pp. 2057-2079.
13. Triwahyono S., Jalil A. A., Mukti R. R., Musthofa M., Razali N.A.M., Aziz M.A.A. Appl. Catal. A: Gen. 2011. No. 407. Pp. 91-99.
14. Logan S. R. Addison Wesley Longman. 1996. Pp. 280.
15. Akhmetov I. V., Gubaidullin I. M., Koledina K. F., Safin R. R. Elektrotekhnicheskie i informatsionnye kompleksy i sistemy. 2015. Vol. 11. No. 2. Pp. 53-58.
16. Akhmetov I. V., Gubaidullin I. M. Obozrenie prikladnoi i promyshlennoi matematiki. 2009. Vol. 16. No. 5. S. 80.
17. Polak L. S., Gol'denberg M. Ya., Levitskii A. A. Vychislitel'nye metody v khimicheskoi kinetike [Computational methods in chemical kinetics]. Moscow: Nauka. 1984. Pp. 280.
18. Arushanyan O. B., Zaletkin S. V. Chislennoe reshenie obyknovennykh differentsial'nykh uravnenii na Fortrane [Numerical solution of ordinary differential equations on Fortran]. Moscow: izd-vo MGU. 1990. Pp. 336.
Received 30.08.2018.
