CC BY
57
УДК 519.622
Л.В. Кнауб, А.Е. Новиков, Е.А. Новиков
МОДЕЛИРОВАНИЕ КИНЕТИКИ ХИМИЧЕСКИХ РЕАКЦИЙ МЕТОДОМ ВТОРОГО ПОРЯДКА ДЛЯ НЕЯВНЫХ СИСТЕМ*
Построен двухстадийный L-устойчивый метод второго порядка для решения неявных систем. Метод отличается от классических схем приближенным нахождением производной решения. Приведены результаты численного моделирования кинетики химической реакции.
Ключевые слова: дифференциально-алгебраическая задача, метод типа Розенброка, химическая кинетика.
L.V. Knaub, A.E. Novikov, E.A. Novikov
MODELING OF THE CHEMICAL REACTION KINETICS BY THE SECOND ORDER METHOD
FOR IMPLICIT SYSTEMS
The two-stage L-stable method of the second order for solving the implicit systems is developed. The method differs from the classical schemes by the approximate finding of the solution derivative. The numerical modeling results of the chemical reaction kinetics are presented.
Key words: differential-algebraic problem, the method of Rosenbrock type, chemical kinetics.
Введение. При моделировании кинетики некоторых химических реакций возникает проблема решения задачи Коши для жестких систем, не разрешенных относительно производной вида [1]
^(х,х',г) = 0, х(г0) = х0, г0 < г < гк, (1)
где x и F - вещественные Л/-мерные вектор-функции; I - независимая переменная. Современные численные методы обычно предполагают задание явной зависимости производной [1-2], то есть
х = /(г, х), х(/о) = Хо, г, < г < гк. (2)
Приведение задачи (1) к виду (2) требует дополнительных затрат на шаг интегрирования, которые связаны с обращением матрицы Fy. Разрешенная задача (2) жесткая, что приводит к необходимости применения ¿-устойчивых методов, которые нуждаются в обращении матрицы Якоби. Наиболее известные алгоритмы решения задачи (1) основаны на многошаговых численных формулах [3]. Однако многие практические задачи описываются жесткими непрерывно-дискретными системами, которые в современной литературе называют гибридными задачами [4]. Для них характерным является непрерывное описание режимов, смена которых происходит дискретно. В такой ситуации многошаговые методы могут быть малоэффективны. При решении жестких задач широкое распространение получили методы типа Розенброка [5]. Здесь построен двухстадийный ¿-устойчивый метод второго порядка точности решения неявных задач. Метод отли-
* Работа поддержана грантом РФФИ (проект №14-01 -00047).
чается от классических схем типа Розенброка приближенным нахождением производной решения. Приведены результаты численного моделирования кинетики химической реакции.
Численная схема. Используя обозначение х'=у, задачу (1) можно переписать в виде системы дифференциально-алгебраических уравнений
у = у, ^ (х, у, г) — 0, ^ < г < гк
(3)
с начальными условиями х(?о)=хо и о)=уо. Дополнительное условие у(?о)=уо можно вычислить, например, решая задачу Р(хо,уДо)=0 методом установления [6]. Ниже будем предполагать существование и единственность решения задачи (3) и невырожденность матрицы 0п=гпу+аЛРт, где а - числовой коэффициент, к - шаг интегрирования, а и - частные производные функции Р(хуД), вычисленные в точке П Тогда т-стадийный метод типа Розенброка для численного решения задачи (3) имеет вид [6]
X 1 — X
п+х п
+ 1 1 ' Уп+1 Уп + 1 1
1—х
Вкх — к¥ •
п I пу
г 1-х N
у •у
1—х
ак1 К
пг
-к¥
Уп + 1 в ук
V
+1 в,к? - Уп+2 в.,к. у'+к •! в
V у—х у—х х—х
1-х
X
кУ —
ак
1 -х
к? - к •
Уп + 1 ву к
V у—х
(4)
где а, р' в - числовые коэффициенты. При применении численных формул (4) для решения задачи (2) получим классические методы типа Розенброка [5]. Так как уо вычисляется приближенно, будем предполагать выполненным соотношение ||Р(хоуоДо)||<СЛР, где ||-|| - некоторая норма в С - константа, не зависящая от величины шага интегрирования, р - порядок точности метода.
Для решения (3) рассмотрим формулу типа Розенброка вида
Хп+х — Хп + РхК + Р2к2Х' Уп+х — Уп + РхМ + Р2к2 ' ОкХх — к • Уп - акГп1 - ^ (Хп, Уп, гп)], кУ — а- х (к? - кУп) / к, (5)
— к¥щ .(Уп + в2х кхУ )-ак2 К - кКп (х + в2х кХ - Уп + вх КЛ + вгх к ),
к( — а
- х
к2Х - к (Уп + в2 х кхУ )
/ к.
Разложим стадии кх и ку в ряды Тейлора по степеням Л и подставим в первую формулу (5). Сравнивая полученное представление приближенного решения с рядом Тейлора для точного решения, получим условия второго порядка точности схемы (5)
Рх + Р2 — х ' в2хР2 — 0-5 - а ■
Разлагая =Г(х^+1 ,уп+1,Гп+1) в ряд Тейлора с учетом (6), получим
= я"2 (а2 - 2а + 0.5)^ + О(И2)
Условие F„+1=O(h2) приводит к уравнению а2-2а+0.5=0. Данное соотношение обеспечивает Ь -устойчивость схемы (5) при решении разрешенной задачи. Уравнение а2-2а+0,5=0 имеет два вещественных корня. Из численных экспериментов следует, что корень а=1-0^2 предпочтительнее, потому что он приводит к более эффективным расчетам. Требование ¿-устойчивости промежуточной численной формулы приводит к соотношению в1=а. В результате получим коэффициенты ¿-устойчивого метода (5) второго порядка
Р1
= в21 = а = 1 - 0.5л/2, р2 = 1 - а = 0.5л/2
Контроль точности. Неравенство для контроля точности построим по аналогии [6-8], то есть на каждом шаге следует проверять неравенство
II ¿2х - К 11<е-
где к 1х и к 2х заданы при описании численной схемы (5), е-требуемая точность расчетов, где ||-|| - некоторая норма в №. В расчетах использовалась норма
II? п\\ = тах
11 п|1 1<г<^
х | + г
п
где г - положительное число. Если |уп'|<г, то по Ай компоненте решения контролируется абсолютная ошибка ге; в противном случае контролируется относительная ошибка е. В отличие от методов типа Розенброка решения разрешенных систем производная в численной формуле (5) вычисляется приближенно, и поэтому при выборе шага дополнительно проверяется неравенство ||Оп ^п||<е.
Тестовый пример. Задача взята из [9], схема реакции имеет вид
1 к1
2БЬВ+—СО,^БЬВТ+Н-О- ZLA+FLB< к2'К >FLBT+ZHU-2 2
к3 1 ¿4
FLB+2ZHU+CO2 ^LB+nitrate- Бт2Ни+-СО2 ^ZLA+H2O -
FLB+ZHU <-»FLB.ZHU-
где к есть константа скоростей стадий. Последнее уравнение реакции описывает равновесие. Значение Кг вычисляется по формуле
К = [ FLB.ZHU ] / ([ FLB] • [ ZHU ])
и используется для оценки остальных параметров. Скорости элементарных стадий вычисляются следующим образом:
Г = к • [FLB]4 • [СО2 - г2 = к2 • [FLBT] • [ZHU], г3 = К- • к2 • [FLB] • [ZLA]-г4 = к4 • [FLB] • [ZHU]2 - г5 = к4 • [FLB.ZHU]2 • [СО2.
Приток углекислого газа на единицу объема обозначается через Fin, где к1А - коэффициент массопе-реноса, Н - константа Генри и р(002) - парциальное давление углекислого газа. В расчетах значение дав-
ления р(С02) предполагается независимым от углекислого газа [СО2], причем кл =18,7, к2=о,58' кз=о,о9, к4=о,42, К=34'4' к/Д=3'3' =115г83г Н=737 и р(С02) =од
Процесс начинается путем смешивания 0,444 моль/литр реагента Р1_В с 0,007 моль/литр элемента 7Ии. В начальный момент времени концентрация углекислого газа полагается равным 0,00123 моль/литр. По переменной уб начальное условие задается формулой уо,б=Квуо,1уол По остальным компонентам начальные условия равны нулю. Теперь, вводя обозначения концентраций реагентов по формулам у 1=[РЬВ]'
у2=[С02]' у3=[РЬВТ]' у4=[7Ии]' У5 — [7ЬЛ] и уб=[РЬВ.7Ии], получим систему уравнений вида
МУ — f (У)' У(0) — Уо' У (0) — У0' У Е R6' 0 < г < х80.
(6)
Ранг матрицы М равен пяти и она имеет вид
гх 0 0 0 0 0Л 0 х 0 0 0 0
м
0 0 х 0 0 0
0 0 0 х 0 0
0 0 0 0 х 0
0 0 0 0 0 0
Функция Ау) записывается следующим образом:
Л
"2гх + Г2
Г3
г4' Л2 —-0-5гх - Г4 - 0-5г5 + ¥т' /э
Г - Г2 + Гэ '
/4
-г2 + г3
■2Г4 ' /5
Г3 + Г5 ' Л — К, • Ух • У4 - Уб ■
Значения переменных г' х < 1 < 5, и ¥п задаются формулами:
Г — кх • Ух4 • У2^ ' Г2 — к2 • Уз • У4' Г3 — V к2 • Ух ^ У5 '
к
Г4 — к3 • Ух • У42' Г5 — к4 • Уб2 • У22'
¥п — мл • [р(СО2)/н - у2] ■
Начальные условия имеют вид
У0 — (0.444, 0.00х23, 0., 0.007, 0., К, • у0Д • У0,4). у0 — /(У0) ■
Результаты расчетов. Расчеты проводились на PC Core(TM) i7-3770S CPU@3.10GHz с двойной точностью. Задаваемая точность £ расчетов полагалась равной 10-2. Моделирование выполняется на интервале времени [оЛ8о]. Решение в конце интервала интегрирования имеет вид
ух — 0.хх5' у2 — 0.х20 40-2' у3 — 0.хбх'
у4 — 0.366 •Ю-3' у5 — 0.х7х • х0-х' уб — 0.487 •Ю-2 ■
Зависимость концентрации CO2 от времени приведена на рисунке.
Зависимость концентрации CO2 от времени
Заключение. Из результатов сравнения построенного алгоритма интегрирования с методом типа Ро-зенброка для разрешенных задач следует, что эффективности алгоритмов отличаются незначительно. Отсюда можно сделать вывод, что построенный алгоритм можно применять для решения как разрешенных, так и неразрешенных относительно производной задач.
Литература
1. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. - М.: Мир, 1999. - 685 с.
2. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. - М.: Мир, 1990. - 512 с.
3. Gear C.W., PetzoldL. ODE methods for solution differential-algebraic systems // SIAM J. Numerical Analysis. - 1984. - Vol. 21, № 4. - P. 716-728.
4. Новиков Е.А., Шорников Ю.В. Компьютерное моделирование жестких гибридных систем. - Новосибирск: Изд-во НГТУ, 2012. - 451 с.
5. Rosenbrock H.H. Some general implicit processes for the numerical solution of differential equations // Computer. - 1963. - № 5. - P. 329-330.
6. Новиков Е.А., Юматова Л.А. Некоторые методы решения обыкновенных дифференциальных уравнений, неразрешенных относительно производной //ДАН СССР. - 1987. - Т. 295, № 4. - С. 809-812.
7. Демидов Г.В., Новиков Е.А. Оценка ошибки одношаговых методов интегрирования обыкновенных дифференциальных уравнений // Численные методы механической сплошной среды. - 1985. - Т. 16, №1. - С. 27-39.
8. Новиков Е.А. Явные методы для жестких систем. - Новосибирск: Наука, 1997. - 197 с.
9. Mazzia F., Magherini C. Test set for initial value problem solvers. - Department of mathematics University of Bari, Italy, 2008. - 295 p.
