CC BY
22
ГОРЮНОВ М. В., ПЕСКОВА Е. Е.
РАЗРАБОТКА ПРОГРАММНОГО МОДУЛЯ ДЛЯ РЕШЕНИЯ УРАВНЕНИЙ ХИМИЧЕСКОЙ КИНЕТИКИ1
Аннотация. Работа посвящена разработке программного модуля для решения уравнений химической кинетики с использованием специализированной явной схемы второго порядка. Были решены системы уравнений, описывающие механизм пиролиза этана. Показано, что реализованный алгоритм дает заметный выигрыш по времени выполнения.
Ключевые слова: химическая кинетика, жесткие системы дифференциальных уравнений, пиролиз этана.
GORYUNOV M. V., PESKOVA E. E.
DEVELOPMENT OF SOFTWARE MODULE FOR SOLVING CHEMICAL KINETICS EQUATIONS
Abstract. The study deals with the development of a specific software module for solving chemical kinetics equations using a specialized explicit second-order scheme. The equations systems describing the ethane pyrolysis mechanism have been solved. It is shown that the implemented algorithm gives a noticeable gain on the execution time.
Keyword: chemical kinetics, hard systems of differential equations, ethane pyrolysis.
Введение. Химия - экспериментальная наука, в которой перед практическим использованием результатов, исследования веществ и реакций должны проверяться на опыте с последующими рекомендациями. На сегодняшний день актуальным подходом является математическое моделирование химических процессов. Оно зарекомендовало себя как эффективный и безопасный метод исследования, не требующий больших материальных затрат. Основные причины этого - технический прогресс, который определяется, прежде всего, фантастическими достижениями компьютерных технологий и успехами развития теоретических представлений о строении веществ.
Моделирование химических реакций в газах имеет огромную область практических приложений, таких как работа ТЭЦ, бензиновых и дизельных моторов, газофазные химические реакторы.
При численном решении задач химической кинетики возникают трудности, связанные с жесткостью решаемой системы дифференциальных уравнений. Жесткими называют задачи, решение которых содержит компоненты с резко различными характерными
1 Работа выполнена при поддержке Минобрнауки РФ (№ 1.6958.2017/8.9) и Российского фонда фундаментальных исследований (проект № 18-31-00102).
1
временами изменения [1]. При решении жестких систем ОДУ используются многошаговые методы Гира и схемы Розенброка. Однако данные схемы обладают высокой трудоемкостью.
Для численного решения задач химической кинетики в работе [2] предложен метод, основанный на специфической форме правых частей, которыми обладают дифференциальные уравнения, построенные на основе кинетических схем. Это явный алгоритм, имеющий малую трудоемкость и существенно превосходящий известные методы по простоте и надежности.
Численный метод. Предположим, что реагирующая смесь состоит из М компонент. Химический процесс можно описать совокупностью N элементарных стадий:
м м
1=1 1=1
Здесь А[ - химические элементы, у[п, - стехиометрические коэффициенты компонента А^ в стадии п.
Система уравнений химической кинетики записывается следующим образом:
N
йС; V"1
аГ = £1Мп-1'1п)пп, 1 = 1.....м, (1)
п=1
где С1 - концентрация компонента А^, №п - скорость элементарной стадии реакции:
м м
= К П(С1У[п - к-п , 1=1 1=1
кп, к-п - константы скоростей прямой и обратной стадии реакции, которые определяются из выражений Аррениуса [3]:
кп=Ап • е( кг).
Здесь Ап - предэкспоненциальный множитель, Еп - энергия активации. Согласно [2], систему (1) можно представить в следующем виде:
йС[
= -С№(С) + С = (С1,С2,...,СМ), 1 = 1.....М.
Здесь с1 > 0, щ(с) > 0, > 0.
Решение данной системы находится простыми итерациями [2]: ci + t^i(cs)(1 + TVi(cs)/2) _s (с + с) л0
С?+1 = -2-> °S = -п-> °0 = С-
1 + TVi(cs) + (rVi(cs)) /2 2
Здесь Ci - решение в исходный момент времени, q - решение в новый момент времени.
При численном решении системы выполняем две итерации, последующие итерации не повышают порядок точности и ухудшают надежность схемы [2].
Программная реализация. Рассмотрим элементарную стадию реакции следующего
вида:
3а + Ъ = с + 2d
Занумеруем компоненты: а — 1, Ъ -2, с — 3, d — 4. Для каждой реакции указываются номера компонент a, b,c,d в том порядке, в котором они записаны в реакции. В случае отсутствия некоторой компоненты записываем 0. Также будем полагать, что 0 не могут быть равны одновременно а и Ъ, с и d.
Следующий этап - это запись коэффициентов, стоящих перед каждой компонентой. Коэффициент перед компонентой а - а, Ъ - Ъ, с —~с, d — d.
Таким образом, каждая реакция определяется набором параметров а, Ь, с, d, a, b,~c, d, а также предэкспоненциальным множителем А и энергией активации Е. Эти 10 величин определяют структуру элементарной стадии реакции:
S = (а, Ь, с, d, а, Ь, с, d, А, Е) (2)
Схема реакции представляется в виде таблицы, число строк в которой равно числу элементарных стадий. Добавление новых стадий сводится к добавлению строк этой таблицы.
Реальные химические задачи описываются большим количеством элементарных стадий с сотнями реагирующих компонент. В данной работе реализован алгоритм, который определяет <pi и автоматически. Этот алгоритм, реализованный в функциях solve_phi и solve_psi, сводится к цепочке последовательных проверок, которые находятся внутри двух циклов: внешнего и внутреннего. Внешний 0...count_substance - позволяет выбрать конкретное вещество. Внутренний 0...count_reaction - позволяет пройти по всем реакциям и найти совпадения, которые впоследствии будут обработаны исходя из условий.
Если <pi - это ^ для i-того вещества; ^j - это ^ для i-того вещества; kj - константа скорости для у-той стадии реакции; са - концентрация вещества, с номером а; сь -концентрация вещества, с номером Ь.
Для phi совпадения ищем среди а и Ь.
Если исходные вещества реакции одинаковые и оба совпадают с /-тым веществом
г ^
(а = Ь), то = + (а + Ь) • kj • ct a+b 1
Если а равно i-тому веществу и Ъ = 0, то tyi = + а • kj • саа-1, если Ъ любое другое вещество, то = + а • kj • саа-1съ.
Если b равно i-тому веществу и а = 0, то tyi = tyi + b • kj • сьъ-1, если а любое другое вещество, то + b • kj • cbb-1ca.
Для psi ищем совпадения в с и d.
Если с равно i-тому веществу, возникает три случая. Если а = 0, то ^j = + с • kj • сьь; если b = 0, то = + с • kj • саа; если a и b разные вещества, то ^j = + с • kj • саа •
cbb.
Если й равно /-тому веществу, возникают также три случая. Если а = 0, то = +
й • к] • сьь; если Ь = 0, то ^ = ^ + й • к^ • саа; если а и Ь разные вещества, тогда ^ = ^ + ^ • к] • са.а • сьь•
Вычислительный эксперимент. Описанный алгоритм был реализован на языке программирования С++. Для проверки работы программного модуля была решена система уравнений, описывающая процесс пиролиза этана. В качестве кинетической схемы принята брутто-реакция пиролиза этана, состоящая из двух элементарных стадий [4; 5].
Таблица 1
Схема и кинетические параметры механизма брутто-реакции пиролиза этана
№ Стадия АиУс V (моль •с) /моль
1 ^2^6 ^ СгН4 + Н2 1.0QE + 16 2.5Е + 5
2 2С2Н6 ^ С2Н4 + 2СН4 3.16Е + 16 2.7Е + 5
Для записи данных о каждой элементарной стадии реакции пронумеруем компоненты смеси: В1 = [С2Н6], В2 = [С2Н4], В3 = [Н2], В4 = [СН4].
Соответственно, имеем следующую запись для двух стадий реакции: (1,0,2,3,1,0,1,1,1.08 • 1016,2.5 • 105) (1,0,2,4,2,0,1,2,3.16 • 1016,2.7 • 105) Задачу решаем со следующими начальными условиями:
Т = 800К,с1(0) = 1,с2(0) = 0,с3(0) = 0, с4(0) = 0
На рисунке 1 представлены результаты работы программы в момент установления. Полученные результаты хорошо согласуются с результатами, полученными другими
авторами [6; 7]. В таблице 2 представлены значения концентраций в установившемся течении.
Таблица 2
Сравнение результатов работы программы с полученными ранее
Концентрации Результаты текущей работы Результаты, полученные аналитически в [6]
0 0
с? 0.939841 0.94
Сз 0.878950 0.88
с4 0.121782 0.12
Этан Этилен Водород Метан
200 400 600 800 1000 1200 1400
0
Рис. 1. Изменение концентраций веществ при температуре 800 К, брутто-схема.
Для более точного описания химических превращений в вычислительной практике, как правило, используются схемы с большим количеством элементарных стадий реакции. Покажем, что разработанный программный код дает выигрыш по времени вычислений в сравнении с явным методом Эйлера и методом Розенброка второго порядка, реализованного в системе МЛТЬЛВ.
Для построения математической модели пиролиза этана принята 15-ти стадийная схема реакции, включающая 12 компонент смеси, представленная в таблице 3 [8].
Таблица 3
Схема и кинетические параметры радикального механизма пиролиза этана
№ Стадия \%А-1,1/с или л/ / (моль • с) р кДж/ Ьь /моль
1 С?Нб ^ СН3' + СН3' 16.0 360.0
2 СНз' + С?Нб ^ СН4 + С?Н5' 10.0 50.0
3 С2Н5' ^ С2Н4 + н- 13.5 170.0
4 н• + С?Нб + С?Н5' 11.0 40.0
5 Н' + С2Н4 ^ С?Н5' 10.4 8.4
6 СНз' + С2Н4 ^п- С3Н7' 10.9 33.0
7 п - С3Н7' ^ СНз' + С2Н4 13.9 137.0
8 С2Н5' + С2Н5' ^ С2Н4 + С?Нб 10.0 8.4
9 п - С3Н7' + С?Н4 ^ С?Н5' + СзНб 7.4 27.6
10 СНз' + С?Щ ^ СН4 + С?Нз' 8.6 35.0
11 СНз' + С?Нз' ^ СН4 + С?Н? 9.95 3.2
12 С?Нз + Н' ^ с?н? + н? 10.0 0.0
13 С2Н4 ^ 'с?Н4 15.8 253.0
14 'С?Н'4 + С?Нб ^ СНз' +п- СзЩ' 14.7 216.0
15 'с?Н4 ^ С2Н4 5.38 0.0
Введем следующие обозначения: В1 = [С2Н6], В? = [С2Н4], Вз = [Н2], В4 = [СН4], В5 = [СН3'1 Вб = [С?Н5'], В7 = [Н'], В8 = [п- С3Н7Ч В9 = [С3НбЪ Вю = [С?Н3'1 В11 = [С?Н?1 В12 = [-С?Н4].
Формируем записи элементарных стадий реакции согласно таблице 3 (табл. 4).
Таблица 4
Запись стадий реакции для дальнейших расчетов
№ Запись стадий реакции
1 (1,0,5,0,1,0,2,0,16,360)
2 (5,1,4,6,1,1,1,1,10,10,50)
3 (6,0,2,7,1,0,1,1,13,170)
4 (7,1,3,6,1,1,1,1,11,40)
5 (7,2,6,0,1,1,1,0,10.4,8.4)
6 (5,2,8,0,1,1,1,0,10.9,33)
7 (8,0,5,2,1,0,1,1,13.9,137)
8 (6,6,2,1,1,1,1,1,10,8.4)
9 (8,2,6,9,1,1,1,1,7.4,27.6)
10 (5,2,4,10,1,1,1,1,8.6,35)
11 (5,10,4,11,1,1,1,1,9.35,3.2)
12 (10,7,11,3,1,1,1,1,15,10,0)
13 (2,0,12,0,1,0,1,0,15.8,253)
14 (12,1,5,8,1,1,1,1,14.7,216)
15 (12,0,2,0,1,0,1,0,5.38,0)
Задачу решаем со следующими начальными условиями:
Г = 900 К,с1(0) = 1, с2(0) = 0, ..., с12(0) = 0 Были проведены расчеты с использованием трех вышеупомянутых методов, получены хорошо согласованные друг с другом результаты. На рис. 2 представлен результат работы разработанного программного модуля.
Рис. 2. Изменение концентраций веществ при температуре 800 К.
В таблице 5 представлено время выполнения каждого метода, расчеты велись до времени t = 1 сек.
Таблица 5
Время расчета методов
Разработанный код на основе метода [2] Разработанный код на основе явного метода Эйлера Метод Розенброка второго порядка в системе МЛТЬЛВ
Время 3 мс 29109 мс 138 мс
Таким образом, разработанный программный код имеет преимущество по времени выполнения по отношению к явному методу Эйлера и широко известному методу Розенброка второго порядка для решения жестких ОДУ, реализованному в системе МЛТЬЛВ.
Заключение. С помощью разработанного программного модуля можно проводить моделирование задач химической кинетики, включающих любое количество элементарных стадий. В дальнейшем планируется включение данного модуля в разрабатываемый пакет программ [9; 10] для расчета динамики многокомпонентного ламинарного газового потока с учетом химических реакций.
ЛИТЕРАТУРА
1. Галанин М. П., Ходжаева С. Р. Методы решения жестких обыкновенных дифференциальных уравнений. Результаты тестовых расчетов // Препринт ИПМ им. М. В. Келдыша. - 2013. - № 98. - 29 с.
2. Белов А. А., Калиткин Н. Н., Кузьмина Л. В. Моделирование химической кинетики в газах // Математическое моделирование. - 2016. - Т. 28, № 8. - С. 46-64.
3. Штиллер В. Уравнение Аррениуса и неравновесная кинетика. - М.: Мир, 2000. -176 с.
4. Мухина Т. Н., Барабанов Н. Л., Бабаш С. Е. и др. Пиролиз углеводородного сырья. -М.: Химия, 1987. - 240 с.
5. Губайдуллин И. М., Пескова Е. Е., Язовцева О. С. Математическая модель динамики многокомпонентного газа на примере брутто-реакции пиролиза этана [Электронный ресурс] // Огарев-online. - 2016. - № 20. - Режим доступа: http://journal.mrsu.ru/arts/matematicheskaya-model-dinamiki-mnogokomponentnogo-gaza-na-primere-brutto-reakcii-piroliza-etana.
6. Язовцева О. С. Локальная покомпонентная асимптотическая эквивалентность и ее применение к исследованию устойчивости по части переменных [Электронный ресурс] // Огарев-online. - 2017. - № 13. - Режим доступа: http://journal.mrsu.ru/arts/lokalnaya-pokomponentnaya-asimptoticheskaya-ekvivalentnost-i-ee-primenenie-k-issledovaniyu-
ustoj chivosti-po-chasti -peremennyx.
7. Назаров В. И., Пескова Е. Е., Язовцева О. С. Численное моделирование жестких систем с использованием (4,2)-метода [Электронный ресурс] // Огарев-online. - 2017. -№ 13. - Режим доступа: http://joumal.mrsu.ru/arts/chislennoe-modelirovanie-zhestkix-sistem-s-ispolzovaniem-42-metoda.
8. Nurislamova L. F, Stoyanovskaya O. P., Stadnichenko O. A., Gubaidullin I. M., Snytnikov V. N., Novichkova A. V. Few-Step Kinetic Model of Gaseous Autocatalytic Ethane Pyrolysis and Its Evaluation by Means of Uncertainty and Sensitivity Analysis // Chemical Product and Process Modeling. - 2014. - 9 (2). - pp. 143-154.
9. Жалнин Р. В., Пескова Е. Е., Стадниченко О. А., Тишкин В. Ф. Математическое моделирование динамики многокомпонентного газа с использованием WENO схем на примере пиролиза этана // Журнал Средневолжского математического общества. - 2016. - Т. 18, № 3. - С. 98-106. 10. Жалнин Р. В., Пескова Е. Е., Стадниченко О. А., Тишкин В. Ф. Моделирование течения многокомпонентного реагирующего газа с использованием алгоритмов высокого порядка точности // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. - 2017. - Т. 27, № 4. - С. 608-617.
